%Exercise series 3
%
%Reproduce the main features of Figure 18.5.

%1) Use the following code to load the time series data

fid = fopen('ts001.bin','r');

%big endian format, ieee-be, 16 bits integer, unsigned, uint16
ts = fread(fid,inf,'uint16','ieee-be');

fclose(fid);

%The number of samples per second is 1200

%2) plot the entire data set (45 mins), a 45 seconds subset and
%a 4.5 seconds subset

%3) plot a histogram of the pixel luminance values. Use 100 bins.
%Plot a histogram of the log10 luminance values; use again 100 bins.


%4) Compute an estimate of the power spectrum of the time series. Read the
%description of the function pwelch to make your task easier. 

%First, subtract the mean value from the time series. Then, split the
%data in 10 segments of 4.5 mins.
%
%Compute the power spectrum over each segment. Use the function pwelch
%with a window of 4096 elements, an overlap of 2048 points and an fft size
%4096. Pass a cut-off frequency of 1200 samples/s to the pwelch function
%and use the 'twosided' option. Use the largest number of data samples that
%is a multiple of 4096 within the 4.5 min data segment. 

%store the power spectrum at positive frequencies and compute its mean and
%SD across the 10 segments. Plot them on a graph.

%Store the base 10 logarithm of the power spectrum at positive frequencies,
%as well as the base 10 logarithm of the positive frequencies. Compute the 
%mean log10(power spectrum) and its SD. Plot them as a function of the base 
%10 logarithm of the positive frequencies.