%=================================================================
% Aplicacion de  promedios moviles sobre senales  sinusoidales
% simples.
%
% Se busca ver el efecto que  tiene el filtro sobre la senal compuesta
% por tonos (Figura 1) y el efecto que tiene el filtro sobre los tonos
% directamente.   Con esto  ultimo se  puede interpretar  la respuesta
% impulsiva del filtro (Figura 2)

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

clear all ; close all % Limpio variables


% =================================================================
% Filtro de promedio movil
% =================================================================

   h = (1/4.0)*[ 1 2 1]; % promedio pesado

% h = (1/3.0)*[ 1 1 1]; % promedio simple

  disp('filtro ='),disp(h)

% =================================================================
%  Parametros y senales de interes
% =================================================================

fM = 50;           % Frecuencia de muestreo, [Hz].
fN = 0.5*fM;        % Frecuencia de Nyquist (A partir de aqui, ALIASING)
dt = 1.0/fM;        % Intervalo de muestreo,   [s].

f1 = 1 ; amp1 = 1.0;  % Frecuencia [Hz] y amplitud de la componente s1

% --------------------------------------------
% Casos  interesantes   para  el   promedio  simple:
%
% a) w=2pi/3 =>  f2=1/(3*dt). Aqui la Rta en frecuencia  de h tiene un
% notch. Esta componente es anulada.
%
% b) w=5pi/6 (a la  mitad entre 2pi/3 y pi). Hay cambio  de fase en pi
% en  la zona  de  reyeccion  y la  amplitud  no  es despreciable.
%

%  f2=1/(3*dt); amp2=1.0; %  caso a) Frecuencia [Hz]  y amplitud de
                         %  la componente s2.
  
  f2=5/(12*dt); amp2=1.0; %  caso b) Frecuencia [Hz] y  amplitud de la
%                            componente  s2. Nota:  Hacer  zoom en  el
%                            grafico correspondiente  para apreciar el
%                            cambio de fase en pi.

% --------------------------------------------

t  = (0:dt:5*1/f1); % Vector de tiempos asociado a dt, [s]
N  = length(t);

s1 = sin(2*pi*f1*t);        % componente s1
s2 = sin(2*pi*f2*t);        % componente s2
s  = amp1 * s1 + amp2 * s2; % senal de interes s = s1 + s2


% =================================================================
% Aplicacion del filtro sobre las senales
% =================================================================


s1f = conv(h,s1)(2:N+1); % senal s1 filtrada, sin efectos de borde
s2f = conv(h,s2)(2:N+1); % senal s1 filtrada, sin efectos de borde
sf  = conv(h,s)(2:N+1);  % senal s filtrada, sin efectos de borde


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

clf;                  % Limpio figuras

figure(1)

subplot(2,1,1)
hold on
plot(t,s,'r;senal = s1+s2;')
plot(t,s1,'k')
hold off
axis([ 0 t(N) -2 +2])
ylabel('amplitud')
subplot(2,1,2)
hold on
plot(t,sf,'g;senal filtrada;')
plot(t,s1,'k')
xlabel('tiempo (s)')         % agrego etiqueta al eje x
ylabel('amplitud')       % agrega etiqueta al eje y
axis([ 0 t(N) -2 +2])
hold off

% Aqui puede  interpretarse el  efecto que tiene  el filtro  sobre las
% tonos  que componen  la senal.   Se  esta viendo  aqui la  respuesta
% impulsiva actuando  sobre dos  sinusoides. El  efecto es  alterar la
% amplitud y la  fase (si el operador tiene fase  distinta de cero) de
% los tonos de entrada s1 y s2.
%
% Nota: se grafica  s1 y s2 con amplitud unitaria  para poder comparar
% el efecto del filtro sobre la amplitud entre las figuras.


figure(2)

subplot(2,2,1)
 plot(t,s1,'r;s1;')
 axis([ 0 t(N) -2 +2])
 ylabel('amplitud')
subplot(2,2,2)
 plot(t,s1f,'g;s1 filtrada;')
 axis([ 0 t(N) -2 +2])
subplot(2,2,3)
plot(t,s2,'r;s2;')
 axis([ 0 t(N) -2 +2])
 xlabel('tiempo (s)')
 ylabel('amplitud')  
subplot(2,2,4)
 plot(t,s2f,'g;s2 filtrada;')
 axis([ 0 t(N) -2 +2])
 xlabel('tiempo (s)')


% descomentar la siguiente linea para guardar grafico EPS)
% print(strcat('filtro_prom.eps'),'-deps')
