FFT algorithm for Delphi 2
David Ullrich <ullrich@math.okstate.edu>
Here's an FFT that handles 256 data points in about 0.008 seconds on a P66 (with 72MB, YMMV). Nothing but Delphi.
This one came out a lot nicer than the one I did a year ago. It's probably not optimal; if we want an optimal FFT we have to buy hardware, what the heck.
But I don't think it's too bad, performance-wise. There's a little bit of recursion involved, but the recursion doesn't involve copying any data, just a few pointers; if we have an array of length N = 2^d then the depth of the recursion is just d. Possibly it could be improved by unwrapping the recursion, it's not clear whether it would be worth the trouble. (But probably a person could get substantial inprovement with relatively little effort by unwrapping the bottom layer or two of the recursion, ie by saying
--------------------------------------------------------------------------------
if Depth < 2 then
{do what needs to be done}
--------------------------------------------------------------------------------
instead of the current 'if Depth = 0 then...' This would eliminate function calls that do nothing but make assignments, a good thing, while OTOH unwrapping all of the resursion would be trickier, and wouldn't seem as productive, since most of the function calls that would be eliminated do much more than just an assignment.)
There's a lookup table used for the sines and cosines; it could be that this is the wrong way to do it for large arrays, seems to work just fine for small to medium arrays.
Probably on a mchine with a lot of RAM a person would use VirtualAlloc(... PAGE_NOCACHE) for Src, Dest, and the lookup tables.
If anybody notices anything stupid about the way something's done not mentioned above please mention it.
What does it do, exactly? There are FFT's and FFT's - this one does the 'complex FT', that being the one I understand and care about. By this I mean that if N = 2^d and Src^ and Dest^ are arrays of N TComplexes, then a call
--------------------------------------------------------------------------------
FFT(d, Src, Dest)
--------------------------------------------------------------------------------
will fill in Dest with the complex FT: after the call Dest^[j] will equal
--------------------------------------------------------------------------------
1/sqrt(N) * Sum(k=0.. N - 1 ; EiT(2*Pi(j*k/N)) * Src^[k])
--------------------------------------------------------------------------------
, where EiT(t) = cos(t) + i sin(t) . Ie, the standard Fourier Transform.
Comes in two versions: In the first version I use a TComplex, with functions to manipulate the complex numbers. In the second version everything's real - instead of arrays Src and Dest of complexes we have arrays SrcR, SrcI, DestR, DestI of reals (for the real and imagionary parts), and all those function calls are written out inline. The first one is much easier for me to make sense of, the second version is much faster. (They both give the 'complex FFT'.) With little programs that test whether it does what it should by checking Plancherel (aka Parseval). It really does work, btw - if it doesn't work for you it's because I garbled something in the process of deleting stupid comments. The complex version:
--------------------------------------------------------------------------------
***
unit cplx;
interface
type
PReal = ^TReal;
TReal = extended;
PComplex = ^TComplex;
TComplex = record
r : TReal;
i : TReal;
end;
function MakeComplex(x, y: TReal): TComplex;
function Sum(x, y: TComplex) : TComplex;
function Difference(x, y: TComplex) : TComplex;
function Product(x, y: TComplex): TComplex;
function TimesReal(x: TComplex; y: TReal): TComplex;
function PlusReal(x: TComplex; y: TReal): TComplex;
function EiT(t: TReal):TComplex;
function ComplexToStr(x: TComplex): string;
function AbsSquared(x: TComplex): TReal;
implementation
uses SysUtils;
function MakeComplex(x, y: TReal): TComplex;
begin
with result do
begin
r:=x;
i:= y;
end;
end;
function Sum(x, y: TComplex) : TComplex;
begin
with result do
begin
r:= x.r + y.r;
i:= x.i + y.i;
end;
end;
function Difference(x, y: TComplex) : TComplex;
begin
with result do
begin
r:= x.r - y.r;
i:= x.i - y.i;
end;
end;
function EiT(t: TReal): TComplex;
begin
with result do
begin
r:= cos(t);
i:= sin(t);
end;
end;
function Product(x, y: TComplex): TComplex;
begin
with result do
begin
r:= x.r * y.r - x.i * y.i;
i:= x.r * y.i + x.i * y.r;
end;
end;
function TimesReal(x: TComplex; y: TReal): TComplex;
begin
with result do
begin
r:= x.r * y;
i:= x.i * y;
end;
end;
function PlusReal(x: TComplex; y: TReal): TComplex;
begin
with result do
begin
r:= x.r + y;
i:= x.i;
end;
end;
function ComplexToStr(x: TComplex): string;
begin
result:= FloatToStr(x.r)
+ ' + '
+ FloatToStr(x.i)
+ 'i';
end;
function AbsSquared(x: TComplex): TReal;
begin
result:= x.r*x.r + x.i*x.i;
end;
end.
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
unit cplxfft1;
interface
uses Cplx;
type
PScalar = ^TScalar;
TScalar = TComplex; {Making conversion to real version easier}
PScalars = ^TScalars;
TScalars = array[0..High(integer) div SizeOf(TScalar) - 1]
of TScalar;
const
TrigTableDepth: word = 0;
TrigTable : PScalars = nil;
procedure InitTrigTable(Depth: word);
procedure FFT(Depth: word;
Src: PScalars;
Dest: PScalars);
{REQUIRES allocating
(integer(1) shl Depth) * SizeOf(TScalar)
bytes for Src and Dest before call!}
implementation
procedure DoFFT(Depth: word;
Src: PScalars;
SrcSpacing: word;
Dest: PScalars);
{the recursive part called by FFT when ready}
var j, N: integer; Temp: TScalar; Shift: word;
begin
if Depth = 0 then
begin
Dest^[0]:= Src^[0];
exit;
end;
N:= integer(1) shl (Depth - 1);
DoFFT(Depth - 1, Src, SrcSpacing * 2, Dest);
DoFFT(Depth - 1, @Src^[SrcSpacing], SrcSpacing * 2, @Dest^[N] );
Shift:= TrigTableDepth - Depth;
for j:= 0 to N - 1 do
begin
Temp:= Product(TrigTable^[j shl Shift],
Dest^[j + N]);
Dest^[j + N]:= Difference(Dest^[j], Temp);
Dest^[j] := Sum(Dest^[j], Temp);
end;
end;
procedure FFT(Depth: word;
Src: PScalars;
Dest: PScalars);
var j, N: integer; Normalizer: extended;
begin
N:= integer(1) shl depth;
if Depth TrigTableDepth then
InitTrigTable(Depth);
DoFFT(Depth, Src, 1, Dest);
Normalizer:= 1 / sqrt(N) ;
for j:=0 to N - 1 do
Dest^[j]:= TimesReal(Dest^[j], Normalizer);
end;
procedure InitTrigTable(Depth: word);
var j, N: integer;
begin
N:= integer(1) shl depth;
ReAllocMem(TrigTable, N * SizeOf(TScalar));
for j:=0 to N - 1 do
TrigTable^[j]:= EiT(-(2*Pi)*j/N);
TrigTableDepth:= Depth;
end;
initialization
;
finalization
ReAllocMem(TrigTable, 0);
end.
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
unit DemoForm;
interface
uses
Windows, Messages, SysUtils, Classes, Graphics, Controls, Forms, Dialogs,
StdCtrls;
type
TForm1 = class(TForm)
Button1: TButton;
Memo1: TMemo;
Edit1: TEdit;
Label1: TLabel;
procedure Button1Click(Sender: TObject);
private
{ Private declarations }
public
{ Public declarations }
end;
var
Form1: TForm1;
implementation
{$R *.DFM}
uses cplx, cplxfft1, MMSystem;
procedure TForm1.Button1Click(Sender: TObject);
var j: integer; s:string;
src, dest: PScalars;
norm: extended;
d,N,count:integer;
st,et: longint;
begin
d:= StrToIntDef(edit1.text, -1) ;
if d <1 then
raise exception.Create('depth must be a positive integer');
N:= integer(1) shl d ;
GetMem(Src, N*Sizeof(TScalar));
GetMem(Dest, N*SizeOf(TScalar));
for j:=0 to N-1 do
begin
src^[j]:= MakeComplex(random, random);
end;
begin
st:= timeGetTime;
FFT(d, Src, dest);
et:= timeGetTime;
end;
Memo1.Lines.Add('N = ' + IntToStr(N));
Memo1.Lines.Add('expected norm: ' +#9+ FloatToStr(N*2/3));
norm:=0;
for j:=0 to N-1 do norm:= norm + AbsSquared(src^[j]);
Memo1.Lines.Add('Data norm: '+#9+FloatToStr(norm));
norm:=0;
for j:=0 to N-1 do norm:= norm + AbsSquared(dest^[j]);
Memo1.Lines.Add('FT norm: '+#9#9+FloatToStr(norm));
Memo1.Lines.Add('Time in FFT routine: '+#9
+ inttostr(et - st)
+ ' ms.');
Memo1.Lines.Add(' ');
FreeMem(Src);
FreeMem(DEst);
end;
end.
--------------------------------------------------------------------------------
**** The real version:
****
--------------------------------------------------------------------------------
unit cplxfft2;
interface
type
PScalar = ^TScalar;
TScalar = extended;
PScalars = ^TScalars;
TScalars = array[0..High(integer) div SizeOf(TScalar) - 1]
of TScalar;
const
TrigTableDepth: word = 0;
CosTable : PScalars = nil;
SinTable : PScalars = nil;
procedure InitTrigTables(Depth: word);
procedure FFT(Depth: word;
SrcR, SrcI: PScalars;
DestR, DestI: PScalars);
{REQUIRES allocating
(integer(1) shl Depth) * SizeOf(TScalar)
bytes for SrcR, SrcI, DestR and DestI before call!}
implementation
procedure DoFFT(Depth: word;
SrcR, SrcI: PScalars;
SrcSpacing: word;
DestR, DestI: PScalars);
{the recursive part called by FFT when ready}
var j, N: integer;
TempR, TempI: TScalar;
Shift: word;
c, s: extended;
begin
if Depth = 0 then
begin
DestR^[0]:= SrcR^[0];
DestI^[0]:= SrcI^[0];
exit;
end;
N:= integer(1) shl (Depth - 1);
DoFFT(Depth - 1, SrcR, SrcI, SrcSpacing * 2, DestR, DestI);
DoFFT(Depth - 1,
@SrcR^[srcSpacing],
@SrcI^[SrcSpacing],
SrcSpacing * 2,
@DestR^[N],
@DestI^[N]);
Shift:= TrigTableDepth - Depth;
for j:= 0 to N - 1 do
begin
c:= CosTable^[j shl Shift];
s:= SinTable^[j shl Shift];
TempR:= c * DestR^[j + N] - s * DestI^[j + N];
TempI:= c * DestI^[j + N] + s * DestR^[j + N];
DestR^[j + N]:= DestR^[j] - TempR;
DestI^[j + N]:= DestI^[j] - TempI;
DestR^[j]:= DestR^[j] + TempR;
DestI^[j]:= DestI^[j] + TempI;
end;
end;
procedure FFT(Depth: word;
SrcR, SrcI: PScalars;
DestR, DestI: PScalars);
var j, N: integer; Normalizer: extended;
begin
N:= integer(1) shl depth;
if Depth TrigTableDepth then
InitTrigTables(Depth);
DoFFT(Depth, SrcR, SrcI, 1, DestR, DestI);
Normalizer:= 1 / sqrt(N) ;
for j:=0 to N - 1 do
begin
DestR^[j]:= DestR^[j] * Normalizer;
DestI^[j]:= DestI^[j] * Normalizer;
end;
end;
procedure InitTrigTables(Depth: word);
var j, N: integer;
begin
N:= integer(1) shl depth;
ReAllocMem(CosTable, N * SizeOf(TScalar));
ReAllocMem(SinTable, N * SizeOf(TScalar));
for j:=0 to N - 1 do
begin
CosTable^[j]:= cos(-(2*Pi)*j/N);
SinTable^[j]:= sin(-(2*Pi)*j/N);
end;
TrigTableDepth:= Depth;
end;
initialization
;
finalization
ReAllocMem(CosTable, 0);
ReAllocMem(SinTable, 0);
end.
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
unit demofrm;
interface
uses
Windows, Messages, SysUtils, Classes, Graphics,
Controls, Forms, Dialogs, cplxfft2, StdCtrls;
type
TForm1 = class(TForm)
Button1: TButton;
Memo1: TMemo;
Edit1: TEdit;
Label1: TLabel;
procedure Button1Click(Sender: TObject);
private
{ Private declarations }
public
{ Public declarations }
end;
var
Form1: TForm1;
implementation
{$R *.DFM}
uses MMSystem;
procedure TForm1.Button1Click(Sender: TObject);
var SR, SI, DR, DI: PScalars;
j,d,N:integer;
st, et: longint;
norm: extended;
begin
d:= StrToIntDef(edit1.text, -1) ;
if d <1 then
raise exception.Create('depth must be a positive integer');
N:= integer(1) shl d;
GetMem(SR, N * SizeOf(TScalar));
GetMem(SI, N * SizeOf(TScalar));
GetMem(DR, N * SizeOf(TScalar));
GetMem(DI, N * SizeOf(TScalar));
for j:=0 to N - 1 do
begin
SR^[j]:=random;
SI^[j]:=random;
end;
st:= timeGetTime;
FFT(d, SR, SI, DR, DI);
et:= timeGetTime;
memo1.Lines.Add('N = '+inttostr(N));
memo1.Lines.Add('expected norm: '+#9+FloatToStr(N*2/3));
norm:=0;
for j:=0 to N - 1 do
norm:= norm + SR^[j]*SR^[j] + SI^[j]*SI^[j];
memo1.Lines.Add('Data norm: '+#9+FloatToStr(norm));
norm:=0;
for j:=0 to N - 1 do
norm:= norm + DR^[j]*DR^[j] + DI^[j]*DI^[j];
memo1.Lines.Add('FT norm: '+#9#9+FloatToStr(norm));
memo1.Lines.Add('Time in FFT routine: '+#9+inttostr(et-st));
memo1.Lines.add('');
(*for j:=0 to N - 1 do
Memo1.Lines.Add(FloatToStr(SR^[j])
+ ' + '
+ FloatToStr(SI^[j])
+ 'i');
for j:=0 to N - 1 do
Memo1.Lines.Add(FloatToStr(DR^[j])
+ ' + '
+ FloatToStr(DI^[j])
+ 'i');*)
FreeMem(SR, N * SizeOf(TScalar));
FreeMem(SI, N * SizeOf(TScalar));
FreeMem(DR, N * SizeOf(TScalar));
FreeMem(DI, N * SizeOf(TScalar));
end;
end.