From df7061ed77b9c7e412e4bb068b2eda81025cff70 Mon Sep 17 00:00:00 2001
From: Erich Eckner <git@eckner.net>
Date: Fri, 21 Jul 2017 10:39:43 +0200
Subject: fftunit: getrenntes Einlesen von Real- und Imaginaerteilen jetzt auch
 von Pointer und fuer doAlleResIms und doAlleResSmi (braucht man fuer
 voll-komplexe fft2d)

---
 fftunit.inc | 64 +++++++++++++++++++++++++----------------------
 fftunit.pas | 82 ++++++++++++++++++++++++++++++++++++++-----------------------
 2 files changed, 85 insertions(+), 61 deletions(-)

diff --git a/fftunit.inc b/fftunit.inc
index 7f5cc34..8d465f6 100644
--- a/fftunit.inc
+++ b/fftunit.inc
@@ -10,7 +10,7 @@ begin
   {$IFDEF fftNormierung}
   lenFak:=1/sqrt(len);
   {$ENDIF}
-  case ino of
+  case inO of
     doRes:
       for i:=0 to len-1 do begin
         {$IFDEF tFFTCooleyTukey}
@@ -32,7 +32,7 @@ begin
         {$ENDIF}
       end;
     doAlleResSmi:
-      for i:=0 to hlen + (len and $1) - 1 do begin
+      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};
@@ -44,33 +44,33 @@ begin
     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
+      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};
+        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};
+        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
+        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;
+        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
+      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};
@@ -85,11 +85,11 @@ begin
       end;
       if not odd(len) then begin
         {$IFDEF tFFTCooleyTukey}
-        res[1]:=(q+hlen*schritt)^ {$IFDEF fftNormierung} * lenFak{$ENDIF}; // perm[hlen]=1
+        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;
+        res[hLen]:=(q+hLen*schritt)^ {$IFDEF fftNormierung} * lenFak{$ENDIF};
+        ims[hLen]:=0;
         {$ENDIF}
       end;
     end;
@@ -102,7 +102,7 @@ end;
 {$ENDIF}
 
 {$IFDEF fftunitGetrenntLaden}
-//procedure tFFTAlgorithmus.laden(invers: boolean; qRe,qIm: tExtendedArray);
+//procedure tFFTAlgorithmus.laden(invers: boolean; qRe,qIm: pSingle; schritt: longint);
 var
   i,inv:  longint;
   {$IFDEF fftNormierung}
@@ -113,15 +113,17 @@ begin
   {$IFDEF fftNormierung}
   lenFak:=1/sqrt(len);
   {$ENDIF}
-  case ino of
-    doGetrennt:
+  case inO of
+    doGetrennt,doAlleResIms,doAlleResSmi:
+      // Real- und Imaginärteile sind _immer_ in gleicher Reihenfolge;
+      // die ggf. verdrehte Dimension ist senkrecht dazu! (benötigt für fft2d)
       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};
+        res[perm[i]]:=(qRe+i*schritt)^ {$IFDEF fftNormierung} * lenFak{$ENDIF};
+        ims[perm[i]]:=(qIm+i*schritt)^ * inv {$IFDEF fftNormierung} * lenFak{$ENDIF};
         {$ELSE}
-        res[i]:=qRe[i] {$IFDEF fftNormierung} * lenFak{$ENDIF};
-        ims[i]:=qIm[i] * inv {$IFDEF fftNormierung} * lenFak{$ENDIF};
+        res[i]:=(qRe+i*schritt)^ {$IFDEF fftNormierung} * lenFak{$ENDIF};
+        ims[i]:=(qIm+i*schritt)^ * inv {$IFDEF fftNormierung} * lenFak{$ENDIF};
         {$ENDIF}
       end;
     else
@@ -136,7 +138,7 @@ var
   i,inv: longint;
 begin
   inv:=1-2*byte(invers);
-  case outo of
+  case outO of
     doRes:
       for i:=0 to len-1 do
         (z+i*schritt)^:=res[i];
@@ -152,21 +154,21 @@ begin
       end;
     doResIms: begin
       z^:=res[0];
