%% =================================================================
%% MÉTODO DE LAS RAÍCES
%%
%% !OBS: el directorio debe contener la subrutina ricker.m
%%       y la ondícula "wavelet".
%% =================================================================

clear all; close;

% ------------------------
% Elijo el dato de entrada
% ------------------------

dato = 1;   % Ondícula wavelet
%dato = 2;   % Ondícula de Ricker

if (dato == 1)
 
  load('wavelet');
  rw = wavelet'/max(wavelet);
  dt = 1;
  N  = length(rw); 

elseif (dato == 2)

  % Parámetros de la ondícula de Ricker
  leng  = 56 ; % (ms) Longitud en tiempo
  dt    = 1  ; % (ms) Intervalo de muestreo en tiempo
  f0    = 30 ; % (Hz) Frecuencia central

  % Llamo a la subrutina ricker y defino la ondícula
  rw = ricker(leng,dt,f0);
  N  = length(rw);

end

% ----------------------------------------
% Calculamos la equivalente de fase mínima
% ----------------------------------------

% Calculamos la autocorrelacion de la ondícula
xrw = xcorr(rw);

% Calculamos las raíces de la autocorrelación
root = roots(xrw);

% Generamos una ondícula de fase mínima con las raíces que 
% estén fuera del círculo unidad
aux(1) = 1;
for i = 1:length(root)
  if abs(root(i)) > 1
    aux = conv(aux,[1,(-1/root(i))]);
  end
end

% Esto porque algunas veces agrega una raíz de más o de menos
minrw      = zeros(1,N);
minrw(1:N) = aux(1:N);

% Escalamos la ondícula de fase mínima para que tenga 
% la misma media que la ondícula de fase cero original
k     = sign(rw(1))*mean(abs(rw))/(sign(minrw(1))*mean(abs(minrw)));
minrw = k*real(minrw);

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

% -------------------
% Ondículas en tiempo
% -------------------

figure(1); clf;

subplot(2,1,1)
stem((0:N-1)*dt,rw); hold on
plot((0:N-1)*dt,rw); hold off
axis([ 0 (N-1)*dt -1.2 1.2 ])
title('Ondicula de entrada')

subplot(2,1,2)
stem((0:N-1)*dt,minrw,'r'); hold on
plot((0:N-1)*dt,minrw,'r'); hold off
axis([0 (N-1)*dt -1.2 1.2])
title('Ondicula equivalente de fase minima')

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

% --------------------------
% Espectros con zero padding
% --------------------------

figure(2); clf;

% Zero padding
Npad        = 1024;
rw(Npad)    = 0;
minrw(Npad) = 0;
dw          = 2*pi/Npad;

% Espectros
RW    = abs(fft(rw));
MINRW = abs(fft(minrw));

% Ticks en frecuencia angular digital
xt      =  0:pi/2:2*pi;
xtlabel = ["0";"{/Symbol p}/2";"{/Symbol p}";"3/2{/Symbol p}";"{/Symbol 2p}"];

subplot(2,1,1)
plot((0:Npad-1)*dw,RW)
title('Espectro de la ondicula de entrada')
set(gca,"xtick",xt,"xticklabel",xtlabel);
axis([0 2*pi])

subplot(2,1,2)
plot((0:Npad-1)*dw,MINRW,'r')
title('Espectro de la ondicula equivalente de fase minima')
set(gca,"xtick",xt,"xticklabel",xtlabel);
axis([0 2*pi])

% Guardo el grafico en un postscript
print('metodo_de_raices_b.eps','-depsc');
