% M-file for Lecture 4.

p1 = 1;
if( p1 == 1 ) 
    % Convolution problem.
    t = [-10:0.01:10]; s=[-3:0.01:3];
    f1=1/6; f2 = 3; 
    gt = 5*cos(2*pi*f1*t+0.5)+cos(2*pi*f2*t+1.1);
    
    % Set sigma and Gaussian Filter
    sig = 0.1; hs=5*exp(-s.*s/(2*sig^2)); % hs=1/(sqrt(2*pi)*sig)*exp(-s.*s/(2*sig^2));
    h1 = figure(1); set(h1,'Position',[100 300 750 180],'PaperPositionMode','auto');
    % Plot functiion and Gassian
    hold off ; plot(t,gt,'k',s,hs,'r');
    % Generate the result by convolution.  Normalize by H(t) area
    gch = conv(gt,hs)/sum(hs); hold on
    % Matlab generates a series equal to sum of two parts so get start
    % index an plot
    is = (length(s)-1)/2+1; ie = is + length(gt)-1;
    plot(t,gch(is:is+length(gt)-1),'b');
    % Compute using theory. s1 and s2 are losses at frequencies.
    s1 = exp(-(sig*2*pi*f1)^2/2); s2 = exp(-(sig*2*pi*f2)^2/2);
    gf = s1*5*cos(2*pi*f1*t+0.5)+s2*cos(2*pi*f2*t+1.1);
    plot(t,gf,'g');  axis tight; xlim([-4 4]);
    title('Smoothed \sigma=0.1');
    print -dpng S2L04_g01.png
    
    % Set sigma and Gaussian Filter
    sig = 0.25; hs=5*exp(-s.*s/(2*sig^2)); % hs=1/(sqrt(2*pi)*sig)*exp(-s.*s/(2*sig^2));
    h2 = figure(2); set(h2,'Position',[120 280 750 180],'PaperPositionMode','auto');
    % Plot functiion and Gassian
    hold off ; plot(t,gt,'k',s,hs,'r');
    % Generate the result by convolution.  Normalize by H(t) area
    gch = conv(gt,hs)/sum(hs); hold on
    % Matlab generates a series equal to sum of two parts so get start
    % index an plot
    is = (length(s)-1)/2+1; ie = is + length(gt)-1;
    plot(t,gch(is:is+length(gt)-1),'b');
    % Compute using theory. s1 and s2 are losses at frequencies.
    s1 = exp(-(sig*2*pi*f1)^2/2); s2 = exp(-(sig*2*pi*f2)^2/2);
    gf = s1*5*cos(2*pi*f1*t+0.5)+s2*cos(2*pi*f2*t+1.1);
    plot(t,gf,'g');  axis tight; xlim([-4 4]);
    title('Smoothed \sigma=0.25');
    print -dpng S2L04_g02.png
    
    % Set sigma and Gaussian Filter
    sig = 0.625; hs=5*exp(-s.*s/(2*sig^2)); % hs=1/(sqrt(2*pi)*sig)*exp(-s.*s/(2*sig^2));
    h3 = figure(3); set(h3,'Position',[140 260 750 180],'PaperPositionMode','auto');
    % Plot functiion and Gassian
    hold off ; plot(t,gt,'k',s,hs,'r');
    % Generate the result by convolution.  Normalize by H(t) area
    gch = conv(gt,hs)/sum(hs); hold on
    % Matlab generates a series equal to sum of two parts so get start
    % index an plot
    is = (length(s)-1)/2+1; ie = is + length(gt)-1;
    plot(t,gch(is:is+length(gt)-1),'b');
    % Compute using theory. s1 and s2 are losses at frequencies.
    s1 = exp(-(sig*2*pi*f1)^2/2); s2 = exp(-(sig*2*pi*f2)^2/2);
    gf = s1*5*cos(2*pi*f1*t+0.5)+s2*cos(2*pi*f2*t+1.1);
    plot(t,gf,'g');  axis tight; xlim([-4 4]);
    title('Smoothed \sigma=0.625');
    print -dpng S2L04_g03.png
end


