clear all
close all
clc
T=0.1;
%%%ws=2*pi/T;
W=2;
tc=0:0.0001:10; %%% simulating the continous domain
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xc=sin(W*tc);
%xs=sin(W*ts)+j*cos(W*ts);
ts=0:T:10;
sig=0.5;
f=sin(W*ts);
noise=sig*randn(1,length(ts));
xs=f+noise;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(3,1,1)
plot(ts,f,'k-')
title('hidden signal')
xlabel('t')
axis([min(ts) max(ts) -3 3])
subplot(3,1,2)
plot(ts,noise,'r-')
title('noise')
xlabel('t')
axis([min(ts) max(ts) -3 3])
subplot(3,1,3)
plot(ts,xs,'b-')
axis([min(ts) max(ts) -3 3])
title('observed signals')
xlabel('t')


typefilter=1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
switch typefilter
    case 1
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%% filter (MA) FIR 1
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       % M=5;
        M=10;
       % M=20;
       % M=50;
        b=1/M*ones(1,M);
        %%%%%%
        for n=M:length(xs) %%% n seria n-2.... preguntar a Luca
            Xv=xs(n:-1:n-M+1);
        y(n)=b*Xv';
        end
        figure
        subplot(3,1,1)
        plot(ts(M:end),f(M:end),'k-')
        title('hidden signal')
        xlabel('t')
        axis([min(ts) max(ts) -3 3])
        subplot(3,1,2)
        plot(ts(M:end),xs(M:end),'b-')
        axis([min(ts) max(ts) -3 3])
        title('observed signals')
        xlabel('t')
        subplot(3,1,3)
        plot(ts(M:end),y(M:end),'m-')
        title('FILTERED SIGNAL')
        xlabel('t')
        axis([min(ts) max(ts) min(y)-1 max(y)+1])
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        figure
        W=-pi:0.01:pi;
        for k=1:M
            FourierElement(k,:)=b(k).*exp(-(k-1)*j*W); %%% cuidado k y k-1....luego explico
        end
        F=sum(FourierElement);
        plot(W,abs(F))
        set(gca,'fontSize',20,'fontWeight','Bold')
        xlabel('\Omega')
         hold on
        plot([-pi -pi],[0 max(abs(F))],'r--','LineWidth',2)
        plot([pi pi],[0 max(abs(F))],'r--','LineWidth',2)
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case 2
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%% filter (MA) FIR 2 
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
       b=[1 -1];
        for n=2:length(xs) %%% n seria n-2.... preguntar a Luca
          y(n)=b(1)*xs(n)*b(2)*xs(n-1);
        end
        figure
        subplot(3,1,1)
        plot(ts(2:end),f(2:end),'k-')
        title('hidden signal')
        xlabel('t')
        axis([min(ts) max(ts) -3 3])
        subplot(3,1,2)
        plot(ts(2:end),xs(2:end),'b-')
        axis([min(ts) max(ts) -3 3])
        title('observed signals')
        xlabel('t')
        subplot(3,1,3)
        plot(ts(2:end),y(2:end),'m-')
        title('FILTERED SIGNAL')
        xlabel('t')
        axis([min(ts) max(ts) min(y)-1 max(y)+1]) 
        figure
         W=-pi:0.01:pi;
        F2=1+exp(-j*W)-2*exp(-2*j*W);
         plot(W,abs(F2))
        set(gca,'fontSize',20,'fontWeight','Bold')
        xlabel('\Omega') 
        hold on
        plot([-pi -pi],[0 max(abs(F2))],'r--','LineWidth',2)
        plot([pi pi],[0 max(abs(F2))],'r--','LineWidth',2)
        title('Respuesta en Frecuencia del filtro')
     case 3
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%% filter (AR) IIR 
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%     
       p=0.9;
       y(1)=0;  %%% condiciones iniciales
        for n=2:length(xs) %%% n seria n-2.... preguntar a Luca
          y(n)=p*y(n-1)+xs(n);
        end
      figure
        subplot(3,1,1)
        plot(ts(2:end),f(2:end),'k-')
        title('hidden signal')
        xlabel('t')
        axis([min(ts) max(ts) -3 3])
        subplot(3,1,2)
        plot(ts(2:end),xs(2:end),'b-')
        axis([min(ts) max(ts) -3 3])
        title('observed signals')
        xlabel('t')
        subplot(3,1,3)
        plot(ts(2:end),y(2:end),'m-')
        title('FILTERED SIGNAL')
        xlabel('t')
        axis([min(ts) max(ts) min(y)-1 max(y)+1])   
        %%%%
        figure   
        W=-pi:0.01:pi;
        F3=1./(1-p*exp(-j*W));
        plot(W,abs(F3),'k')
        title('Respuesta en Frecuencia del filtro')
        hold on
         set(gca,'fontSize',20,'fontWeight','Bold')
         xlabel('\Omega') 
        plot([-pi -pi],[0 max(abs(F3))],'r--','LineWidth',2)
        plot([pi pi],[0 max(abs(F3))],'r--','LineWidth',2)
        case 4
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%% filter (AR) IIR 
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%     
       %p=[0.5 0.2];
       %p=[0.2 0.4];
       % p=[0.5 -0.2];
       %  p=[0.5 -0.5];
        p=[0.9 -0.9]; 
       disp('los polos son')
       polos=roots([1 -p(1) -p(2)])
       modulos_polos=abs(polos)
       y(1)=0;  %%% condiciones iniciales
       y(2)=0;  %%% condiciones iniciales
        for n=3:length(xs) %%% n seria n-2.... preguntar a Luca
          y(n)=p(1)*y(n-1)+p(2)*y(n-2)+xs(n);
        end
      figure
        subplot(3,1,1)
        plot(ts(2:end),f(2:end),'k-')
        title('hidden signal')
        xlabel('t')
        axis([min(ts) max(ts) -3 3])
        subplot(3,1,2)
        plot(ts(2:end),xs(2:end),'b-')
        axis([min(ts) max(ts) -3 3])
        title('observed signals')
        xlabel('t')
        subplot(3,1,3)
        plot(ts(2:end),y(2:end),'m-')
        title('FILTERED SIGNAL')
        xlabel('t')
        axis([min(ts) max(ts) min(y)-1 max(y)+1])   
        %%%%
        figure
        
        W=-pi:0.01:pi;
        F4=1./(1-p(1)*exp(-j*W)-p(2)*exp(-2*j*W));
        plot(W,abs(F4),'k')
        hold on
        set(gca,'fontSize',20,'fontWeight','Bold')
        title('Respuesta en Frecuencia del filtro')
         xlabel('\Omega') 
        plot([-pi -pi],[0 max(abs(F4))],'r--','LineWidth',2)
        plot([pi pi],[0 max(abs(F4))],'r--','LineWidth',2)
end

