Mercurial > hg > qm-dsp

annotate tests/TestFFT.cpp @ 384:279991b6ebe7

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Complex fft tests
author Chris Cannam <c.cannam@qmul.ac.uk>
date Fri, 01 Nov 2013 17:11:50 +0000
parents 3c7338aff6a8
children 8c86761a5533
rev   line source
c@343 1 /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */
c@335 2
c@335 3 #include "dsp/transforms/FFT.h"
c@335 4
c@335 5 #define BOOST_TEST_DYN_LINK
c@335 6 #define BOOST_TEST_MAIN
c@335 7
c@335 8 #include <boost/test/unit_test.hpp>
c@335 9
c@355 10 #include <stdexcept>
c@355 11
c@335 12 BOOST_AUTO_TEST_SUITE(TestFFT)
c@335 13
c@335 14 #define COMPARE_CONST(a, n) \
c@335 15 for (int cmp_i = 0; cmp_i < (int)(sizeof(a)/sizeof(a[0])); ++cmp_i) { \
c@335 16 BOOST_CHECK_SMALL(a[cmp_i] - n, 1e-14); \
c@335 17 }
c@335 18
c@335 19 #define COMPARE_ARRAY(a, b) \
c@335 20 for (int cmp_i = 0; cmp_i < (int)(sizeof(a)/sizeof(a[0])); ++cmp_i) { \
c@335 21 BOOST_CHECK_SMALL(a[cmp_i] - b[cmp_i], 1e-14); \
c@335 22 }
c@335 23
c@339 24 //!!! need at least one test with complex time-domain signal
c@335 25
c@335 26 BOOST_AUTO_TEST_CASE(forwardArrayBounds)
c@335 27 {
c@335 28 // initialise bins to something recognisable, so we can tell
c@335 29 // if they haven't been written
c@335 30 double in[] = { 1, 1, -1, -1 };
c@335 31 double re[] = { 999, 999, 999, 999, 999, 999 };
c@335 32 double im[] = { 999, 999, 999, 999, 999, 999 };
c@335 33 FFT(4).process(false, in, 0, re+1, im+1);
c@335 34 // And check we haven't overrun the arrays
c@335 35 BOOST_CHECK_EQUAL(re[0], 999.0);
c@335 36 BOOST_CHECK_EQUAL(im[0], 999.0);
c@335 37 BOOST_CHECK_EQUAL(re[5], 999.0);
c@335 38 BOOST_CHECK_EQUAL(im[5], 999.0);
c@335 39 }
c@335 40
c@339 41 BOOST_AUTO_TEST_CASE(r_forwardArrayBounds)
c@339 42 {
c@339 43 // initialise bins to something recognisable, so we can tell
c@339 44 // if they haven't been written
c@339 45 double in[] = { 1, 1, -1, -1 };
c@339 46 double re[] = { 999, 999, 999, 999, 999, 999 };
c@339 47 double im[] = { 999, 999, 999, 999, 999, 999 };
c@339 48 FFTReal(4).forward(in, re+1, im+1);
c@339 49 // And check we haven't overrun the arrays
c@339 50 BOOST_CHECK_EQUAL(re[0], 999.0);
c@339 51 BOOST_CHECK_EQUAL(im[0], 999.0);
c@339 52 BOOST_CHECK_EQUAL(re[5], 999.0);
c@339 53 BOOST_CHECK_EQUAL(im[5], 999.0);
c@339 54 }
c@339 55
c@335 56 BOOST_AUTO_TEST_CASE(inverseArrayBounds)
c@335 57 {
c@335 58 // initialise bins to something recognisable, so we can tell
c@335 59 // if they haven't been written
c@339 60 double re[] = { 0, 1, 0, 1 };
c@339 61 double im[] = { 0, -2, 0, 2 };
c@335 62 double outre[] = { 999, 999, 999, 999, 999, 999 };
c@335 63 double outim[] = { 999, 999, 999, 999, 999, 999 };
c@339 64 FFT(4).process(true, re, im, outre+1, outim+1);
c@335 65 // And check we haven't overrun the arrays
c@335 66 BOOST_CHECK_EQUAL(outre[0], 999.0);
c@335 67 BOOST_CHECK_EQUAL(outim[0], 999.0);
c@335 68 BOOST_CHECK_EQUAL(outre[5], 999.0);
c@335 69 BOOST_CHECK_EQUAL(outim[5], 999.0);
c@335 70 }
c@335 71
c@339 72 BOOST_AUTO_TEST_CASE(r_inverseArrayBounds)
c@339 73 {
c@339 74 // initialise bins to something recognisable, so we can tell
c@339 75 // if they haven't been written
c@339 76 double re[] = { 0, 1, 0 };
c@339 77 double im[] = { 0, -2, 0 };
c@339 78 double outre[] = { 999, 999, 999, 999, 999, 999 };
c@339 79 FFTReal(4).inverse(re, im, outre+1);
c@339 80 // And check we haven't overrun the arrays
c@339 81 BOOST_CHECK_EQUAL(outre[0], 999.0);
c@339 82 BOOST_CHECK_EQUAL(outre[5], 999.0);
c@339 83 }
c@339 84
c@339 85 BOOST_AUTO_TEST_CASE(dc)
c@339 86 {
c@339 87 // DC-only signal. The DC bin is purely real
c@339 88 double in[] = { 1, 1, 1, 1 };
c@339 89 double re[] = { 999, 999, 999, 999 };
c@339 90 double im[] = { 999, 999, 999, 999 };
c@339 91 FFT(4).process(false, in, 0, re, im);
c@339 92 BOOST_CHECK_EQUAL(re[0], 4.0);
c@339 93 BOOST_CHECK_EQUAL(re[1], 0.0);
c@339 94 BOOST_CHECK_EQUAL(re[2], 0.0);
c@339 95 BOOST_CHECK_EQUAL(re[3], 0.0);
c@339 96 COMPARE_CONST(im, 0.0);
c@339 97 double back[4];
c@339 98 double backim[4];
c@339 99 FFT(4).process(true, re, im, back, backim);
c@339 100 COMPARE_ARRAY(back, in);
c@339 101 COMPARE_CONST(backim, 0.0);
c@339 102 }
c@339 103
c@339 104 BOOST_AUTO_TEST_CASE(r_dc)
c@339 105 {
c@339 106 // DC-only signal. The DC bin is purely real
c@339 107 double in[] = { 1, 1, 1, 1 };
c@339 108 double re[] = { 999, 999, 999, 999 };
c@339 109 double im[] = { 999, 999, 999, 999 };
c@339 110 FFTReal(4).forward(in, re, im);
c@339 111 BOOST_CHECK_EQUAL(re[0], 4.0);
c@339 112 BOOST_CHECK_EQUAL(re[1], 0.0);
c@339 113 BOOST_CHECK_EQUAL(re[2], 0.0);
c@339 114 BOOST_CHECK_EQUAL(re[3], 0.0);
c@339 115 COMPARE_CONST(im, 0.0);
c@339 116 double back[4];
c@339 117 // check conjugates are reconstructed
c@339 118 re[3] = 999;
c@339 119 im[3] = 999;
c@339 120 FFTReal(4).inverse(re, im, back);
c@339 121 COMPARE_ARRAY(back, in);
c@339 122 }
c@339 123
c@384 124 BOOST_AUTO_TEST_CASE(c_dc)
c@384 125 {
c@384 126 // DC-only signal. The DC bin is purely real
c@384 127 double rin[] = { 1, 1, 1, 1 };
c@384 128 double iin[] = { 1, 1, 1, 1 };
c@384 129 double re[] = { 999, 999, 999, 999 };
c@384 130 double im[] = { 999, 999, 999, 999 };
c@384 131 FFT(4).process(false, rin, iin, re, im);
c@384 132 BOOST_CHECK_EQUAL(re[0], 4.0);
c@384 133 BOOST_CHECK_EQUAL(re[1], 0.0);
c@384 134 BOOST_CHECK_EQUAL(re[2], 0.0);
c@384 135 BOOST_CHECK_EQUAL(re[3], 0.0);
c@384 136 BOOST_CHECK_EQUAL(im[0], 4.0);
c@384 137 BOOST_CHECK_EQUAL(im[1], 0.0);
c@384 138 BOOST_CHECK_EQUAL(im[2], 0.0);
c@384 139 BOOST_CHECK_EQUAL(im[3], 0.0);
c@384 140 double back[4];
c@384 141 double backim[4];
c@384 142 FFT(4).process(true, re, im, back, backim);
c@384 143 COMPARE_ARRAY(back, rin);
c@384 144 COMPARE_ARRAY(backim, iin);
c@384 145 }
c@384 146
c@339 147 BOOST_AUTO_TEST_CASE(sine)
c@339 148 {
c@339 149 // Sine. Output is purely imaginary
c@339 150 double in[] = { 0, 1, 0, -1 };
c@339 151 double re[] = { 999, 999, 999, 999 };
c@339 152 double im[] = { 999, 999, 999, 999 };
c@339 153 FFT(4).process(false, in, 0, re, im);
c@339 154 COMPARE_CONST(re, 0.0);
c@339 155 BOOST_CHECK_EQUAL(im[0], 0.0);
c@339 156 BOOST_CHECK_EQUAL(im[1], -2.0);
c@339 157 BOOST_CHECK_EQUAL(im[2], 0.0);
c@339 158 BOOST_CHECK_EQUAL(im[3], 2.0);
c@339 159 double back[4];
c@339 160 double backim[4];
c@339 161 FFT(4).process(true, re, im, back, backim);
c@339 162 COMPARE_ARRAY(back, in);
c@339 163 COMPARE_CONST(backim, 0.0);
c@339 164 }
c@339 165
c@339 166 BOOST_AUTO_TEST_CASE(r_sine)
c@339 167 {
c@339 168 // Sine. Output is purely imaginary
c@339 169 double in[] = { 0, 1, 0, -1 };
c@339 170 double re[] = { 999, 999, 999, 999 };
c@339 171 double im[] = { 999, 999, 999, 999 };
c@339 172 FFTReal(4).forward(in, re, im);
c@339 173 COMPARE_CONST(re, 0.0);
c@339 174 BOOST_CHECK_EQUAL(im[0], 0.0);
c@339 175 BOOST_CHECK_EQUAL(im[1], -2.0);
c@339 176 BOOST_CHECK_EQUAL(im[2], 0.0);
c@339 177 BOOST_CHECK_EQUAL(im[3], 2.0);
c@339 178 double back[4];
c@339 179 // check conjugates are reconstructed
c@339 180 re[3] = 999;
c@339 181 im[3] = 999;
c@339 182 FFTReal(4).inverse(re, im, back);
c@339 183 COMPARE_ARRAY(back, in);
c@339 184 }
c@339 185
c@339 186 BOOST_AUTO_TEST_CASE(cosine)
c@339 187 {
c@339 188 // Cosine. Output is purely real
c@339 189 double in[] = { 1, 0, -1, 0 };
c@339 190 double re[] = { 999, 999, 999, 999 };
c@339 191 double im[] = { 999, 999, 999, 999 };
c@339 192 FFT(4).process(false, in, 0, re, im);
c@339 193 BOOST_CHECK_EQUAL(re[0], 0.0);
c@339 194 BOOST_CHECK_EQUAL(re[1], 2.0);
c@339 195 BOOST_CHECK_EQUAL(re[2], 0.0);
c@339 196 BOOST_CHECK_EQUAL(re[3], 2.0);
c@339 197 COMPARE_CONST(im, 0.0);
c@339 198 double back[4];
c@339 199 double backim[4];
c@339 200 FFT(4).process(true, re, im, back, backim);
c@339 201 COMPARE_ARRAY(back, in);
c@339 202 COMPARE_CONST(backim, 0.0);
c@339 203 }
c@339 204
c@339 205 BOOST_AUTO_TEST_CASE(r_cosine)
c@339 206 {
c@339 207 // Cosine. Output is purely real
c@339 208 double in[] = { 1, 0, -1, 0 };
c@339 209 double re[] = { 999, 999, 999, 999 };
c@339 210 double im[] = { 999, 999, 999, 999 };
c@339 211 FFTReal(4).forward(in, re, im);
c@339 212 BOOST_CHECK_EQUAL(re[0], 0.0);
c@339 213 BOOST_CHECK_EQUAL(re[1], 2.0);
c@339 214 BOOST_CHECK_EQUAL(re[2], 0.0);
c@339 215 BOOST_CHECK_EQUAL(re[3], 2.0);
c@339 216 COMPARE_CONST(im, 0.0);
c@339 217 double back[4];
c@339 218 // check conjugates are reconstructed
c@339 219 re[3] = 999;
c@339 220 im[3] = 999;
c@339 221 FFTReal(4).inverse(re, im, back);
c@339 222 COMPARE_ARRAY(back, in);
c@339 223 }
c@384 224
c@384 225 BOOST_AUTO_TEST_CASE(c_cosine)
c@384 226 {
c@384 227 // Cosine. Output is purely real
c@384 228 double rin[] = { 1, 0, -1, 0 };
c@384 229 double iin[] = { 1, 0, -1, 0 };
c@384 230 double re[] = { 999, 999, 999, 999 };
c@384 231 double im[] = { 999, 999, 999, 999 };
c@384 232 FFT(4).process(false, rin, iin, re, im);
c@384 233 BOOST_CHECK_EQUAL(re[0], 0.0);
c@384 234 BOOST_CHECK_EQUAL(re[1], 2.0);
c@384 235 BOOST_CHECK_EQUAL(re[2], 0.0);
c@384 236 BOOST_CHECK_EQUAL(re[3], 2.0);
c@384 237 BOOST_CHECK_EQUAL(im[0], 0.0);
c@384 238 BOOST_CHECK_EQUAL(im[1], 2.0);
c@384 239 BOOST_CHECK_EQUAL(im[2], 0.0);
c@384 240 BOOST_CHECK_EQUAL(im[3], 2.0);
c@384 241 double back[4];
c@384 242 double backim[4];
c@384 243 FFT(4).process(true, re, im, back, backim);
c@384 244 COMPARE_ARRAY(back, rin);
c@384 245 COMPARE_ARRAY(backim, iin);
c@384 246 }
c@339 247
c@339 248 BOOST_AUTO_TEST_CASE(sineCosine)
c@339 249 {
c@339 250 // Sine and cosine mixed
c@339 251 double in[] = { 0.5, 1, -0.5, -1 };
c@339 252 double re[] = { 999, 999, 999, 999 };
c@339 253 double im[] = { 999, 999, 999, 999 };
c@339 254 FFT(4).process(false, in, 0, re, im);
c@339 255 BOOST_CHECK_EQUAL(re[0], 0.0);
c@339 256 BOOST_CHECK_CLOSE(re[1], 1.0, 1e-12);
c@339 257 BOOST_CHECK_EQUAL(re[2], 0.0);
c@339 258 BOOST_CHECK_CLOSE(re[3], 1.0, 1e-12);
c@339 259 BOOST_CHECK_EQUAL(im[0], 0.0);
c@339 260 BOOST_CHECK_CLOSE(im[1], -2.0, 1e-12);
c@339 261 BOOST_CHECK_EQUAL(im[2], 0.0);
c@339 262 BOOST_CHECK_CLOSE(im[3], 2.0, 1e-12);
c@339 263 double back[4];
c@339 264 double backim[4];
c@339 265 FFT(4).process(true, re, im, back, backim);
c@339 266 COMPARE_ARRAY(back, in);
c@339 267 COMPARE_CONST(backim, 0.0);
c@339 268 }
c@339 269
c@339 270 BOOST_AUTO_TEST_CASE(r_sineCosine)
c@339 271 {
c@339 272 // Sine and cosine mixed
c@339 273 double in[] = { 0.5, 1, -0.5, -1 };
c@339 274 double re[] = { 999, 999, 999, 999 };
c@339 275 double im[] = { 999, 999, 999, 999 };
c@339 276 FFTReal(4).forward(in, re, im);
c@339 277 BOOST_CHECK_EQUAL(re[0], 0.0);
c@339 278 BOOST_CHECK_CLOSE(re[1], 1.0, 1e-12);
c@339 279 BOOST_CHECK_EQUAL(re[2], 0.0);
c@339 280 BOOST_CHECK_CLOSE(re[3], 1.0, 1e-12);
c@339 281 BOOST_CHECK_EQUAL(im[0], 0.0);
c@339 282 BOOST_CHECK_CLOSE(im[1], -2.0, 1e-12);
c@339 283 BOOST_CHECK_EQUAL(im[2], 0.0);
c@339 284 BOOST_CHECK_CLOSE(im[3], 2.0, 1e-12);
c@339 285 double back[4];
c@339 286 // check conjugates are reconstructed
c@339 287 re[3] = 999;
c@339 288 im[3] = 999;
c@339 289 FFTReal(4).inverse(re, im, back);
c@339 290 COMPARE_ARRAY(back, in);
c@339 291 }
c@384 292
c@384 293 BOOST_AUTO_TEST_CASE(c_sineCosine)
c@384 294 {
c@384 295 double rin[] = { 1, 0, -1, 0 };
c@384 296 double iin[] = { 0, 1, 0, -1 };
c@384 297 double re[] = { 999, 999, 999, 999 };
c@384 298 double im[] = { 999, 999, 999, 999 };
c@384 299 FFT(4).process(false, rin, iin, re, im);
c@384 300 BOOST_CHECK_EQUAL(re[0], 0.0);
c@384 301 BOOST_CHECK_EQUAL(re[1], 4.0);
c@384 302 BOOST_CHECK_EQUAL(re[2], 0.0);
c@384 303 BOOST_CHECK_EQUAL(re[3], 0.0);
c@384 304 COMPARE_CONST(im, 0.0);
c@384 305 double back[4];
c@384 306 double backim[4];
c@384 307 FFT(4).process(true, re, im, back, backim);
c@384 308 COMPARE_ARRAY(back, rin);
c@384 309 COMPARE_ARRAY(backim, iin);
c@384 310 }
c@339 311
c@339 312 BOOST_AUTO_TEST_CASE(nyquist)
c@339 313 {
c@339 314 double in[] = { 1, -1, 1, -1 };
c@339 315 double re[] = { 999, 999, 999, 999 };
c@339 316 double im[] = { 999, 999, 999, 999 };
c@339 317 FFT(4).process(false, in, 0, re, im);
c@339 318 BOOST_CHECK_EQUAL(re[0], 0.0);
c@339 319 BOOST_CHECK_EQUAL(re[1], 0.0);
c@339 320 BOOST_CHECK_EQUAL(re[2], 4.0);
c@339 321 BOOST_CHECK_EQUAL(re[3], 0.0);
c@339 322 COMPARE_CONST(im, 0.0);
c@339 323 double back[4];
c@339 324 double backim[4];
c@339 325 FFT(4).process(true, re, im, back, backim);
c@339 326 COMPARE_ARRAY(back, in);
c@339 327 COMPARE_CONST(backim, 0.0);
c@339 328 }
c@339 329
c@339 330 BOOST_AUTO_TEST_CASE(r_nyquist)
c@339 331 {
c@339 332 double in[] = { 1, -1, 1, -1 };
c@339 333 double re[] = { 999, 999, 999, 999 };
c@339 334 double im[] = { 999, 999, 999, 999 };
c@339 335 FFTReal(4).forward(in, re, im);
c@339 336 BOOST_CHECK_EQUAL(re[0], 0.0);
c@339 337 BOOST_CHECK_EQUAL(re[1], 0.0);
c@339 338 BOOST_CHECK_EQUAL(re[2], 4.0);
c@339 339 BOOST_CHECK_EQUAL(re[3], 0.0);
c@339 340 COMPARE_CONST(im, 0.0);
c@339 341 double back[4];
c@339 342 // check conjugates are reconstructed
c@339 343 re[3] = 999;
c@339 344 im[3] = 999;
c@339 345 FFTReal(4).inverse(re, im, back);
c@339 346 COMPARE_ARRAY(back, in);
c@339 347 }
c@339 348
c@339 349 BOOST_AUTO_TEST_CASE(dirac)
c@339 350 {
c@339 351 double in[] = { 1, 0, 0, 0 };
c@339 352 double re[] = { 999, 999, 999, 999 };
c@339 353 double im[] = { 999, 999, 999, 999 };
c@339 354 FFT(4).process(false, in, 0, re, im);
c@339 355 BOOST_CHECK_EQUAL(re[0], 1.0);
c@339 356 BOOST_CHECK_EQUAL(re[1], 1.0);
c@339 357 BOOST_CHECK_EQUAL(re[2], 1.0);
c@339 358 BOOST_CHECK_EQUAL(re[3], 1.0);
c@339 359 COMPARE_CONST(im, 0.0);
c@339 360 double back[4];
c@339 361 double backim[4];
c@339 362 FFT(4).process(true, re, im, back, backim);
c@339 363 COMPARE_ARRAY(back, in);
c@339 364 COMPARE_CONST(backim, 0.0);
c@339 365 }
c@339 366
c@339 367 BOOST_AUTO_TEST_CASE(r_dirac)
c@339 368 {
c@339 369 double in[] = { 1, 0, 0, 0 };
c@339 370 double re[] = { 999, 999, 999, 999 };
c@339 371 double im[] = { 999, 999, 999, 999 };
c@339 372 FFTReal(4).forward(in, re, im);
c@339 373 BOOST_CHECK_EQUAL(re[0], 1.0);
c@339 374 BOOST_CHECK_EQUAL(re[1], 1.0);
c@339 375 BOOST_CHECK_EQUAL(re[2], 1.0);
c@339 376 BOOST_CHECK_EQUAL(re[3], 1.0);
c@339 377 COMPARE_CONST(im, 0.0);
c@339 378 double back[4];
c@339 379 // check conjugates are reconstructed
c@339 380 re[3] = 999;
c@339 381 im[3] = 999;
c@339 382 FFTReal(4).inverse(re, im, back);
c@339 383 COMPARE_ARRAY(back, in);
c@339 384 }
c@339 385
c@355 386 BOOST_AUTO_TEST_CASE(sizes)
c@355 387 {
c@355 388 // Complex supports any size. A single test with an odd size
c@355 389 // will do here, without getting too much into our expectations
c@355 390 // about supported butterflies etc
c@355 391
c@355 392 double in[] = { 1, 1, 1 };
c@355 393 double re[] = { 999, 999, 999 };
c@355 394 double im[] = { 999, 999, 999 };
c@355 395 FFT(3).process(false, in, 0, re, im);
c@355 396 BOOST_CHECK_EQUAL(re[0], 3.0);
c@355 397 BOOST_CHECK_EQUAL(re[1], 0.0);
c@355 398 BOOST_CHECK_EQUAL(re[2], 0.0);
c@355 399 COMPARE_CONST(im, 0.0);
c@355 400 double back[3];
c@355 401 double backim[3];
c@355 402 FFT(3).process(true, re, im, back, backim);
c@355 403 COMPARE_ARRAY(back, in);
c@355 404 COMPARE_CONST(backim, 0.0);
c@355 405 }
c@355 406
c@355 407 BOOST_AUTO_TEST_CASE(r_sizes)
c@355 408 {
c@355 409 // Real supports any even size, but not odd ones
c@355 410
c@355 411 BOOST_CHECK_THROW(FFTReal r(3), std::invalid_argument);
c@355 412
c@355 413 double in[] = { 1, 1, 1, 1, 1, 1 };
c@355 414 double re[] = { 999, 999, 999, 999, 999, 999 };
c@355 415 double im[] = { 999, 999, 999, 999, 999, 999 };
c@355 416 FFTReal(6).forward(in, re, im);
c@355 417 BOOST_CHECK_EQUAL(re[0], 6.0);
c@355 418 BOOST_CHECK_EQUAL(re[1], 0.0);
c@355 419 BOOST_CHECK_EQUAL(re[2], 0.0);
c@355 420 BOOST_CHECK_EQUAL(re[3], 0.0);
c@355 421 BOOST_CHECK_EQUAL(re[4], 0.0);
c@355 422 BOOST_CHECK_EQUAL(re[5], 0.0);
c@355 423 COMPARE_CONST(im, 0.0);
c@355 424 double back[6];
c@355 425 FFTReal(6).inverse(re, im, back);
c@355 426 COMPARE_ARRAY(back, in);
c@355 427 }
c@355 428
c@335 429 BOOST_AUTO_TEST_SUITE_END()
c@335 430