fftunit neu

1 files changed, 326 insertions, 0 deletions

diff --git a/fftunit.pas b/fftunit.pas
new file mode 100644
index 0000000..f78935f
--- /dev/null
+++ b/fftunit.pas
@@ -0,0 +1,326 @@
+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;
+  public
+    procedure dumpResIms;
+    constructor create(laenge: longint; inputOrdnung,outputOrdnung: tFFTDatenordnung);
+    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;
+  end;
+
+  tFFTCooleyTukey = class(tFFTAlgorithmus)
+  private
+    wre,wim: tExtendedArray;
+    perm:    tLongintArray;
+  public
+    constructor create(laenge: longint; inputOrdnung,outputOrdnung: tFFTDatenordnung);
+    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);
+    destructor destroy; override;
+    procedure ausfuehren; override;
+  end;
+
+  function createFFTAlgorithmus(laenge: longint; inputOrdnung,outputOrdnung: tFFTDatenordnung): tFFTAlgorithmus;
+
+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;
+
+// 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;
+
+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}
+
+{$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;
+
+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;
+
+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.
+