% M-file for Lecture 8.

p1 = 0;
if( p1 == 1 ); 
    % Generate AR process.  First AR(2) Xt,2=0.75X(t-1)-0.5X(t-2)+e(t)
    rand('seed',120);
    X2 = zeros(1126,1); 
    for k = 1:1124  % Run the process for 100 steps to get past start up.  
        % Alternative would be to generate the first two values with
        % correct covariance matrix
        X2(k+2)=-0.75*X2(k+1)-0.5*X2(k)+randn(1);
    end
    h1 = figure(6); set(h1,'Position',[100 300 750 200],'PaperPositionMode','auto');
    plot([1:1024],X2([103:1126]),'k'); axis tight; ylim([-5 5]); 
    title('AR(2) = -0.75X(t-1)-0.5X(t-2)+\epsilon');
    % print -dpng S2L02_AR2.png
    
    % 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
    h1 = figure(7); set(h1,'Position',[120 280 750 200],'PaperPositionMode','auto');
    plot([1:1024],X4([105:1128]),'k'); axis tight; ylim([-100 100]); 
    title('AR(4) = 2.7607*X4(k+3))-3.8106*X4(k+2)+2.65235*X4(k+1)-0.9238*K4(k)+\epsilon');
    % print -dpng S2L02_AR4.png

end

p2 = 0;
if ( p2 == 1 )
    % Show the Fejer's Kernel.
    f=[-0.5:0.001:0.5]; Dt=1;
    N=4; F4 = Dt*(sin(N*pi*f*Dt).*sin(N*pi*f*Dt))./(N*sin(pi*f*Dt).*sin(pi*f*Dt));
    dB_F4=10*log10(F4);
    h1 = figure(1); set(h1,'Position',[100 300 750 200],'PaperPositionMode','auto');
    plot(f,dB_F4,'k'); axis tight; ylim([-40 20]); title('Fejers Kernel N=4 '); ylabel('dB');
    print -dpng S2L08_F04.png
    
    N=16; F16 = Dt*(sin(N*pi*f*Dt).*sin(N*pi*f*Dt))./(N*sin(pi*f*Dt).*sin(pi*f*Dt));
    dB_F16=10*log10(F16);
    h1 = figure(2); set(h1,'Position',[120 280 750 200],'PaperPositionMode','auto');
    plot(f,dB_F16,'k'); axis tight; ylim([-40 20]); title('Fejers Kernel N=16 '); ylabel('dB');
    print -dpng S2L08_F16.png
    
    N=64; F64 = Dt*(sin(N*pi*f*Dt).*sin(N*pi*f*Dt))./(N*sin(pi*f*Dt).*sin(pi*f*Dt));
    dB_F64=10*log10(F64);
    h1 = figure(3); set(h1,'Position',[140 260 750 200],'PaperPositionMode','auto');
    plot(f,dB_F64,'k'); axis tight; ylim([-40 20]); title('Fejers Kernel N=64 '); ylabel('dB');
    xlabel('Frequency');
    print -dpng S2L08_F64.png
end

