% M-file for Lecture 6.

p1 = 1;
if( p1 == 1 ) 
    % Generate the spectrum for Fig 145
    f = [0:0.005:4]; lf = length(f); 
    SXa = 1.9*exp(-2*f.^2) + exp(-6*(f-1.25).^2);
    %
    % Now generate sub-sampled spectrum. Notice the way the spectrum is
    % folded.  If were continuing the next block wouuld be added in
    % increasing order, and then reversed again.
    SX1 = zeros(1,lf);
    SX1([1:200]) = SXa([1:200])+SXa([400:-1:201]);
    % Now generate for f=2
    SX2 = zeros(1,lf);
    SX2([1:400]) = SXa([1:400])+SXa([800:-1:401]);
    
   
    h1 = figure(1); set(h1,'Position',[100 200 750 500],'PaperPositionMode','auto');
    subplot(2,1,1); hold off
    plot(f([1:601]),SXa([1:601]),'r'); hold on;
    plot(f([1:601]),SX1([1:601]),'k');
    title('Folded Spectrum for Nyquist f = 1 '); xlabel('Frequency'), ylabel('Spectral Density');
    subplot(2,1,2); hold off
    plot(f([1:601]),SXa([1:601]),'r'); hold on;
    plot(f([1:601]),SX2([1:601]),'k');
    title('Folded Spectrum for Nyquist f = 2 '); xlabel('Frequency'), ylabel('Spectral Density');
    
    print -dpng S2L06_F145.png 
    
end

p2 =1 ; 
if( p2 == 1 )
    randn('seed',120);
    % Generate AR(4) X(t) =
    % 2.7607X(t-1)-3.8106X(t-2)+2.65235X(t-3)-0.9238*X(t-4)+e

    f=[0:0.001:.5];
    denom = 1-2.7607*exp(-i*2*pi*f)+3.8106*exp(-i*2*pi*2*f)-2.65235*exp(-i*2*pi*3*f)+0.9238*exp(-i*2*pi*4*f);
    Sf = 1./(conj(denom).*denom); dB_Sf = 10*log10(Sf);
    h1 = figure(2); set(h1,'Position',[100 200 750 500],'PaperPositionMode','auto'); hold off; 
    subplot(2,1,1); 
    plot(f,dB_Sf,'k'); axis tight;
    title('SDF for AR(4) model'); xlabel('Frequency'); ylabel('sdf (dB)');
    % Now compute the acvs
    t=[0:50]; nt = length(t); df = 0.001;
    for n = 1:nt
        avcs(n)=sum(Sf.*exp(i*2*pi*f*t(n)))*df;
    end    
    subplot(2,1,2); axis tight;
    plot(t,2*avcs,'k-+');title('ACVS for AR(4)'); xlabel('Time'); ylabel('ACVS');
    print -dpng S2L06_P148a.png 
        
    % Part two: Now add any oscillating term at f = 0.35
    Sfa = Sf + 1*exp(-(f-0.35).^2/(0.01^2)); dB_Sfa=10*log10(Sfa);
    h1 = figure(3); set(h1,'Position',[110 190 750 500],'PaperPositionMode','auto'); hold off; 
    subplot(2,1,1); 
    plot(f,dB_Sfa,'k'); axis tight;
    title('SDF for AR(4) model + unit 0.35 signal'); xlabel('Frequency'); ylabel('sdf (dB)');
    % Now compute the acvs
    t=[0:50]; nt = length(t); df = 0.001;
    for n = 1:nt
        avcs(n)=sum(Sfa.*exp(i*2*pi*f*t(n)))*df;
    end
    
    subplot(2,1,2); axis tight;
    plot(t,2*avcs,'k-+');title('ACVS for AR(4)'); xlabel('Time'); ylabel('ACVS');
    print -dpng S2L06_P148b.png 
    
    % Part 3: Now add any oscillating term at f = 0.35 with larger
    % amplitude
    Sfa = Sf + 1e3*exp(-(f-0.35).^2/(0.01^2)); dB_Sfa=10*log10(Sfa);
    h1 = figure(4); set(h1,'Position',[110 190 750 500],'PaperPositionMode','auto'); hold off; 
    subplot(2,1,1); 
    plot(f,dB_Sfa,'k'); axis tight;
    title('SDF for AR(4) model + 1e3 0.35 signal'); xlabel('Frequency'); ylabel('sdf (dB)');
    % Now compute the acvs
    t=[0:50]; nt = length(t); df = 0.001;
    for n = 1:nt
        avcs(n)=sum(Sfa.*exp(i*2*pi*f*t(n)))*df;
    end
    
    subplot(2,1,2); axis tight;
    plot(t,2*avcs,'k-+');title('ACVS for AR(4) + 1e3 signal'); xlabel('Time'); ylabel('ACVS');
    print -dpng S2L06_P148c.png 

    
end

p3 = 1
if( p3 == 1 )
    % Generate AR(4) X(t) =
    % 2.7607X(t-1)-3.8106X(t-2)+2.65235X(t-3)-0.9238*X(t-4)+e
    X4 = zeros(1128,1);
    
    for k = 1:1128
        X4(k+4)=2.7607*X4(k+3)-3.8106*X4(k+2)+2.65235*X4(k+1)-0.9238*X4(k)+randn(1);
    end
    
    % Now compute the ACVS from sample
    for k = 0:256
        acvs_ts(k+1) = (X4(100:1024-k)'*X4(100+k:1024))/256/8;
    end
    h1 = figure(5); set(h1,'Position',[100 180 750 400],'PaperPositionMode','auto'); hold off; 
    
    hold off; plot([0:256],acvs_ts,'r+-');
    f=[0:0.001:.5];
    denom = 1-2.7607*exp(-i*2*pi*f)+3.8106*exp(-i*2*pi*2*f)-2.65235*exp(-i*2*pi*3*f)+0.9238*exp(-i*2*pi*4*f);
    Sf = 1./(conj(denom).*denom); dB_Sf = 10*log10(Sf);
    t=[0:256]; nt = length(t); df = 0.001;
    for n = 1:nt
        acvs(n)=sum(Sf.*exp(i*2*pi*f*t(n)))*df;
    end  
    hold on; plot([0:256],acvs,'k'); axis tight; hold off;
    title('ACVS from sample and theory'); xlabel('Lag'); ylabel('ACVS');
    print -dpng S2L06_ACVS.png 
    
    Gp = fft(X4([100:1124]))/16; dB_SG = 10*log10(conj(Gp).*Gp);
    h1 = figure(6); set(h1,'Position',[110 190 750 400],'PaperPositionMode','auto'); hold off;
    plot([0:512]/1024,dB_SG(1:513),'k');
    hold on; 
    plot(f,dB_Sf,'r'); 
    title('Spectrum from Sample');  xlabel('Frequency'); ylabel('sdf (dB)');
    
    print -dpng S2L06_SDF.png 
    
end