p2 = 1;
if( p2 == 1 ) 
    % Convolution problem.  This time we will convolve with a box car
    % +-delta in width
    t = [-10:0.01:10]; 
    f1=1/6; f2 = 3; 
    gt = 5*cos(2*pi*f1*t+0.5)+cos(2*pi*f2*t+1.1);
    
    % Set sigma and box car filter width -delta to delta
    delta = 0.17; s=[-delta-0.01:0.01:delta+0.01]; hs=ones(1,length(s))/(2*delta); hs(1) = 0; hs(length(s))=0;
    h1 = figure(4); set(h1,'Position',[100 300 750 200],'PaperPositionMode','auto');
    % Plot functiion and Gassian
    hold off ; plot(t,gt,'k',s,hs,'r');
    % Generate the result by convolution.  Normalize by H(t) area
    gch = conv(gt,hs)/sum(hs); hold on
    % Matlab generates a series equal to sum of two parts so get start
    % index an plot
    is = (length(s)-1)/2+1; ie = is + length(gt)-1;
    plot(t,gch(is:is+length(gt)-1),'b');
    % Compute using theory. s1 and s2 are losses at frequencies.
    s1 = sinc(2*f1*delta); s2 = sinc(2*f2*delta);
    gf = s1*5*cos(2*pi*f1*t+0.5)+s2*cos(2*pi*f2*t+1.1);
    plot(t,gf,'g');  axis tight; xlim([-4 4]);
    title('Smoothed Box Car \delta =0.17');
    print -dpng S2L04_b01.png
    
    % Set sigma and box car filter width -delta to delta
    delta = 0.25; s=[-delta-0.01:0.01:delta+0.01]; hs=ones(1,length(s))/(2*delta); hs(1) = 0; hs(length(s))=0;
    h2 = figure(5); set(h2,'Position',[120 280 750 200],'PaperPositionMode','auto');
    % Plot functiion and Gassian
    hold off ; plot(t,gt,'k',s,hs,'r');
    % Generate the result by convolution.  Normalize by H(t) area
    gch = conv(gt,hs)/sum(hs); hold on
    % Matlab generates a series equal to sum of two parts so get start
    % index an plot
    is = (length(s)-1)/2+1; ie = is + length(gt)-1;
    plot(t,gch(is:is+length(gt)-1),'b');
    % Compute using theory. s1 and s2 are losses at frequencies.
    s1 = sinc(2*f1*delta); s2 = sinc(2*f2*delta);
    gf = s1*5*cos(2*pi*f1*t+0.5)+s2*cos(2*pi*f2*t+1.1);
    plot(t,gf,'g');  axis tight; xlim([-4 4]);
    title('Smoothed Box Car \delta =0.25');
    print -dpng S2L04_b02.png
    
end

p3 = 1;
if( p3 == 1 ) 
    % Plots of Dirichlet's kernel.
    f = [-0.5:0.002:0.5]; 
    % Generate each of the terms in the sum
    m = 4; D4f = zeros(1,length(f));
    for t = -m:m
        D4f = D4f + exp(i*2*pi*f*t)/(2*m+1);
    end
    % The following is an analytic formula for the Dirichlet's Kernel.
    % (f=0 poses a problem).
    D4a = exp(-i*2*pi*m*f).*(exp(i*2*pi*(2*m+1)*f)-1)./(exp(i*2*pi*f)-1)/(2*m+1);
    h1 = figure(6); set(h1,'Position',[100 300 750 200],'PaperPositionMode','auto');
    % Plot functiion and Gassian
    plot(f,real(D4f),'k'); axis tight; 
    title('Dirichlets Kernal m= 4'); xlabel('Frequency \DeltaT=1'); ylabel('D_{2m+1}');
    print -dpng S2L04_d04.png
    
    m = 16; D16f = zeros(1,length(f));
    for t = -m:m
        D16f = D16f + exp(i*2*pi*f*t)/(2*m+1);
    end
    h1 = figure(7); set(h1,'Position',[120 280 750 200],'PaperPositionMode','auto');
    % Plot functiion and Gassian
    plot(f,real(D16f),'k'); axis tight; xlabel('Frequency \DeltaT=1'); ylabel('D_{2m+1}');
    title('Dirichlets Kernal m=16');
    print -dpng S2L04_d16.png

    m = 64; D64f = zeros(1,length(f));
    for t = -m:m
        D64f = D64f + exp(i*2*pi*f*t)/(2*m+1);
    end
    h1 = figure(8); set(h1,'Position',[140 260 750 200],'PaperPositionMode','auto');
    % Plot functiion and Gassian
    plot(f,real(D64f),'k'); axis tight; xlabel('Frequency \DeltaT=1'); ylabel('D_{2m+1}');
    title('Dirichlets Kernal m=64');
    print -dpng S2L04_d64.png

   
