%% =================================================================
%% DECONVOLUCIÓN DE ONDÍCULA MAL CONDICIONADA POR CEROS EN EL ESPECTRO
%% 
%% Agregando ceros en el espectro de amplitud de una ondícula de fase
%% mínima: ONDICULA MAL CONDICIONADA.
%% La matriz de autocorrelación de la ondícula resultante estará
%% mal condicionada.
%% =================================================================

clear all; close;

%% ---------------------------------------------------------------==
%% Ondícula de fase mínima con ventana de Hamming
%% ---------------------------------------------------------------==

N  = 40 ;   % muestras de la ondícula
w0 = pi/3 ; % frecuencia central de la ondícula de fase mínima a generar
dW = w0/N ; % paso en frecuencia

% Generamos una ondícula de fase mínima y le aplicamos un taper de Hamming:
for n = 0:N-1
  hamm(n+1) = 0.54 - 0.46*cos(pi*(N-1-n)/(N-1));  % taper de Hamming
  a(n+1)    = hamm(n+1)*(0.8^n)*cos((w0-dW*n)*n); % ondicula de fase mín.
end

% -----------------------------------------------------------
% Agrego un dipolo con un cero muy cercano el círculo unidad
% y agrego también el conjungado del dipolo para que ondícula
% la ondícula continue siendo real.
% -----------------------------------------------------------

w0 = (7/8)*pi; % frecuencia del cero a agregar
dipo = [1,-0.999*exp(i*w0)];
a0 = conv(a, dipo);
a0 = conv(a0, conj(dipo));
a0 = real(a0);
N0 = length(a0); % tamaño de la nueva ondícula con el cero

%% ---------------------------------------------------------------==
%% Filtro inverso Wiener
%% ---------------------------------------------------------------==

Nf  = 15;                    % longitud del filtro

% ------------------------
% Filtro inverso sin ceros
% ------------------------

Nd = length(a) + Nf - 1 ;    % dimensión de la salida
d  = zeros(Nd,1); d(1) = 1 ; % salida deseada
A  = convmtx(a,Nf).' ;       % matriz de convolución
R  = A.'*A ;                 % matriz de autocorrelación
g  = A.'*d ;                 % correlación entre entrada y salida
fw = inv(R)*g;               % filtro Wiener
ad = conv(a,fw);             % deconvolución

% ------------------------
% Filtro inverso con ceros
% ------------------------

Nd0 = length(a0) + Nf - 1 ;     % dimensión de la convolucion
d0  = zeros(Nd0,1); d0(1) = 1 ; % salida deseada
A0  = convmtx(a0,Nf)' ;         % matriz de convolución    
R0  = A0.'*A0 ;                 % matriz de autocorrelación
g0  = A0.'*d0 ;                 % correlación entre entrada y salida
fw0 = inv(R0)*g0 ;              % filtro Wiener
a0d = conv(a0,fw0) ;            % deconvolución

%% --------------------------------------------------------------
%% Zero padding y espectros de amplitud
%% --------------------------------------------------------------

N2  = 1024 ;
M2  = N2/2 + 1 ;
WN2 = 2*pi/N2 ;

% ---------
% Sin ceros
% ---------

a2     = a;
a2(N2) = 0;
A2     = abs(fft(a2));  % espectro de amplitud sin ceros

ad2     = ad;
ad2(N2) = 0;
AD2     = abs(fft(ad2));  % espectro de amplitud ondicula deconv

% ---------
% Con ceros
% ---------

a02     = a0;
a02(N2) = 0;
A02     = abs(fft(a02)); % espectro de amplitud con ceros

a0d2     = a0d;
a0d2(N2) = 0;
A0D2     = abs(fft(a0d2));  % espectro de amplitud ondicula deconv

%%--------------------------------------------------------------
%% GRAFICOS 
%% -------------------------------------------------------------

 dB = 0; % Grafico los espectros en escala lineal
