% M-file for Lecture 5.

p1 = 0;
if( p1 == 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.
    clear all;
    f= [0:0.5/1024:0.5]; t = [0:31];
    [E,V] = dpss(32,4,32);  % Bandwidth is 4/32.
    nf = length(f); Gp0 = zeros(1,nf);
    for k = 1:nf
        Gp0(k) = sum(E(:,1).*exp(-i*2*pi*f(k)*t)');
    end    
    dB_Gp0 = 10*log10(conj(Gp0).*Gp0); 
   
    h1 = figure(1); set(h1,'Position',[100 200 750 500],'PaperPositionMode','auto');
    subplot(2,1,1);
    plot(t,E(:,1),'k.'); axis tight; ylim([-0.4 0.4]); 
    title('DPSS for k=0, N=32, W=4/N ');

    subplot(2,1,2);
    plot(f,dB_Gp0,'k'); axis tight; ylim([-140 20]); 
    hold on; plot([1/8 1/8],[-130 10],'r'); hold off
    
    title('Spectrum of DPSS for k=0, N=32, W=4/N ');
    print -dpng S2L05_dpss_32_0.png 
    % Next eigen value
    for k = 1:nf
        Gp1(k) = sum(E(:,2).*exp(-i*2*pi*f(k)*t)');
    end    
    dB_Gp1 = 10*log10(conj(Gp1).*Gp1); 
   
    h2 = figure(2); set(h2,'Position',[120 180 750 500],'PaperPositionMode','auto');
    subplot(2,1,1);
    plot(t,E(:,2),'k.'); axis tight; ylim([-0.4 0.4]); 
    title('DPSS for k=1, N=32, W=4/N ');

    subplot(2,1,2);
    plot(f,dB_Gp1,'k'); axis tight; ylim([-140 20]); 
    hold on; plot([1/8 1/8],[-130 10],'r'); hold off
    
    title('Spectrum of DPSS for k=1, N=32, W=4/N ');
    print -dpng S2L05_dpss_32_1.png 
    
    % Next eigen value 2
    for k = 1:nf
        Gp2(k) = sum(E(:,3).*exp(-i*2*pi*f(k)*t)');
    end    
    dB_Gp2 = 10*log10(conj(Gp2).*Gp2); 
   
    h3 = figure(3); set(h3,'Position',[140 160 750 500],'PaperPositionMode','auto');
    subplot(2,1,1);
    plot(t,E(:,3),'k.'); axis tight; ylim([-0.4 0.4]); 
    title('DPSS for k=2, N=32, W=4/N ');

    subplot(2,1,2);
    plot(f,dB_Gp2,'k'); axis tight; ylim([-140 20]); 
    hold on; plot([1/8 1/8],[-130 10],'r'); hold off
    
    title('Spectrum of DPSS for k=2, N=32, W=4/N ');
    print -dpng S2L05_dpss_32_2.png 
    % Next eigen value 3
    for k = 1:nf
        Gp3(k) = sum(E(:,4).*exp(-i*2*pi*f(k)*t)');
    end    
    dB_Gp3 = 10*log10(conj(Gp3).*Gp3); 
   
    h4 = figure(4); set(h4,'Position',[160 140 750 500],'PaperPositionMode','auto');
    subplot(2,1,1);
    plot(t,E(:,4),'k.'); axis tight; ylim([-0.4 0.4]); 
    title('DPSS for k=3, N=32, W=4/N ');

    subplot(2,1,2);
    plot(f,dB_Gp3,'k'); axis tight; ylim([-140 20]); 
    hold on; plot([1/8 1/8],[-130 10],'r'); hold off
    
    title('Spectrum of DPSS for k=3, N=32, W=4/N ');
    print -dpng S2L05_dpss_32_3.png 
    
    
end

p2 = 0;
if( p2 == 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.
    clear all;
    f= [0:0.5/1024:0.5]; t = [0:98];
    [E,V] = dpss(99,4,99);  % Bandwidth is 4/32.
    nf = length(f); Gp0 = zeros(1,nf);
    for k = 1:nf
        Gp0(k) = sum(E(:,1).*exp(-i*2*pi*f(k)*t)');
    end    
    dB_Gp0 = 10*log10(conj(Gp0).*Gp0); 
   
    h1 = figure(5); set(h1,'Position',[100 200 750 500],'PaperPositionMode','auto');
    subplot(2,1,1);
    plot(t,E(:,1),'k.'); axis tight; ylim([-0.4 0.4]); 
    title('DPSS for k=0, N=99, W=4/N ');

    subplot(2,1,2);
    plot(f,dB_Gp0,'k'); axis tight; ylim([-140 20]); 
    hold on; plot([4/99 4/99],[-130 10],'r'); hold off
    
    title('Spectrum of DPSS for k=0, N=99, W=4/N ');
    print -dpng S2L05_dpss_99_0.png 
    % Next eigen value
    for k = 1:nf
        Gp1(k) = sum(E(:,2).*exp(-i*2*pi*f(k)*t)');
    end    
    dB_Gp1 = 10*log10(conj(Gp1).*Gp1); 
   
    h2 = figure(6); set(h2,'Position',[120 180 750 500],'PaperPositionMode','auto');
    subplot(2,1,1);
    plot(t,E(:,2),'k.'); axis tight; ylim([-0.4 0.4]); 
    title('DPSS for k=1, N=99, W=4/N ');

    subplot(2,1,2);
    plot(f,dB_Gp1,'k'); axis tight; ylim([-140 20]); 
    hold on; plot([4/99 4/99],[-130 10],'r'); hold off
    
    title('Spectrum of DPSS for k=1, N=99, W=4/N ');
    print -dpng S2L05_dpss_99_1.png 
    
    % Next eigen value 2
    for k = 1:nf
        Gp2(k) = sum(E(:,3).*exp(-i*2*pi*f(k)*t)');
    end    
    dB_Gp2 = 10*log10(conj(Gp2).*Gp2); 
   
    h3 = figure(7); set(h3,'Position',[140 160 750 500],'PaperPositionMode','auto');
    subplot(2,1,1);
    plot(t,E(:,3),'k.'); axis tight; ylim([-0.4 0.4]); 
    title('DPSS for k=2, N=99, W=4/N ');

    subplot(2,1,2);
    plot(f,dB_Gp2,'k'); axis tight; ylim([-140 20]); 
    hold on; plot([4/99 4/99],[-130 10],'r'); hold off
    
    title('Spectrum of DPSS for k=2, N=99, W=4/N ');
    print -dpng S2L05_dpss_99_2.png 
    % Next eigen value 3
    for k = 1:nf
        Gp3(k) = sum(E(:,4).*exp(-i*2*pi*f(k)*t)');
    end    
    dB_Gp3 = 10*log10(conj(Gp3).*Gp3); 
   
    h4 = figure(8); set(h4,'Position',[160 140 750 500],'PaperPositionMode','auto');
    subplot(2,1,1);
    plot(t,E(:,4),'k.'); axis tight; ylim([-0.4 0.4]); 
    title('DPSS for k=3, N=99, W=4/N ');

    subplot(2,1,2);
    plot(f,dB_Gp3,'k'); axis tight; ylim([-140 20]); 
    hold on; plot([4/99 4/99],[-130 10],'r'); hold off
    
    title('Spectrum of DPSS for k=3, N=99, W=4/N ');
    print -dpng S2L05_dpss_99_3.png 
    
    
end

p3 = 0;
if( p3 == 1 ) 
    % Concentration problem: Plot the eigen values
    clear all;
    t0 = [0:31]; t1 = [0:98]; 
    [E,V0] = dpss(32,4,32);  % Bandwidth is 4/99.
    [E,V1] = dpss(99,4,99);  % Bandwidth is 4/99.
       
    h1 = figure(9); set(h1,'Position',[100 200 750 500],'PaperPositionMode','auto');
    subplot(2,1,1);
    plot(t0,V0,'k.'); axis tight; ylim([-0.1 1.1]); 
    title('Eigenvalues for DPSS N=32, W=4/N ');
    hold on; plot([8 8],[0 1],'r'); hold off

    subplot(2,1,2);
    plot(t1,V1,'k.'); axis tight; ylim([-0.1 1.1]); 
    hold on; plot([8 8],[0 1],'r'); hold off
    
    title('Eigenvalues for DPSS N=99, W=4/N ');
    print -dpng S2L05_dpss_eig.png 
end

p4 = 1;
if( p4 == 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.
    clear all;
    f= [0:0.5/1024:0.5]; t = [0:127];
    [E,V] = dpss(128,4,128);  % Bandwidth is 4/32.
    nf = length(f); Gp0 = zeros(1,nf);
    for k = 1:nf
        Gp0(k) = sum(E(:,1).*exp(-i*2*pi*f(k)*t)');
    end    
    dB_Gp0 = 10*log10(conj(Gp0).*Gp0); 
    % Now do an fft at nominal length
    Gp0_fftn = fft(E(:,1)); fn = [0:0.5/64:0.5];
    dB_Gp0_fftn = 10*log10(conj(Gp0_fftn).*Gp0_fftn); 
    % FFT at twice length (padded with zeros)
    Gp0_fftl = fft(E(:,1),256);fl = [0:0.5/128:0.5];
    dB_Gp0_fftl = 10*log10(conj(Gp0_fftl).*Gp0_fftl); 
   
    h1 = figure(10); set(h1,'Position',[100 200 750 500],'PaperPositionMode','auto');
    subplot(2,1,1); hold off; 
    plot(f,dB_Gp0,'r'); axis tight; ylim([-140 20]); 
    title('Spectrum DPSS for k=0, N=128 FFT 128, W=4/N '); hold on;
    plot(fn,dB_Gp0_fftn([1:65]),'k-+');

    subplot(2,1,2); hold off;
    plot(f,dB_Gp0,'r'); axis tight; ylim([-140 20]); 
    hold on; plot(fl,dB_Gp0_fftl([1:129]),'b-+'); 
    
    title('Spectrum DPSS for k=0, N=128 FFT 256, W=4/N ');
    print -dpng S2L05_dpss_128_fft_0.png 
    
    %%% Now repeat for k=8
    nf = length(f); Gp8 = zeros(1,nf);
    for k = 1:nf
        Gp8(k) = sum(E(:,9).*exp(-i*2*pi*f(k)*t)');
    end    
    dB_Gp8 = 10*log10(conj(Gp8).*Gp8); 
    % Now do an fft at nominal length
    Gp8_fftn = fft(E(:,9)); fn = [0:0.5/64:0.5];
    dB_Gp8_fftn = 10*log10(conj(Gp8_fftn).*Gp8_fftn); 
    % FFT at twice length (padded with zeros)
    Gp8_fftl = fft(E(:,9),256);fl = [0:0.5/128:0.5];
    dB_Gp8_fftl = 10*log10(conj(Gp8_fftl).*Gp8_fftl); 
   
    h2 = figure(11); set(h2,'Position',[120 180 750 500],'PaperPositionMode','auto');
    subplot(2,1,1); hold off; 
    plot(f,dB_Gp8,'r'); axis tight; ylim([-140 20]); 
    title('Spectrum DPSS for k=8, N=128 FFT 128, W=4/N '); hold on;
    plot(fn,dB_Gp8_fftn([1:65]),'k-+');

    subplot(2,1,2); hold off;
    plot(f,dB_Gp8,'r'); axis tight; ylim([-140 20]); 
    hold on; plot(fl,dB_Gp8_fftl([1:129]),'b-+'); 
    
    title('Spectrum DPSS for k=8, N=128 FFT 256, W=4/N ');
    print -dpng S2L05_dpss_128_fft_8.png 

end