% SPECTRUM.M

% CALCULA ESPECTOS DE SERIES TEMPORALES DISCRETAS DE ACUERDO LA MÉTODO
% DE BLACKMAN-TUKEY (LIBRO DE OTNES, PAG: 270-272), CON EL CORRESPONDIENTE
% CONTINUO NULO Y BANDAS DE SIGNIFICANCIA AL 5% O AL 10% (CLIMATIC CHANGE,
% NOTA TÉCNICA Nº 79, PAG: 36-42).
% LA MATRIZ DE ENTRADA DE DATOS DEBE PRESENTAR LAS SERIES A LAS QUE SE LES
% VA A CALCULAR EL ESPECTRO EN COLUMNA. TODAS LAS SERIES DEBEN TENER EL 
% MISMO NÚMERO DE REGISTROS (LA MISMA LONGITUD).
% AL COMENZAR EL PROGRAMA REMOVERÁ LAS MEDIAS DE CADA SERIE Y APLICARÁ UNA
% VENTANA COSENO EN LOS EXTREMOS. LUEGO APLICARÁ UNA VENTANA HANN, HAMMING
% O PARZEN A ELECCIÓN (SE ELIGE DESDE PANTALLA).
% LUEGO GRAFICARÁ LOS ESPECTROS DE CADA SERIE Y GUARDARÁ LOS DATOS EN UNA
% MATRIZ CON EL NOMBRE "Resul" CUYAS COLUMNAS SON:
% col1: frecuencias
% col2: espectro de la serie 1
% col3: limite inferior del intervalo de confianza (serie 1)
% col4: continuo nulo (serie 1)
% col5: limite superior del intervalo de confianza (serie 1)
% col6: espectro de la serie 2
% col7: limite inferior del intervalo de confianza (serie 2)
% col8: continuo nulo (serie 2)
% col9: limite superior del intervalo de confianza (serie 2)
% Y ASI SIGUIENDO CON TODAS LAS SERIES PRESENTES EN LA MATRIZ DE ENTRADA.
% LA MATRIZ "Resul" PUEDE LUEGO EXPORTARSE PARA GRAFICAR LOS ESPECTROS EN
% ORIGIN Y ALLI TRABAJARLOS.
%--------------------------------------------------------------------------


% Entrada de datos

Matriz=xlsread('rain9506.xls');
m=input('Máximo lag m = ');
disp('Elija la Ventana: (1)Hann  (2)Hamming  (3)Parzen')
ventana=input('Ventana nº: ');
disp('Elija intervalos de confianza: (1)95% (2)90%')
intervalo=input('Intervalo nº: ');
[n,q]=size(Matriz);


% Calcula las frecuencias y las guarda

for k=0:m
 f(k+1)=k/(2*m);
 Resul(k+1,1)=f(k+1);
end


% Calcula intervalos de confianza

disp('Grados de libertad')
gl=(2*n-(m/2))/m

if intervalo==1
 chiinf=(.06414329*gl^1.157371-.09347972)/(1.783669*gl^(-1.319044)+.1701124)
 chisup=(.9295456*gl^.5908365+3.231091)/(.8309965*gl^(-.3875086)-.006014)
end

if intervalo==2
 chiinf=(.3559888*gl^2.408943-.3692226)/(.6930017*gl^1.299689+4.887728)
 chisup=(3.23178*gl^.9261549+19.64045)/(6.454432*gl^(-.8282872)+1.975337)
end


% Cálculo de la Ventana

if ventana==1
 for r=0:m
  u(r+1)=.5*(1+cos(pi*r/m));
 end
end

if ventana==2
 for r=0:m
  u(r+1)=.54+.46*cos(pi*r/m);
 end
end

if ventana==3
 for r=0:m
  if r<=(m/2)
   u(r+1)=1-6*((r/m)^2)*(1-(r/m));
  end
  if r>(m/2)
   u(r+1)=2*(1-(r/m))^3;
  end
 end
end


for j=1:q
x=Matriz(:,j);

 
  % Remoción de la media de la serie

  xm=mean(x);
  for i=1:n
   xd(i)=x(i)-xm;
  end


  % Ventana coseno en los extremos

  l=floor(.1*n);
  for i=1:l
   c=.5*(1-cos((i-1)*pi/(l-1)));
   xd(i)=c*xd(i);
   xd(n-i+1)=c*xd(n-i+1);
  end


  % Función de autocovarianzas

  for r=0:m
   s=0;
   for i=1:n-r
    s=s+xd(i)*xd(i+r);
   end
   Cov(r+1)=s/(n-r);
  end


  % Aplicación de la ventana

  R=Cov.*u;


  % Cálculo de espectro y su promedio

  Gsum=0;
  for k=0:m
   s=0;
   for r=1:m-1
    s=s+R(r+1)*cos(pi*r*k/m);
   end
   G(k+1)=2*(R(1)+2*s+R(m+1)*cos(pi*k));
   Gsum=Gsum+G(k+1);
  end
  Gmed=Gsum/(m+1);


  % Continuo nulo

  cor=corrcoef(x(2:n),x(1:n-1));
  r1=cor(1,2);
  for k=0:m
   Cn(k+1)=Gmed*((1-r1^2)/(1+r1^2-2*r1*cos(pi*k/m)));
   Cnsup(k+1)=Cn(k+1)*chisup/gl;
   Cninf(k+1)=Cn(k+1)*chiinf/gl;
  end


  % Gráfica

  figure
  plot(f,G,'-k',f,Cninf,'.r',f,Cn,'--g',f,Cnsup,'.r')
  

  % Guarda datos

  for k=0:m
   Resul(k+1,4*j-2)=G(k+1);
   Resul(k+1,4*j-1)=Cninf(k+1);
   Resul(k+1,4*j)=Cn(k+1);
   Resul(k+1,4*j+1)=Cnsup(k+1);
  end

  disp('frec  densidadespectral  cnsup  contnulo   cninf')
  disp(Resul)
  

  % Valores significativos

  disp('Valores significativos: f, 1/f, Cnsup, G')
  disp('Serie nº: '),j

  for k=1:m
   if G(k+1)>Cnsup(k+1)
    disp([f(k+1),1/f(k+1),Cnsup(k+1),G(k+1)])
   end
  end


  clear Cn Cninf Cnsup G Gmed Gsum R c cor
  clear i k l r r1 s xd xm x Cov

end

disp('Listo!')
