%% =================================================================
%% OPERADOR DE PROMEDIO PESADO APLICADO A ONDICULA DE RICKER 
%% CONTAMINADA CON RUIDO ALEATORIO
%%
%% !OBS: el directorio debe contener la subrutina ricker.m
%% =================================================================

clear all;
close all;

% =================================================================
% (1) Ondicula de Ricker

% Parametros de la ondicula de Ricker

leng = 100 ;  % (ms) Longitud en tiempo   (ms == 1E-03 s)
dt     = 1   ;   %  (ms)  Intervalo de muestreo en tiempo
f0     = 50  ;  %  (Hz)  Frecuencia "central"

rw = ricker(leng,dt,f0);
N = length(rw);        % N es impar por cómo es la función ricker()
M = (N - 1)/2;

% =================================================================
% (2) Ruido aleatorio 

amp_noise = 1./5.;   % amplitud del ruido / amplitud de la ondicula

noise = amp_noise * ( 2*rand(1,N) - 1 );

% =================================================================
% (3) Ondicula con ruido 

s = rw + noise;

% =================================================================
% (4) Filtro 

h   = [0.25 0.5 0.25];  % operador promedio pesado
%                 ^   
%                t=0

sd = conv(s, h); % aplicacion del filtro, dim = dim(h) + dim(rw2) - 1
sd = sd(2:N+1);  % quito los efectos de borde, dim = dim(rw2)

% =================================================================
% Dominio transformado
% ================================================================

% ================================================================
% (1b) Espectro de amplitud de la Ondicula de Ricker (normalizado)

RW  = abs(fft(rw));
maxrw = max(RW);
RW  = RW/maxrw;

% =================================================================
% (2b) Espectro de amplitud del ruido (normalizado)

NOISE = abs(fft(noise));
NOISE = NOISE/maxrw;

% =================================================================
% (3b) Espectro de amplitud de la ondicula con ruido (normalizado)

S = abs(fft(s));
S = S/maxrw;

% =================================================================
% (4b) Espectro de amplitud de la ondicula filtrada  (normalizado)

SD = abs(fft(sd));
SD = SD/maxrw;

% =================================================================
% GRAFICOS EN TIEMPO 
% =================================================================

figure(1); clf;

% =================================================================
% Ondicula
  
subplot(2,2,1)
plot((-M:M).*dt,rw,'g')
title ('Ondicula de Ricker')
axis([-M,M,-0.6,1.1])

% =================================================================
% Ruido
  
subplot(2,2,3)
plot((-M:M).*dt,noise,'r')
title ('Ruido aleatorio')
axis([-M,M,-0.6,1.1])

% =================================================================
% Ondicula + ruido
  
subplot(2,2,2)
plot((-M:M).*dt,s)
title ('Ondicula con ruido')
axis([-M,M,-0.6,1.1])
  
% =================================================================
% Ondicula + ruido filtrada
  
subplot(2,2,4)
plot((-M:M).*dt,sd,'g')
title ('Ondicula filtrada')
axis([-M,M,-0.6,1.1])

% =================================================================
% Guardo un .eps con el gráfico generado:

print('promedio_aplicado_a.eps','-deps');


% =================================================================
% GRAFICOS EN FRECUENCIA 
% =================================================================

figure(2); clf;

dw = 2*pi/N;

% =================================================================
% |H(w)| de la ondicula

subplot(2,2,1);
plot((0:N-1)*dw,RW,'g-')
title ('{Espectro de amplitud de ondicula}')
axis([0,2*pi,-0.1,1.1])
xt=0:pi/2:2*pi;
xtlabel=["0";"{/Symbol p}/2";"{/Symbol p}";"3/2{/Symbol p}";"{/Symbol 2p}"];
set(gca,"xtick",xt,"xticklabel",xtlabel);

% =================================================================
% |H(w)| del ruido

subplot(2,2,3);
plot((0:N-1)*dw,NOISE,'r-')
title ('Espectro de amplitud de ruido')
axis([0,2*pi,-0.1,1.1])
xt=0:pi/2:2*pi;
xtlabel=["0";"{/Symbol p}/2";"{/Symbol p}";"3/2{/Symbol p}";"{/Symbol 2p}"];
set(gca,"xtick",xt,"xticklabel",xtlabel);

% =================================================================
% |H(w)| de la ondicula con ruido

subplot(2,2,2);
plot((0:N-1)*dw,S,'-')
% grafico del espetro de amplitud del filtro
%hold on
%h=[h zeros(1,length(RW1)-length(h))];
%plot((0:N-1)*dw,abs(fft(h)),"5")
title ('Espectro de amp. de odicula con ruido')
axis([0,2*pi,-0.1,1.1])
xt=0:pi/2:2*pi;
xtlabel=["0";"{/Symbol p}/2";"{/Symbol p}";"3/2{/Symbol p}";"{/Symbol 2p}"];
set(gca,"xtick",xt,"xticklabel",xtlabel);

% =================================================================
% |H(w)| de la ondicula filtrada

subplot(2,2,4);
plot((0:N-1)*dw,SD,'g-')
title ('Espectro de amp. de ondicula filtrada')
axis([0,2*pi,-0.1,1.1])
xt=0:pi/2:2*pi;
xtlabel=["0";"{/Symbol p}/2";"{/Symbol p}";"3/2{/Symbol p}";"{/Symbol 2p}"];
set(gca,"xtick",xt,"xticklabel",xtlabel);


% =================================================================
% Guardo un .eps con el gráfico generado:

print('promedio_aplicado_b.eps','-deps');

% =============================== X ===============================