% This code confirms 4-season periodicity of the the rainfall data for 
% Miles (see Ex. 22.32, p.412) through the use of the Fourier transform.
close all     % closes all windows opened during the previous runs
clear all     % clears the memory from the results of previous runs
clc           % clears screen
%Example 22.19
l=1;   % SOI            
l=2;   % rainfall for Dalby            
%l=3;   % rainfall for Miles            
soirain
sd=sdm(:,l)-mean(sdm(:,l));
n=40
m=length(sd)-n+1
[i,j]=meshgrid(1:n,0:m-1);
w=sd(i+j);
tot_var=sum(w(:).^2)/m
[u,s,v]=svd(w,0);

plot(j(1:n,1:1),w(1,1:n))
hold on
plot(j(1:n,1:1),w(41,1:n)+30)
plot(j(1:n,1:1),w(81,1:n)+60)
plot(j(1:n,1:1),w(121,1:n)+90)
hold off
pause
subplot(1,2,1)
plot(diag(s),'.'),xlabel('j'),ylabel('\sigma_j')
subplot(1,2,2)
plot(cumsum(diag(s).^2)/sum(diag(s).^2),'.')
axis([0 n 0 1]),xlabel('k')
ylabel('\Sigma_{j=1}^k\lambda_j/\Sigma_{j=1}^n\lambda_j') 
pause
subplot(3,1,1)
hold on
plot(v(:,1))
plot(v(:,2),'--'),axis([0 n 1.05*min(v(:,2)) 1.05*max(v(:,2))])
subplot(3,1,2)
hold on
plot(v(:,3))
plot(v(:,4),'--'),axis([0 n 1.05*min(v(:,4)) 1.05*max(v(:,4))])
subplot(3,1,3)
hold on
plot(v(:,5))
plot(v(:,6),'--'),axis([0 n 1.05*min(v(:,6)) 1.05*max(v(:,6))])
xlabel('Seasons') 
pause

%Fourier transform of the same data:
n=length(sd);
w=(0:n-1)/n;               % frequency in cycles per season
y=fft(sd);                 % Fourier transform of the rainfall data
pwr=abs(y).^2;             % Power spectrum 
subplot(2,1,1)
plot(w(1:fix(n/2)),pwr(1:fix(n/2)))     % Power spectrum plot. Note that
xlabel('Frequency, cycles per season')  % only half of the frequency 
                                        % range is plotted because the 
                                        % frequencies higher than 0.5
                                        % cycles per season cannot be 
                                        % resolved (distinguished from 
                                        % lower frequencies)
subplot(2,1,2)
plot(1./w(2:fix(n/2)),pwr(2:fix(n/2))),axis([0 30 0 1.05*max(pwr(2:fix(n/2)))])
xlabel('Period, seasons per cycle')     % The most prominant cycle has
                                        % the largest power. Period T is
                                        % T=1/w.