close all;clc;clear
tic
rand('state',3)

load pcm_B05_12_hpt_0851_monthly.dat
t=1:1188;
t=t';
y=pcm_B05_12_hpt_0851_monthly(:,1);
y0=pcm_B05_12_hpt_0851_monthly(:,1);

for i=1:length(y)-1
    step(i)=y(i+1)-y(i);
end
SDy=std(y);%% Normalize to variance of 1
y=y/SDy;
LL=2^(16);
Y = fft(y,LL);
Pyy = Y.* conj(Y)/(LL);
f=[0:1:round(0.1*LL)]/LL;
%
plot(f,Pyy(1:round(0.1*LL)+1),'k')
hold
%
% RANDOM STEP TRIAL
numtest=10
chknum=numtest/10
for itest = 1:numtest
    if chknum*fix(itest/chknum)==itest
        itest
    end
    yc(1,1) = y(1);% Permute
    RP=randperm(length(y)-1);
    for istep = 1:length(y)-1
        yc(istep+1,1) = yc(istep,1) + step(RP(istep));
    end
    z=yc;
    SDz=std(z);
    z=z/SDz;% Normalize to variance of 1
    Y = fft(z,LL);
    Pzz = Y.* conj(Y)/(LL);
    for isa=1:round(0.1*LL)+1
        Qsave(itest,isa)=Pzz(isa);
    end
end
for isa=1:round(0.1*LL)+1
    ave(isa)=mean(Qsave(:,isa));
    sd(isa)=std(Qsave(:,isa));
end
plot(f(1:round(0.1*LL)+1),ave,'r')

% RANDOM BLOCK TRIAL
chknum=numtest/10
for itest = 1:numtest
    if chknum*fix(itest/chknum)==itest
        itest
    end
    RP=randperm(12);% Permute
    is=0;
    for ib=1:12
        if RP(ib)==12
             for j = 1:41
                 is = is+1;
                 yc(is,1)=y0((RP(ib)-1)*41+j);
             end
        end
    end
    z=yc;
    SDz=std(z);
    z=z/SDz;% Normalize to variance of 1
    Y = fft(z,LL);
    Pzz = Y.* conj(Y)/(LL);
    for isa=1:round(0.1*LL)+1
        Qsave(itest,isa)=Pzz(isa);
    end
end
plot(f(1:round(0.1*LL)+1),ave,'b')

ylabel('Spectral Power')
xlabel('Frequency (cycles/month)')
title('DFTPS')
axis([.0001 0.01 0 12])
print  -depsc dftps.eps

[im re]=max(Pyy(801:60000));
fprintf('The FT imaginary index: %d\n', 800+re )
fprintf('The amplitude of the largest peak is: %0.2f\n', Pyy(re+800) )
fprintf('The frequency of the largest peak is: %0.4f\n', f(800+re) )
fprintf('The period of the largest peak (in years) is: %0.1f\n', (1/f(800+re))/12 )

toc