-      for i:=1 to hlen + (len and $1) - 1 do begin
+      for i:=1 to hLen + (len and $1) - 1 do begin
         (z+i*schritt)^:=res[i];
-        (z+(i+hlen)*schritt)^:=ims[i] * inv;
+        (z+(i+hLen)*schritt)^:=ims[i] * inv;
       end;
       if not odd(len) then
-        (z+hlen*schritt)^:=res[hlen];
+        (z+hLen*schritt)^:=res[hLen];
     end;
     doResSmi: begin
       z^:=res[0];
-      for i:=1 to hlen + (len and $1) - 1 do begin
+      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];
+        (z+hLen*schritt)^:=res[hLen];
     end;
     doBetr:
       for i:=0 to len-1 do
@@ -181,16 +183,18 @@ end;
 {$ENDIF}
 
 {$IFDEF fftunitGetrenntSpeichern}
-//procedure tFFTAlgorithmus.speichern(invers: boolean; zRe,zIm: tExtendedArray);
+//procedure tFFTAlgorithmus.speichern(invers: boolean; zRe,zIm: pSingle; schritt: longint);
 var
   i,inv: longint;
 begin
   inv:=1-2*byte(invers);
-  case outo of
-    doGetrennt:
+  case outO of
+    doGetrennt,doAlleResIms,doAlleResSmi:
+      // Real- und Imaginärteile sind _immer_ in gleicher Reihenfolge;
+      // die ggf. verdrehte Dimension ist senkrecht dazu! (benötigt für fft2d)
       for i:=0 to len-1 do begin
-        zRe[i]:=res[i];
-        zIm[i]:=ims[i] * inv;
+        (zRe+i*schritt)^:=res[i];
+        (zIm+i*schritt)^:=ims[i] * inv;
       end;
     else
       fehler('Diese Outputordnung ist in dieser Methode nicht implementiert!');;
diff --git a/fftunit.pas b/fftunit.pas
index b4b03bf..7b5995f 100644
--- a/fftunit.pas
+++ b/fftunit.pas
@@ -14,9 +14,9 @@ type
 
   tFFTAlgorithmus = class
   private
-    len,hlen:      longint;
+    len,hLen:      longint;
     res,ims:       tExtendedArray;
-    ino,outo:      tFFTDatenordnung;
+    inO,outO:      tFFTDatenordnung;
     procedure dumpResIms;
   public
     constructor create(laenge: longint; inputOrdnung,outputOrdnung: tFFTDatenordnung); overload;
@@ -25,14 +25,20 @@ type
     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 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: tExtendedArray); 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)
@@ -46,7 +52,9 @@ type
     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 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;
 
@@ -157,10 +165,10 @@ end;
 constructor tFFTAlgorithmus.create(laenge: longint; inputOrdnung,outputOrdnung: tFFTDatenordnung);
 begin
   inherited create;
-  ino:=inputOrdnung;
-  outo:=outputOrdnung;
+  inO:=inputOrdnung;
+  outO:=outputOrdnung;
   len:=laenge;
