{+doc
+NOMBRE: Autoco.
+CREACION: 8.12.1990.
+AUTORES: Ruben Chaer.
+REGISTRO:
+TIPO: Unidad Pascal.
+PROPOSITO: Funcion autocorrelación discreta.
+PROYECTO: General.

+REVISION:
+AUTOR:
+DESCRIPCION:
-doc}

unit AutoCo;

interface

uses
  xMatDefs, SysUtils, matreal,
  WRFFTF01, WRFFTB01, WRFFTI01;

procedure AutoCorrelacion(n: longint; var dat);

// Retorna el espectro de potencia en el mismo vector.
// En la primer posición el término de continua, y luego armónico, 1, 2, et.
procedure h2f(n: longint; var dat; afft_base: TWRFFT_base = nil); overload;
procedure h2f(dat: TDAOfNReal; afft_base: TWRFFT_base = nil); overload;
procedure h2f(dat: TVectR; afft_base: TWRFFT_base = nil); overload;

implementation

const
  ErrMsg: string = 'AUTOCO: No hay suficiente Memoria';

procedure AutoCorrelacion(n: longint; var dat);
type
  VR = array[1..1024 * 1024] of NReal;
var
  pd, pa, pb: ^VR;
  a0: NReal;
  k: longint;
  fft_base: TWRFFT_base;
begin
  { Inicilalizamos el módulo de FFTs para trabajar con
  vectores de largo  2n }
  fft_base := TWRFFT_base.Create(2 * n);

  try
    GetMem(pd, n * 2 * sizeOf(NReal));
    GetMem(pa, fft_base.NBarrasLaterales * sizeOf(NReal));
    GetMem(pb, fft_base.NBarrasLaterales * sizeOf(NReal));
  except
    raise Exception.Create(ErrMsg);
  end;


  { C lculo del témino de DC }
  a0 := 0;
  for k := 1 to n do
  begin
    a0 := a0 + VR(dat)[k];
  end;
  a0 := a0 / n;

  { Quitamos el trmino de DC a los datos y duplicamos el
  vector de datos (en longitud), llenando con 0 los nuevos valores. }
  for k := 1 to n do
  begin
    pd^[k] := VR(dat)[k]; // - a0;
    pd^[k + n] := 0;
  end;


  { C lculamos los coeficientes de la serie discreta de Fourier
  correspondiente al vector pd con 2n datos.
    a0 : coeficiente de continua. (debe ser 0 pues se lo restamos antes)
    pa[k] :coeficiente del coseno del arm¢nico k, k:1..n
    pb[k] :coeficiente del seno del arm¢nico k, k:1..n  }
  FFTF(fft_base, pd^, a0, pa^, pb^);

  { Ahora, acomodamos a0, pa y pb de manera que representen los
  coeficientes de la transformada de la autocorrenlaci¢n de los
  datos. En otras palabras, calculamos el espectro de potencia. }
  a0 := a0 * a0;
  for k := 1 to fft_base.NBarrasLaterales do
  begin
    pa^[k] := (sqr(pa^[k]) + sqr(pb^[k])) / 2;
    pb^[k] := 0;
  end;
  //  pa^[n] := 2 * pa^[n]; {???, fundamentar esto}

  FFTB(fft_base, pd^, a0, pa^, pb^);
  fft_base.Free;

  for k := 1 to n do
    VR(dat)[k] := 2 * pd^[k];

  FreeMem(pd, n * 2 * SizeOf(NReal));
  FreeMem(pa, fft_base.NBarrasLaterales * SizeOf(NReal));
  FreeMem(pb, fft_base.NBarrasLaterales * sizeOf(NReal));

end;


procedure h2f(n: longint; var dat; afft_base: TWRFFT_base = nil);
type
  VR = array[1..1024 * 1024] of NReal;
var
  pd, pa, pb: ^VR;
  a0: NReal;
  k: longint;
  fft_base: TWRFFT_base;
  flg_freebase: boolean;
begin

  if afft_base <> nil then
  begin
    flg_freebase := False;
    fft_base := afft_base;
    fft_base.ReInicializar(n);
  end
  else
  begin
    fft_base := TWRFFT_base.Create(n);
    flg_freebase := True;
  end;
  pd := @dat;

  try
    GetMem(pa, fft_base.NBarrasLaterales * sizeOf(NReal));
    GetMem(pb, fft_base.NBarrasLaterales * SizeOf(NReal));
  except
    raise Exception.Create(ErrMsg);
  end;

  { C lculamos los coeficientes de la serie discreta de Fourier
  correspondiente al vector pd con 2n datos.
    a0 : coeficiente de continua. (debe ser 0 pues se lo restamos antes)
    pa[k] :coeficiente del coseno del arm¢nico k, k:1..n
    pb[k] :coeficiente del seno del arm¢nico k, k:1..n  }
  FFTF(fft_base, pd^, a0, pa^, pb^);

  { Ahora, acomodamos a0, pa y pb de manera que representen los
  coeficientes de la transformada de la autocorrenlaci¢n de los
  datos. En otras palabras, calculamos el espectro de potencia. }
  VR(dat)[1] := a0 * a0;
  for k := 1 to fft_base.NBarrasLaterales do
    VR(dat)[1 + k] := (sqr(pa^[k]) + sqr(pb^[k])) / 2;


  if flg_freebase then
    fft_base.Free;

  FreeMem(pb, fft_base.NBarrasLaterales * SizeOf(NReal));
  FreeMem(pa, fft_base.NBarrasLaterales * SizeOf(NReal));
end;

procedure h2f(dat: TDAOfNReal; afft_base: TWRFFT_base = nil);
begin
  h2f(length(dat), dat[0], afft_base);
end;

procedure h2f(dat: TVectR; afft_base: TWRFFT_base = nil);
begin
  h2f(dat.n, dat.pv[1], afft_base);
end;

end.
