% Open file 'Tsup.mat' with Lake Superior air temperatures
t=Tsup(:,1);             % copy the series into separate variables
y=Tsup(:,2);

L = max(size(y));                     % Length of signal

plot(t/24,y);
% plot series (time divided by 24 to plot in days rather than hours)
title('Lake Superior air temperature','FontSize',20);
xlabel('time (days)','FontSize',20);
ylabel('Temperature (deg C)','FontSize',20);

%%

NFFT = 2^nextpow2(L); % Next power of 2 from length of y -- this is the length of series that fft analyzes

Y = fft(y);          % run fft to obtain the amplitudes of Fourier components

% Look at Y in the Array Editor. The elements of Y are *complex numbers*.
% They have a real and an imaginary part, which correspond to the Fourier
% coefficients Am and Bm, respectively.
% The first element is A0. It is always real and describes the time average of the series. 

f = 1/2*linspace(0,1,NFFT/2);     % define a frequency vector that corresponds to data in Y.

% Plot Fourier spectrum. The 'fft' function generates output Y, which
% contains two symmetric spectra, so we only need the first half of Y.

plot(f(1:NFFT/2)*365*24,abs(Y(1:NFFT/2)))

title('Fourier power spectrum','FontSize',20)
xlabel('Frequency (1/yr)','FontSize',20)
ylabel('Fourier power','FontSize',20)

% The 'abs' operator takes the absolute value of Y, i.e. sqrt(Am^2+Bm^2).
% The frequency axis is multiplied by 365*24 to plot it in inverse years,
% rather than inverse hours.

% You can see a lot of low-frequency components. These result from
% low-frequency trends in the data and just contaminate the results. 
% This is why it is important to first remove linear trends! 

% Zoom in on the plot and investigate the components. Find the component,
% corresponding to the annual signal. You will need to zoom in strongly!


%%
% plot the Fourier spectrum against time

plot(1./f(1:NFFT/2)/24,abs(Y(1:NFFT/2))) 
xlabel('period (days)')

% find the important components and check their corresponding periods.