fft in matheunit.pas eingefuegt

1 files changed, 73 insertions, 0 deletions

diff --git a/matheunit.pas b/matheunit.pas
index 8e6426c..710a46b 100644
--- a/matheunit.pas
+++ b/matheunit.pas
@@ -45,6 +45,7 @@ procedure copyArray(i: tExtendedArray; out o: tExtendedArray); overload;
 function nullfunktion(x: extended): extended;
 function exprtofloat(st: boolean; s: string; kvs: tKnownValues; cbgv: tCallBackGetValue): extended;
 function knownValue(nam: string; val: extended): tKnownValue;
+procedure fft(var res,ims: tExtendedArray; inv: boolean); inline;
 
 implementation
 
@@ -502,5 +503,77 @@ 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.