close all; clear all;

load('wavelet');
a =  wavelet'/max(wavelet);

% ------------------------------------------------------------------
% Primer filtro inverso de Wiener 
% ------------------------------------------------------------------

N  = length(a);              % tamaño de entrada   
Nf = 200;                    % longitud del operador (filtro "largo") 
Nd = N + Nf - 1 ;            % dimensión de la salida
d  = [1 ; zeros(Nd-1,1) ];   % salida deseada d = (1,0,....,0)
A  = convmtx(a, Nf).';       % matriz de convolución
R  = A.'*A;                  % matriz de autocorrelación
g  = A.'*d;                  % correlación entre entrada y salida
f  = inv(R)*g;               % filtro inverso Wiener

% ------------------------------------------------------------------
% Segundo filtro inverso de Wiener (ond. equivalente de fase mínima)
% ------------------------------------------------------------------

Nf = length(f);              % tamaño de entrada
Ny = N;                      % longitud del operador (igual a ondicula) 
Nd = Nf + Ny - 1;            % longitud de salida deseada
d  = [1 ; zeros(Nd-1,1) ];   % salida deseada d = (1,0,....,0)
A2 = convmtx(f, Ny);         % matriz de convolución
R2 = A2.'*A2;                % matriz de autocorrelación
g2 = A2.'*d;                 % correlación entre entrada y salida
y  = inv(R2)*g2;             % filtro inverso == ondicula eqv fase minima

% Escalamos la ondícula de fase mínima para que tenga autocorrelacion a lag 0
% (o sea, energia) que la ondicula original
k  = sqrt(R(1,1)/sumsq(y)); 
y = k*y;

figure (1); clf;

subplot(211)
plot((0:N-1), a, 'b-o;wavelet;')
axis([0 N-1 -1.2 1.2])

subplot(212)
plot((0:N-1), y, 'r-o;fase minima;');
axis([0 N-1 -1.2 1.2])
