% M-file for Sec2 Lec 01.
% Load the Perival and Walden Data Sets
ld = 1 ;
if( ld == 1 ) 
%
%    WIND SPEED TIME SERIES (UPPER PLOT OF FIGURE 2)
%     SOURCE: PMEL/NOAA SEATTLE (BERNIE WALTER)
%     DELTA T: 0.025 SECOND
%     SAMPLE SIZE: 128
    WS = load('WindS.dat');
%     WILLAMETTE RIVER TIME SERIES (TOP PLOT OF FIGURE 505)
%     SOURCE: U. S. GEOLOGICAL SURVEY
%     DELTA T: 1 MONTH
%     SAMPLE SIZE: 395
%     NOTE: FIRST 128 VALUES ALSO PLOTTED IN LOWER PLOT OF FIGURE 2
    WRAll = load('WillametteR.dat');
    WR = WRAll(1:128);
%     ATOMIC CLOCK TIME SERIES (UPPER PLOT OF FIGURE 3)
%     SOURCE: U. S. NAVAL OBSERVATORY (DON PERCIVAL)
%     DELTA T: 1 DAY
%     SAMPLE SIZE: 128
    AC = load('AtomC.dat');
%     OCEAN NOISE TIME SERIES (LOWER PLOT OF FIGURE 3)
%     SOURCE: APPLIED PHYSICS LABORATORY
%     DELTA T: 1 SECOND
%     SAMPLE SIZE: 128
    OD = load('OceanN.dat');
end
% Now plot the data set
set(0,'DefaultAxesFontSize',12);

p1 = 1;
p2 = 1;

if ( p1  == 1 ) 
    h1 = figure(1); set(h1,'Position',[100 300 750 250],'PaperPositionMode','auto');
    plot([1:128],WS,'k'); axis tight
    title('Wind Speed');
    print -dpng WS_ts.png

    h2 = figure(2); set(h2,'Position',[120 280 750 250],'PaperPositionMode','auto');
    plot([1:128],WR,'k'); axis tight
    title('WillametteR');
    print -dpng WR_ts.png

    h3 = figure(3); set(h3,'Position',[140 260 750 250],'PaperPositionMode','auto');
    plot([1:128],AC,'k'); axis tight
    title('Atomic Clock');
    print -dpng AC_ts.png

    h4 = figure(4); set(h4,'Position',[160 240 750 250],'PaperPositionMode','auto');
    plot([1:128],OD,'k'); axis tight
    title('Ocean Noise');
    print -dpng OD_ts.png
end

if ( p2 == 1 ) 
    lag = 1;
    h1 = figure(5); set(h1,'Position',[100 300 250 250],'PaperPositionMode','auto');
    plot(WS(1:128-lag),WS(1+lag:128),'k.'); axis tight; daspect([1 1 1]);
    title('Wind Speed: Lag 1');
    print -dpng WS_lag1.png
    
    h2 = figure(6); set(h2,'Position',[120 280 250 250],'PaperPositionMode','auto');
    plot(WR(1:128-lag),WR(1+lag:128),'k.'); axis tight; daspect([1 1 1]);
    title('WillametteR: Lag 1');
    print -dpng WR_lag1.png
    
    h3 = figure(7); set(h3,'Position',[140 260 250 250],'PaperPositionMode','auto');
    plot(AC(1:128-lag),AC(1+lag:128),'k.'); axis tight; daspect([1 1 1]);
    title('Atomic Clock: Lag 1');
    print -dpng AC_lag1.png
    
    h4 = figure(8); set(h4,'Position',[160 240 250 250],'PaperPositionMode','auto');
    plot(OD(1:128-lag),OD(1+lag:128),'k.'); axis tight; daspect([1 1 1]);
    title('Ocean Noise: Lag 1');
    print -dpng OD_lag1.png

end

p3 = 1;

if( p3 == 1) 
 
