%% ================================================================
%% INTERPRETANDO EL ESPECTRO DE AMPLITUD Y SUS ARM�NICOS
%%
%% Dado el espectro de una ond�cula, se va generando la ond�cula en 
%% tiempo tomando paulatinamente los valores de amplitud del 
%% espectro ("arm�nicos). Observamos as� c�mo la sumas parciales 
%% convergen a la ond�cula a medida que aumenta la cantidad de 
%% puntos tomados en el dominio de la frecuencia.
%%
%% !OBS: el directorio debe contener la subrutina ricker.m
%% ================================================================

clear; clf; close;

% =================================================================
% Ond�cula de Ricker
% =================================================================

leng  = 50  ; % (ms) Longitud en tiempo
dt    = 1   ; % (ms) Intervalo de muestreo en tiempo
f0    = 80  ; % (Hz) Frecuencia central

% Llamo a la subrutina ricker y defino la ond�cula
s = ricker(leng,dt,f0);
N = length(s);
M     = (N - 1)/2;
t     = (-M:M)*dt;

% =================================================================
% Espectro de amplitud
% =================================================================

S  = abs(fft(s)) ;
dw = 2*pi/N      ; % paso en frencuencia

% =================================================================
% Picando valores de la TDF
% =================================================================

% Comenzamos con un vector de ceros. 
% Es decir, todas las amplitudes son nulas, al sumar las funciones 
% de base obtendriamos la funcion s_n = 0, para todo n.

SP = zeros(1,N); 

% Se van tomando paulatinamente los valores dictados por el espectro.
% De esta forma a medida que se van tomando cada vez mas valores de la
% transformada, en tiempo estaremos generando una se�al cada vez mas
% parecida a la ondicula original.

for i = 1:13 % bucle para ir agregando para arm�nicos

  SP(i) = S(i);
  if i > 1
    SP(N-(i-2)) = S(i); % herm�tico para que sp sea real.
  end

  % ===============================================================
  % Arm�nico asociado al valor del espectro tomando
  % ===============================================================

  w = dw * (i-1);               % frecuencia del arm�nico
  a = S(i) * cos(w*t) / max(S); % arm�nico

  % ===============================================================
  % Suma parcial en tiempo
  % ===============================================================

  sp = real(ifft(SP)); 
  sp = [sp(M+2:N) sp(1:M+1)]; 

  % ===============================================================
  % GR�FICOS 
  % ===============================================================

  figure(1)

  % Ond�cula
  
  subplot(5,1,1)
  plot((-M:M)*dt,s,'go')
  axis([-M,M,-0.6,1.1])
  ylabel('s_n')
  title('Interpretando el Espectro de Amplitud','fontsize',12)

  % |H(w)| de la ond�cula

  subplot(5,1,2);
  stem((0:N-1)*dw,S/max(S),'g-')
  ylabel('|S_k|')
  axis([0,2*pi,-0.1,1.1])
  xt      = 0:pi/2:2*pi;
  xtlabel = ["0";"{/Symbol p}/2";"{/Symbol p}";
             "3/2{/Symbol p}";"{/Symbol 2p}"];
  set(gca,"xtick",xt,"xticklabel",xtlabel);

  % |H(w)| de la suma parcial

  subplot(5,1,3);
  stem((0:N-1)*dw,SP/max(S),'r-')
  ylabel('|SP_k|')
  axis([0,2*pi,-0.1,1.1])
  xtlabel = ["0";"{/Symbol p}/2";"{/Symbol p}";
             "3/2{/Symbol p}";"{/Symbol 2p}"];
  set(gca,"xtick",xt,"xticklabel",xtlabel);

  % Arm�nico

  subplot(5,1,4)
  plot((-M:M)*dt,a,'b')
  axis([-M,M,-1.1,1.1])
  ylabel('a_n')

  % Suma parcial en tiempo

  subplot(5,1,5)
  plot((-M:M)*dt,sp,'ro')
  axis([-M,M,-0.6,1.1])
  ylabel('sp_n')

  printf('Click para continuar tomando puntos del Espectro ... \n');
  pause

endfor % finalizo bucle para los arm�nicos

% Arm�nicos

a0  = S(1)  * cos(0*dw*t)  / max(S);
a1  = S(2)  * cos(1*dw*t)  / max(S);
a2  = S(3)  * cos(2*dw*t)  / max(S);
a3  = S(4)  * cos(3*dw*t)  / max(S);
a4  = S(5)  * cos(4*dw*t)  / max(S);
a5  = S(6)  * cos(5*dw*t)  / max(S);
a6  = S(7)  * cos(6*dw*t)  / max(S);
a7  = S(8)  * cos(7*dw*t)  / max(S);
a8  = S(9)  * cos(8*dw*t)  / max(S);
a9  = S(10) * cos(9*dw*t)  / max(S);
a10 = S(11) * cos(10*dw*t) / max(S);
a11 = S(12) * cos(11*dw*t) / max(S);
a12 = S(13) * cos(12*dw*t) / max(S);

suma = a0+a1+a2+a3+a4+a5+a6+a7+a8+a9+a10+a11+a12;

% Grafico arm�nicos

figure(2)

subplot(14,1,1)
plot((-M:M)*dt,a0,'b')
title('Los armonicos y su suma','fontsize',12)
ylabel('a_{0}')
axis([-M,M,-1,1])
axis("off")

subplot(14,1,2)
plot((-M:M)*dt,a1,'b')
ylabel('a_{1}')
axis([-M,M,-1,1])
axis("off")

subplot(14,1,3)
plot((-M:M)*dt,a2,'b')
ylabel('a_{2}')
axis([-M,M,-1,1])
axis("off")

subplot(14,1,4)
plot((-M:M)*dt,a3,'b')
ylabel('a_{3}')
axis([-M,M,-1,1])
axis("off")

subplot(14,1,5)
plot((-M:M)*dt,a4,'b')
ylabel('a_{4}')
axis([-M,M,-1,1])
axis("off")
  
subplot(14,1,6)
plot((-M:M)*dt,a5,'b')
ylabel('a_{5}')
axis([-M,M,-1,1])
axis("off")

subplot(14,1,7)
plot((-M:M)*dt,a6,'b')
ylabel('a_{6}')
axis([-M,M,-1.1,1.1])
axis("off")

subplot(14,1,8)
plot((-M:M)*dt,a7,'b')
ylabel('a_{7}')
axis([-M,M,-1,1])
axis("off")

subplot(14,1,9)
plot((-M:M)*dt,a8,'b')
ylabel('a_{8}')
axis([-M,M,-1,1])
axis("off")

subplot(14,1,10)
plot((-M:M)*dt,a9,'b')
ylabel('a_{9}')
axis([-M,M,-1,1])
axis("off")

subplot(14,1,11)
plot((-M:M)*dt,a10,'b')
ylabel('a_{10}')
axis([-M,M,-1,1])
axis("off")

subplot(14,1,12)
plot((-M:M)*dt,a11,'b')
ylabel('a_{11}')
axis([-M,M,-1,1])
axis("off")

subplot(14,1,13)
plot((-M:M)*dt,a12,'b')
ylabel('a_{12}')
axis([-M,M,-1,1])
axis("off")

subplot(14,1,14)
plot((-M:M)*dt,suma,'g')
ylabel('suma')
axis([-M,M,min(suma),max(suma)])
axis("off")

% =============================== X ===============================
