close all;clc;clear
tic
rand('state',7)
format short g
%
load pcm_B05_12_hpt_0851_monthly.dat
%%%
%%% Find the significance of the 12 month cycle in various 48 
%%% month segments of the time series.
%%%
elements=48
t=1:elements; t=t';
%
%Soil Moisture Data: series 1-3
data(1:elements,1)=pcm_B05_12_hpt_0851_monthly(600:647,2);
data(1:elements,2)=pcm_B05_12_hpt_0851_monthly(648:695,2);
data(1:elements,3)=pcm_B05_12_hpt_0851_monthly(696:743,2);
%Soil Moisture Data: series 4-7
data(1:elements,4)=pcm_B05_12_hpt_0851_monthly(809:856,1);
data(1:elements,5)=pcm_B05_12_hpt_0851_monthly(905:952,1);
data(1:elements,6)=pcm_B05_12_hpt_0851_monthly(1001:1048,1);
num_series=6;
fprintf('Searching %d time series of length %d\n', num_series, elements)
%
xp=zeros(size(t));
l=2^(16);
f=[0:1:round(0.1*l)]/l;
iterations=10000
%
for series=1:num_series
    series
    x1=data(:,series);
    x1 = x1/std(x1);
    x2=fft(x1,l);
    p=x2.* conj(x2)/(l);
    %%Find the max peak near 12 months
    [a b]=max(p(4501:6000));  
    ifreq1=4500+b;
    %
    x=data(:,series).';
    x0=x;
    %
    for i=1:length(x)-1
        step(i)=x(i+1)-x(i);
    end
    rw_count=0;
    rb_count=0;
    for itest = 1:iterations
    %%% Monte Carlo Random Walk Permutation
        xp(1,1) = x(1);
        rp=randperm(length(x)-1);
        for istep = 1:length(x)-1
            xp(istep+1,1) = xp(istep,1) + step(rp(istep));
        end
        z=xp/std(xp);
        xmc=fft(z,l);
        pz=xmc.* conj(xmc)/(l);
    
        if pz(ifreq1) >= p(ifreq1)
            rw_count = rw_count + 1;
        end
    %%% Monte Carlo Random Block Permutation, set for 4 blocks 
    %%% of length 12 (months)
        rp=randperm(12);
        is=0;
        for ib=1:12
            if rp(ib)==12
                for j = 1:4
                    is = is+1;
                    xp(is)=x0((rp(ib)-1)*4+j);
                end
            end
        end
        z=xp/std(xp);
        xmc=fft(z,l);
        pz=xmc.* conj(xmc)/(l);
        if pz(ifreq1) >= p(ifreq1)
            rb_count = rb_count + 1;
        end
    end
    fprintf('\n')
    fprintf('The subset no. is: %d\n', series )
    fprintf('FT imaginary index: %d, Peak amplitude: %0.5f\n', 4500+b, a )
    fprintf('Frequency: %0.4f, Period (in months): %0.1f\n', f(4500+b), 1/f(4500+b) )
    fprintf('Num. times MCRW exceeded peak amp: %d, MCRB: %d\n', rw_count, rb_count )
    fprintf('MCRW significance: %0.3f, MCRB significance: %0.3f\n', rw_count/iterations, rb_count/iterations )
end
%%%
%%% Find the significance of the 12 month cycle anywhere  
%%% in the time series.
%%%
elements=1188
t=1:elements; t=t';
%Soil Moisture Data: series 1
data(1:elements,1)=pcm_B05_12_hpt_0851_monthly(:,2);
%Soil Moisture Data: series 2
data(1:elements,2)=pcm_B05_12_hpt_0851_monthly(:,1);
num_series=2;
%
xp=zeros(size(t));
l=2^(16);
f=[0:1:round(0.1*l)]/l;
%iterations=10
%
for series=1:num_series
    series
    x1=data(:,series);
    x1 = x1/std(x1);
    x2=fft(x1,l);
    p=x2.* conj(x2)/(l);
    %%Find the max peak near 12 months
    [a b]=max(p(4501:6000));  
    ifreq1=4500+b;
    %
    x=data(:,series).';
    x0=x;
    %
    for i=1:length(x)-1
        step(i)=x(i+1)-x(i);
    end
    rw_count=0;
    rb_count=0;
    for itest = 1:iterations
    %%% Monte Carlo Random Walk Permutation
        xp(1,1) = x(1);
        rp=randperm(length(x)-1);
        for istep = 1:length(x)-1
            xp(istep+1,1) = xp(istep,1) + step(rp(istep));
        end
        z=xp/std(xp);
        xmc=fft(z,l);
        pz=xmc.* conj(xmc)/(l);
    
        if pz(ifreq1) >= p(ifreq1)
            rw_count = rw_count + 1;
        end
    %%% Monte Carlo Random Block Permutation, set for 99 blocks 
    %%% of length 12 (months)
        rp=randperm(12);
        is=0;
        for ib=1:12
            if rp(ib)==12
                for j = 1:99
                    is = is+1;
                    xp(is)=x0((rp(ib)-1)*99+j);
                end
            end
        end
        z=xp/std(xp);
        xmc=fft(z,l);
        pz=xmc.* conj(xmc)/(l);
        if pz(ifreq1) >= p(ifreq1)
            rb_count = rb_count + 1;
        end
    end
    fprintf('\n')
    fprintf('The subset no. is: %d\n', series )
    fprintf('FT imaginary index: %d, Peak amplitude: %0.5f\n', 4500+b, a )
    fprintf('Frequency: %0.4f, Period (in months): %0.1f\n', f(4500+b), 1/f(4500+b) )
    fprintf('Num. times MCRW exceeded peak amp: %d, MCRB: %d\n', rw_count, rb_count )
    fprintf('MCRW significance: %0.3f, MCRB significance: %0.3f\n', rw_count/iterations, rb_count/iterations )
end
toc