end

p4 = 1;
if( p4 == 1 ) 
    % Plots of Dirichlet's kernel.
    f = [-0.5:0.002:0.5]; 
    % Now generate function from Figure 92, Page 92, PW
    f1 = 1/4-1/50; f2 = 1/4+1/50; S=10000;
    gpf = exp(-S*(f-f1).*(f-f1))+exp(-S*(f+f1).*(f+f1))+exp(-S*(f-f2).*(f-f2))+exp(-S*(f+f2).*(f+f2));

    % Generate each of the terms in the sum
    m = 4; D4f = zeros(1,length(f));
    for t = -m:m
        D4f = D4f + exp(i*2*pi*f*t)/(2*m+1);
    end
    % Now convolute D4f and gpf
    gp4f = conv(gpf,real(D4f))/sum(real(D4f)); 
    is = (length(f)-1)/2+1; ie = is + length(f)-1;
    h1 = figure(9); set(h1,'Position',[100 300 750 200],'PaperPositionMode','auto');
    % Plot functiion and Gassian
    plot(f,gpf,'k',f,gp4f(is:ie),'r'); axis tight; 
    title('Concolution with Dirichlets Kernal m= 4 '); xlabel('Frequency \DeltaT=1'); 
    print -dpng S2L04_R1d04.png
    
    % Generate each of the terms in the sum
    m = 16; D16f = zeros(1,length(f));
    for t = -m:m
        D16f = D16f + exp(i*2*pi*f*t)/(2*m+1);
    end
    % Now convolute D4f and gpf
    gp16f = conv(gpf,real(D16f))/sum(real(D16f)); 
    is = (length(f)-1)/2+1; ie = is + length(f)-1;
    h2 = figure(10); set(h2,'Position',[120 280 750 200],'PaperPositionMode','auto');
    % Plot functiion and Gassian
    plot(f,gpf,'k',f,gp16f(is:ie),'r'); axis tight; 
    title('Concolution with Dirichlets Kernal m=16 '); xlabel('Frequency \DeltaT=1'); 
    print -dpng S2L04_R1d16.png

    % Generate each of the terms in the sum
    m = 64; D64f = zeros(1,length(f));
    for t = -m:m
        D64f = D64f + exp(i*2*pi*f*t)/(2*m+1);
    end
    % Now convolute D4f and gpf
    gp64f = conv(gpf,real(D64f))/sum(real(D64f)); 
    is = (length(f)-1)/2+1; ie = is + length(f)-1;
    h3 = figure(11); set(h3,'Position',[140 260 750 200],'PaperPositionMode','auto');
    % Plot functiion and Gassian
    plot(f,gpf,'k',f,gp64f(is:ie),'r'); axis tight; 
    title('Concolution with Dirichlets Kernal m=64 '); xlabel('Frequency \DeltaT=1'); 
    print -dpng S2L04_R1d64.png
    
end

t5 = 0;
if( t5 == 1 )
    
    f = [-0.5:0.002:0.5]; 
    % Now generate function from Figure 92, Page 92, PW
    f1 = 1/4-1/50; f2 = 1/4+1/50; S=10000;
    gpf = exp(-S*(f-f1).*(f-f1))+exp(-S*(f+f1).*(f+f1))+exp(-S*(f-f2).*(f-f2))+exp(-S*(f+f2).*(f+f2));
    m = 64; 
    gt = zeros(1,2*m+1);
    for n = 0:m
        gt(n+1) = sum(gpf.*exp(i*2*pi*f*n));
        gt(2*m+1-n) = sum(gpf.*exp(-i*2*pi*f*n));
    end
    gt = gt /(length(f)); gt(1)
    gt(1) = 0;
    % Now compute power spectrum
    fs = exp(-i*2*pi*f);
    gp04e = zeros(1,length(f)); 
    for n = 0:m
        gp04e = gp04e + gt(n+1)*exp(-i*2*pi*f*n) + gt(length(gt)-n)*exp(i*2*pi*f*n);
    end
    %gp04e = gp04e - m/2/(length(f)-1);
    hold off, plot(f,real(gp04e),'b',f,imag(gp04e),'r'); hold on
    plot(f,gp64f(is:ie),'k');
