% M-file for Lecture 10.
%%
p1 = 0;
if( p1 == 1 )
    %  Barrlett Window
    m = 30; wtm = zeros(64,1);
    wtm([1:31]) = 1 - [0:30]/m;
    % Wmf = fft(wtm,2048); db_Wmf = 10*log10(real(Wmf));
    f = [0:0.001:0.5];
    Wmf = sin(m*pi*f).^2./(m*(sin(pi*f).^2)); db_Wmf = 10*log10(real(Wmf));
    
    h1 = figure(1); set(h1,'Position',[100 300 750 600],'PaperPositionMode','auto');
    subplot(2,2,1);
    plot([0:63],wtm,'.');  axis tight; ylim([-0.4 1.2]); xlabel('\tau'); ylabel('W_{\tau,m}');
    subplot(2,2,2);
    %plot(0.5*[0:1023]/1024,db_Wmf([1:1024]),'k');
    plot(f, db_Wmf,'k'); axis tight, ylim([-80 20]); xlabel('frequency'); ylabel('W_m(f) (dB)');
    
    % Rectangular taper
    hr = ones(64,1)/sqrt(64);
    % Now use 244a:
    lf = length(f); Umr = zeros(1,lf); 
    for tau = -63:63
        sum = 0;
        for t = 1:64-abs(tau)
            sum = sum + hr(t)*hr(t+abs(tau));
        end
        in = abs(tau)+1; 
        Umr = Umr + wtm(in)*sum*exp(-i*2*pi*f*tau);
    end
    dB_Umr = 10*log10(real(Umr));
    subplot(2,2,3);
    plot(f, dB_Umr,'k'); axis tight, ylim([-80 20]); xlabel('frequency'); ylabel('U_m(f) Rectangular (dB)');
    
    % Discrete Prolate Spheriodal sequence NW = 4;
    [E,V] = dpss(64,4);  % Bandwidth is 4/64.
    hd = E(:,1);
    % Now use 244a:
    lf = length(f); Umd = zeros(1,lf); 
    for tau = -63:63
        sum = 0;
        for t = 1:64-abs(tau)
            sum = sum + hd(t)*hd(t+abs(tau));
        end
        in = abs(tau)+1; 
        Umd = Umd + wtm(in)*sum*exp(-i*2*pi*f*tau);
    end
    dB_Umd = 10*log10(real(Umd));
    subplot(2,2,4);
    plot(f, dB_Umd,'k'); axis tight, ylim([-80 20]); xlabel('frequency'); ylabel('U_m(f) DPSS NW=4 (dB)');
                
    print -dpng S2L10_Barlett.png
end
 
%%
p2 = 0;
if( p2 == 1 )
    %  Daniell window m=20, N 64
    m = 20; wtm = zeros(64,1); tw = [0:63]; 
    wtm = sinc(tw/m);
    
    % Wmf = fft(wtm,2048); db_Wmf = 10*log10(real(Wmf));
    f = [0:0.001:0.5]; lf = length(f); Wmf = zeros(1,lf);
    for tau = -63:63
        in = abs(tau)+1;
        Wmf = Wmf + sinc(tau/m)*exp(-i*2*pi*f*tau);
    end
    
    db_Wmf = 10*log10(abs(real(Wmf)));
    
    h1 = figure(2); set(h1,'Position',[100 300 750 600],'PaperPositionMode','auto');
    subplot(2,2,1);
    plot([0:63],wtm,'.');  axis tight; ylim([-0.4 1.2]); xlabel('\tau'); ylabel('W_{\tau,m}');
    subplot(2,2,2);
    %plot(0.5*[0:1023]/1024,db_Wmf([1:1024]),'k');
    plot(f, db_Wmf,'k'); axis tight, ylim([-80 20]); xlabel('frequency'); ylabel('W_m(f) (dB)');
    
    % Rectangular taper
    hr = ones(64,1)/sqrt(64);
    % Now use 244a:
    Umr = zeros(1,lf); 
    for tau = -63:63
        sum = 0;
        for t = 1:64-abs(tau)
            sum = sum + hr(t)*hr(t+abs(tau));
        end
        in = abs(tau)+1; 
        Umr = Umr + wtm(in)*sum*exp(-i*2*pi*f*tau);
    end
    dB_Umr = 10*log10(real(Umr));
    subplot(2,2,3);
    plot(f, dB_Umr,'k'); axis tight, ylim([-80 20]); xlabel('frequency'); ylabel('U_m(f) Rectangular (dB)');
    
    % Discrete Prolate Spheriodal sequence NW = 4;
    [E,V] = dpss(64,4);  % Bandwidth is 4/64.
    hd = E(:,1);
    % Now use 244a:
    lf = length(f); Umd = zeros(1,lf); 
    for tau = -63:63
        sum = 0;
        for t = 1:64-abs(tau)
            sum = sum + hd(t)*hd(t+abs(tau));
        end
        in = abs(tau)+1; 
        Umd = Umd + wtm(in)*sum*exp(-i*2*pi*f*tau);
    end
    dB_Umd = 10*log10(real(Umd));
    subplot(2,2,4);
    plot(f, dB_Umd,'k'); axis tight, ylim([-80 20]); xlabel('frequency'); ylabel('U_m(f) DPSS NW=4 (dB)');
                
    print -dpng S2L10_Daniell.png