%db = 1; % Grafico los espectros en decibeles

figure(1); clf;

% ------------------------------------
% Gráficos sin ceros
% ------------------------------------

% Ondícula

subplot(4,2,1)
plot((0:N-1),a,';ondicula;b-')
title('Sin ceros')
axis([-1 N -1.2 1.2])

% Espectro de amplitud

subplot(4,2,3)
if dB==0 % (escala lineal)
  plot((0:M2-1).*WN2,A2(1:M2)./max(A2(1:M2)),';espectro;r-')
else % (decibeles)
  plot((0:M2-1).*WN2,20.*log10(A2(1:M2)./max(A2)),';espectro;r-')
end
axis([0 pi -0.2 1.2])
xt      = [0;pi/4;pi/2;3*pi/4;pi];
xtlabel = ["0";"{/Symbol p}/4";"{/Symbol p/2}";"3/4{/Symbol p}";"{/Symbol p}"];
set(gca,"xtick",xt,"xticklabel",xtlabel);

% Ondicula deconvolucionada

subplot(4,2,5)
plot((0:N-1),ad(1:N)/max(ad),';decon;g-')
axis([-1 N -0.2 1.2])

% Espectro de amplitud de ondicula deconvolucionada

subplot(4,2,7)
if dB==0 % (escala lineal)
  plot((0:M2-1).*WN2,AD2(1:M2)./max(AD2(1:M2)),';esp. decon;c-')
else % (decibeles)
  plot((0:M2-1).*WN2,20.*log10(AD2(1:M2)./max(AD2)),';esp. decon;c-')
end
legend("location","southwest")   
axis([0 pi -0.2 1.2])
xt      = [0;pi/4;pi/2;3*pi/4;pi];
xtlabel = ["0";"{/Symbol p}/4";"{/Symbol p/2}";"3/4{/Symbol p}";"{/Symbol p}"];
set(gca,"xtick",xt,"xticklabel",xtlabel);


% ------------------------------------
% Gráficos con ceros
% ------------------------------------

% Ondícula

subplot(4,2,2)
plot((0:N0-1),a0./max(a0),';ondicula;b-')
title('Con ceros')
axis([-1 N0 -1.2 1.2])

% Espectro de amplitud

subplot(4,2,4)
if dB==0 % (escala lineal)
  plot((0:M2-1).*WN2,A02(1:M2)./max(A02(1:M2)),';espectro;r-')
  axis([0,pi,-0.2,1.2])
else    % (decibeles)
  plot((0:M2-1).*WN2,20.*log10(A02(1:M2)./max(A02)),';espectro;r-')
end
xt      = [0;pi/4;pi/2;3*pi/4;pi];
xtlabel = ["0";"{/Symbol p}/4";"{/Symbol p/2}";"3/4{/Symbol p}";"{/Symbol p}"];
set(gca,"xtick",xt,"xticklabel",xtlabel);

% Deconvolución

subplot(4,2,6)
plot((0:N0-1),a0d(1:N0)/max(a0d),';decon;g-')
axis([-1 N0 -0.2 1.2])

% Espectro de amplitud de ondicula deconvolucionada

subplot(4,2,8)
if dB==0 % (escala lineal)
  plot((0:M2-1).*WN2,A0D2(1:M2)./max(A0D2(1:M2)),';esp. decon;c-')
else % (decibeles)
  plot((0:M2-1).*WN2,20.*log10(A0D2(1:M2)./max(A0D2)),';esp. decon;c-')
end
legend("location","southwest")   
axis([0 pi -0.1 1.1])
xt      = [0;pi/4;pi/2;3*pi/4;pi];
xtlabel = ["0";"{/Symbol p}/4";"{/Symbol p/2}";"3/4{/Symbol p}";"{/Symbol p}"];
set(gca,"xtick",xt,"xticklabel",xtlabel);

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

print('add_freq_zeros.eps','-depsc')
