% ============================================================
% INTERPOLACION LINEAL
%
% IN     h:   secuencia de entrada
%        dt1: intervalo de muestreo de entrada
%        dt2: intervalo de muestreo de salida
%
% OUT    hi:  secuencia interpolada
%
% Obs: El directorio  de trabajo debe contener el  archivo "wavelet" y
% la funcion "lin_interp.m".

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

clear all ;
clear all ;
close all ;
clf ;

% ============================================================
% Defino los intervalos de muestreo de entrada y salida
% ============================================================

dt1 = 2 ; % dt1 > dt2
dt2 = 1 ;
if dt2 > dt1
  printf('ERROR: Para interpolar dt2 debe ser menor que dt1')
  break
end

% ============================================================
% Relación de interpolación
% ============================================================

k = dt1 / dt2; 

% ============================================================
% Cargo secuencia y  normalizo:
% ============================================================

h = load("wavelet");
h = h/max(abs(h));   
N = length(h) ; % Número de muestras de la secuencia original

% ============================================================
% Secuencia interpolada linealmente hi
% ============================================================

 hi = lin_interp(h,k-1);
 
 Ni=length(hi);
 
% ============================================================
% Espectro de la secuencia original (con zero padding)
% ============================================================

% Zero padding
h2     = h ;
NN     = 1024 ;
h2(NN) =    0 ;

% Espectro original
H = abs(fft(h2)) ;

% Observar  los  lobulos secundarios  que  aparecen  en Hi  debido  al
% aliasing. Esto  es debido a que  el operador lineal no  es nulo para
% frecuencias por encima del Nyquist original. Esto conduce a
% llevar dentro del espectro interpolado la parte de amplitudes
% originales por encima de Nyquist.

% Repito el espectro k veces para graficar 
Hrep = H;
for i = 2:k
  Hrep =[Hrep ; H]; 
end

% ============================================================
% Espectro de la secuencia interpolada (con zero padding)
% ============================================================

% Zero padding
h2i     = hi ;
h2i(NN) = 0;

% Espectro interpolado
Hi = abs(fft(h2i)) ;

% ============================================================
% GRÁFICOS
% ============================================================

% ------
% TIEMPO
% ------

subplot(3,1,1) ;


hold on
plot((0:N-1)*dt1,h,'bo;Original;',"markersize",3);
plot((0:Ni-1)*dt2,hi,'ro;Interpolada;',"markersize",5);
stem((0:Ni-1)*dt2,hi,'color','r',"markersize",5);
stem((0:N-1)*dt1,h,'color','b',"markersize",4,"filled");
hold off
xlabel('t [mseg]');
ylabel('h_n')

% ------------------
% FRECUENCIA ANGULAR
% ------------------

% Intervalos de muestreo en frecuencia angular (en Khz)

Ws   = 2*pi/dt1;  % frecuencia de muestreo de secuencia
dW0  = Ws/NN;     % para la secuencia original repetida
Wsi  = 2*pi/dt2;  % frecuencia de muestreo de secuencia interpolada
dW0i = Wsi/NN;    % para la secuencia interpolada

% Ticks para el eje de frecuencia angular

xt      = 0:Wsi/4:Wsi;
xtlabel = ["0";"{/Symbol W}_s/2";"{/Symbol W}_s";
           "3/2{/Symbol W}_s";"2{/Symbol W}_s"];

subplot(3,1,2) ;
hold on
plot((0:k*NN-1)*dW0,Hrep,';Original;b.-')
plot((0:NN-1)*dW0i,Hi,';Interpolada;r-') ;
hold off

xlabel('{/Symbol W}');
ylabel('|H_k({/Symbol W})|')
set(gca,"xtick",xt,"xticklabel",xtlabel);
axis([0,Wsi,0,k*max(abs(H))]);

% --------------------------
% FRECUENCIA ANGULAR DIGITAL
% --------------------------

% Intervalos de muestreo en frecuencia angular digital

dw0  = 2*pi/NN ; % para ambas secuencias con ZERO PADDING

% Ticks para el eje de frecuencia angular digital

xt      = 0:pi/2:2*pi;
xtlabel = ["0";"{/Symbol p}/2";"{/Symbol p}";"3/2{/Symbol p}";"{/Symbol 2p}"];

subplot(3,1,3) ;
hold on
plot((0:NN-1)*dw0,H,';Original;b.-')
plot((0:NN-1)*dw0,Hi,';Interpolada;r-') ;
hold off
xlabel('{/Symbol w}');
ylabel('|H_k({/Symbol w})|')
set(gca,"xtick",xt,"xticklabel",xtlabel);
axis([0,2*pi,0,k*max(abs(H))]);

% ============================================================
% SALIDA
% ============================================================

% print('sinc_interp.eps','-depsc');
