dario mambro | cb97184 | 2020-03-28 00:22:33 +0100 | [diff] [blame] | 1 | /* |
hayati ayguen | 63794b2 | 2020-03-29 03:14:43 +0200 | [diff] [blame] | 2 | Copyright (c) 2013 Julien Pommier ( pommier@modartt.com ) |
| 3 | Copyright (c) 2020 Dario Mambro ( dario.mambro@gmail.com ) |
| 4 | Copyright (c) 2020 Hayati Ayguen ( h_ayguen@web.de ) |
dario mambro | cb97184 | 2020-03-28 00:22:33 +0100 | [diff] [blame] | 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 | |
hayati ayguen | c974c1d | 2020-03-29 03:39:30 +0200 | [diff] [blame] | 44 | /* define own constants required to turn off g++ extensions .. */ |
| 45 | #ifndef M_PI |
| 46 | #define M_PI 3.14159265358979323846 /* pi */ |
| 47 | #endif |
| 48 | |
dario mambro | cb97184 | 2020-03-28 00:22:33 +0100 | [diff] [blame] | 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 |
| 61 | Ttest(int N, bool useOrdered) |
| 62 | { |
hayati ayguen | 5553119 | 2020-04-02 23:11:14 +0200 | [diff] [blame] | 63 | typedef pffft::Fft<T> Fft; |
| 64 | typedef typename pffft::Fft<T>::Scalar FftScalar; |
| 65 | typedef typename Fft::Complex FftComplex; |
dario mambro | cb97184 | 2020-03-28 00:22:33 +0100 | [diff] [blame] | 66 | |
hayati ayguen | 5553119 | 2020-04-02 23:11:14 +0200 | [diff] [blame] | 67 | const bool cplx = pffft::Fft<T>::isComplexTransform(); |
hayati ayguen | 61ec6da | 2020-03-29 08:48:01 +0200 | [diff] [blame] | 68 | const double EXPECTED_DYN_RANGE = Fft::isDoubleScalar() ? 215.0 : 140.0; |
dario mambro | cb97184 | 2020-03-28 00:22:33 +0100 | [diff] [blame] | 69 | |
hayati ayguen | ca8c5f8 | 2020-04-13 01:39:14 +0200 | [diff] [blame] | 70 | assert(Fft::isPowerOfTwo(N)); |
dario mambro | cb97184 | 2020-03-28 00:22:33 +0100 | [diff] [blame] | 71 | |
hayati ayguen | 61ec6da | 2020-03-29 08:48:01 +0200 | [diff] [blame] | 72 | Fft fft = Fft(N); // instantiate and prepareLength() for length N |
dario mambro | cb97184 | 2020-03-28 00:22:33 +0100 | [diff] [blame] | 73 | |
dario mambro | a1cfad4 | 2020-03-30 05:50:26 +0200 | [diff] [blame] | 74 | #if __cplusplus >= 201103L || (defined(_MSC_VER) && _MSC_VER >= 1900) |
| 75 | |
hayati ayguen | 5553119 | 2020-04-02 23:11:14 +0200 | [diff] [blame] | 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) |
dario mambro | a1cfad4 | 2020-03-30 05:50:26 +0200 | [diff] [blame] | 83 | #else |
| 84 | |
hayati ayguen | 5553119 | 2020-04-02 23:11:14 +0200 | [diff] [blame] | 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) |
dario mambro | a1cfad4 | 2020-03-30 05:50:26 +0200 | [diff] [blame] | 91 | #endif |
hayati ayguen | 61ec6da | 2020-03-29 08:48:01 +0200 | [diff] [blame] | 92 | |
| 93 | // work with complex - without the capabilities of a higher c++ standard |
hayati ayguen | 5553119 | 2020-04-02 23:11:14 +0200 | [diff] [blame] | 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; |
dario mambro | cb97184 | 2020-03-28 00:22:33 +0100 | [diff] [blame] | 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) { |
hayati ayguen | ee17cb0 | 2020-08-28 04:47:22 +0200 | [diff] [blame] | 123 | Xs[2 * j] = (FftScalar)( amp * cos(phi) ); /* real part */ |
| 124 | Xs[2 * j + 1] = (FftScalar)( amp * sin(phi) ); /* imag part */ |
dario mambro | cb97184 | 2020-03-28 00:22:33 +0100 | [diff] [blame] | 125 | } else |
hayati ayguen | ee17cb0 | 2020-08-28 04:47:22 +0200 | [diff] [blame] | 126 | Xs[j] = (FftScalar)( amp * cos(phi) ); /* only real part */ |
dario mambro | cb97184 | 2020-03-28 00:22:33 +0100 | [diff] [blame] | 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 { |
hayati ayguen | 1c193e9 | 2020-03-29 16:49:05 +0200 | [diff] [blame] | 138 | fft.forwardToInternalLayout(X, R); /* use R for reordering */ |
| 139 | fft.reorderSpectrum(R, Y); /* have canonical order in Y[] for power calculations */ |
dario mambro | cb97184 | 2020-03-28 00:22:33 +0100 | [diff] [blame] | 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 */ |
hayati ayguen | 61ec6da | 2020-03-29 08:48:01 +0200 | [diff] [blame] | 221 | if (useOrdered) |
| 222 | fft.inverse(Y, Z); |
| 223 | else |
hayati ayguen | 1c193e9 | 2020-03-29 16:49:05 +0200 | [diff] [blame] | 224 | fft.inverseFromInternalLayout(R, Z); /* inverse() from internal Layout */ |
dario mambro | cb97184 | 2020-03-28 00:22:33 +0100 | [diff] [blame] | 225 | |
| 226 | errSum = 0.0; |
| 227 | for (j = 0; j < (cplx ? (2 * N) : N); ++j) { |
| 228 | /* scale back */ |
hayati ayguen | 61ec6da | 2020-03-29 08:48:01 +0200 | [diff] [blame] | 229 | Zs[j] /= N; |
dario mambro | cb97184 | 2020-03-28 00:22:33 +0100 | [diff] [blame] | 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 | } |
hayati ayguen | 1c193e9 | 2020-03-29 16:49:05 +0200 | [diff] [blame] | 249 | |
| 250 | // using the std::vector<> base classes .. no need for alignedFree() for X, Y, Z and R |
dario mambro | cb97184 | 2020-03-28 00:22:33 +0100 | [diff] [blame] | 251 | |
| 252 | return retError; |
| 253 | } |
| 254 | |
| 255 | bool |
| 256 | test(int N, bool useComplex, bool useOrdered) |
| 257 | { |
| 258 | if (useComplex) { |
hayati ayguen | ca11241 | 2020-04-13 00:19:40 +0200 | [diff] [blame] | 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 | ; |
dario mambro | cb97184 | 2020-03-28 00:22:33 +0100 | [diff] [blame] | 270 | } else { |
hayati ayguen | ca11241 | 2020-04-13 00:19:40 +0200 | [diff] [blame] | 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 | ; |
dario mambro | cb97184 | 2020-03-28 00:22:33 +0100 | [diff] [blame] | 282 | } |
| 283 | } |
| 284 | |
| 285 | int |
| 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) { |
hayati ayguen | ca11241 | 2020-04-13 00:19:40 +0200 | [diff] [blame] | 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 |
dario mambro | cb97184 | 2020-03-28 00:22:33 +0100 | [diff] [blame] | 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 | |
hayati ayguen | ca11241 | 2020-04-13 00:19:40 +0200 | [diff] [blame] | 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 |
dario mambro | cb97184 | 2020-03-28 00:22:33 +0100 | [diff] [blame] | 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) |
hayati ayguen | ca11241 | 2020-04-13 00:19:40 +0200 | [diff] [blame] | 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"); |
dario mambro | cb97184 | 2020-03-28 00:22:33 +0100 | [diff] [blame] | 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 | } |