% =================================================================
% DECONVOLUCIÓN DE ONDÍCULAS EQUIVALENTES DE FASE MÍNIMA, LINEAL y 
% MÁXIMA.
%
% A partir de tres dipolos elementales de fase mínima y de sus 
% dipolos equivalentes de fase máxima, se generan tres ondículas 
% equivalentes de tamaño 7 de fase mínima, lineal, y máxima.
%
% Se deconvolucionan las tres ondículas generadas mediante un 
% filtro inverso de Wiener para el cual la salida deseada está dada 
% por un impulso unitario de  posisición variable. La posición de 
% dicho impulso unitario se determina mediante la variable "delay".
% =================================================================

clear all; close all;

% ==================================================================
% Posición del spike en la salida deseada
% DESCOMENTAR SEGUN CORRESPONDA
% ==================================================================

 delay =  1;  % salida deseada de fase mínima
% delay = 11;  % salida deseada de fase lineal (mixta)
% delay = 21;  % salida deseada de fase máxima

% ==================================================================
% Genero las ondículas
% ==================================================================

% Dipolos de fase mínima

a1 = [1.0; 0.2];
b1 = [1.0;-0.4];
c1 = [1.0; 0.6];

% Dipolos de fase máxima

a2 = [ 0.2;1.0];
b2 = [-0.4;1.0];
c2 = [ 0.6;1.0];

% Ondículas equivalentes

wmin = conv(conv(conv(a1,a1),conv(b1,b1)),conv(c1,c1)); % Fase mínima
wlin = conv(conv(conv(a1,a2),conv(b1,b2)),conv(c1,c2)); % Fase lineal
wmax = conv(conv(conv(a2,a2),conv(b2,b2)),conv(c2,c2)); % Fase máxima

% ==================================================================
% Filtro Wiener
% ==================================================================

Nf = 15;           % tamaño del filtro
Nw = length(wmin); % longitud de ondículas
Nd = Nw + Nf - 1;  % longitud de salida deseada

% --------------
% Salida deseada
% --------------

d        = zeros(Nd,1);  % armo vector de ceros del largo Nd
d(delay) = 1;            % agrego el impulso en la posición delay

% --------------
% Calculo del filtro Wiener
% --------------
for k = [1 2 3];

    if(k==1)    
        a = wmin;
    elseif(k==2)    
        a = wlin;
    else
        a = wmax;
    end
    
A = convmtx(a,Nf);
f = inv(A'*A)*(A'*d);

% -------
% Grafico
% -------
Na=length(a);
figure(k); clf;

subplot(2,2,1)
if(k==1)
    stem((0:Nw-1),a,'b');
    title('Entrada de fase minima')
    axis([-1 Nw])
elseif(k==2)
    stem((0:Nw-1),a,'b');
    title('Entrada de fase lineal')
    axis([-1 Nw])
else
    stem((0:Nw-1),a,'b');
    title('Entrada de fase maxima')
    axis([-1 Nw])
end

subplot(2,2,2)
stem((0:Nd-1),d,'g');
title('Salida deseada')
axis([-1 Nd])

subplot(2,2,3)
stem((0:Nf-1),f,'r');
title('Filtro obtenido')
axis([-1 Nf])

subplot(2,2,4)
stem((0:Nd-1),conv(a,f),'m');
title('Salida del filtro')
axis([-1 Nd])

% -----------------------
% Guardo en un postscript
% -----------------------

% print(strcat('dcon_',num2str(k),'.eps'),'-depsc')

end

% ==================================================================
% ==================================================================
