unit uDFT;
interface
uses Math;
//
// Autor Matze - siehe: https://www.delphipraxis.net/597828-post1.html
//
type
TComplex =
record
re, im: Extended;
end;
TComplexArray =
array of TComplex;
procedure DFT(
var a: TComplexArray);
implementation
function AddC(a, b: TComplex): TComplex;
begin
Result.re := a.re + b.re;
Result.im := a.im + b.im;
end;
function SubC(a, b: TComplex): TComplex;
begin
Result.re := a.re - b.re;
Result.im := a.im - b.im;
end;
function MulC(a, b: TComplex): TComplex;
begin
Result.re := a.re * b.re - a.im * b.im;
Result.im := a.re * b.im + a.im * b.re;
end;
function MakeC(re, im: extended): TComplex;
begin
Result.re := re;
Result.im := im;
end;
procedure shuffle(
var a: TComplexArray; n, lo: Integer);
var I, m: Integer;
b: TComplexArray;
begin
m := n
shr 1;
setlength(b, m);
for I := 0
to m - 1
do
b[i] := a[lo + i];
for I := 0
to m - 1
do
a[lo + i + i + 1] := a[lo + i + m];
for I := 0
to m - 1
do
a[lo + i + i] := b[i];
end;
procedure DoFFT(
var a: TComplexArray; n, lo: Integer; w: TComplex);
var I, m: Integer;
z, v, h: TComplex;
begin
if n
and (n - 1) = 0
then
begin
if n > 1
then
begin
m := n
shr 1;
z := MakeC(1, 0);
for I := lo
to lo + m - 1
do
begin
h := SubC(a[i], a[i + m]);
a[i] := AddC(a[i], a[i + m]);
a[i + m] := MulC(h, z);
z:=MulC(z,w);
end;
v := MulC(w, w);
DoFFT(a, m, lo, v);
DoFFT(a, m, lo + m, v);
shuffle(a, n, lo);
end;
end;
end;
procedure DFT(
var a: TComplexArray);
begin
DoFFT(a, length(a), 0, MakeC(cos(2 * Pi / length(a)),
sin(2 * Pi / length(a))));
end;
end.