%% ==================================================================
%% DISEÑO DE FILTRO PASABAJOS DE BUTTERWORTH
%% ==================================================================

clear all; close;

% -------------------------------------------------------------------
% Parámetros del filtro
% -------------------------------------------------------------------

N  = 51 ;   % número de muestras del filtro
No = 10 ;   % orden 
wc = pi/2 ; % frecuencia de corte

% -------------------------------------------------------------------
% Defino el filtro en frecuencia
% -------------------------------------------------------------------

M  = (N-1)/2 ; % Hasta Nyquist (pi)
dw =  2*pi/N ;

for k = 0:M
    wk=k*dw;
    BW(k+1) = sqrt( 1./(1+(wk/wc)^(2*No)) );
end

BW = [ BW(1:M+1) BW(M+1:-1:2) ] ; % filtro simétrico

% -------------------------------------------------------------------
% Calculo el filtro en tiempo
% -------------------------------------------------------------------

bw = real(ifft(BW)) ;
bw = [bw((M+2):N) bw(1:M+1)]; % filtro simetrico en tiempo

% -------------------------------------------------------------------
% Zero padding
% -------------------------------------------------------------------

NN      = 1024 ;
bw2     = bw ;
bw2(NN) = 0 ; 
BW2     = abs(fft(bw2)) ;

% -------------------------------------------------------------------
% Expresión analítica de |H(w)|
% -------------------------------------------------------------------

MM  = NN/2 ;    % Hasta Nyquist (pi)
dw2 = 2*pi/NN ;

for k = 0:MM-1
    H(k+1) = sqrt(1./(1+(k*dw2/wc)^(2*No)));
end

H = [ H(1:MM) H(MM:-1:1) ] ; % Hermiticidad

% ===================================================================
% GRAFICOS
% ===================================================================

figure(1); clf;

% Operador en tiempo

subplot(3,1,1)
stem((-M:M),bw,'b-o') ;
axis([-M,M,-0.2,0.6]);
ylabel('bw_n')
title('Filtro pasabajos Butterworth','fontsize',14)

% Respuesta en frecuencia

xt      = 0:pi/2:pi;
xtlabel = ["0";"{/Symbol p}/2";"{/Symbol p}"];
yt      = [0 1];

subplot(3,1,2)
hold on
plot((0:N-1)*dw,BW,';BW({/Symbol w});bo');
plot((0:NN-1)*dw2,BW2,';BW_{pad}({/Symbol w});r-');
plot((0:NN-1)*dw2,H,';H({/Symbol w});g-');
hold off
ylabel('H(\omega)')
set(gca,"xtick",xt,"xticklabel",xtlabel,"ytick",yt);
axis([0,pi,-0.2,1.2]);

% Respuesta en frecuencia en decibeles

subplot(3,1,3)
plot((0:NN-1).*dw2,20.*log10(BW2),';BW_{pad}({/Symbol w});r-')
axis([0,pi,-120,20]);
ylabel('H(\omega)[dB]')
set(gca,"xtick",xt,"xticklabel",xtlabel);

% -------------------------------------------------------------------
% Salida
% -------------------------------------------------------------------

print('butterworth.eps','-depsc');