% The ACF computation is based on Box, Jenkins, Reinsel, pages 30-34, 188.
    meanWS = mean(WS);
    FWS = fft(WS-meanWS,256);  % 256 = 2^(nextpow2(length(Series)) + 1);
    FWS = FWS .* conj(FWS);
    ACSWS = ifft(FWS);
    ACSWS = real(ACSWS/ACSWS(1)); % Normalize and take real part
    % Now compute the using Pearson Directly
    VarWS = (WS-meanWS)'*(WS-meanWS);
    for k = 0:32
        rhoWS(k+1) = (WS(1:128-k)-meanWS)'*(WS(1+k:128)-meanWS)/VarWS;
    end
    h4 = figure(9); set(h4,'Position',[100 300 350 250],'PaperPositionMode','auto');
    plot([0:32],ACSWS([1:33]),'k-',[0:32],rhoWS([1:33]),'r.'); axis tight; grid on;
    title('Wind Speed: ACS'); xlabel('Lag');
    print -dpng WS_ACS.png

    meanWR = mean(WR);
    FWR = fft(WR-meanWR,256);  % 256 = 2^(nextpow2(length(Series)) + 1);
    FWR = FWR .* conj(FWR);
    ACSWR = ifft(FWR);
    ACSWR = real(ACSWR/ACSWR(1)); % Normalize and take real part
    % Now compute the using Pearson Directly
    VarWR = (WR-meanWR)'*(WR-meanWR);
    for k = 0:32
        rhoWR(k+1) = (WR(1:128-k)-meanWR)'*(WR(1+k:128)-meanWR)/VarWR;
    end
    h4 = figure(10); set(h4,'Position',[120 280 350 250],'PaperPositionMode','auto');
    plot([0:32],ACSWR([1:33]),'k-',[0:32],rhoWR([1:33]),'r.'); axis tight; grid on;
    title('WillametteR: ACS'); xlabel('Lag');
    print -dpng WR_ACS.png

    meanAC = mean(AC);
    FAC = fft(AC-meanAC,256);  % 256 = 2^(nextpow2(length(Series)) + 1);
    FAC = FAC .* conj(FAC);
    ACSAC = ifft(FAC);
    ACSAC = real(ACSAC/ACSAC(1)); % Normalize and take real part
    % Now compute the using Pearson Directly
    VarAC = (AC-meanAC)'*(AC-meanAC);
    for k = 0:32
        rhoAC(k+1) = (AC(1:128-k)-meanAC)'*(AC(1+k:128)-meanAC)/VarAC;
    end
    h4 = figure(11); set(h4,'Position',[140 260 350 250],'PaperPositionMode','auto');
    plot([0:32],ACSAC([1:33]),'k-',[0:32],rhoAC([1:33]),'r.'); axis tight; grid on;
    title('Atomic Clock: ACS'); xlabel('Lag');
    print -dpng AC_ACS.png

    meanOD = mean(OD);
    FOD = fft(OD-meanOD,256);  % 256 = 2^(nextpow2(length(Series)) + 1);
    FOD = FOD .* conj(FOD);
    ACSOD = ifft(FOD);
    ACSOD = real(ACSOD/ACSOD(1)); % Normalize and take real part
    % Now compute the using Pearson Directly
    VarOD = (OD-meanOD)'*(OD-meanOD);
    for k = 0:32
        rhoOD(k+1) = (OD(1:128-k)-meanOD)'*(OD(1+k:128)-meanOD)/VarOD;
    end
    h4 = figure(12); set(h4,'Position',[160 240 350 250],'PaperPositionMode','auto');
    plot([0:32],ACSOD([1:33]),'k-',[0:32],rhoOD([1:33]),'r.'); axis tight; grid on;
    title('Ocean Noise: ACS'); xlabel('Lag');
    print -dpng OD_ACS.png
    
end
   
p4 = 1;
if( p4 == 1 )
    
    % Spectra
    FWS = fft(WS)*(2/128);  % Matlab fft does not multiply by 2/N
    FWS = FWS .* conj(FWS);
    dB_FWS = 10*log10(FWS);
    h1 = figure(13); set(h1,'Position',[100 300 700 250],'PaperPositionMode','auto');
    plot([1:64]/128,dB_FWS(2:65)); xlabel('Frequency cycles/sample time'); axis tight; grid on;
    ylabel('dB 10*log10(Sj)'); title('Wind Speed: Power Spectrum');
    print -dpng WS_PS.png
    
    % Spectra
    FWR = fft(WR)*(2/128);  % Matlab fft does not multiply by 2/N
    FWR = FWR .* conj(FWR);
    dB_FWR = 10*log10(FWR);
    h2 = figure(14); set(h2,'Position',[120 280 700 250],'PaperPositionMode','auto');
    plot([1:64]/128,dB_FWR(2:65)); xlabel('Frequency cycles/sample time'); axis tight; grid on;
    ylabel('dB 10*log10(Sj)'); title('WillametteR: Power Spectrum');
    print -dpng WR_PS.png

    % Spectra
    FAC = fft(AC)*(2/128);  % Matlab fft does not multiply by 2/N
    FAC = FAC .* conj(FAC);
    dB_FAC = 10*log10(FAC);
    h3 = figure(15); set(h3,'Position',[140 260 700 250],'PaperPositionMode','auto');
    plot([1:64]/128,dB_FAC(2:65)); xlabel('Frequency cycles/sample time'); axis tight; grid on;
    ylabel('dB 10*log10(Sj)'); title('Atomic Clock: Power Spectrum');
    print -dpng AC_PS.png

    % Spectra
    FOD = fft(OD)*(2/128);  % Matlab fft does not multiply by 2/N
    FOD = FOD .* conj(FOD);
    dB_FOD = 10*log10(FOD);
    h4 = figure(16); set(h4,'Position',[160 240 700 250],'PaperPositionMode','auto');
    plot([1:64]/128,dB_FOD(2:65)); xlabel('Frequency cycles/sample time'); axis tight; grid on;
    ylabel('dB 10*log10(Sj)'); title('Ocean Noise: Power Spectrum');
    print -dpng OD_PS.png
    
end