1 /*
2 Copyright (c) 2013 Julien Pommier ( [email protected] )
3 Copyright (c) 2020 Dario Mambro ( [email protected] )
4 Copyright (c) 2020 Hayati Ayguen ( [email protected] )
5
6 Small test & bench for PFFFT, comparing its performance with the scalar
7 FFTPACK, FFTW, and Apple vDSP
8
9 How to build:
10
11 on linux, with fftw3:
12 gcc -o test_pffft -DHAVE_FFTW -msse -mfpmath=sse -O3 -Wall -W pffft.c
13 test_pffft.c fftpack.c -L/usr/local/lib -I/usr/local/include/ -lfftw3f -lm
14
15 on macos, without fftw3:
16 clang -o test_pffft -DHAVE_VECLIB -O3 -Wall -W pffft.c test_pffft.c fftpack.c
17 -L/usr/local/lib -I/usr/local/include/ -framework Accelerate
18
19 on macos, with fftw3:
20 clang -o test_pffft -DHAVE_FFTW -DHAVE_VECLIB -O3 -Wall -W pffft.c
21 test_pffft.c fftpack.c -L/usr/local/lib -I/usr/local/include/ -lfftw3f
22 -framework Accelerate
23
24 as alternative: replace clang by gcc.
25
26 on windows, with visual c++:
27 cl /Ox -D_USE_MATH_DEFINES /arch:SSE test_pffft.c pffft.c fftpack.c
28
29 build without SIMD instructions:
30 gcc -o test_pffft -DPFFFT_SIMD_DISABLE -O3 -Wall -W pffft.c test_pffft.c
31 fftpack.c -lm
32
33 */
34
35 #include "pffft.hpp"
36
37 #include <assert.h>
38 #include <math.h>
39 #include <stdio.h>
40 #include <stdlib.h>
41 #include <string.h>
42 #include <time.h>
43
44 /* define own constants required to turn off g++ extensions .. */
45 #ifndef M_PI
46 #define M_PI 3.14159265358979323846 /* pi */
47 #endif
48
49 /* maximum allowed phase error in degree */
50 #define DEG_ERR_LIMIT 1E-4
51
52 /* maximum allowed magnitude error in amplitude (of 1.0 or 1.1) */
53 #define MAG_ERR_LIMIT 1E-6
54
55 #define PRINT_SPEC 0
56
57 #define PWR2LOG(PWR) ((PWR) < 1E-30 ? 10.0 * log10(1E-30) : 10.0 * log10(PWR))
58
59 template<typename T>
60 bool
Ttest(int N,bool useOrdered)61 Ttest(int N, bool useOrdered)
62 {
63 typedef pffft::Fft<T> Fft;
64 typedef typename pffft::Fft<T>::Scalar FftScalar;
65 typedef typename Fft::Complex FftComplex;
66
67 const bool cplx = pffft::Fft<T>::isComplexTransform();
68 const double EXPECTED_DYN_RANGE = Fft::isDoubleScalar() ? 215.0 : 140.0;
69
70 assert(Fft::isPowerOfTwo(N));
71
72 Fft fft = Fft(N); // instantiate and prepareLength() for length N
73
74 #if __cplusplus >= 201103L || (defined(_MSC_VER) && _MSC_VER >= 1900)
75
76 // possible ways to declare/instatiate aligned vectors with C++11
77 // some lines require a typedef of above
78 auto X = fft.valueVector(); // for X = input vector
79 pffft::AlignedVector<typename Fft::Complex> Y = fft.spectrumVector(); // for Y = forward(X)
80 pffft::AlignedVector<FftScalar> R = fft.internalLayoutVector(); // for R = forwardInternalLayout(X)
81 pffft::AlignedVector<T> Z = fft.valueVector(); // for Z = inverse(Y) = inverse( forward(X) )
82 // or Z = inverseInternalLayout(R)
83 #else
84
85 // possible ways to declare/instatiate aligned vectors with C++98
86 pffft::AlignedVector<T> X = fft.valueVector(); // for X = input vector
87 pffft::AlignedVector<FftComplex> Y = fft.spectrumVector(); // for Y = forward(X)
88 pffft::AlignedVector<typename Fft::Scalar> R = fft.internalLayoutVector(); // for R = forwardInternalLayout(X)
89 pffft::AlignedVector<T> Z = fft.valueVector(); // for Z = inverse(Y) = inverse( forward(X) )
90 // or Z = inverseInternalLayout(R)
91 #endif
92
93 // work with complex - without the capabilities of a higher c++ standard
94 FftScalar* Xs = reinterpret_cast<FftScalar*>(X.data()); // for X = input vector
95 FftScalar* Ys = reinterpret_cast<FftScalar*>(Y.data()); // for Y = forward(X)
96 FftScalar* Zs = reinterpret_cast<FftScalar*>(Z.data()); // for Z = inverse(Y) = inverse( forward(X) )
97
98 int k, j, m, iter, kmaxOther;
99 bool retError = false;
100 double freq, dPhi, phi, phi0;
101 double pwr, pwrCar, pwrOther, err, errSum, mag, expextedMag;
102 double amp = 1.0;
103
104 for (k = m = 0; k < (cplx ? N : (1 + N / 2)); k += N / 16, ++m) {
105 amp = ((m % 3) == 0) ? 1.0F : 1.1F;
106 freq = (k < N / 2) ? ((double)k / N) : ((double)(k - N) / N);
107 dPhi = 2.0 * M_PI * freq;
108 if (dPhi < 0.0)
109 dPhi += 2.0 * M_PI;
110
111 iter = -1;
112 while (1) {
113 ++iter;
114
115 if (iter)
116 printf("bin %d: dphi = %f for freq %f\n", k, dPhi, freq);
117
118 /* generate cosine carrier as time signal - start at defined phase phi0 */
119 phi = phi0 =
120 (m % 4) * 0.125 * M_PI; /* have phi0 < 90 deg to be normalized */
121 for (j = 0; j < N; ++j) {
122 if (cplx) {
123 Xs[2 * j] = (FftScalar)( amp * cos(phi) ); /* real part */
124 Xs[2 * j + 1] = (FftScalar)( amp * sin(phi) ); /* imag part */
125 } else
126 Xs[j] = (FftScalar)( amp * cos(phi) ); /* only real part */
127
128 /* phase increment .. stay normalized - cos()/sin() might degrade! */
129 phi += dPhi;
130 if (phi >= M_PI)
131 phi -= 2.0 * M_PI;
132 }
133
134 /* forward transform from X --> Y .. using work buffer W */
135 if (useOrdered)
136 fft.forward(X, Y);
137 else {
138 fft.forwardToInternalLayout(X, R); /* use R for reordering */
139 fft.reorderSpectrum(R, Y); /* have canonical order in Y[] for power calculations */
140 }
141
142 pwrOther = -1.0;
143 pwrCar = 0;
144
145 /* for positive frequencies: 0 to 0.5 * samplerate */
146 /* and also for negative frequencies: -0.5 * samplerate to 0 */
147 for (j = 0; j < (cplx ? N : (1 + N / 2)); ++j) {
148 if (!cplx && !j) /* special treatment for DC for real input */
149 pwr = Ys[j] * Ys[j];
150 else if (!cplx && j == N / 2) /* treat 0.5 * samplerate */
151 pwr = Ys[1] *
152 Ys[1]; /* despite j (for freq calculation) we have index 1 */
153 else
154 pwr = Ys[2 * j] * Ys[2 * j] + Ys[2 * j + 1] * Ys[2 * j + 1];
155 if (iter || PRINT_SPEC)
156 printf("%s fft %d: pwr[j = %d] = %g == %f dB\n",
157 (cplx ? "cplx" : "real"),
158 N,
159 j,
160 pwr,
161 PWR2LOG(pwr));
162 if (k == j)
163 pwrCar = pwr;
164 else if (pwr > pwrOther) {
165 pwrOther = pwr;
166 kmaxOther = j;
167 }
168 }
169
170 if (PWR2LOG(pwrCar) - PWR2LOG(pwrOther) < EXPECTED_DYN_RANGE) {
171 printf("%s fft %d amp %f iter %d:\n",
172 (cplx ? "cplx" : "real"),
173 N,
174 amp,
175 iter);
176 printf(" carrier power at bin %d: %g == %f dB\n",
177 k,
178 pwrCar,
179 PWR2LOG(pwrCar));
180 printf(" carrier mag || at bin %d: %g\n", k, sqrt(pwrCar));
181 printf(" max other pwr at bin %d: %g == %f dB\n",
182 kmaxOther,
183 pwrOther,
184 PWR2LOG(pwrOther));
185 printf(" dynamic range: %f dB\n\n",
186 PWR2LOG(pwrCar) - PWR2LOG(pwrOther));
187 retError = true;
188 if (iter == 0)
189 continue;
190 }
191
192 if (k > 0 && k != N / 2) {
193 phi = atan2(Ys[2 * k + 1], Ys[2 * k]);
194 if (fabs(phi - phi0) > DEG_ERR_LIMIT * M_PI / 180.0) {
195 retError = true;
196 printf("%s fft %d bin %d amp %f : phase mismatch! phase = %f deg "
197 "expected = %f deg\n",
198 (cplx ? "cplx" : "real"),
199 N,
200 k,
201 amp,
202 phi * 180.0 / M_PI,
203 phi0 * 180.0 / M_PI);
204 }
205 }
206
207 expextedMag = cplx ? amp : ((k == 0 || k == N / 2) ? amp : (amp / 2));
208 mag = sqrt(pwrCar) / N;
209 if (fabs(mag - expextedMag) > MAG_ERR_LIMIT) {
210 retError = true;
211 printf("%s fft %d bin %d amp %f : mag = %g expected = %g\n",
212 (cplx ? "cplx" : "real"),
213 N,
214 k,
215 amp,
216 mag,
217 expextedMag);
218 }
219
220 /* now convert spectrum back */
221 if (useOrdered)
222 fft.inverse(Y, Z);
223 else
224 fft.inverseFromInternalLayout(R, Z); /* inverse() from internal Layout */
225
226 errSum = 0.0;
227 for (j = 0; j < (cplx ? (2 * N) : N); ++j) {
228 /* scale back */
229 Zs[j] /= N;
230 /* square sum errors over real (and imag parts) */
231 err = (Xs[j] - Zs[j]) * (Xs[j] - Zs[j]);
232 errSum += err;
233 }
234
235 if (errSum > N * 1E-7) {
236 retError = true;
237 printf("%s fft %d bin %d : inverse FFT doesn't match original signal! "
238 "errSum = %g ; mean err = %g\n",
239 (cplx ? "cplx" : "real"),
240 N,
241 k,
242 errSum,
243 errSum / N);
244 }
245
246 break;
247 }
248 }
249
250 // using the std::vector<> base classes .. no need for alignedFree() for X, Y, Z and R
251
252 return retError;
253 }
254
255 bool
test(int N,bool useComplex,bool useOrdered)256 test(int N, bool useComplex, bool useOrdered)
257 {
258 if (useComplex) {
259 return
260 #ifdef PFFFT_ENABLE_FLOAT
261 Ttest< std::complex<float> >(N, useOrdered)
262 #endif
263 #if defined(PFFFT_ENABLE_FLOAT) && defined(PFFFT_ENABLE_DOUBLE)
264 &&
265 #endif
266 #ifdef PFFFT_ENABLE_DOUBLE
267 Ttest< std::complex<double> >(N, useOrdered)
268 #endif
269 ;
270 } else {
271 return
272 #ifdef PFFFT_ENABLE_FLOAT
273 Ttest<float>(N, useOrdered)
274 #endif
275 #if defined(PFFFT_ENABLE_FLOAT) && defined(PFFFT_ENABLE_DOUBLE)
276 &&
277 #endif
278 #ifdef PFFFT_ENABLE_DOUBLE
279 Ttest<double>(N, useOrdered)
280 #endif
281 ;
282 }
283 }
284
285 int
main(int argc,char ** argv)286 main(int argc, char** argv)
287 {
288 int N, result, resN, resAll, k, resNextPw2, resIsPw2, resFFT;
289
290 int inp_power_of_two[] = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 511, 512, 513 };
291 int ref_power_of_two[] = { 1, 2, 4, 4, 8, 8, 8, 8, 16, 512, 512, 1024 };
292
293 resNextPw2 = 0;
294 resIsPw2 = 0;
295 for (k = 0; k < (sizeof(inp_power_of_two) / sizeof(inp_power_of_two[0]));
296 ++k) {
297 #ifdef PFFFT_ENABLE_FLOAT
298 N = pffft::Fft<float>::nextPowerOfTwo(inp_power_of_two[k]);
299 #else
300 N = pffft::Fft<double>::nextPowerOfTwo(inp_power_of_two[k]);
301 #endif
302 if (N != ref_power_of_two[k]) {
303 resNextPw2 = 1;
304 printf("pffft_next_power_of_two(%d) does deliver %d, which is not "
305 "reference result %d!\n",
306 inp_power_of_two[k],
307 N,
308 ref_power_of_two[k]);
309 }
310
311 #ifdef PFFFT_ENABLE_FLOAT
312 result = pffft::Fft<float>::isPowerOfTwo(inp_power_of_two[k]);
313 #else
314 result = pffft::Fft<double>::isPowerOfTwo(inp_power_of_two[k]);
315 #endif
316 if (inp_power_of_two[k] == ref_power_of_two[k]) {
317 if (!result) {
318 resIsPw2 = 1;
319 printf("pffft_is_power_of_two(%d) delivers false; expected true!\n",
320 inp_power_of_two[k]);
321 }
322 } else {
323 if (result) {
324 resIsPw2 = 1;
325 printf("pffft_is_power_of_two(%d) delivers true; expected false!\n",
326 inp_power_of_two[k]);
327 }
328 }
329 }
330 if (!resNextPw2)
331 printf("tests for pffft_next_power_of_two() succeeded successfully.\n");
332 if (!resIsPw2)
333 printf("tests for pffft_is_power_of_two() succeeded successfully.\n");
334
335 resFFT = 0;
336 for (N = 32; N <= 65536; N *= 2) {
337 result = test(N, 1 /* cplx fft */, 1 /* useOrdered */);
338 resN = result;
339 resFFT |= result;
340
341 result = test(N, 0 /* cplx fft */, 1 /* useOrdered */);
342 resN |= result;
343 resFFT |= result;
344
345 result = test(N, 1 /* cplx fft */, 0 /* useOrdered */);
346 resN |= result;
347 resFFT |= result;
348
349 result = test(N, 0 /* cplx fft */, 0 /* useOrdered */);
350 resN |= result;
351 resFFT |= result;
352
353 if (!resN)
354 printf("tests for size %d succeeded successfully.\n", N);
355 }
356
357 if (!resFFT)
358 printf("all pffft transform tests (FORWARD/BACKWARD, REAL/COMPLEX, "
359 #ifdef PFFFT_ENABLE_FLOAT
360 "float"
361 #endif
362 #if defined(PFFFT_ENABLE_FLOAT) && defined(PFFFT_ENABLE_DOUBLE)
363 "/"
364 #endif
365 #ifdef PFFFT_ENABLE_DOUBLE
366 "double"
367 #endif
368 ") succeeded successfully.\n");
369
370 resAll = resNextPw2 | resIsPw2 | resFFT;
371 if (!resAll)
372 printf("all tests succeeded successfully.\n");
373 else
374 printf("there are failed tests!\n");
375
376 return resAll;
377 }
378