end

p5 = 1;
if( p5 == 1 ) 
    % Look at Gibbs phenomena
    f = [-0.5:0.002:0.5]; 
    % Now generate function from Figure 93, Page 93, PW
    f0 = 1/4;
    gpf = zeros(1,length(f));
    ind = find(abs(f)<=f0);
    gpf(ind) = 1.0;

    % Generate each of the terms in the sum
    m = 4; D4f = zeros(1,length(f));
    for t = -m:m
        D4f = D4f + exp(i*2*pi*f*t)/(2*m+1);
    end
    % Now convolute D4f and gpf
    gp4f = conv(gpf,real(D4f))/sum(real(D4f)); 
    is = (length(f)-1)/2+1; ie = is + length(f)-1;
    h1 = figure(12); set(h1,'Position',[100 300 750 200],'PaperPositionMode','auto');
    % Plot functiion and Gassian
    plot(f,gpf,'k',f,gp4f(is:ie),'r'); axis tight; 
    title('Box Car Concolution with Dirichlets Kernal m= 4 '); xlabel('Frequency \DeltaT=1'); 
    rms = std(real(gpf-gp4f(is:ie))); maxerr = max(abs(gpf-gp4f(is:ie)));
    fprintf(' For m %3d RMS %5.2f  Max %5.2f\n',m, rms, maxerr)
    print -dpng S2L04_R2d04.png
    
    % Generate each of the terms in the sum
    m = 16; D16f = zeros(1,length(f));
    for t = -m:m
        D16f = D16f + exp(i*2*pi*f*t)/(2*m+1);
    end
    % Now convolute D4f and gpf
    gp16f = conv(gpf,real(D16f))/sum(real(D16f)); 
    is = (length(f)-1)/2+1; ie = is + length(f)-1;
    h2 = figure(13); set(h2,'Position',[120 280 750 200],'PaperPositionMode','auto');
    % Plot functiion and Gassian
    plot(f,gpf,'k',f,gp16f(is:ie),'r'); axis tight; 
    title('Box Car Concolution with Dirichlets Kernal m=16 '); xlabel('Frequency \DeltaT=1'); 
    rms = std(real(gpf-gp16f(is:ie))); maxerr = max(abs(gpf-gp16f(is:ie)));
    fprintf(' For m %3d RMS %5.2f  Max %5.2f\n',m, rms, maxerr)
    print -dpng S2L04_R2d16.png

    % Generate each of the terms in the sum
    m = 64; D64f = zeros(1,length(f));
    for t = -m:m
        D64f = D64f + exp(i*2*pi*f*t)/(2*m+1);
    end
    % Now convolute D4f and gpf
    gp64f = conv(gpf,real(D64f))/sum(real(D64f)); 
    is = (length(f)-1)/2+1; ie = is + length(f)-1;
    h3 = figure(14); set(h3,'Position',[140 260 750 200],'PaperPositionMode','auto');
    % Plot functiion and Gassian
    plot(f,gpf,'k',f,gp64f(is:ie),'r'); axis tight; 
    title('Box Car Concolution with Dirichlets Kernal m=64 '); xlabel('Frequency \DeltaT=1'); 
    rms = std(real(gpf-gp64f(is:ie))); maxerr = max(abs(gpf-gp64f(is:ie)));
    fprintf(' For m %3d RMS %5.2f  Max %5.2f\n',m, rms, maxerr)
    print -dpng S2L04_R2d64.png
    
end
    