end
%%
p3 = 0;
if( p3 == 1 )
    %  Parzen window m=37, N 64
    m = 37; wtm = zeros(64,1); tw = [0:63]; 
    wtm([1:floor(m/2)+1]) = 1-6*([0:floor(m/2)]/m).^2+6*([0:floor(m/2)]/m).^3;
    wtm([floor(m/2)+2:m]) = 2*(1-[floor(m/2)+1:m-1]/m).^3;
    
    % Wmf = fft(wtm,2048); db_Wmf = 10*log10(real(Wmf));
    f = [0:0.001:0.5]; lf = length(f); Wmf = zeros(1,lf);
    %Wmf = 4*(3-2*(sin(pi*f).^2)).*(sin(m*pi*f/2).^4)./(m^3*(sin(pi*f).^4));
    for t = 0:63
        Wmf = Wmf + wtm(t+1)*exp(-i*2*pi*f*t) + wtm(t+1)*exp(i*2*pi*f*t);
    end
        
    
    db_Wmf = 10*log10(abs(real(Wmf-1)));
    
    h1 = figure(3); set(h1,'Position',[100 300 750 600],'PaperPositionMode','auto');
    subplot(2,2,1);
    plot([0:63],wtm,'.');  axis tight; ylim([-0.4 1.2]); xlabel('\tau'); ylabel('W_{\tau,m}');
    subplot(2,2,2);
    %plot(0.5*[0:1023]/1024,db_Wmf([1:1024]),'k');
    plot(f, db_Wmf,'k'); axis tight, ylim([-80 20]); xlabel('frequency'); ylabel('W_m(f) (dB)');
    
    % Rectangular taper
    hr = ones(64,1)/sqrt(64);
    % Now use 244a:
    Umr = zeros(1,lf); 
    for tau = -63:63
        sum = 0;
        for t = 1:64-abs(tau)
            sum = sum + hr(t)*hr(t+abs(tau));
        end
        in = abs(tau)+1; 
        Umr = Umr + wtm(in)*sum*exp(-i*2*pi*f*tau);
    end
    dB_Umr = 10*log10(real(Umr));
    subplot(2,2,3);
    plot(f, dB_Umr,'k'); axis tight, ylim([-80 20]); xlabel('frequency'); ylabel('U_m(f) Rectangular (dB)');
    
    % Discrete Prolate Spheriodal sequence NW = 4;
    [E,V] = dpss(64,4);  % Bandwidth is 4/64.
    hd = E(:,1);
    % Now use 244a:
    lf = length(f); Umd = zeros(1,lf); 
    for tau = -63:63
        sum = 0;
        for t = 1:64-abs(tau)
            sum = sum + hd(t)*hd(t+abs(tau));
        end
        in = abs(tau)+1; 
        Umd = Umd + wtm(in)*sum*exp(-i*2*pi*f*tau);
    end
    dB_Umd = 10*log10(real(Umd));
    subplot(2,2,4);
    plot(f, dB_Umd,'k'); axis tight, ylim([-80 20]); xlabel('frequency'); ylabel('U_m(f) DPSS NW=4 (dB)');
                
    print -dpng S2L10_Parzen.png
