% M-file for Lecture 9.

p1 = 0;
if( p1 == 1 ); 
    % Generate sdf of white noise
    randn('seed',101);
    wn_all = randn(2048,1);
    wn_64 = wn_all([1:32:2048]);
    fw64 = fft(wn_64); Sw64 = conj(fw64).*fw64/64; dB_Sw64 = 10*log10(Sw64); f = [0:0.5/64:0.5-1e-6];
    fw64f = fft(wn_64,2048); Sw64f = conj(fw64f).*fw64f/64; dB_Sw64f = 10*log10(Sw64f); ff = [0:0.5/2048:0.5-1e-6];
    
    h1 = figure(1); set(h1,'Position',[100 300 750 200],'PaperPositionMode','auto');
    plot(f,dB_Sw64,'r+',ff,dB_Sw64f,'k'); axis tight; ylim([-40 20]); 
    
    title('White Noise N=64 '); xlabel('Frequency'); ylabel('dB');
    print -dpng S2L09_WN0064.png
    
    % Run with 128 samples
    wn_128 = wn_all([1:16:2048]);
    fw128 = fft(wn_128); Sw128 = conj(fw128).*fw128/128; dB_Sw128 = 10*log10(Sw128); f = [0:0.5/128:0.5-1e-6];
    fw128f = fft(wn_128,2048); Sw128f = conj(fw128f).*fw128f/128; dB_Sw128f = 10*log10(Sw128f); ff = [0:0.5/2048:0.5-1e-6];
    
    h2 = figure(2); set(h2,'Position',[120 280 750 200],'PaperPositionMode','auto');
    plot(f,dB_Sw128,'r+',ff,dB_Sw128f,'k'); axis tight; ylim([-40 20]); 
    
    title('White Noise N=128 '); xlabel('Frequency'); ylabel('dB');
    print -dpng S2L09_WN0128.png
    
    % Run with 1024 samples
    wn_1024 = wn_all([1:2:2048]);
    fw1024 = fft(wn_1024); Sw1024 = conj(fw1024).*fw1024/1024; dB_Sw1024 = 10*log10(Sw1024); f = [0:0.5/1024:0.5-1e-6];
    fw1024f = fft(wn_1024,2048); Sw1024f = conj(fw1024f).*fw1024f/1024; dB_Sw1024f = 10*log10(Sw1024f); ff = [0:0.5/2048:0.5-1e-6];
    
    h3 = figure(3); set(h3,'Position',[140 260 750 200],'PaperPositionMode','auto');
    plot(f,dB_Sw1024,'r+',ff,dB_Sw1024f,'k'); axis tight; ylim([-40 20]); 
    
    title('White Noise N=1024 '); xlabel('Frequency'); ylabel('dB');
    print -dpng S2L09_WN1024.png
        
end

p2 = 0;
if( p2 == 1 )
    % Generate AR(4) X(t) =
    % 2.7607X(t-1)-3.8106X(t-2)+2.65235X(t-3)-0.9238*X(t-4)+e
    randn('seed',100);  % 102 generates an unbiased sample.  Depending on sample different S(f).
    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
    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);
    
    Gp = fft(X4([100:1124])); dB_SG = 10*log10(conj(Gp).*Gp/1024);
    h1 = figure(4); set(h1,'Position',[100 200 750 200],'PaperPositionMode','auto'); hold off;
    plot([0:512]/1024,dB_SG(1:513),'k');
    hold on; 
    plot(f,dB_Sf,'r'); ylim([-30 50]);
    title('Spectrum from Sample');  xlabel('Frequency'); ylabel('sdf (dB)'); hold off;
    
    print -dpng S2L09_AR4_P.png 
    
    % Now apply a cosine taper (Hanning Window).
    Ch = [cos([0:512]*pi/512) cos([511:-1:0]*pi/512)];
    hth = (1-Ch')/2; %Maghth = sum(hth.*hth); hth = hth/sqrt(Maghth);
    X4h = X4([100:1124]).*hth;
    Gph = fft(X4h); dB_SGh = 10*log10(conj(Gph).*Gph/1024);
    h2 = figure(5); set(h2,'Position',[110 190 750 200],'PaperPositionMode','auto'); hold off;
    plot([0:512]/1024,dB_SGh(1:513),'k');
    hold on; 
    plot(f,dB_Sf,'r'); ylim([-30 50]);
    title('Spectrum from Sample with Hanning Taper');  xlabel('Frequency'); ylabel('sdf (dB)'); hold off;
    
    print -dpng S2L09_AR4_H.png
    

    % Now apply a cosine taper (Hanning Window).
    [E,V] = dpss(1024,4);  % Bandwidth is 4/1024.
    htd = E(:,1);
    X4d = X4([101:1124]).*htd*sum(htd);
    Gpd = fft(X4d); dB_SGd = 10*log10(conj(Gpd).*Gpd/1024);
    h3 = figure(6); set(h3,'Position',[120 180 750 200],'PaperPositionMode','auto'); hold off;
    plot([0:512]/1024,dB_SGd(1:513),'k');
    hold on; 
    plot(f,dB_Sf,'r'); ylim([-30 50]);
    title('Spectrum from Sample with DPSS NW=4');  xlabel('Frequency'); ylabel('sdf (dB)'); hold off;
    
    print -dpng S2L09_AR4_DPSSNW4.png
end

p3=1;
if( p3 == 1 )
    % Chi squared plots
    x = [0.01:0.05:10]; r =2 ; % r is degrees of freedom
    chi2p = 1-gammainc(x/2,r/2);
    m05 = find(chi2p < 0.05); x05 = x(m05(1)); dB_x05=10*log10(x05);
    m95 = find(chi2p < 0.95); x95 = x(m95(1)); dB_x95=10*log10(x95);
    
    chi2d = exp(-x/2).*x.^(r/2-1)./(2^r/2*gamma(r/2));
    lx = [-20:0.1:20];
    for k = 1:length(lx)-1
        x1 = 10^(lx(k)/10); x2 = 10^(lx(k+1)/10);
        chi1 = 1-gammainc(x1/2,r/2);chi2 = 1-gammainc(x2/2,r/2);
        chi2log(k) = (chi1-chi2)/(lx(k+1)-lx(k));
        chi2arg(k) = (lx(k)+lx(k+1))/2;
    end
    
    h1 = figure(7); set(h1,'Position',[120 180 750 500],'PaperPositionMode','auto'); hold off;
    subplot(2,1,1);
    plot(x,chi2d,'k'); title('\chi^2 probabolity density 2-degrees of freedom '); xlabel('u'); ylabel('PDF');  
    hold on, plot([x05 x05 NaN x95 x95],[0 0.4 NaN 0 0.5],'r'); 
    plot([2 2],[0 0.4],'b');hold off;
    subplot(2,1,2);
    plot(chi2arg,chi2log,'k'); xlabel('10log10(u) dB'); ylabel('PDF');
    hold on, plot([dB_x05 dB_x05 NaN dB_x95 dB_x95],[0 0.09 NaN 0 0.09],'r'); 
    plot([10*log10(2/exp(0.577215665)) 10*log10(2/exp(0.577215665))],[0 0.09],'b');hold off;
    
    print -dpng S2L09_Chi2.png
    
end