unit fftunit;

{$mode objfpc}{$H+}

interface

uses
  Classes, SysUtils, lowlevelunit;

{$DEFINE fftNormierung}

type
  tFFTDatenordnung = (doResIms,doResSmi,doAlleResIms,doAlleResSmi,doRes,doBetr,doBetrQdr,doGetrennt);

  tFFTAlgorithmus = class
  private
    len,hlen:      longint;
    res,ims:       tExtendedArray;
    ino,outo:      tFFTDatenordnung;
    procedure dumpResIms;
  public
    constructor create(laenge: longint; inputOrdnung,outputOrdnung: tFFTDatenordnung); overload;
    constructor create(original: tFFTAlgorithmus); overload;
    destructor destroy; override;
    procedure laden(invers: boolean; q: pSingle; schritt: longint); overload; virtual;
    procedure laden(invers: boolean; q: pDouble; schritt: longint); overload; virtual;
    procedure laden(invers: boolean; q: pExtended; schritt: longint); overload; virtual;
    procedure laden(invers: boolean; qRe,qIm: tExtendedArray); overload; virtual;
    procedure ausfuehren; virtual; abstract;
    procedure speichern(invers: boolean; z: pSingle; schritt: longint); overload; virtual;
    procedure speichern(invers: boolean; z: pDouble; schritt: longint); overload; virtual;
    procedure speichern(invers: boolean; z: pExtended; schritt: longint); overload; virtual;
    procedure speichern(invers: boolean; zRe,zIm: tExtendedArray); overload; virtual;
    procedure summen(amAnfang: boolean; var energie: extended; var nurNullen: boolean);
    function dumpParams: string;
  end;

  tFFTCooleyTukey = class(tFFTAlgorithmus)
  private
    wre,wim: tExtendedArray;
    perm:    tLongintArray;
  public
    constructor create(laenge: longint; inputOrdnung,outputOrdnung: tFFTDatenordnung); overload;
    constructor create(original: tFFTCooleyTukey); overload;
    destructor destroy; override;
    procedure laden(invers: boolean; q: pSingle; schritt: longint); overload; override;
    procedure laden(invers: boolean; q: pDouble; schritt: longint); overload; override;
    procedure laden(invers: boolean; q: pExtended; schritt: longint); overload; override;
    procedure laden(invers: boolean; qRe,qIm: tExtendedArray); overload; override;
    procedure ausfuehren; override;
  end;

  tFFTBluestein = class(tFFTAlgorithmus)
  private
    rRe,rIm,sRe,sIm: array of extended;
    subFFT:          tFFTCooleyTukey;
  public
    constructor create(laenge: longint; inputOrdnung,outputOrdnung: tFFTDatenordnung); overload;
    constructor create(original: tFFTBluestein); overload;
    destructor destroy; override;
    procedure ausfuehren; override;
  end;

  function createFFTAlgorithmus(laenge: longint; inputOrdnung,outputOrdnung: tFFTDatenordnung): tFFTAlgorithmus; overload;
  function createFFTAlgorithmus(original: tFFTAlgorithmus): tFFTAlgorithmus; overload;

implementation

uses
  math;

// allgemeine Funktionen *******************************************************

function createFFTAlgorithmus(laenge: longint; inputOrdnung,outputOrdnung: tFFTDatenordnung): tFFTAlgorithmus;
begin
  if inputOrdnung in [doBetr,doBetrQdr] then
    fehler('Ich brauche mehr als Beträge oder Betragsquadrate um eine FFT auszurechnen!');

  if round(power(2,round(ln(laenge)/ln(2))))=laenge then
    result:=tFFTCooleyTukey.create(laenge,inputOrdnung,outputOrdnung)
  else
    result:=tFFTBluestein.create(laenge,inputOrdnung,outputOrdnung);
end;

