FFT algorithm for Delphi 2

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;

function Sum(x, y: TComplex) : TComplex;

begin

with result do

begin

r:= x.r + y.r;

i:= x.i + y.i;

end;

function Difference(x, y: TComplex) : TComplex;

begin

with result do

begin

r:= x.r - y.r;

i:= x.i - y.i;

end;

function EiT(t: TReal): TComplex;

begin

with result do

begin

r:= cos(t);

i:= sin(t);

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;

function TimesReal(x: TComplex; y: TReal): TComplex;

begin

with result do

begin

r:= x.r * y;

i:= x.i * y;

end;

function PlusReal(x: TComplex; y: TReal): TComplex;

begin

with result do

begin

r:= x.r + y;

i:= x.i;

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;

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;

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;

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.