fftunit neu

3 files changed, 525 insertions, 73 deletions

diff --git a/fftunit.inc b/fftunit.inc
new file mode 100644
index 0000000..7f5cc34
--- /dev/null
+++ b/fftunit.inc
@@ -0,0 +1,199 @@
+{$IFDEF fftunitLaden}
+//procedure tFFTAlgorithmus.laden(invers: boolean; const q: single; schritt: longint);
+var
+  i,inv:  longint;
+  {$IFDEF fftNormierung}
+  lenFak: extended;
+  {$ENDIF}
+begin
+  inv:=1-2*byte(invers);
+  {$IFDEF fftNormierung}
+  lenFak:=1/sqrt(len);
+  {$ENDIF}
+  case ino of
+    doRes:
+      for i:=0 to len-1 do begin
+        {$IFDEF tFFTCooleyTukey}
+        res[perm[i]]:=(q+i*schritt)^ {$IFDEF fftNormierung} * lenFak{$ENDIF};
+        ims[perm[i]]:=0;
+        {$ELSE}
+        res[i]:=(q+i*schritt)^ {$IFDEF fftNormierung} * lenFak{$ENDIF};
+        ims[i]:=0;
+        {$ENDIF}
+      end;
+    doAlleResIms:
+      for i:=0 to len-1 do begin
+        {$IFDEF tFFTCooleyTukey}
+        res[perm[i]]:=(q+i*schritt)^ {$IFDEF fftNormierung} * lenFak{$ENDIF};
+        ims[perm[i]]:=(q+(i+len)*schritt)^ * inv {$IFDEF fftNormierung} * lenFak{$ENDIF};
+        {$ELSE}
+        res[i]:=(q+i*schritt)^ {$IFDEF fftNormierung} * lenFak{$ENDIF};
+        ims[i]:=(q+(i+len)*schritt)^ * inv {$IFDEF fftNormierung} * lenFak{$ENDIF};
+        {$ENDIF}
+      end;
+    doAlleResSmi:
+      for i:=0 to hlen + (len and $1) - 1 do begin
+        {$IFDEF tFFTCooleyTukey}
+        res[perm[i]]:=(q+i*schritt)^ {$IFDEF fftNormierung} * lenFak{$ENDIF};
+        ims[perm[i]]:=(q+(2*len-1-i)*schritt)^ * inv {$IFDEF fftNormierung} * lenFak{$ENDIF};
+        {$ELSE}
+        res[i]:=(q+i*schritt)^ {$IFDEF fftNormierung} * lenFak{$ENDIF};
+        ims[i]:=(q+(2*len-1-i)*schritt)^ * inv {$IFDEF fftNormierung} * lenFak{$ENDIF};
+        {$ENDIF}
+      end;
+    doResIms: begin
+      res[0]:=q^ {$IFDEF fftNormierung} * lenFak{$ENDIF};                // perm[0]=0
+      ims[0]:=0;
+      for i:=1 to hlen + (len and $1) - 1 do begin
+        {$IFDEF tFFTCooleyTukey}
+        res[perm[i]]:=(q+i*schritt)^ {$IFDEF fftNormierung} * lenFak{$ENDIF};
+        ims[perm[i]]:=(q+(i+hlen)*schritt)^ * inv {$IFDEF fftNormierung} * lenFak{$ENDIF};
+        res[perm[len-i]]:=res[perm[i]];
+        ims[perm[len-i]]:=-ims[perm[i]];
+        {$ELSE}
+        res[i]:=(q+i*schritt)^ {$IFDEF fftNormierung} * lenFak{$ENDIF};
+        ims[i]:=(q+(i+hlen)*schritt)^ * inv {$IFDEF fftNormierung} * lenFak{$ENDIF};
+        res[len-i]:=res[i];
+        ims[len-i]:=-ims[i];
+        {$ENDIF}
+      end;
+      if not odd(len) then begin
+        {$IFDEF tFFTCooleyTukey}
+        res[1]:=(q+hlen*schritt)^ {$IFDEF fftNormierung} * lenFak{$ENDIF}; // perm[hlen]=1
+        ims[1]:=0;
+        {$ELSE}
+        res[hlen]:=(q+hlen*schritt)^ {$IFDEF fftNormierung} * lenFak{$ENDIF};
+        ims[hlen]:=0;
+        {$ENDIF}
+      end;
+    end;
+    doResSmi: begin
+      res[0]:=q^ {$IFDEF fftNormierung} * lenFak{$ENDIF};                // perm[0]=0
+      ims[0]:=0;
+      for i:=1 to hlen + (len and $1) - 1 do begin
+        {$IFDEF tFFTCooleyTukey}
+        res[perm[i]]:=(q+i*schritt)^ {$IFDEF fftNormierung} * lenFak{$ENDIF};
+        ims[perm[i]]:=(q+(len-i)*schritt)^ * inv {$IFDEF fftNormierung} * lenFak{$ENDIF};
+        res[perm[len-i]]:=res[perm[i]];
+        ims[perm[len-i]]:=-ims[perm[i]];
+        {$ELSE}
+        res[i]:=(q+i*schritt)^ {$IFDEF fftNormierung} * lenFak{$ENDIF};
+        ims[i]:=(q+(len-i)*schritt)^ * inv {$IFDEF fftNormierung} * lenFak{$ENDIF};
+        res[len-i]:=res[i];
+        ims[len-i]:=-ims[i];
+        {$ENDIF}
+      end;
+      if not odd(len) then begin
+        {$IFDEF tFFTCooleyTukey}
+        res[1]:=(q+hlen*schritt)^ {$IFDEF fftNormierung} * lenFak{$ENDIF}; // perm[hlen]=1
+        ims[1]:=0;
+        {$ELSE}
+        res[hlen]:=(q+hlen*schritt)^ {$IFDEF fftNormierung} * lenFak{$ENDIF};
+        ims[hlen]:=0;
+        {$ENDIF}
+      end;
+    end;
+    doBetr,doBetrQdr:
+      fehler('Ich brauche mehr als Beträge oder Betragsquadrate um eine FFT auszurechnen!');
+    else
+      fehler('Diese Inputordnung ist in dieser Methode nicht implementiert!');;
+  end{of case};
+end;
+{$ENDIF}
+
+{$IFDEF fftunitGetrenntLaden}
+//procedure tFFTAlgorithmus.laden(invers: boolean; qRe,qIm: tExtendedArray);
+var
+  i,inv:  longint;
+  {$IFDEF fftNormierung}
+  lenFak: extended;
+  {$ENDIF}
+begin
+  inv:=1-2*byte(invers);
+  {$IFDEF fftNormierung}
+  lenFak:=1/sqrt(len);
+  {$ENDIF}
+  case ino of
+    doGetrennt:
+      for i:=0 to len-1 do begin
+        {$IFDEF tFFTCooleyTukey}
+        res[perm[i]]:=qRe[i] {$IFDEF fftNormierung} * lenFak{$ENDIF};
+        ims[perm[i]]:=qIm[i] * inv {$IFDEF fftNormierung} * lenFak{$ENDIF};
+        {$ELSE}
+        res[i]:=qRe[i] {$IFDEF fftNormierung} * lenFak{$ENDIF};
+        ims[i]:=qIm[i] * inv {$IFDEF fftNormierung} * lenFak{$ENDIF};
+        {$ENDIF}
+      end;
+    else
+      fehler('Diese Inputordnung ist in dieser Methode nicht implementiert!');;
+  end{of case};
+end;
+{$ENDIF}
+
+{$IFDEF fftunitSpeichern}
+//procedure tFFTAlgorithmus.speichern(invers: boolean; var z: single; schritt: longint);
+var
+  i,inv: longint;
+begin
+  inv:=1-2*byte(invers);
+  case outo of
+    doRes:
+      for i:=0 to len-1 do
+        (z+i*schritt)^:=res[i];
+    doAlleResIms:
+      for i:=0 to len-1 do begin
+        (z+i*schritt)^:=res[i];
+        (z+(i+len)*schritt)^:=ims[i] * inv;
+      end;
+    doAlleResSmi:
+      for i:=0 to len-1 do begin
+        (z+i*schritt)^:=res[i];
+        (z+(2*len-1-i)*schritt)^:=ims[i] * inv;
+      end;
+    doResIms: begin
+      z^:=res[0];
+      for i:=1 to hlen + (len and $1) - 1 do begin
+        (z+i*schritt)^:=res[i];
+        (z+(i+hlen)*schritt)^:=ims[i] * inv;
+      end;
+      if not odd(len) then
+        (z+hlen*schritt)^:=res[hlen];
+    end;
+    doResSmi: begin
+      z^:=res[0];
+      for i:=1 to hlen + (len and $1) - 1 do begin
+        (z+i*schritt)^:=res[i];
+        (z+(len-i)*schritt)^:=ims[i] * inv;
+      end;
+      if not odd(len) then
+        (z+hlen*schritt)^:=res[hlen];
+    end;
+    doBetr:
+      for i:=0 to len-1 do
+        (z+i*schritt)^:=sqrt(sqr(res[i])+sqr(ims[i]));
+    doBetrQdr:
+      for i:=0 to len-1 do
+        (z+i*schritt)^:=sqr(res[i])+sqr(ims[i]);
+    else
+      fehler('Diese Outputordnung ist in dieser Methode nicht implementiert!');;
+  end{of case};
+end;
+{$ENDIF}
+
+{$IFDEF fftunitGetrenntSpeichern}
+//procedure tFFTAlgorithmus.speichern(invers: boolean; zRe,zIm: tExtendedArray);
+var
+  i,inv: longint;
+begin
+  inv:=1-2*byte(invers);
+  case outo of
+    doGetrennt:
+      for i:=0 to len-1 do begin
+        zRe[i]:=res[i];
+        zIm[i]:=ims[i] * inv;
+      end;
+    else
+      fehler('Diese Outputordnung ist in dieser Methode nicht implementiert!');;
+  end{of case};
+end;
+{$ENDIF}
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.
+
diff --git a/matheunit.pas b/matheunit.pas
index 431345c..61c1eea 100644
--- a/matheunit.pas
+++ b/matheunit.pas
@@ -49,7 +49,6 @@ function nullfunktion(x: extended): extended;
 function exprtofloat(st: boolean; s: string; kvs: tKnownValues; cbgv: tCallBackGetValue): extended;
 function formelnAuswerten(s: string; kvs: tKnownValues; cbgv: tCallBackGetValue): string;
 function knownValue(nam: string; val: extended): tKnownValue;