function createFFTAlgorithmus(original: tFFTAlgorithmus): tFFTAlgorithmus;
begin
  if original is tFFTCooleyTukey then
    result:=tFFTCooleyTukey.create(original as tFFTCooleyTukey)
  else if original is tFFTBluestein then
    result:=tFFTBluestein.create(original as tFFTBluestein)
  else
    fehler('Unbekannter FFT-Algorithmus-Typ '''+original.className+'!');
end;

// tFFTAlgorithmus *************************************************************
{$DEFINE tFFTAlgorithmus}

constructor tFFTAlgorithmus.create(laenge: longint; inputOrdnung,outputOrdnung: tFFTDatenordnung);
begin
  inherited create;
  ino:=inputOrdnung;
  outo:=outputOrdnung;
  len:=laenge;
  hlen:=laenge div 2;
  fillchar(res,sizeof(res),#0);
  setlength(res,0);
  fillchar(ims,sizeof(ims),#0);
  setlength(ims,0);
end;

constructor tFFTAlgorithmus.create(original: tFFTAlgorithmus);
begin
  inherited create;
  ino:=original.ino;
  outo:=original.outo;
  len:=original.len;
  hlen:=original.hlen;
  fillchar(res,sizeof(res),#0);
  setlength(res,0);
  fillchar(ims,sizeof(ims),#0);
  setlength(ims,0);
end;

destructor tFFTAlgorithmus.destroy;
begin
  setlength(res,0);
  setlength(ims,0);
  inherited destroy;
end;

procedure tFFTAlgorithmus.dumpResIms;
var
  i: longint;
begin
  writeln('***');
  for i:=0 to length(res)-1 do
    writeln(res[i],' ',ims[i]);
  writeln('***');
end;

{$DEFINE fftunitLaden}
procedure tFFTAlgorithmus.laden(invers: boolean; q: pSingle; schritt: longint);
{$INCLUDE fftunit.inc}
procedure tFFTAlgorithmus.laden(invers: boolean; q: pDouble; schritt: longint);
{$INCLUDE fftunit.inc}
procedure tFFTAlgorithmus.laden(invers: boolean; q: pExtended; schritt: longint);
{$INCLUDE fftunit.inc}
{$UNDEF fftunitLaden}

{$DEFINE fftunitGetrenntLaden}
procedure tFFTAlgorithmus.laden(invers: boolean; qRe,qIm: tExtendedArray);
{$INCLUDE fftunit.inc}
{$UNDEF fftunitGetrenntLaden}

{$DEFINE fftunitSpeichern}
procedure tFFTAlgorithmus.speichern(invers: boolean; z: pSingle; schritt: longint);
{$INCLUDE fftunit.inc}
procedure tFFTAlgorithmus.speichern(invers: boolean; z: pDouble; schritt: longint);
{$INCLUDE fftunit.inc}
procedure tFFTAlgorithmus.speichern(invers: boolean; z: pExtended; schritt: longint);
{$INCLUDE fftunit.inc}
{$UNDEF fftunitSpeichern}

{$DEFINE fftunitGetrenntSpeichern}
procedure tFFTAlgorithmus.speichern(invers: boolean; zRe,zIm: tExtendedArray);
{$INCLUDE fftunit.inc}
{$UNDEF fftunitGetrenntSpeichern}

procedure tFFTAlgorithmus.summen(amAnfang: boolean; var energie: extended; var nurNullen: boolean);
var
  i:    longint;
  tmpE: extended;
begin
  tmpE:=0;
  for i:=0 to len-1 do begin
    tmpE:=tmpE + sqr(res[i]) + sqr(ims[i]);
    nurNullen:=nurNullen and (res[i]=0) and (ims[i]=0);
  end;
  {$IFDEF fftNormierung}
  if amAnfang then
    energie:=energie + tmpE*len
  else
  {$ENDIF}
    energie:=energie + tmpE;
end;

function tFFTAlgorithmus.dumpParams: string;
begin
  result:=className+': '+inttostr(len)+' ('+inttostr(hlen)+')';
end;

{$UNDEF tFFTAlgorithmus}

// tFFTCooleyTukey *************************************************************
{$DEFINE tFFTCooleyTukey}

constructor tFFTCooleyTukey.create(laenge: longint; inputOrdnung,outputOrdnung: tFFTDatenordnung);
var
  i,j,loglen: longint;
begin
  inherited create(laenge,inputOrdnung,outputOrdnung);

  setlength(res,len);
  setlength(ims,len);

  loglen:=round(ln(len)/ln(2));
  setlength(wre,hlen);
  setlength(wim,hlen);
  setlength(perm,len);
  for i:=0 to len-1 do begin
    res[i]:=nan;
    ims[i]:=nan;
    perm[i]:=0;
    for j:=0 to loglen-1 do
      perm[i]:=perm[i] or (byte(odd(i shr j)) shl (loglen-1-j));
  end;
  for i:=0 to hlen-1 do begin
    wre[i]:=cos(2*pi*i/len);
    wim[i]:=sin(2*pi*i/len);
  end;
end;

constructor tFFTCooleyTukey.create(original: tFFTCooleyTukey);
begin
  inherited create(original);
  setlength(res,length(original.res));
  setlength(ims,length(original.ims));
  setlength(wre,length(original.wre));
  move(original.wre[0],wre[0],sizeof(wre[0])*length(wre));
  setlength(wim,length(original.wim));
  move(original.wim[0],wim[0],sizeof(wim[0])*length(wim));
  setlength(perm,length(original.perm));
  move(original.perm[0],perm[0],sizeof(perm[0])*length(perm));
end;

destructor tFFTCooleyTukey.destroy;
begin
  setlength(wre,0);
  setlength(wim,0);
  setlength(perm,0);
  inherited destroy;
end;

{$DEFINE fftunitLaden}
procedure tFFTCooleyTukey.laden(invers: boolean; q: pSingle; schritt: longint);
{$INCLUDE fftunit.inc}
procedure tFFTCooleyTukey.laden(invers: boolean; q: pDouble; schritt: longint);
{$INCLUDE fftunit.inc}
procedure tFFTCooleyTukey.laden(invers: boolean; q: pExtended; schritt: longint);
{$INCLUDE fftunit.inc}
{$UNDEF fftunitLaden}

{$DEFINE fftunitGetrenntLaden}
procedure tFFTCooleyTukey.laden(invers: boolean; qRe,qIm: tExtendedArray);
{$INCLUDE fftunit.inc}
{$UNDEF fftunitGetrenntLaden}

procedure tFFTCooleyTukey.ausfuehren;
var
  i,dist,absch,wnum,wstep: longint;
  t1,t2:                   extended;
begin
  dist:=1;
  wstep:=hlen;
  while dist<len do begin
    absch:=0;
    while absch<len do begin
      wnum:=0;
      for i:=0 to dist-1 do begin
        // x_links:  [absch+j]
        // x_rechts: [absch+j+dist]
        t1:=wRe[wnum]*res[absch+i+dist] - wIm[wnum]*ims[absch+i+dist];
        t2:=wRe[wnum]*ims[absch+i+dist] + wIm[wnum]*res[absch+i+dist];

        res[absch+i+dist]:=res[absch+i]-t1;
        ims[absch+i+dist]:=ims[absch+i]-t2;
        res[absch+i]:=res[absch+i]+t1;
        ims[absch+i]:=ims[absch+i]+t2;
        wnum:=wnum+wstep;
      end;
      absch:=absch+2*dist;
    end;
    wstep:=wstep div 2;
    dist:=dist*2;
  end;
end;

{$UNDEF tFFTCooleyTukey}

// tFFTBluestein ***************************************************************

constructor tFFTBluestein.create(laenge: longint; inputOrdnung,outputOrdnung: tFFTDatenordnung);
var
  preFFT: tFFTCooleyTukey;
  e:      extended;
  i:      longint;
begin
  inherited create(laenge,inputOrdnung,outputOrdnung);
  setlength(res,round(power(2,ceil(ln(2*len-1)/ln(2)))));
  setlength(ims,length(res));
  setlength(rRe,length(res));
  setlength(rIm,length(res));
  setlength(sRe,length(res));
  setlength(sIm,length(res));
  subFFT:=tFFTCooleyTukey.create(length(res),doGetrennt,doGetrennt);

  for i:=0 to subFFT.len-1 do begin
    e:=sqr(min(i,subFFT.len-i))*pi/len;
    rRe[i]:=Cos(e);
    rIm[i]:=Sin(e);
  end;

  preFFT:=tFFTCooleyTukey.create(length(rRe),doGetrennt,doGetrennt);
  preFFT.laden(false,rRe,rIm);
  preFFT.ausfuehren;
  preFFT.speichern(false,sRe,sIm);
  preFFT.free;

  {$IFDEF fftNormierung}
  for i:=0 to subFFT.len-1 do begin
    sRe[i]:=sRe[i]*sqrt(subFFT.len);
    sIm[i]:=sIm[i]*sqrt(subFFT.len);
  end;
  {$ENDIF}
end;

constructor tFFTBluestein.create(original: tFFTBluestein);
begin
  inherited create(original);
  setlength(res,length(original.res));
  setlength(ims,length(original.ims));
  setlength(rRe,length(original.rRe));
  move(original.rRe[0],rRe[0],sizeof(rRe[0])*length(rRe));
  setlength(rIm,length(original.rIm));
  move(original.rIm[0],rIm[0],sizeof(rIm[0])*length(rIm));
  setlength(sRe,length(original.sRe));
  move(original.sRe[0],sRe[0],sizeof(sRe[0])*length(sRe));
  setlength(sIm,length(original.sIm));
  move(original.sIm[0],sIm[0],sizeof(sIm[0])*length(sIm));

  subFFT:=tFFTCooleyTukey.create(original.subFFT);
end;

destructor tFFTBluestein.destroy;
begin
  setlength(rRe,0);
  setlength(rIm,0);
  setlength(sRe,0);
  setlength(sIm,0);
  subFFT.free;
  res:=nil;
  ims:=nil;
  inherited destroy;
end;

procedure tFFTBluestein.ausfuehren;
var
  i:   longint;
  tmp: extended;
begin
  // a[j] := werte[j] * CC(raumFaktoren[j])
  for i:=0 to len-1 do begin
    tmp:=res[i];
    res[i]:=  tmp*rRe[i] + ims[i]*rIm[i];
    ims[i]:=- tmp*rIm[i] + ims[i]*rRe[i];
  end;
  for i:=len to length(res)-1 do begin
    res[i]:=0;
    ims[i]:=0;
  end;

  // a[j] := iFFT(a[j])
  subFFT.laden(true,res,ims);
  subFFT.ausfuehren;
  subFFT.speichern(true,res,ims);

  for i:=0 to length(res)-1 do begin
    res[i]:=res[i];
    ims[i]:=ims[i];
  end;

  // a[j] := a[j] * spektralFaktoren[j]
  for i:=0 to length(res)-1 do begin
    tmp:=res[i];
    res[i]:=tmp*sRe[i] - ims[i]*sIm[i];
    ims[i]:=tmp*sIm[i] + ims[i]*sRe[i];
  end;

  // a[j] := iFFT(a[j])
  subFFT.laden(false,res,ims);
  subFFT.ausfuehren;
  subFFT.speichern(false,res,ims);

  for i:=0 to length(res)-1 do begin
    res[i]:=res[i];
    ims[i]:=ims[i];
  end;

  // werte[j] := a[j] * CC ( raumFaktoren[j] )
  for i:=0 to len-1 do begin
    tmp:=res[i];
    res[i]:=  tmp*rRe[i] + ims[i]*rIm[i];
    ims[i]:=- tmp*rIm[i] + ims[i]*rRe[i];
  end;
end;

end.