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: pSingle; schritt: longint); overload; virtual;
    procedure laden(invers: boolean; qRe,qIm: pDouble; schritt: longint); overload; virtual;
    procedure laden(invers: boolean; qRe,qIm: pExtended; schritt: longint); 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: pSingle; schritt: longint); overload; virtual;
    procedure speichern(invers: boolean; zRe,zIm: pDouble; schritt: longint); overload; virtual;
    procedure speichern(invers: boolean; zRe,zIm: pExtended; schritt: longint); overload; virtual;
    procedure summen(amAnfang: boolean; var energie: extended; var nurNullen: boolean);
    function dumpParams: string;
    property inOrdnung: tFFTDatenordnung read inO;
    property outOrdnung: tFFTDatenordnung read outO;
  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: pSingle; schritt: longint); overload; override;
    procedure laden(invers: boolean; qRe,qIm: pDouble; schritt: longint); overload; override;
    procedure laden(invers: boolean; qRe,qIm: pExtended; schritt: longint); 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;

  function strToFftDo(out fftDo: tFFTDatenordnung; s: string): boolean;
  function fftDoToStr(fftDo: tFFTDatenordnung): string;

implementation

uses
  math, myStringListUnit;

// 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;

function strToFftDo(out fftDo: tFFTDatenordnung; s: string): boolean;
var
  bekannteBefehle: tMyStringList;
begin
  fftDo:=low(tFFTDatenordnung);
  bekannteBefehle:=tMyStringList.create;
  repeat
    if bekannteBefehle.count>0 then inc(fftDo);
    bekannteBefehle.add(fftDoToStr(fftDo));
  until (s=bekannteBefehle.last) or (fftDo=high(tFFTDatenordnung));
  if fftDo<high(tFFTDatenordnung) then begin
    bekannteBefehle.free;
    result:=true;
    exit;
  end;
  result:=false;
  bekannteBefehle.sort;
  gibAus('Kenne FFT-Datenordnung '''+s+''' nicht!'#10'Ich kenne:'#10+bekannteBefehle.text,3);
  bekannteBefehle.free;
  result:=false;
end;

function fftDoToStr(fftDo: tFFTDatenordnung): string;
begin
  case fftDo of
    doResIms:
      result:='Realteile:Imaginärteile';
    doResSmi:
      result:='Realteile:Imaginärteile umgedreht';
    doAlleResIms:
      result:='alle Realteile:Imaginärteile';
    doAlleResSmi:
      result:='alle Realteile:Imaginärteile umgedreht';
    doRes:
      result:='Realteile';
    doBetr:
      result:='Beträge';
    doBetrQdr:
      result:='Betragsquadrate';
    doGetrennt:
      result:='Realteile und Imaginärteile getrennt';
    else
      result:='unbekannte fftDatenOrdnung';
  end{of case};
end;

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

constructor tFFTAlgorithmus.create(laenge: longint; inputOrdnung,outputOrdnung: tFFTDatenordnung);
begin
  inherited create;
  inO:=inputOrdnung;
  outO:=outputOrdnung;
  if (inO=doAlleResSmi) or (outO=doAlleResSmi) then
    raise exception.create('Die Datenordnung doAlleResSmi ist momentan noch verbugt!');
  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;
  if (inO=doAlleResSmi) or (outO=doAlleResSmi) then
    raise exception.create('Die Datenordnung doAlleResSmi ist momentan noch verbugt!');
  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: pSingle; schritt: longint);
{$INCLUDE fftunit.inc}
procedure tFFTAlgorithmus.laden(invers: boolean; qRe,qIm: pDouble; schritt: longint);
{$INCLUDE fftunit.inc}
procedure tFFTAlgorithmus.laden(invers: boolean; qRe,qIm: pExtended; schritt: longint);
{$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: pSingle; schritt: longint);
{$INCLUDE fftunit.inc}
procedure tFFTAlgorithmus.speichern(invers: boolean; zRe,zIm: pDouble; schritt: longint);
{$INCLUDE fftunit.inc}
procedure tFFTAlgorithmus.speichern(invers: boolean; zRe,zIm: pExtended; schritt: longint);
{$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: pSingle; schritt: longint);
{$INCLUDE fftunit.inc}
procedure tFFTCooleyTukey.laden(invers: boolean; qRe,qIm: pDouble; schritt: longint);
{$INCLUDE fftunit.inc}
procedure tFFTCooleyTukey.laden(invers: boolean; qRe,qIm: pExtended; schritt: longint);
{$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,pExtended(@(rRe[0])),pExtended(@(rIm[0])),1);
  preFFT.ausfuehren;
  preFFT.speichern(false,pExtended(@(sRe[0])),pExtended(@(sIm[0])),1);
  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,pExtended(@(res[0])),pExtended(@(ims[0])),1);
  subFFT.ausfuehren;
  subFFT.speichern(true,pExtended(@(res[0])),pExtended(@(ims[0])),1);

  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,pExtended(@(res[0])),pExtended(@(ims[0])),1);
  subFFT.ausfuehren;
  subFFT.speichern(false,pExtended(@(res[0])),pExtended(@(ims[0])),1);

  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.