ezfit:fft_algorithm

FFT Algorithm

Source: Vitalie Eşanu, FFT Algorithm for Delphi 2, http://www.esanu.name/delphi/Algorithms/Maths/FFT%20algorithm%20for%20Delphi%202.html, Oct 2008.

FFT Algorithm in Delphi2

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 recursion 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 machine 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.