% ==================================================================
% INTERPOLACION POR ZERO INTERLACING
%
% OBS: El directorio de trabajo debe contener el archivo wavelet
%
% NOTA: ver zero-padding en frecuencia como alternativa
% ==================================================================

clear ; clf ;

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

L  = 2;  % relación  de interpolación  (L-1 ceros  entre las  muestras
       % originales)

% ===================================================================
% Cargo la secuencia y la normalizo
% ===================================================================

h = load('wavelet');
h = h/max(abs(h));

N  = length(h) ; % Número de muestras de la secuencia original
Ni = L*N ;       % Número de muestras de la secuencia interpolada

% ===================================================================
% Calculo la TDF:
% ===================================================================

H  = fft(h) ;
w0 = 2*pi/N ;

% ===================================================================
% Intercalo ceros en la secuencia original en lo que serán 
% las muestras interpoladas (zero interlacing)
% ===================================================================

h0 = zeros(Ni,1) ; % genero vector de Ni ceros
for i = 1:N
  h0(L*(i-1)+1) = h(i) ;
end

% ===================================================================
% Calculo la TDF la nueva secuencia con ceros intercalados
% ===================================================================
% OBS: al agregar ceros duplico el espectro

H0 = fft(h0) ;

% ===================================================================
% Genero la secuencia interpolada en el dominio de las frecuencias,
% teniendo en cuenta el resultado del ejercicio anterior,
% hago cero las frecuencias entre pi/L y 2*pi-pi/L.
%
% wn  = pi   (Frecuencia angular digital de Nyquist)
% wn' = wn/L 
% ===================================================================

% anulo las frecuencias entre pi/L y 2pi-pi/L

w0i = 2*pi/Ni ;
for i = 1:Ni
  w = (i-1)*w0i ;
  if w < (pi/L)
    Hi(i) = H0(i) ;
  elseif w > (2*pi-pi/L)
    Hi(i) = H0(i) ;
  else  
    Hi(i) = 0 ;
  endif
end

% Conservo las amplitudes originales. la TDf inversa  va a dividir
% por NL.

Hi =  L * Hi; %  Ahora la  TDF inversa al dividir por NL resulta en el
              %  factor 1/N, lo mismo que para la secuencia original.

% ===================================================================
% Calculo la secuencia interpolada como la TDF inversa de Hi
% ===================================================================

 hi = real(ifft(Hi)) ;

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

% Graficos en tiempo

dt1=L;
dt2=1;

subplot(3,2,1);
stem(0:dt1:Ni-1,h,'b')
title('En tiempo')
legend('Secuencia original')
axis([0 Ni -1.4 1.4])

subplot(3,2,3);
stem(0:dt2:Ni-1,h0,'g')
legend('Entrelazado de ceros')
axis([0 Ni -1.4 1.4])

subplot(3,2,5);
stem(0:dt2:Ni-1,hi,'r')
legend('Secuencia interpolada')
axis([0 Ni -1.4 1.4])

% Grafico en frecuencia angular digital

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

subplot(3,2,2);
plot((0:N-1)*w0,abs(H),'b;Secuencia original;')
title('En frecuencia')
set(gca,"xtick",xt,"xticklabel",xtlabel);
axis([0,2*pi,-1 10]);

subplot(3,2,4);
plot((0:Ni-1)*w0i,abs(H0),'g;Entrelazado de ceros;')
set(gca,"xtick",xt,"xticklabel",xtlabel);
axis([0,2*pi,-1 10]);

% Observe  que  este resultado  es  identico  a haber  realizado  zero
% padding en frecuencia.

subplot(3,2,6);
plot((0:Ni-1)*w0i,abs(Hi),'r;Secuencia interpolada;') ;
set(gca,"xtick",xt,"xticklabel",xtlabel);
axis([0,2*pi,-1 10]);

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

% Guardo el gráfico en un postscript
%print('zero_interl.eps','-depsc')
