% =================================================================
% FILTRO PREDICTIVO DEL ERROR (PEF: PREDICTION ERROR FILTER)
%

% Si en el filtro de Wiener se  pide que la salida deseada sea igual a
% la entrada  retrasada una muestra,
%
%  f_0 * x_{n-1} + f_1 * x_{n-2} + f_2 * x_{n-3} + ... + f_{N-1} * x_{n-(N-1)} = x_{n+1},
%
% dicho filtro (prediction filter, pf) predecirá la parte NO aleatoria
% de la entrada.  Por lo tanto, podrá usarse para  remover el ruido de
% carácter aleatorio de la señal.
%
% Del mismo modo, la diferencia entre  la entrada y lo predicho por el
% filtro brinda una estimación de  error aleatorio de la señal. Luego,
% si  pf es  el filtro  predictor,  entonces el  vector [1  ; -pf]  se
% denomina filtro predictor del  error (prediction error filter, pef).
% =================================================================

clear all; close all;


Nf = 9;      % Longitud del filtro predictivo

% =================================================================
% CONSTRUCCION DE LA SENAL SINTETICA
% =================================================================

% Generamos una señal sinusoidal (completamente determinística)
N  = 200;    % Longitud de la senal
Nw = 20;     % Longitud de bordes de la ventana para suavizar
w1 = pi/12;  % Frecuencia de la senal

amp_senal=1.0;
amp_ruido=0.8;

for n = 0:N-1
  signal(n+1) = amp_senal*sin(w1*n);
end

% Ventana con bordes de Hamming para truncar el ruido
for n = 0:N-1
  if n < Nw
    win(n+1) = 0.54 - 0.46*cos(pi*n/Nw);
  elseif n >= N-Nw
    win(n+1) = 0.54 - 0.46*cos(pi*(n+1-N)/Nw);
  else
    win(n+1) = 1 ;
  end
end


% Generamos ruido aleatorio
noise = amp_ruido*(2*rand(1,N)-1);
%noise = noise.*win; % ventana

% Sumamos senal y ruido
data = signal + noise;
% data=signal; % senal sin ruido


% =================================================================
% CONSTRUCCION DEL FILTRO PREDICTIVO
% =================================================================

  acor = xcorr(data); % correlacion de la senal == conv(fliplr(data),(data))
  a = acor(N:N+Nf-1);      % primeros Nf lags de la autocorrelacion

%  correlacion entre  la senal y la  salida deseada. En este  caso por
%  ser predictor este  vector son los lags de  la autocorrelacion de
%  la señal desplazados una unidad.
    
  c(:,1) = acor((N)+1:(N+Nf-1)+1); % (:,1) : obligo a que sea vector columna

  % filtro predictor:
  
    % el sistema tiene estructura Toeplitz,
    
     pf  = inv(toeplitz(a))*c;
     pef = [ 1 ; -pf ];       %  prediction error filter (pef)
   
   % Como el sistema tiene estructura Toeplitz, puede usarse el algoritmo
   % de Levinson  para resolver.   La subrutina "levinson"  obtiene el
   % "prediction error filter"  (pef) por lo cual hay  que arreglar la
   % salida para obtener el filtro predictor (pf).
   %
   %   pef  = levinson(acor(N:N+Nf),Nf+1); % pef
   %   pf  = -pef(2:end)'; % pf
   %
   % OBV: la unica informacion utilizada son los lags de la
   % autocorrelacion del dato de entrada.

% =================================================================
% FILTRADO DE LA SENAL
% =================================================================


% Filtramos la senal con el operador pf
signal2 = conv(data,pf);
%signal2 = conv(data,pf)(1:N); % sin efectos de borde

signal2  = [  0 signal2(1:N-1)  ]; %  considero x0,  para alinear
                                   %  las prediciones  a   partir
                                   %  de  x1  en adelante  respecto
                                   %  del vector original "signal"


% Calculamos el error que predijo el filtro
noise2 = conv(data,pef);
noise2 = noise2(1:N);
noise2 = noise2.*win;


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

figure(1); clf;

ymin=1.2*min(data);
ymax=1.2*max(data);

subplot(3,2,1)
plot((0:N-1),signal,'k')
axis([ 0 N ymin ymax ])
title('Senal')

subplot(3,2,3)
plot((0:N-1),noise,'r')
axis([ 0 N ymin ymax])
title('Ruido')

subplot(3,2,5)
plot((0:N-1),data,'b')
axis([ 0 N ymin ymax ])
title('Senal+Ruido')

subplot(3,2,2)
plot((0:N-1),signal,'k',(0:N-1),signal2,'g:')
axis([ 0 N ymin ymax])
title('Senal estimada con pf')

subplot(3,2,4)
plot(1:N,noise,'r',1:N,noise2,'g:')
axis([ 0 N ymin ymax])
title('Ruido estimado con  pef')

figure(2)
subplot(2,1,1)
plot((0:N-1),signal,'k',(0:N-1),data,'b')
axis([ 0 N ymin ymax ])
title('Senal+Ruido')

subplot(2,1,2)
plot((0:N-1),signal,'k',(0:N-1),signal2,'g:')
axis([ 0 N ymin ymax])
title('Senal estimada con filtro predictivo')

% ==================================================================
% Salida
% ==================================================================

% print('filtro_predictivo.eps','-depsc')
