%%%%%% FUNÇÃO DO FILTRO %%%%%%

function [Ho p]=analog_butter(dp,dr,or,op)
  delta1=sqrt((2*dp-dp^2)/(1-dp)^2);
  delta2=sqrt((1-dr^2)/(dr^2));
  N=ceil((log((delta2/delta1)))/log(or/op))

  c=((delta1/op^N)+(delta2/or^N))/2
  o=0:0.01:1;

  H2=1./(1+(c^2)*(o.^(2*(N-1))));
  H=sqrt(H2);


  for k=1:2*(N-1)
    sk(k)=j*(1/(c)^(1/(N-1)))*exp(j*((pi/2*(N-1))+(pi*(k-1)/(N-1))));
  end
  
  V=(real(sk)~=0);
  V=(real(sk)<0);
  idx=0;
  for k=1:2*(N-1)
      if(V(k)==1)
         idx=idx+1;
         p(idx)=sk(k); 
      end
  end

  prod=1;
  for i=1:idx
      prod=prod*(-p(i)); 
  end

  Ho=prod;
end;
%%%%%%%%%%%%%%%%%%

%%%%%% EX 2 %%%%%%
wp=0.1*pi; op=1; wr=0.5*pi;
% normalização
or=1/(tan(wp/2));
%
dr=0.01; dp=0.05;

[Ho p]=analog_butter(dp,dr,or,op);

[a m]=size(p);
prodt=1; prod=1; prods=1;
pd=zeros(1,m); num=zeros(1,m); numt=[1]; dent=[1];

for i=1:m
  prodt=prodt*((-p(i))/(1-p(i)));
  pd(i)=((1+p(i))/(1-p(i)));
  num(i)=(-1);
  numt=conv(numt,[1 -1]);
  dent=conv(dent,[1 pd(i)]);
end
Ho=prodt;
num=poly(num);

%%%%%% ITEM A %%%%%%
num=Ho*num; den=poly(pd);

[Hx w]=freqz(num,den);
figure(1);
hold on
%desnormalização
plot(w*(tan(wp/2))/pi,abs(Hx))
%
title('RF PB');

%%%%%% ITEM B %%%%%%
numt=Ho*numt;
[Ht w]=freqz(numt,dent);
figure(2);
hold on
plot(w/pi*(tan(wp/2)),abs(Ht))
title('RF PA');

%%%%%% ITEM C %%%%%%
num=real(num); den=real(den); ts=0.5;
sys = tf(num,den,ts)
numt=real(numt); dent=real(dent); ts=0.5;
sys = tf(numt,dent,ts)

Wp=0.1;  Ws=0.5; Rp= -20*log10(1-0.05); Rs= 40;
[n,Wn] = buttord(Wp,Ws,Rp,Rs);
[z,p,k] = butter(n,Wn);
sos = zp2sos(z,p,k)

figure(3)
[Hd,w]=freqz(sos,500,2);
plot(w,abs(Hd));
title('Questão 3 - Comparação');
%%%%%%%%%%%%%%%%%%

%%%%%% EX 3 %%%%%%
Wp=0.1; Ws=0.5; Rp= -20*log10(1-0.05); Rs= 40;
[n,Wn] = buttord(Wp,Ws,Rp,Rs);
[z,p,k] = butter(n,Wn);
sos = zp2sos(z,p,k);
[Hd,w]=freqz(sos,500,2);

figure(3)
plot(w,abs(Hd),'k');
hold on
%%%%%%%%%%%%%%%%%%

%%%%%% EX 4 %%%%%%
Ap=0.05; Ar=10^(Rs/(-20));
n1=(((-10*log10(Ap*Ar))-13)/(2.324*(Ws-Wp)*pi))+1; n1=fix(n1); n1=n1+2;
[h err]= firpm([n1], [0 Wp Ws 1], [1 1 0 0], [1 Ap/Ar])

figure(4)
hold on
[Hd1,w]=freqz(h);
plot(w/pi,abs(Hd1),'r');
title('Questão 4 - Comparação');

dev = [(10^(Rp/20)-1)/(10^(Rp/20)+1)  10^(-Rs/20)];
a=[1 0]; f=[Wp Ws];
[n,fo,ao,w] = firpmord(f,a,dev);
b = firpm(n+3,fo,ao,w);
[Hd2,w]=freqz(b);

plot(w/pi,abs(Hd2)); plot(w/pi,abs(Hx))
legend('Chebyshev + Kaiser','Chebyshev + firpmod','IIR projetado');
%%%%%%%%%%%%%%%%%%

%%%%%% EX 5 %%%%%%
figure(5); plot(n1,s); title('sinal limpo');

figure(6); plot(n1,sr); title('sinal original');

%%%%%% ITEM A %%%%%%
N1=1001; n1=0:N1-1; s=sin(((pi/20)*n1))+cos(((19*pi)/20*n1));
v=.04*randn(1001,1); v=v'; sr=s+v;

filtered_pb=filter (num,den, sr);
figure(7)
plot(n1,filtered_pb)
title('sinal filtrado');
    
figure(8)
hold on 
plot(n1,sr); plot(n1,filtered_pb); plot(n1,s);
title('Comparação');
legend('Sinal limpo','Sinal filtrado','Sinal limpo');

%%%%%% ITEM B %%%%%%
filtered_pa=filter (numt,dent, sr);
figure(9)
plot(n1,filtered_pa)
title('Sinal filtrado');
    
figure(12)
hold on 
plot(n1,sr); plot(n1,filtered_pa); plot(n1,s)
title('Comparação');
legend('Sinal limpo','Sinal filtrado','Sinal original');
%%%%%%%%%%%%%%%%%%