p3 = 0;
if ( p3 == 1 )
    % Biases in Periodograms for AR2 and AR4 processes.
    eps = 1.e-10;
    f=[-0.5+eps:0.0005:0.5]; Dt=1;
    N=16; F16 = Dt*(sin(N*pi*f*Dt).*sin(N*pi*f*Dt))./(N*sin(pi*f*Dt).*sin(pi*f*Dt));
    N=64; F64 = Dt*(sin(N*pi*f*Dt).*sin(N*pi*f*Dt))./(N*sin(pi*f*Dt).*sin(pi*f*Dt));
    N=256; F256 = Dt*(sin(N*pi*f*Dt).*sin(N*pi*f*Dt))./(N*sin(pi*f*Dt).*sin(pi*f*Dt));
    N=1024; F1024 = Dt*(sin(N*pi*f*Dt).*sin(N*pi*f*Dt))./(N*sin(pi*f*Dt).*sin(pi*f*Dt));
    
    % AR2: X2(k+2)=0.75*X2(k+1)-0.5*X2(k)+randn(1);
    denom = 1-0.75*exp(-i*2*pi*f)+0.5*exp(-i*2*pi*2*f);
    SfAR2 = 1./(conj(denom).*denom); dB_SfAR2 = 10*log10(SfAR2);
    lf = length(f);
    
    % Now convolve to get reponse
    SfAR2_16 = conv(F16,SfAR2)/(length(F16));
    dB_SfAR2_16=10*log10(SfAR2_16);
    h1 = figure(4); set(h1,'Position',[100 300 750 600],'PaperPositionMode','auto');
    subplot(2,1,1); 
    plot(f([lf/2:lf]),dB_SfAR2([lf/2:lf]),'k',f([lf/2:lf]),dB_SfAR2_16([lf:lf+lf/2]),'r');
    axis tight; ylim([-10 10]); title('AR(2), N=16 '); ylabel('dB'); xlabel('Frequency');
    
    SfAR2_64 = conv(F64,SfAR2)/(length(F64));
    dB_SfAR2_64=10*log10(SfAR2_64);
    subplot(2,1,2); 
    plot(f([lf/2:lf]),dB_SfAR2([lf/2:lf]),'k',f([lf/2:lf]),dB_SfAR2_64([lf:lf+lf/2]),'r');
    axis tight; ylim([-10 10]);title('AR(2), N=64 '); ylabel('dB'); xlabel('Frequency');
    
     print -dpng S2L08_AR2.png
          
    % AR(4) Process(t) =
    % A$4 = 2.7607X(t-1)-3.8106X(t-2)+2.65235X(t-3)-0.9238*X(t-4)+e
    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);
    SfAR4 = 1./(conj(denom).*denom); dB_SfAR4 = 10*log10(SfAR4);
    
    % Now convolve to get reponse
    SfAR4_16 = conv(F16,SfAR4)/(length(F16));
    dB_SfAR4_16=10*log10(SfAR4_16);
    h2 = figure(5); set(h2,'Position',[120 280 750 600],'PaperPositionMode','auto');
    subplot(2,1,1); 
    plot(f([lf/2:lf]),dB_SfAR4([lf/2:lf]),'k',f([lf/2:lf]),dB_SfAR4_16([lf:lf+lf/2]),'r');
    axis tight; ylim([-40 60]);title('AR(2), N=16 '); ylabel('dB'); xlabel('Frequency');
    
    SfAR4_64 = conv(F64,SfAR4)/(length(F64));
    dB_SfAR4_64=10*log10(SfAR4_64);
    subplot(2,1,2); 
    plot(f([lf/2:lf]),dB_SfAR4([lf/2:lf]),'k',f([lf/2:lf]),dB_SfAR4_64([lf:lf+lf/2]),'r');
    axis tight; ylim([-40 60]); title('AR(2), N=64 '); ylabel('dB'); xlabel('Frequency');
    
    print -dpng S2L08_AR4a.png
    
    SfAR4_256 = conv(F256,SfAR4)/(length(F256));
    dB_SfAR4_256=10*log10(SfAR4_256);
    h2 = figure(6); set(h2,'Position',[140 260 750 600],'PaperPositionMode','auto');
    subplot(2,1,1); 
    plot(f([lf/2:lf]),dB_SfAR4([lf/2:lf]),'k',f([lf/2:lf]),dB_SfAR4_256([lf:lf+lf/2]),'r');
    axis tight; ylim([-40 60]);title('AR(2), N=256 '); ylabel('dB'); xlabel('Frequency');
    
    SfAR4_1024 = conv(F1024,SfAR4)/(length(F1024));
    dB_SfAR4_1024=10*log10(SfAR4_1024);
    subplot(2,1,2); 
    plot(f([lf/2:lf]),dB_SfAR4([lf/2:lf]),'k',f([lf/2:lf]),dB_SfAR4_1024([lf:lf+lf/2]),'r');
    axis tight; ylim([-40 60]); title('AR(2), N=1024 '); ylabel('dB'); xlabel('Frequency');

    print -dpng S2L08_AR4b.png

end

