% ================================================================
% DECIMACION Y REINTERPOLACION DE UNA ONDICULA DE RICKER.
%
% A medida que se aumenta la frecuencia central de la ondicula, el dato
% decimado sufre de aliasing en frecuencia. Una manera de disminuir los
% efectos de aliasing consiste en aplicar un filtro pasabajos como paso
% previo a la decimacion.
%
% Una vez  obtenida la senal  DECIMADA, se emplea  la interpolacion
% sinc para volver  al intervalo de muestreo original.  Al hacer esto,
% los  efectos  del   aliasing  sobre  la  senal   decimada  se  hacen
% evidentes. Esto  sirve para  comprender el ALIASING  y el  metodo de
% INTERPOLACION SINC, que recupera la  version de banda limitada de la
% senal de entrada, no pudiendo generar nueva informacion mas alla del
% Nyquist  original.  Es decir,  si  bien  los puntos  decimados  caen
% exactamente sobre la ondicula de Ricker, si hay efectos de aliasing,
% al  reinterpolar se  obtiene  una senal  claramente  diferente a  la
% ondicula.  Esta senal  necesariamente  tiene un  contenido menor  de
% frecuencias, como puede observarse.
%
% OBS: el directorio debe contener la subrutina ricker.m
% ================================================================


clear all; close;

%  Ondícula de Ricker  muestreada a 2 ms

leng = 101  ; % (ms) Longitud en tiempo
dt   = 2    ; % (ms) Intervalo de muestreo en tiempo
f0   = 40   ; % (Hz) Frecuencia central

% Llamo a la subrutina ricker y defino la ondícula
[s,Ns] = ricker(leng,dt,f0);
Ms     = (Ns - 1)/2;
t      = (-Ms:Ms)*dt;


% =====================================
% Decimacion tomando muestras cada k puntos de la senal original
% =====================================


% Si k=2 el aliasing es despreciable  y puedo recuperar la misma senal
% de entrada al  deshacer la decimacion por medio  de la interpolacion
% sinc.

k = 8;  

sd = s(1:k:Ns); % se toma muestras cada k puntos de la senal original

Nd = length(sd);


% =====================================
% Interpolacion de la senal decimada
% =====================================

% Puede obtenerse la senal original (s)  de la decimada (sd) siempre y
% cuando el efecto del aliasing sea despreciable al pasar de s a sd.
    
si  = sinc_interp(sd,k-1);

Ni  = length(si)  % numero de puntos interpolados


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

N2 = 1024; % zero padding

% Espectro de la ondícula con zero padding
s2     = s;
s2(N2) = 0;
S      = fft(s2);

% Espectro de la ondícula decimada con zero padding
sd2     = sd;
sd2(N2) = 0;
SD      = k*fft(sd2);

dw = 2*pi/N2;

subplot(3,1,1);
hold on
plot((0:Ns-1),s,'k')
plot((0:Ns-1),s,'bo;Original;',"markersize",4)
plot((0:Nd-1)*k,sd,'ro;Decimada;');
hold off     
xlabel('t [mseg]');
ylabel('h_n')
axis([ 0 Ns-1 1.2*min(s) 1.2*max(s) ]);

subplot(3,1,2);
hold on
plot((0:N2-1).*dw,abs(S),'b-;Original;')
plot((0:N2-1).*dw,abs(SD),'r-;Decimada;')
hold off     
xlabel('{/Symbol w}');
ylabel('|H_k({/Symbol w})|')
xt      = 0:pi/2:2*pi;
xtlabel = ["0";"{/Symbol p}/2";"{/Symbol p}";"3/2{/Symbol p}";"2{/Symbol p}"];
yt      = [0 1];
set(gca,"xtick",xt,"xticklabel",xtlabel);
axis([0,2*pi,-0.5,5+max(abs(S))]);

subplot(3,1,3);
hold on
plot((0:Ns-1),s,'k');
plot((0:Ni-1),si,'go;Interpolado;',"markersize",4);
plot((0:Nd-1)*k,sd,'ro');
hold off
xlabel('t [mseg]');
ylabel('h_n')
axis([ 0 Ns-1 1.2*min(s) 1.2*max(s) ]);



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

print('decimacion_interpolacion.eps','-depsc');