-  hlen:=laenge div 2;
+  hLen:=laenge div 2;
   fillchar(res,sizeOf(res),#0);
   setLength(res,0);
   fillchar(ims,sizeOf(ims),#0);
@@ -170,10 +178,10 @@ end;
 constructor tFFTAlgorithmus.create(original: tFFTAlgorithmus);
 begin
   inherited create;
-  ino:=original.ino;
-  outo:=original.outo;
+  inO:=original.inO;
+  outO:=original.outO;
   len:=original.len;
-  hlen:=original.hlen;
+  hLen:=original.hLen;
   fillchar(res,sizeOf(res),#0);
   setLength(res,0);
   fillchar(ims,sizeOf(ims),#0);
@@ -207,7 +215,11 @@ procedure tFFTAlgorithmus.laden(invers: boolean; q: pExtended; schritt: longint)
 {$UNDEF fftunitLaden}
 
 {$DEFINE fftunitGetrenntLaden}
-procedure tFFTAlgorithmus.laden(invers: boolean; qRe,qIm: tExtendedArray);
+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}
 
@@ -221,7 +233,11 @@ procedure tFFTAlgorithmus.speichern(invers: boolean; z: pExtended; schritt: long
 {$UNDEF fftunitSpeichern}
 
 {$DEFINE fftunitGetrenntSpeichern}
-procedure tFFTAlgorithmus.speichern(invers: boolean; zRe,zIm: tExtendedArray);
+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}
 
@@ -245,7 +261,7 @@ end;
 
 function tFFTAlgorithmus.dumpParams: string;
 begin
-  result:=className+': '+intToStr(len)+' ('+intToStr(hlen)+')';
+  result:=className+': '+intToStr(len)+' ('+intToStr(hLen)+')';
 end;
 
 {$UNDEF tFFTAlgorithmus}
@@ -261,8 +277,8 @@ begin
   setLength(res,len);
   setLength(ims,len);
   loglen:=round(ln(len)/ln(2));
-  setLength(wRe,hlen);
-  setLength(wIm,hlen);
+  setLength(wRe,hLen);
+  setLength(wIm,hLen);
   setLength(perm,len);
   for i:=0 to len-1 do begin
     res[i]:=nan;
@@ -271,7 +287,7 @@ begin
     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
+  for i:=0 to hLen-1 do begin
     wRe[i]:=cos(2*pi*i/len);
     wIm[i]:=sin(2*pi*i/len);
   end;
@@ -308,36 +324,40 @@ procedure tFFTCooleyTukey.laden(invers: boolean; q: pExtended; schritt: longint)
 {$UNDEF fftunitLaden}
 
 {$DEFINE fftunitGetrenntLaden}
-procedure tFFTCooleyTukey.laden(invers: boolean; qRe,qIm: tExtendedArray);
+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;
+  i,dist,absch,wNum,wStep: longint;
   t1,t2:                   extended;
 begin
   dist:=1;
-  wstep:=hlen;
+  wStep:=hLen;
   while dist<len do begin
     absch:=0;
     while absch<len do begin
-      wnum:=0;
+      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];
+        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;
+        wNum:=wNum+wStep;
       end;
       absch:=absch+2*dist;
     end;
-    wstep:=wstep div 2;
+    wStep:=wStep div 2;
     dist:=dist*2;
   end;
 end;
@@ -368,9 +388,9 @@ begin
   end;
 
   preFFT:=tFFTCooleyTukey.create(length(rRe),doGetrennt,doGetrennt);
-  preFFT.laden(false,rRe,rIm);
+  preFFT.laden(false,pExtended(@(rRe[0])),pExtended(@(rIm[0])),1);
   preFFT.ausfuehren;
-  preFFT.speichern(false,sRe,sIm);
+  preFFT.speichern(false,pExtended(@(sRe[0])),pExtended(@(sIm[0])),1);
   preFFT.free;
 
   {$IFDEF fftNormierung}
@@ -427,9 +447,9 @@ begin
   end;
 
   // a[j] := iFFT(a[j])
-  subFFT.laden(true,res,ims);
+  subFFT.laden(true,pExtended(@(res[0])),pExtended(@(ims[0])),1);
   subFFT.ausfuehren;
-  subFFT.speichern(true,res,ims);
+  subFFT.speichern(true,pExtended(@(res[0])),pExtended(@(ims[0])),1);
 
   for i:=0 to length(res)-1 do begin
     res[i]:=res[i];
@@ -444,9 +464,9 @@ begin
   end;
 
   // a[j] := iFFT(a[j])
-  subFFT.laden(false,res,ims);
+  subFFT.laden(false,pExtended(@(res[0])),pExtended(@(ims[0])),1);
   subFFT.ausfuehren;
-  subFFT.speichern(false,res,ims);
+  subFFT.speichern(false,pExtended(@(res[0])),pExtended(@(ims[0])),1);
 
   for i:=0 to length(res)-1 do begin
     res[i]:=res[i];
-- 
cgit v1.2.3-70-g09d2