% =================================================================
% INVERSA CALCULADA MEDIANTE:
%
% a) SERIE GEOMETRICA
% b) FILTRO INVERSO WIENER (SPIKE DECONVOLUTION)
%
% OBS: Este ejercicio debe realizarse primero "en papel". 
%      Este codigo es para controlar la solucion.
% =================================================================

clear all; close all;

% =================================================================
% Parametros
% =================================================================

a  = [1,0.6];   % secuencia de fase minima
Na = length(a); % longitud de la secuencia
Nf =  4;        % longitud del filtro

% Salida deseada:   d = (1,0,....,0)

Nd   = Na + Nf - 1;  % dimension de la convolucion
d    = zeros(Nd,1);
d(1) = 1.0;   

% =================================================================
%  a) Inversa a partir de la serie geometrica:
% =================================================================

for i = 1:Nf
  fg(i) = (-a(2))^(i-1);
end

fg                     % inversa por serie geometrica

% =================================================================
% Verificacion del resultado
% =================================================================

cg      = conv(a,fg.') % convolucion con la entrada
e1      = cg - d       % vector de error
sum2_e1 = e1'*e1       % error cuadratico medio
disp('')

% =================================================================
% b) Inversa con filtro Wiener
% =================================================================
% A f = d => A'A f = A'd ==> f = inv(A'A) (A'd)

A       = convmtx(a,Nf).' % matriz de convolucion
R       = A'*A           % matriz de autocorrelacion. Calcule
                          % xcorr(a) para verificar las entradas de
                          % esta matriz.
Inv     = inv(R)          % inversa de R
g       = A'*d           % correlacion entre entrada y salida deseada
fw      = inv(R)*g        % filtro Wiener inverso

% =================================================================
% Verificacion del resultado
% =================================================================

cw      = conv(a,fw)     % convolucion con la entrada
e2      = cw - d         % vector de error
sum2_e2 = e2'*e2         % error cuadratico medio
disp('')

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

Na=length(a);

clf;

figure(1)
subplot(2,3,1)
stem(a,'b')
axis([0,Na+1,-1.2,1.2])
title('entrada')
subplot(2,3,2)
stem(fg,'r')
title('inversa geometrica')
axis([0,Nf+1,-1.2,1.2])
subplot(2,3,3)
stem(cg,'g')
axis([0,Na+Nf,-1.2,1.2])
title('convolucion')

subplot(2,3,4)
stem(a,'b')
axis([0,Na+1,-1.2,1.2])
title('entrada')
subplot(2,3,5)
stem(fw,'m')
axis([0,Nf+1,-1.2,1.2])
title('inversa Wiener')
subplot(2,3,6)
stem(cw,'g')
axis([0,Na+Nf,-1.2,1.2])
title('convolucion')