end
%%
p4 = 0;
if( p4 == 1 )
    %  Papoulis Window m = 34 N=64 
    m = 34; wtm = zeros(64,1); tw = [0:m-1]; 
    wtm([1:m])=1/pi*abs(sin(pi*tw/m))+(1-tw/m).*cos(pi*tw/m);
    
    % Wmf = fft(wtm,2048); db_Wmf = 10*log10(real(Wmf));
    f = [0:0.001:0.5]; lf = length(f); Wmf = zeros(1,lf);
    %Wmf = 4*(3-2*(sin(pi*f).^2)).*(sin(m*pi*f/2).^4)./(m^3*(sin(pi*f).^4));
    for t = 0:63
        Wmf = Wmf + wtm(t+1)*exp(-i*2*pi*f*t) + wtm(t+1)*exp(i*2*pi*f*t);
    end
        
    
    db_Wmf = 10*log10(abs(real(Wmf-1)));
    
    h1 = figure(4); set(h1,'Position',[100 300 750 600],'PaperPositionMode','auto');
    subplot(2,2,1);
    plot([0:63],wtm,'.');  axis tight; ylim([-0.4 1.2]); xlabel('\tau'); ylabel('W_{\tau,m}');
    subplot(2,2,2);
    %plot(0.5*[0:1023]/1024,db_Wmf([1:1024]),'k');
    plot(f, db_Wmf,'k'); axis tight, ylim([-80 20]); xlabel('frequency'); ylabel('W_m(f) (dB)');
    
    % Rectangular taper
    hr = ones(64,1)/sqrt(64);
    % Now use 244a:
    Umr = zeros(1,lf); 
    for tau = -63:63
        sum = 0;
        for t = 1:64-abs(tau)
            sum = sum + hr(t)*hr(t+abs(tau));
        end
        in = abs(tau)+1; 
        Umr = Umr + wtm(in)*sum*exp(-i*2*pi*f*tau);
    end
    dB_Umr = 10*log10(real(Umr));
    subplot(2,2,3);
    plot(f, dB_Umr,'k'); axis tight, ylim([-80 20]); xlabel('frequency'); ylabel('U_m(f) Rectangular (dB)');
    
    % Discrete Prolate Spheriodal sequence NW = 4;
    [E,V] = dpss(64,4);  % Bandwidth is 4/64.
    hd = E(:,1);
    % Now use 244a:
    lf = length(f); Umd = zeros(1,lf); 
    for tau = -63:63
        sum = 0;
        for t = 1:64-abs(tau)
            sum = sum + hd(t)*hd(t+abs(tau));
        end
        in = abs(tau)+1; 
        Umd = Umd + wtm(in)*sum*exp(-i*2*pi*f*tau);
    end
    dB_Umd = 10*log10(real(Umd));
    subplot(2,2,4);
    plot(f, dB_Umd,'k'); axis tight, ylim([-80 20]); xlabel('frequency'); ylabel('U_m(f) DPSS NW=4 (dB)');
                
    print -dpng S2L10_Papoulis.png
end

%%
p5 = 0;
if( p5 == 1 )
    % Generate the effective degrees of freedom for the wosa estminator
    N= 1024 ; Ns = 64; 
    ht = sqrt(2/3/(Ns+1))*(1-cos(2*pi*[1:Ns]/(Ns+1)));
    nuf64 = zeros(Ns,1); ov64=zeros(Ns,1);
    for n = 1:Ns   % Loop over shift factor
        % Compute number of blocks
        Nb = (N-Ns)/n+1; ov64(n) = (1-n/Ns)*100;
        s1 = 0;
        for m=1:Nb-1
            s2 = 0;
            for t = 1:Ns
                if( t+m*n <= Ns ) s2 = s2+ht(t)*ht(t+m*n); end
            end
            s1 = s1 + (1-m/Nb)*s2^2;
        end
        nuf64(n) = 2*Nb/(1+2*s1);
    end
    
    % Now repeat of 128
    N= 1024 ; Ns = 128; 
    ht = sqrt(2/3/(Ns+1))*(1-cos(2*pi*[1:Ns]/(Ns+1)));
    nuf128 = zeros(Ns,1); ov128=zeros(Ns,1);
    for n = 1:Ns   % Loop over shift factor
        % Compute number of blocks
        Nb = (N-Ns)/n+1; ov128(n) = (1-n/Ns)*100;
        s1 = 0;
        for m=1:Nb-1
            s2 = 0;
            for t = 1:Ns
                if( t+m*n <= Ns ) s2 = s2+ht(t)*ht(t+m*n); end
            end
            s1 = s1 + (1-m/Nb)*s2^2;
        end
        nuf128(n) = 2*Nb/(1+2*s1);
    end

    % Now repeat of 256
    N= 1024 ; Ns = 256; 
    ht = sqrt(2/3/(Ns+1))*(1-cos(2*pi*[1:Ns]/(Ns+1)));
    nuf256 = zeros(Ns,1); ov256=zeros(Ns,1);
    for n = 1:Ns   % Loop over shift factor
        % Compute number of blocks
        Nb = (N-Ns)/n+1; ov256(n) = (1-n/Ns)*100;
        s1 = 0;
        for m=1:Nb-1
            s2 = 0;
            for t = 1:Ns
                if( t+m*n <= Ns ) s2 = s2+ht(t)*ht(t+m*n); end
            end
            s1 = s1 + (1-m/Nb)*s2^2;
        end
        nuf256(n) = 2*Nb/(1+2*s1);
    end
    
    h1 = figure(5); set(h1,'Position',[100 300 750 500],'PaperPositionMode','auto');
    plot(ov64,nuf64,'r');hold on; 
    plot(ov128,nuf128,'b');
    plot(ov256,nuf256,'k');
    legend('N_s = 64','N_s = 128 ', 'N_s = 256','Location','NorthWest');
    title('Effective degrees of freedom for COSA for sample size 1024 '); xlabel('Percent Overlap');
    ylabel('Effective degrees of freedom');
    
    print -dpng S2L10_cosa_dof.png
    
    