-procedure fft(var res,ims: tExtendedArray; inv: boolean); inline;
 
 implementation
 
@@ -593,77 +592,5 @@ begin
   result.value:=val;
 end;
 
-procedure fft(var res,ims: tExtendedArray; inv: boolean);
-var
-  haL,i,j,k,n,dist,wstep,absch,wnum: longint;
-  umsortierung:                      tLongintArray;
-  wRe,wIm:                           tExtendedArray;
-  t1,t2:                             extended;
-begin
-  haL:=length(res) div 2;
-  n:=round(ln(length(res))/ln(2));
-
-  if round(power(2,n))<>length(res) then
-    raise exception.create(inttostr(length(res))+' ist keine Zweierpotenz!');
-
-  setlength(umsortierung,2*haL);
-  for j:=0 to 2*haL-1 do begin
-    k:=0;
-    for i:=0 to n-1 do             // Bits umkehren
-      k:=2*k + byte(odd(j shr i));
-    umsortierung[j]:=k;
-  end;
-  setlength(wRe,haL);
-  setlength(wIm,haL);
-  for i:=0 to haL-1 do begin
-    wRe[i]:=cos(pi*i/haL);
-    wIm[i]:=sin(pi*i/haL);
-  end;
-
-  for j:=0 to 2*haL-1 do begin
-    if umsortierung[j]>j then begin
-      t1:=res[j];
-      res[j]:=res[umsortierung[j]];
-      res[umsortierung[j]]:=t1;
-      t1:=ims[j];
-      ims[j]:=ims[umsortierung[j]] * (1-2*byte(inv));
-      ims[umsortierung[j]]:=t1 * (1-2*byte(inv));
-    end;
-    if (umsortierung[j]=j) and inv then
-      ims[j]:=-ims[j];
-  end;
-
-  dist:=1;
-  wstep:=haL;
-
-  while dist<2*haL do begin
-    absch:=0;
-    while absch<2*haL do begin
-      wnum:=0;
-      for j:=0 to dist-1 do begin
-        // x_links:  [absch+j]
-        // x_rechts: [absch+j+dist]
-        t1:=wRe[wnum]*res[absch+j+dist] - wIm[wnum]*ims[absch+j+dist];
-        t2:=wRe[wnum]*ims[absch+j+dist] + wIm[wnum]*res[absch+j+dist];
-
-        res[absch+j+dist]:=res[absch+j]-t1;
-        ims[absch+j+dist]:=ims[absch+j]-t2;
-        res[absch+j]:=res[absch+j]+t1;
-        ims[absch+j]:=ims[absch+j]+t2;
-        wnum:=wnum+wstep;
-      end;
-      absch:=absch+2*dist;
-    end;
-    wstep:=wstep div 2;
-    dist:=dist*2;
-  end;
-
-  t1:=1/sqrt(length(res));
-  for i:=0 to 2*haL-1 do begin
-    res[i]:=res[i]*t1;
-    ims[i]:=ims[i]*t1;
-  end;
-end;
-
 end.