1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
|
program testapp;
{$APPTYPE CONSOLE}
uses
SysUtils,
fftreal in 'fftreal.pas',
Math,
Windows;
var
nbr_points : longint;
x, f : pflt_array;
fft : TFFTReal;
i : longint;
PI : double;
areal, img : double;
f_abs : double;
buffer_size : longint;
nbr_tests : longint;
time0, time1, time2 : int64;
timereso : int64;
offset : longint;
t0, t1 : double;
nbr_s_chn : longint;
tempp1, tempp2 : pflt_array;
begin
(*______________________________________________
*
* Exactness test
*______________________________________________
*)
WriteLn('Accuracy test:');
WriteLn;
nbr_points := 16; // Power of 2
GetMem(x, nbr_points * sizeof_flt);
GetMem(f, nbr_points * sizeof_flt);
fft := TFFTReal.Create(nbr_points); // FFT object initialized here
// Test signal
PI := ArcTan(1) * 4;
for i := 0 to nbr_points-1 do
begin
x^[i] := -1 + sin (3*2*PI*i/nbr_points)
+ cos (5*2*PI*i/nbr_points) * 2
- sin (7*2*PI*i/nbr_points) * 3
+ cos (8*2*PI*i/nbr_points) * 5;
end;
// Compute FFT and IFFT
fft.do_fft(f, x);
fft.do_ifft(f, x);
fft.rescale(x);
// Display the result
WriteLn('FFT:');
for i := 0 to nbr_points div 2 do
begin
areal := f^[i];
if (i > 0) and (i < nbr_points div 2) then
img := f^[i + nbr_points div 2]
else
img := 0;
f_abs := Sqrt(areal * areal + img * img);
WriteLn(Format('%5d: %12.6f %12.6f (%12.6f)', [i, areal, img, f_abs]));
end;
WriteLn;
WriteLn('IFFT:');
for i := 0 to nbr_points-1 do
WriteLn(Format('%5d: %f', [i, x^[i]]));
WriteLn;
FreeMem(x);
FreeMem(f);
fft.Free;
(*______________________________________________
*
* Speed test
*______________________________________________
*)
WriteLn('Speed test:');
WriteLn('Please wait...');
WriteLn;
nbr_points := 1024; // Power of 2
buffer_size := 256*nbr_points; // Number of flt_t (float or double)
nbr_tests := 10000;
assert(nbr_points <= buffer_size);
GetMem(x, buffer_size * sizeof_flt);
GetMem(f, buffer_size * sizeof_flt);
fft := TFFTReal.Create(nbr_points); // FFT object initialized here
// Test signal: noise
for i := 0 to nbr_points-1 do
x^[i] := Random($7fff) - ($7fff shr 1);
// timing
QueryPerformanceFrequency(timereso);
QueryPerformanceCounter(time0);
for i := 0 to nbr_tests-1 do
begin
offset := (i * nbr_points) and (buffer_size - 1);
tempp1 := f;
inc(tempp1, offset);
tempp2 := x;
inc(tempp2, offset);
fft.do_fft(tempp1, tempp2);
end;
QueryPerformanceCounter(time1);
for i := 0 to nbr_tests-1 do
begin
offset := (i * nbr_points) and (buffer_size - 1);
tempp1 := f;
inc(tempp1, offset);
tempp2 := x;
inc(tempp2, offset);
fft.do_ifft(tempp1, tempp2);
fft.rescale(x);
end;
QueryPerformanceCounter(time2);
t0 := ((time1-time0) / timereso) / nbr_tests;
t1 := ((time2-time1) / timereso) / nbr_tests;
WriteLn(Format('%d-points FFT : %.0f us.', [nbr_points, t0 * 1000000]));
WriteLn(Format('%d-points IFFT + scaling: %.0f us.', [nbr_points, t1 * 1000000]));
nbr_s_chn := Floor(nbr_points / ((t0 + t1) * 44100 * 2));
WriteLn(Format('Peak performance: FFT+IFFT on %d mono channels at 44.1 KHz (with overlapping)', [nbr_s_chn]));
WriteLn;
FreeMem(x);
FreeMem(f);
fft.Free;
WriteLn('Press [Return] key to terminate...');
ReadLn;
end.
|
