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