end

%%
p6 = 0;
if( p6 == 1 )
    % Generate the effective degrees of freedom for the wosa estminator
    N= 1024 ; Ns = 64; 
    [E,V] = dpss(Ns,4); ht = E(:,1);
    nuf64 = zeros(Ns,1); ov64=zeros(Ns,1);
    for n = 1:Ns   % Loop over shift factor
        % Compute number of blocks
        Nb = (N-Ns)/n+1; ov64(n) = (1-n/Ns)*100;
        s1 = 0;
        for m=1:Nb-1
            s2 = 0;
            for t = 1:Ns
                if( t+m*n <= Ns ) s2 = s2+ht(t)*ht(t+m*n); end
            end
            s1 = s1 + (1-m/Nb)*s2^2;
        end
        nuf64(n) = 2*Nb/(1+2*s1);
    end
    
    % Now repeat of 128
    N= 1024 ; Ns = 128; 
    [E,V] = dpss(Ns,4); ht = E(:,1);
    nuf128 = zeros(Ns,1); ov128=zeros(Ns,1);
    for n = 1:Ns   % Loop over shift factor
        % Compute number of blocks
        Nb = (N-Ns)/n+1; ov128(n) = (1-n/Ns)*100;
        s1 = 0;
        for m=1:Nb-1
            s2 = 0;
            for t = 1:Ns
                if( t+m*n <= Ns ) s2 = s2+ht(t)*ht(t+m*n); end
            end
            s1 = s1 + (1-m/Nb)*s2^2;
        end
        nuf128(n) = 2*Nb/(1+2*s1);
    end

    % Now repeat of 256
    N= 1024 ; Ns = 256; 
    [E,V] = dpss(Ns,4); ht = E(:,1);
    nuf256 = zeros(Ns,1); ov256=zeros(Ns,1);
    for n = 1:Ns   % Loop over shift factor
        % Compute number of blocks
        Nb = (N-Ns)/n+1; ov256(n) = (1-n/Ns)*100;
        s1 = 0;
        for m=1:Nb-1
            s2 = 0;
            for t = 1:Ns
                if( t+m*n <= Ns ) s2 = s2+ht(t)*ht(t+m*n); end
            end
            s1 = s1 + (1-m/Nb)*s2^2;
        end
        nuf256(n) = 2*Nb/(1+2*s1);
    end
    
    h1 = figure(6); set(h1,'Position',[100 300 750 500],'PaperPositionMode','auto');
    plot(ov64,nuf64,'r');hold on; 
    plot(ov128,nuf128,'b');
    plot(ov256,nuf256,'k');
    legend('N_s = 64','N_s = 128 ', 'N_s = 256','Location','NorthWest');
    title('Effective degrees of freedom for COSA for sample size 1024 dpss NW4'); xlabel('Percent Overlap');
    ylabel('Effective degrees of freedom');
    
    print -dpng S2L10_cosa_dpss.png
    
    
end

