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.