% =================================================================
% FILTRO BUTTERWORTH APLICADO EN EL DOMINIO DE LA TDF a "wavelet"
%
% OBS: el directorio debe contener la ondícula "wavelet"
% =================================================================

clear all; close all

% -----------------------
% Lectura del dato
% -----------------------

a = load('wavelet') ;
N = length(a) ;

% -----------------------
% Parametros del filtro Butterworth
% -----------------------
wc = (4/5)*(pi/2); % frecuencia de corte
No = 10;  % orden del filtro (bajo orden)
% No = 30;  % orden del filtro (alto orden)

% ==================================================================
% AGREGO CEROS (descomentar para agregar)
% ==================================================================

% a = [ zeros(N,1) ; a ] ; % agrego N ceros adelante 
% a = [ a ; zeros(N,1) ] ; % agrego N ceros atrás

N = length(a); % recalculo longitud del dato

dw = 2*pi/N;


% Genero el Butterworth en frecuencia

for i = 1:N
    w = dw*(i-1) ;
    if w <= pi
        BW(i) = sqrt(1/(1+(w/wc)^(2*No))) ;
    else 
        BW(i) = sqrt(1/(1+((2*pi-w)/wc)^(2*No))) ;
    end
end

BW=BW';   % traspongo para que quede una columna como el vector "a"

% -----------------------
% Filtro en frecuencia
% -----------------------

% Calculo la respuesta en frecuencia de la ondícula
A = fft(a) ;

% Filtro en frecuencias y anti-transformo
AF = A.*BW ;
af = real(ifft(AF)) ;

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

% Zero padding para la ondícula y transformo
N2     = 1024 ;
dw2    = 2*pi/N2 ;
a2     = a ;
a2(N2) = 0 ;
A2     = fft(a2);

% Zero padding para ondícula filtrada y transformo
af2     = af ;
af2(N2) = 0 ;
AF2     = fft(af2) ;

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

% Rta en frecuencia del filtro BW

Npad = length(BW);
Npad = 1014;
M    = floor(N/2);
%bw   = real(fftshift(ifft(BW)));
bw= real(ifft(BW)); % operador circular (zero padding en el centro)

bw   = [bw(M+1:N)' bw(1:M)' ]; % operador lineal (zero padding al final)

bwF  = abs(fft(bw,Npad)); % zero padding si el operador es lineal
wpad =(0:Npad-1)*(2*pi/Npad);
w =(0:N-1)*(2*pi/N);


figure(1)
% Grafico en tiempo
subplot(2,1,1);
hold on
plot((0:N-1),a,'g;Ondicula;');
plot((0:N-1),af,'r;Ondicula filtrada;');
hold off
% Grafico en frecuencia
subplot(2,1,2);
hold on
plot((0:N2-1)*dw2,abs(A2),'g-;Ondicula;');
plot((0:N-1)*dw,max(abs(A2)).*BW,'bo; BW_k;',"markersize",3);
plot(wpad,max(abs(A2))*bwF,'m;Filtro BW;');
plot((0:N2-1)*dw2,abs(AF2),'r-;Ondicula filtrada;');
hold off
xt      = 0:pi/2:pi;
xtlabel = ["0";"{/Symbol p}/2";"{/Symbol p}"];
set(gca,"xtick",xt,"xticklabel",xtlabel);
axis([0,pi,-0.5,5+max(abs(A2))]);

% Guardo el grafico en un postscript
print('bw_wavelet.eps','-depsc');

figure(2)
hold on
plot(wpad,20*log10(bwF),'m;Filtro BW;');
plot(w,20*log10(BW),'b;Respuesta BW deseada ;');
hold off
axis([0 pi, ])
xlabel('\omega')
ylabel('|BW|, dB')

figure(3)
hold on
plot((0:N-1),a,'g;Ondicula;');
plot((0:N-1),af,'r;Ondicula filtrada;');
hold off
