% M-file for Sec2 Lec 07.
% 
p1 = 0 ;
if( p1 == 1 ) 
    % ERP File is index, MJD, Xpole (mas) Ypole (mas) LOD (ms).  We noise
    % up the xpole position for this plot
    ERP_all = load('mit_060422.erp');
    t = ERP_all(:,1); Xp = ERP_all(:,3)-mean(ERP_all(:,3)); 
    randn('seed',0);
    Xr = Xp + randn(100,1);
    g1=[1/4 1/2 1/4];
    X1s = conv(Xr,g1);
    LenXr = length(Xr);
    X1r = Xr([2:LenXr-1])-X1s(3:LenXr); tr=t([2:LenXr-1]);
    f=[0:0.01:0.5];
    Gf1 = cos(pi*f).^2; Gr1=sin(pi*f).^2;
    
    h1 = figure(1); set(h1,'Position',[100 200 750 500],'PaperPositionMode','auto');
    subplot(2,1,1); hold off;
    plot(f,Gf1,'k',f,Gr1,'r');title('Transferr function G1  '); xlabel('Frequency (cpd)'); ylabel('G(f)');
    axis tight; ylim([0 4]);
    subplot(2,1,2);
    plot(t,Xr,'k+',t([2:LenXr-1]),X1s([3:LenXr]),'b-');  hold on;
    plot(tr,X1r,'ro'); 
    title('Time Series  '); xlabel('Time (days)'); ylabel('Xpole and Residual (mas)'); 
    hold off;
    
    print -dpng S2L07_F172.png 

    %%%% Now repeat with the g2 filter: This us G1 multiplied by 2.
    g2=[1/2 1 1/2];
    X2s = conv(Xr,g2);
    LenXr = length(Xr);
    X2r = Xr([2:LenXr-1])-X2s(3:LenXr); tr=t([2:LenXr-1]);
    f=[0:0.01:0.5];
    Gf2 = conj(0.5*exp(i*2*pi*f)+1+0.5*exp(-i*2*pi*f)).*(0.5*exp(i*2*pi*f)+1+0.5*exp(-i*2*pi*f))
    Gr2 = conj((exp(i*2*pi*f)+exp(-i*2*pi*f))/2).*(exp(i*2*pi*f)+exp(-i*2*pi*f))/2;
    
    h2 = figure(2); set(h2,'Position',[120 180 750 500],'PaperPositionMode','auto');
    subplot(2,1,1); hold off;
    plot(f,Gf2,'k',f,Gr2,'r');title('Transferr function G2 = 2*G1 '); xlabel('Frequency (cpd)  '); ylabel('G(f)  ');
    axis tight; ylim([0 4]);
    subplot(2,1,2);
    plot(t,Xr,'k+',t([2:LenXr-1]),X2s([3:LenXr]),'b-');  hold on;
    plot(tr,X2r,'ro'); 
    title('Time Series  '); xlabel('Time (days)  '); ylabel('Xpole and Residual (mas)'); 
    hold off;
    
    print -dpng S2L07_F173.png 
    
end

p2 = 1 ;
if( p2 == 1 ) 
    % LSQ band pass filter K=65, W=0.1
    f=[0:0.001:0.5]; u=[0:32]; W=0.1;
    % Generate postive side of gu
    gu = zeros(1,33);
    gu(1)=2*W; gu([2:33])=sin(2*pi*W*u([2:33]))./(pi*u([2:33]));
    GK = zeros(length(f),1);
    for n = 1:length(f);
       GK(n) = gu(1)+sum(gu([2:33]).*exp(-i*2*pi*f(n)*u([2:33])))+sum(gu([2:33]).*exp(i*2*pi*f(n)*u([2:33])));
    end
    db_GK2 = 10*log10(conj(GK).*GK);
    h1 = figure(3); set(h1,'Position',[100 200 750 400],'PaperPositionMode','auto');
    plot(f,db_GK2,'k'); title('Least squares low pass filter; K=65 ');
    xlabel('Frequency'); ylabel('|G_K(f)|^2 (dB)');
    hold on; ylim([-70 10]); plot([W W],[-65 5],'r'); hold off
    
    print -dpng S2L07_F177.png 
    
    % Use a triangular version of impulse respomse
    gu = zeros(1,33); K=65;
    cu = 1-2*u/(K+1); 
    gu(1)=2*W; gu([2:33])=sin(2*pi*W*u([2:33]))./(pi*u([2:33]));
    gu = gu.*cu/(sum(gu.*cu));
    GK = zeros(length(f),1);
    for n = 1:length(f);
       GK(n) = gu(1)+sum(gu([2:33]).*exp(-i*2*pi*f(n)*u([2:33])))+sum(gu([2:33]).*exp(i*2*pi*f(n)*u([2:33])));
    end
    db_GK2 = 10*log10(conj(GK).*GK);
    h2 = figure(4); set(h2,'Position',[110 180 750 400],'PaperPositionMode','auto');
    plot(f,db_GK2,'k'); title('Triangular factor low pass filter; K=65 ');
    xlabel('Frequency'); ylabel('|G_K(f)|^2 (dB)');
    hold on; ylim([-70 10]); plot([W W],[-65 5],'r'); hold off

    print -dpng S2L07_F178.png 
    
end