p6 = 1;
if( p6 == 1 ) 
    % Look at Gibbs phenomena
    f = [-0.5:0.002:0.5]; 
    % Now generate function from Figure 93, Page 93, PW
    f0 = 1/4;
    gpf = zeros(1,length(f));
    ind = find(abs(f)<=f0);
    gpf(ind) = 1.0;

    % Generate each of the terms in the sum
    m = 4; D4f2 = zeros(1,length(f));
    for t = -m:m
        D4f2 = D4f2 + (1-abs(t)/m)*exp(i*2*pi*f*t)/m^2;
    end
    % Now convolute D4f and gpf
    gp4f = conv(gpf,real(D4f2))/sum(real(D4f2)); 
    is = (length(f)-1)/2+1; ie = is + length(f)-1;
    h1 = figure(12); set(h1,'Position',[100 300 750 200],'PaperPositionMode','auto');
    % Plot functiion and Gassian
    plot(f,gpf,'k',f,gp4f(is:ie),'r'); axis tight; 
    title('Box Car Concolution with Fejer Kernal m= 4 '); xlabel('Frequency \DeltaT=1'); 
    rms = std(real(gpf-gp4f(is:ie))); maxerr = max(abs(gpf-gp4f(is:ie)));
    fprintf(' For m %3d RMS %5.2f  Max %5.2f\n',m, rms, maxerr)
    print -dpng S2L04_R2f04.png
    
    % Generate each of the terms in the sum
    m = 16; D16f2 = zeros(1,length(f));
    for t = -m:m
        D16f2 = D16f2 + (1-abs(t)/m)*exp(i*2*pi*f*t)/m^2;
    end
    % Now convolute D4f and gpf
    gp16f = conv(gpf,real(D16f2))/sum(real(D16f2)); 
    is = (length(f)-1)/2+1; ie = is + length(f)-1;
    h2 = figure(13); set(h2,'Position',[120 280 750 200],'PaperPositionMode','auto');
    % Plot functiion and Gassian
    plot(f,gpf,'k',f,gp16f(is:ie),'r'); axis tight; 
    title('Box Car Concolution with Fejer Kernal m=16 '); xlabel('Frequency \DeltaT=1'); 
    rms = std(real(gpf-gp16f(is:ie))); maxerr = max(abs(gpf-gp16f(is:ie)));
    fprintf(' For m %3d RMS %5.2f  Max %5.2f\n',m, rms, maxerr)
    print -dpng S2L04_R2f16.png

    % Generate each of the terms in the sum
    m = 64; D64f = zeros(1,length(f));
    for t = -m:m
        D64f2 = D64f2 + (1-abs(t)/m)*exp(i*2*pi*f*t)/m^2;
    end
    % Now convolute D4f and gpf
    gp64f = conv(gpf,real(D64f2))/sum(real(D64f2)); 
    is = (length(f)-1)/2+1; ie = is + length(f)-1;
    h3 = figure(14); set(h3,'Position',[140 260 750 200],'PaperPositionMode','auto');
    % Plot functiion and Gassian
    plot(f,gpf,'k',f,gp64f(is:ie),'r'); axis tight; 
    title('Box Car Concolution with Fejer Kernal m=64 '); xlabel('Frequency \DeltaT=1'); 
    rms = std(real(gpf-gp64f(is:ie))); maxerr = max(abs(gpf-gp64f(is:ie)));
    fprintf(' For m %3d RMS %5.2f  Max %5.2f\n',m, rms, maxerr)
    print -dpng S2L04_R2f64.png
    
end

p7 = 1;
if( p7 == 1)
    % Figures showing aliasing
    t=[0:0.01:8]; ts=[0:8];
    % 
    TS = cos(ts);
    k = 1; TA = cos((1+2*pi*k)*t);
    h1 = figure(15); set(h1,'Position',[100 300 750 200],'PaperPositionMode','auto');
    plot(t,TA,'k'); hold on; 
    hs = plot(ts,TS,'rs-'); set(hs,'MarkerSize',6,'MarkerFaceColor','r');
    axis tight; title('Alias k=1: cos((1+2*\pi*k)*t)')
    print -dpng S2L04_Alias1.png
    
    k = 3; TA = cos((1+2*pi*k)*t);
    h2 = figure(16); set(h2,'Position',[120 280 750 200],'PaperPositionMode','auto');
    plot(t,TA,'k'); hold on; 
    hs = plot(ts,TS,'rs-'); set(hs,'MarkerSize',6,'MarkerFaceColor','r');
    axis tight; title('Alias k=3: cos((1+2*\pi*k)*t)')
    print -dpng S2L04_Alias3.png

    
end