%%
p7 = 0;
if( p7 == 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(1132,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
    X4n = X4([105:1128]);
    % AR(4) Process(t) =
    % A$4 = 2.7607X(t-1)-3.8106X(t-2)+2.65235X(t-3)-0.9238*X(t-4)+e
    f=[0:0.5/1024:0.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);
    SfAR4 = 1./(conj(denom).*denom); dB_SfAR4 = 10*log10(SfAR4);
    
    [SfWelch,w] = pwelch(X4n,[],[],[],1); % By default generates 8 segements with 50% overlap, Hanning window
    dB_SfWelch=10*log10(SfWelch/2);
    h1 = figure(6); set(h1,'Position',[100 300 750 500],'PaperPositionMode','auto');
    
    plot(f,dB_SfAR4,'r',w,dB_SfWelch,'k'); title('AR(4) generated with pwelch: 8 segements with 50% overlap, Hanning Taper');
    xlabel('frequency'); ylabel('dB');
    
    print -dpng S2L10_pwelchAR4.png
end
    
%% Now look at Multitaper methods
% Here we use DPSS functions as tapers which are orthognal to each other.
% First we compute with each of the orthogonal tapers, saving each one, and
% then we average them togther).
p8 = 1;
if p8 == 1 
    % Generate an AR(4) time series
    % Generate AR(4) X(t) =
    % 2.7607X(t-1)-3.8106X(t-2)+2.65235X(t-3)-0.9238*X(t-4)+e
    N = 1024 ; % Number of samples to be used
    randn('seed',120);
    % Generate longer series so that we get this one started
    X4 = zeros(N+64,1);
    for k = 1:N+64
        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
    X4 = X4(65:N+64);
    
    % OK: Genetrate time series convolved with DPSS NW 4 1
    % Concentration problem: E are the prolate prolate spheroidal sequences
    % (Slepian sequences) function, E(:,1) is k=0, E(:,2) is k=1 etc, V are
    % the eigenvalues.  Call below generates all 32 eigenfunctions.
    [E,V] = dpss(N,4,8);  % Bandwidth is 4/1024.
    
    % So of the first set of DPSS functions
    f=[0:N/2-1]/N;
    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);
    psdX4sv = zeros(N,m);

    for m = 1:8
        X4m = X4.*E(:,m); % DPSS window applied to time series
        % Now get the PSD
        ftX4m = fft(X4m); psdX4sv(:,m) = conj(ftX4m).*ftX4m;
        dB_ftX4m = 10*log10(conj(ftX4m).*ftX4m);
        ftEm  = fft(E(:,m),10240);
        fE =[0:5119]/(2*5120);
        dB_ftEm = 10*log10(conj(ftEm).*ftEm);
        
        h1 = figure(6+m); set(h1,'Position',[100 200 750 500],'PaperPositionMode','auto'); hold off; 
        subplot(2,2,1);
        ptitle = sprintf('DPSS k = %d Sample %d',m-1,N);
        plot(1:N,E(:,m),'k'); xlabel('Time'); ylabel('DPSS function'); title(ptitle); axis tight; 
        subplot(2,2,2)
        plot(1:N,X4m,'k'); xlabel('Time'); ylabel('AR(4) * DPSS function'); axis tight;
        subplot(2,2,3);
        plot(fE(1:200),dB_ftEm(1:200),'k',[4/N 4/N],[-70 30],'r'); axis tight; ylim([-80 40]); xlabel('Frequency'); ylabel('Spectral Window'); 
        subplot(2,2,4);
        plot(f,dB_ftX4m(1:N/2),'k',f,dB_Sf,'r');xlabel('Frequency'); ylabel('PSD'); axis tight;
        
        pfn = sprintf('S2L10_Fig%2.2d_MTComp.png',m+6);
        print('-dpng',pfn);
    end

    % Now averge power spectra
    h1 = figure(15); set(h1,'Position',[100 200 750 500],'PaperPositionMode','auto'); hold off;
    for m = 1:4
        if m == 1 
            psav = psdX4sv(:,1); 
        else
            psav = mean(psdX4sv(:,[1:m])'); 
        end
        dB_psav = 10*log10(psav);
 
        subplot(4,1,m);
        ylab = sprintf('Average PSD to k=%d ',m-1); 
        plot(f,dB_psav(1:N/2),'k',f,dB_Sf,'r');xlabel('Frequency'); ylabel(ylab); axis tight;
    end
    pfn = 'S2L10_Fig15_MTAver.png';
    print('-dpng',pfn);
    h1 = figure(16); set(h1,'Position',[100 200 750 500],'PaperPositionMode','auto'); hold off;
    for m = 5:8
        psav = mean(psdX4sv(:,[1:m])'); dB_psav = 10*log10(psav);
 
        subplot(4,1,m-4);
        ylab = sprintf('Average PSD to k=%d ',m-1); 
        plot(f,dB_psav(1:N/2),'k',f,dB_Sf,'r');xlabel('Frequency'); ylabel(ylab); axis tight;
    end
    pfn = 'S2L10_Fig16_MTAver.png';
    print('-dpng',pfn);

end