program myfft;
{$mode objfpc}{$I-}{$R-}{$R-}{$C-}
uses
windows, sysutils, math;
type
TComplex = packed record
Re, Im: Double;
end;
TComplexArray = array of TComplex;
function Complex(const Re, Im: Double): TComplex;inline;
begin
Result.Re := Re;
Result.Im := Im;
end;
operator + (const a,b:TComplex):TComplex;inline;
begin
Result.Re := A.Re + B.Re;
Result.Im := A.Im + B.Im;
end;
operator - (const A,B:TComplex):TComplex;inline;
begin
Result.Re := A.Re - B.Re;
Result.Im := A.Im - B.Im;
end;
operator * (const A,B:TComplex):TComplex;inline;
begin
Result.Re := A.Re * B.Re - A.Im * B.Im;
Result.Im := A.Re * B.Im + A.Im * B.Re;
end;
procedure FFT(var Data: TComplexArray; Inverse: Boolean);inline;
var
i, j, k, n, m: Integer;
w, wm, t, u: TComplex;
theta, sign: Double;
begin
n := Length(Data);
if n <= 1 then Exit;
{ Bit-reversal permutation }
j := 0;
for i := 0 to n - 1 do
begin
if i < j then
begin
t := Data[i];
Data[i] := Data[j];
Data[j] := t;
end;
k := n shr 1;
while (j >= k) and (k > 0) do
begin
Dec(j, k);
k := k shr 1;
end;
Inc(j, k);
end;
{ Set sign for forward/inverse FFT }
if Inverse then
sign := -1.0
else
sign := 1.0;
{ Iterative FFT }
m := 2;
while m <= n do
begin
theta := sign * 2.0 * Pi / m;
wm := Complex(Cos(theta), Sin(theta));
for k := 0 to (n div m) - 1 do
begin
w := Complex(1.0, 0.0);
for j := 0 to (m div 2) - 1 do
begin
u := Data[k*m + j];
t := w * Data[k*m + j + m div 2];
Data[k*m+j]:= u + t;
Data[k*m + j + m div 2] := u - t;
w := w * wm;
end;
end;
m := m shl 1;
end;
{ Scale for inverse FFT }
if Inverse then
for i := 0 to n - 1 do
begin
Data[i].Re := Data[i].Re / n;
Data[i].Im := Data[i].Im / n;
end;
end;
type
{ must be power of two }
range = 0..16383;
var
Data: TComplexArray = nil;
i: Integer;
time0, time1, time2 : int64;
timereso : int64;
t0, t1 : double;
begin
SetLength(Data, high(range)+1);
(*
{ Test signal: noise }
for i in range do
Data[i].re := Random($7fff) - ($7fff shr 1);
{ Test signal }
for i in range do
begin
Data[i].re := -1 + sin (3*2*PI*i/ (high(range)+1))
+ cos (5*2*PI*i/(high(range)+1)) * 2
- sin (7*2*PI*i/(high(range)+1)) * 3
+ cos (8*2*PI*i/(high(range)+1)) * 5;
end;
*)
{ Test signal (sine wave + noise) }
for i in range do
data[i].re:=sin(2*Pi*i/(high(range)+1))+((Random(100)/100.0)-0.5);
(*
for i in range do
WriteLn(Format('%.2f + %.2fi', [Data[i].Re, Data[i].Im]));
writeln;
*)
// timing
QueryPerformanceFrequency(timereso);
QueryPerformanceCounter(time0);
{ forward FFT }
FFT(Data,false);
QueryPerformanceCounter(time1);
{ inverse FFT }
FFT(Data,true);
QueryPerformanceCounter(time2);
t0 := ((time1-time0) / timereso) / (high(range) + 1);
t1 := ((time2-time1) / timereso) / (high(range) + 1);
WriteLn(Format('%d-points FFT in nanoseconds: %.0f ns.', [high(range)+1, t0 * 1000000000]));
WriteLn(Format('%d-points IFFT + scaling : %.0f ns.', [high(range)+1, t1 * 1000000000]));
(*
{ accuracy check }
for i in range do
WriteLn(Format('%.2f + %.2fi', [Data[i].Re, Data[i].Im]));
*)
end.
Recent
