%% =================================================================
%% MÉTODO DE KOLMOGOROFF
%%
%% !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 por Kolmogoroff
% --------------------------------------------------------

% Calculamos la TF de la ondícula de entrada con zero padding
Npad     = 1024 ;
rw(Npad) = 0 ;
H        = fft(rw) ;

% Calculamos el espectro de potencia, le tomamos el logaritmo,
% y la calculamos su TF inversa
u = ifft(log(H.*conj(H))) ;

% Anulamos los tiempos negativos y dividimos por 2 la componente 
% continua (u0) para generar el u+
u((Npad/2+2):Npad) = 0 ;      % hacemos cero los negativos
u(1)               = u(1)/2 ; % este es el u_0/2

% Calculamos la TF de la parte causal (del u+)
U = fft(u) ;

% Calculamos la exponencial y la TF inversa
minrw = real(ifft(exp(U))) ;

% 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(1:N)); hold on
plot((0:N-1)*dt,rw(1:N)); 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(1:N),'r'); hold on
plot((0:N-1)*dt,minrw(1:N),'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('kolmogoroff_a.eps','-depsc')

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

figure(2); clf;

% Obs: ya se hizo el zer padding para el método
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('kolmogoroff_b.eps','-depsc');