p4 = 1;
if ( p4 == 1 )
    % Look at cosine tapers.
    N  = 64;
    f = [0:0.5/1023:0.5]; lf = length(f);
    
    % Rectangular taper
    ht = ones(64,1)/sqrt(64);
    Hf = fft(ht,2*lf);
    H = conj(Hf).*Hf; db_H = 10*log10(H);
    h1 = figure(7); set(h1,'Position',[100 100 750 550],'PaperPositionMode','auto');
    subplot(2,1,1);
    plot([0:63],ht,'k.'); axis tight; title('Rectangular Taper '); xlabel('Time'); ylabel('h_t');
    ylim([0 0.25]);
    subplot(2,1,2);
    plot(f,db_H([1:lf]),'k'); axis tight; ylim([-100 20]); xlabel('Frequency'); ylabel('|H(f)|^2 dB')
    
    print -dpng S2L08_Taper00.png
    
    % Try as 20% cosine
    p = 0.2 ; Np = floor(p*64/2);
    ht = ones(64,1);
    for k = 1:Np
        ht(k) = (1-cos(2*pi*(k-1)/(2*Np+1)))/2;
        ht(65-k)= (1-cos(2*pi*(k-1)/(2*Np+1)))/2;
    end
    maght = sum(ht.*ht); ht = ht/sqrt(maght);
    Hf = fft(ht,2*lf);
    H = conj(Hf).*Hf; db_H = 10*log10(H);
    h2 = figure(8); set(h2,'Position',[120 80 750 550],'PaperPositionMode','auto');
    subplot(2,1,1); 
    plot([0:63],ht,'k.'); axis tight; title('Cosine Taper p=0.2 '); xlabel('Time'); ylabel('h_t');
    ylim([0 0.25]);
    subplot(2,1,2);
    plot(f,db_H([1:lf]),'k'); axis tight; ylim([-100 20]); xlabel('Frequency'); ylabel('|H(f)|^2 dB')
    print -dpng S2L08_Taper20.png
    
    % Try with 50%
    p = 0.5 ; Np = floor(p*64/2);
    ht = ones(64,1);
    for k = 1:Np
        ht(k) = (1-cos(2*pi*(k-1)/(2*Np+1)))/2;
        ht(65-k)= (1-cos(2*pi*(k-1)/(2*Np+1)))/2;
    end
    maght = sum(ht.*ht); ht = ht/sqrt(maght);
    Hf = fft(ht,2*lf);
    H = conj(Hf).*Hf; db_H = 10*log10(H);
    h3 = figure(9); set(h3,'Position',[140 60 750 550],'PaperPositionMode','auto');
    subplot(2,1,1); 
    plot([0:63],ht,'k.'); axis tight; title('Cosine Taper p=0.5 '); xlabel('Time'); ylabel('h_t');
    ylim([0 0.25]);
    subplot(2,1,2);
    plot(f,db_H([1:lf]),'k'); axis tight; ylim([-100 20]); xlabel('Frequency'); ylabel('|H(f)|^2 dB')
    print -dpng S2L08_Taper50.png
    
    % Try with 50%
    p = 1.0 ; Np = floor(p*64/2);
    ht = ones(64,1);
    for k = 1:Np
        ht(k) = (1-cos(2*pi*(k-1)/(2*Np+1)))/2;
        ht(65-k)= (1-cos(2*pi*(k-1)/(2*Np+1)))/2;
    end
    maght = sum(ht.*ht); ht = ht/sqrt(maght);
    Hf = fft(ht,2*lf);
    H = conj(Hf).*Hf; db_H = 10*log10(H);
    h4 = figure(10); set(h4,'Position',[140 60 750 550],'PaperPositionMode','auto');
    subplot(2,1,1); 
    plot([0:63],ht,'k.'); axis tight; title('Cosine Taper p=1.0 Hanning Taper '); xlabel('Time'); ylabel('h_t');
    ylim([0 0.25]);
    subplot(2,1,2);
    plot(f,db_H([1:lf]),'k'); axis tight; ylim([-100 20]); xlabel('Frequency'); ylabel('|H(f)|^2 dB')
    print -dpng S2L08_Taper100.png
end