Chris@366: /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */ Chris@366: Chris@366: #include "dsp/FFT.h" Chris@366: Chris@366: #define BOOST_TEST_DYN_LINK Chris@366: #define BOOST_TEST_MAIN Chris@366: Chris@366: #include Chris@366: Chris@366: #include Chris@366: Chris@366: BOOST_AUTO_TEST_SUITE(TestFFT) Chris@366: Chris@366: #define COMPARE_CONST(a, n) \ Chris@366: for (int cmp_i = 0; cmp_i < (int)(sizeof(a)/sizeof(a[0])); ++cmp_i) { \ Chris@366: BOOST_CHECK_SMALL(a[cmp_i] - n, 1e-14); \ Chris@366: } Chris@366: Chris@366: #define COMPARE_ARRAY(a, b) \ Chris@366: for (int cmp_i = 0; cmp_i < (int)(sizeof(a)/sizeof(a[0])); ++cmp_i) { \ Chris@366: BOOST_CHECK_SMALL(a[cmp_i] - b[cmp_i], 1e-14); \ Chris@366: } Chris@366: Chris@366: //!!! need at least one test with complex time-domain signal Chris@366: Chris@366: BOOST_AUTO_TEST_CASE(forwardArrayBounds) Chris@366: { Chris@366: // initialise bins to something recognisable, so we can tell if Chris@366: // they haven't been written; and allocate the inputs on the heap Chris@366: // so that, if running under valgrind, we get warnings about Chris@366: // overruns Chris@366: double *in = new double[4]; Chris@366: in[0] = 1; Chris@366: in[1] = 1; Chris@366: in[2] = -1; Chris@366: in[3] = -1; Chris@366: double re[] = { 999, 999, 999, 999, 999, 999 }; Chris@366: double im[] = { 999, 999, 999, 999, 999, 999 }; Chris@366: FFT(4).process(false, in, 0, re+1, im+1); Chris@366: // And check we haven't overrun the arrays Chris@366: BOOST_CHECK_EQUAL(re[0], 999.0); Chris@366: BOOST_CHECK_EQUAL(im[0], 999.0); Chris@366: BOOST_CHECK_EQUAL(re[5], 999.0); Chris@366: BOOST_CHECK_EQUAL(im[5], 999.0); Chris@366: delete[] in; Chris@366: } Chris@366: Chris@366: BOOST_AUTO_TEST_CASE(r_forwardArrayBounds) Chris@366: { Chris@366: // initialise bins to something recognisable, so we can tell if Chris@366: // they haven't been written; and allocate the inputs on the heap Chris@366: // so that, if running under valgrind, we get warnings about Chris@366: // overruns Chris@366: double *in = new double[4]; Chris@366: in[0] = 1; Chris@366: in[1] = 1; Chris@366: in[2] = -1; Chris@366: in[3] = -1; Chris@366: double re[] = { 999, 999, 999, 999, 999, 999 }; Chris@366: double im[] = { 999, 999, 999, 999, 999, 999 }; Chris@366: FFTReal(4).forward(in, re+1, im+1); Chris@366: // And check we haven't overrun the arrays Chris@366: BOOST_CHECK_EQUAL(re[0], 999.0); Chris@366: BOOST_CHECK_EQUAL(im[0], 999.0); Chris@366: BOOST_CHECK_EQUAL(re[5], 999.0); Chris@366: BOOST_CHECK_EQUAL(im[5], 999.0); Chris@366: delete[] in; Chris@366: } Chris@366: Chris@366: BOOST_AUTO_TEST_CASE(inverseArrayBounds) Chris@366: { Chris@366: // initialise bins to something recognisable, so we can tell if Chris@366: // they haven't been written; and allocate the inputs on the heap Chris@366: // so that, if running under valgrind, we get warnings about Chris@366: // overruns Chris@366: double *re = new double[4]; Chris@366: double *im = new double[4]; Chris@366: re[0] = 0; Chris@366: re[1] = 1; Chris@366: re[2] = 0; Chris@366: re[3] = 1; Chris@366: im[0] = 0; Chris@366: im[1] = -2; Chris@366: im[2] = 0; Chris@366: im[3] = 2; Chris@366: double outre[] = { 999, 999, 999, 999, 999, 999 }; Chris@366: double outim[] = { 999, 999, 999, 999, 999, 999 }; Chris@366: FFT(4).process(true, re, im, outre+1, outim+1); Chris@366: // And check we haven't overrun the arrays Chris@366: BOOST_CHECK_EQUAL(outre[0], 999.0); Chris@366: BOOST_CHECK_EQUAL(outim[0], 999.0); Chris@366: BOOST_CHECK_EQUAL(outre[5], 999.0); Chris@366: BOOST_CHECK_EQUAL(outim[5], 999.0); Chris@366: delete[] re; Chris@366: delete[] im; Chris@366: } Chris@366: Chris@366: BOOST_AUTO_TEST_CASE(r_inverseArrayBounds) Chris@366: { Chris@366: // initialise bins to something recognisable, so we can tell if Chris@366: // they haven't been written; and allocate the inputs on the heap Chris@366: // so that, if running under valgrind, we get warnings about Chris@366: // overruns Chris@366: double *re = new double[3]; Chris@366: double *im = new double[3]; Chris@366: re[0] = 0; Chris@366: re[1] = 1; Chris@366: re[2] = 0; Chris@366: im[0] = 0; Chris@366: im[1] = -2; Chris@366: im[2] = 0; Chris@366: double outre[] = { 999, 999, 999, 999, 999, 999 }; Chris@366: FFTReal(4).inverse(re, im, outre+1); Chris@366: // And check we haven't overrun the arrays Chris@366: BOOST_CHECK_EQUAL(outre[0], 999.0); Chris@366: BOOST_CHECK_EQUAL(outre[5], 999.0); Chris@366: delete[] re; Chris@366: delete[] im; Chris@366: } Chris@366: Chris@366: BOOST_AUTO_TEST_CASE(dc) Chris@366: { Chris@366: // DC-only signal. The DC bin is purely real Chris@366: double in[] = { 1, 1, 1, 1 }; Chris@366: double re[] = { 999, 999, 999, 999 }; Chris@366: double im[] = { 999, 999, 999, 999 }; Chris@366: FFT(4).process(false, in, 0, re, im); Chris@366: BOOST_CHECK_EQUAL(re[0], 4.0); Chris@366: BOOST_CHECK_EQUAL(re[1], 0.0); Chris@366: BOOST_CHECK_EQUAL(re[2], 0.0); Chris@366: BOOST_CHECK_EQUAL(re[3], 0.0); Chris@366: COMPARE_CONST(im, 0.0); Chris@366: double back[4]; Chris@366: double backim[4]; Chris@366: FFT(4).process(true, re, im, back, backim); Chris@366: COMPARE_ARRAY(back, in); Chris@366: COMPARE_CONST(backim, 0.0); Chris@366: } Chris@366: Chris@366: BOOST_AUTO_TEST_CASE(r_dc) Chris@366: { Chris@366: // DC-only signal. The DC bin is purely real Chris@366: double in[] = { 1, 1, 1, 1 }; Chris@366: double re[] = { 999, 999, 999, 999 }; Chris@366: double im[] = { 999, 999, 999, 999 }; Chris@366: FFTReal(4).forward(in, re, im); Chris@366: BOOST_CHECK_EQUAL(re[0], 4.0); Chris@366: BOOST_CHECK_EQUAL(re[1], 0.0); Chris@366: BOOST_CHECK_EQUAL(re[2], 0.0); Chris@366: BOOST_CHECK_EQUAL(re[3], 0.0); Chris@366: COMPARE_CONST(im, 0.0); Chris@366: double back[4]; Chris@366: // check conjugates are reconstructed Chris@366: re[3] = 999; Chris@366: im[3] = 999; Chris@366: FFTReal(4).inverse(re, im, back); Chris@366: COMPARE_ARRAY(back, in); Chris@366: } Chris@366: Chris@366: BOOST_AUTO_TEST_CASE(c_dc) Chris@366: { Chris@366: // DC-only signal. The DC bin is purely real Chris@366: double rin[] = { 1, 1, 1, 1 }; Chris@366: double iin[] = { 1, 1, 1, 1 }; Chris@366: double re[] = { 999, 999, 999, 999 }; Chris@366: double im[] = { 999, 999, 999, 999 }; Chris@366: FFT(4).process(false, rin, iin, re, im); Chris@366: BOOST_CHECK_EQUAL(re[0], 4.0); Chris@366: BOOST_CHECK_EQUAL(re[1], 0.0); Chris@366: BOOST_CHECK_EQUAL(re[2], 0.0); Chris@366: BOOST_CHECK_EQUAL(re[3], 0.0); Chris@366: BOOST_CHECK_EQUAL(im[0], 4.0); Chris@366: BOOST_CHECK_EQUAL(im[1], 0.0); Chris@366: BOOST_CHECK_EQUAL(im[2], 0.0); Chris@366: BOOST_CHECK_EQUAL(im[3], 0.0); Chris@366: double back[4]; Chris@366: double backim[4]; Chris@366: FFT(4).process(true, re, im, back, backim); Chris@366: COMPARE_ARRAY(back, rin); Chris@366: COMPARE_ARRAY(backim, iin); Chris@366: } Chris@366: Chris@366: BOOST_AUTO_TEST_CASE(sine) Chris@366: { Chris@366: // Sine. Output is purely imaginary Chris@366: double in[] = { 0, 1, 0, -1 }; Chris@366: double re[] = { 999, 999, 999, 999 }; Chris@366: double im[] = { 999, 999, 999, 999 }; Chris@366: FFT(4).process(false, in, 0, re, im); Chris@366: COMPARE_CONST(re, 0.0); Chris@366: BOOST_CHECK_EQUAL(im[0], 0.0); Chris@366: BOOST_CHECK_EQUAL(im[1], -2.0); Chris@366: BOOST_CHECK_EQUAL(im[2], 0.0); Chris@366: BOOST_CHECK_EQUAL(im[3], 2.0); Chris@366: double back[4]; Chris@366: double backim[4]; Chris@366: FFT(4).process(true, re, im, back, backim); Chris@366: COMPARE_ARRAY(back, in); Chris@366: COMPARE_CONST(backim, 0.0); Chris@366: } Chris@366: Chris@366: BOOST_AUTO_TEST_CASE(r_sine) Chris@366: { Chris@366: // Sine. Output is purely imaginary Chris@366: double in[] = { 0, 1, 0, -1 }; Chris@366: double re[] = { 999, 999, 999, 999 }; Chris@366: double im[] = { 999, 999, 999, 999 }; Chris@366: FFTReal(4).forward(in, re, im); Chris@366: COMPARE_CONST(re, 0.0); Chris@366: BOOST_CHECK_EQUAL(im[0], 0.0); Chris@366: BOOST_CHECK_EQUAL(im[1], -2.0); Chris@366: BOOST_CHECK_EQUAL(im[2], 0.0); Chris@366: BOOST_CHECK_EQUAL(im[3], 2.0); Chris@366: double back[4]; Chris@366: // check conjugates are reconstructed Chris@366: re[3] = 999; Chris@366: im[3] = 999; Chris@366: FFTReal(4).inverse(re, im, back); Chris@366: COMPARE_ARRAY(back, in); Chris@366: } Chris@366: Chris@366: BOOST_AUTO_TEST_CASE(cosine) Chris@366: { Chris@366: // Cosine. Output is purely real Chris@366: double in[] = { 1, 0, -1, 0 }; Chris@366: double re[] = { 999, 999, 999, 999 }; Chris@366: double im[] = { 999, 999, 999, 999 }; Chris@366: FFT(4).process(false, in, 0, re, im); Chris@366: BOOST_CHECK_EQUAL(re[0], 0.0); Chris@366: BOOST_CHECK_EQUAL(re[1], 2.0); Chris@366: BOOST_CHECK_EQUAL(re[2], 0.0); Chris@366: BOOST_CHECK_EQUAL(re[3], 2.0); Chris@366: COMPARE_CONST(im, 0.0); Chris@366: double back[4]; Chris@366: double backim[4]; Chris@366: FFT(4).process(true, re, im, back, backim); Chris@366: COMPARE_ARRAY(back, in); Chris@366: COMPARE_CONST(backim, 0.0); Chris@366: } Chris@366: Chris@366: BOOST_AUTO_TEST_CASE(r_cosine) Chris@366: { Chris@366: // Cosine. Output is purely real Chris@366: double in[] = { 1, 0, -1, 0 }; Chris@366: double re[] = { 999, 999, 999, 999 }; Chris@366: double im[] = { 999, 999, 999, 999 }; Chris@366: FFTReal(4).forward(in, re, im); Chris@366: BOOST_CHECK_EQUAL(re[0], 0.0); Chris@366: BOOST_CHECK_EQUAL(re[1], 2.0); Chris@366: BOOST_CHECK_EQUAL(re[2], 0.0); Chris@366: BOOST_CHECK_EQUAL(re[3], 2.0); Chris@366: COMPARE_CONST(im, 0.0); Chris@366: double back[4]; Chris@366: // check conjugates are reconstructed Chris@366: re[3] = 999; Chris@366: im[3] = 999; Chris@366: FFTReal(4).inverse(re, im, back); Chris@366: COMPARE_ARRAY(back, in); Chris@366: } Chris@366: Chris@366: BOOST_AUTO_TEST_CASE(c_cosine) Chris@366: { Chris@366: // Cosine. Output is purely real Chris@366: double rin[] = { 1, 0, -1, 0 }; Chris@366: double iin[] = { 1, 0, -1, 0 }; Chris@366: double re[] = { 999, 999, 999, 999 }; Chris@366: double im[] = { 999, 999, 999, 999 }; Chris@366: FFT(4).process(false, rin, iin, re, im); Chris@366: BOOST_CHECK_EQUAL(re[0], 0.0); Chris@366: BOOST_CHECK_EQUAL(re[1], 2.0); Chris@366: BOOST_CHECK_EQUAL(re[2], 0.0); Chris@366: BOOST_CHECK_EQUAL(re[3], 2.0); Chris@366: BOOST_CHECK_EQUAL(im[0], 0.0); Chris@366: BOOST_CHECK_EQUAL(im[1], 2.0); Chris@366: BOOST_CHECK_EQUAL(im[2], 0.0); Chris@366: BOOST_CHECK_EQUAL(im[3], 2.0); Chris@366: double back[4]; Chris@366: double backim[4]; Chris@366: FFT(4).process(true, re, im, back, backim); Chris@366: COMPARE_ARRAY(back, rin); Chris@366: COMPARE_ARRAY(backim, iin); Chris@366: } Chris@366: Chris@366: BOOST_AUTO_TEST_CASE(sineCosine) Chris@366: { Chris@366: // Sine and cosine mixed Chris@366: double in[] = { 0.5, 1, -0.5, -1 }; Chris@366: double re[] = { 999, 999, 999, 999 }; Chris@366: double im[] = { 999, 999, 999, 999 }; Chris@366: FFT(4).process(false, in, 0, re, im); Chris@366: BOOST_CHECK_EQUAL(re[0], 0.0); Chris@366: BOOST_CHECK_CLOSE(re[1], 1.0, 1e-12); Chris@366: BOOST_CHECK_EQUAL(re[2], 0.0); Chris@366: BOOST_CHECK_CLOSE(re[3], 1.0, 1e-12); Chris@366: BOOST_CHECK_EQUAL(im[0], 0.0); Chris@366: BOOST_CHECK_CLOSE(im[1], -2.0, 1e-12); Chris@366: BOOST_CHECK_EQUAL(im[2], 0.0); Chris@366: BOOST_CHECK_CLOSE(im[3], 2.0, 1e-12); Chris@366: double back[4]; Chris@366: double backim[4]; Chris@366: FFT(4).process(true, re, im, back, backim); Chris@366: COMPARE_ARRAY(back, in); Chris@366: COMPARE_CONST(backim, 0.0); Chris@366: } Chris@366: Chris@366: BOOST_AUTO_TEST_CASE(r_sineCosine) Chris@366: { Chris@366: // Sine and cosine mixed Chris@366: double in[] = { 0.5, 1, -0.5, -1 }; Chris@366: double re[] = { 999, 999, 999, 999 }; Chris@366: double im[] = { 999, 999, 999, 999 }; Chris@366: FFTReal(4).forward(in, re, im); Chris@366: BOOST_CHECK_EQUAL(re[0], 0.0); Chris@366: BOOST_CHECK_CLOSE(re[1], 1.0, 1e-12); Chris@366: BOOST_CHECK_EQUAL(re[2], 0.0); Chris@366: BOOST_CHECK_CLOSE(re[3], 1.0, 1e-12); Chris@366: BOOST_CHECK_EQUAL(im[0], 0.0); Chris@366: BOOST_CHECK_CLOSE(im[1], -2.0, 1e-12); Chris@366: BOOST_CHECK_EQUAL(im[2], 0.0); Chris@366: BOOST_CHECK_CLOSE(im[3], 2.0, 1e-12); Chris@366: double back[4]; Chris@366: // check conjugates are reconstructed Chris@366: re[3] = 999; Chris@366: im[3] = 999; Chris@366: FFTReal(4).inverse(re, im, back); Chris@366: COMPARE_ARRAY(back, in); Chris@366: } Chris@366: Chris@366: BOOST_AUTO_TEST_CASE(c_sineCosine) Chris@366: { Chris@366: double rin[] = { 1, 0, -1, 0 }; Chris@366: double iin[] = { 0, 1, 0, -1 }; Chris@366: double re[] = { 999, 999, 999, 999 }; Chris@366: double im[] = { 999, 999, 999, 999 }; Chris@366: FFT(4).process(false, rin, iin, re, im); Chris@366: BOOST_CHECK_EQUAL(re[0], 0.0); Chris@366: BOOST_CHECK_EQUAL(re[1], 4.0); Chris@366: BOOST_CHECK_EQUAL(re[2], 0.0); Chris@366: BOOST_CHECK_EQUAL(re[3], 0.0); Chris@366: COMPARE_CONST(im, 0.0); Chris@366: double back[4]; Chris@366: double backim[4]; Chris@366: FFT(4).process(true, re, im, back, backim); Chris@366: COMPARE_ARRAY(back, rin); Chris@366: COMPARE_ARRAY(backim, iin); Chris@366: } Chris@366: Chris@366: BOOST_AUTO_TEST_CASE(nyquist) Chris@366: { Chris@366: double in[] = { 1, -1, 1, -1 }; Chris@366: double re[] = { 999, 999, 999, 999 }; Chris@366: double im[] = { 999, 999, 999, 999 }; Chris@366: FFT(4).process(false, in, 0, re, im); Chris@366: BOOST_CHECK_EQUAL(re[0], 0.0); Chris@366: BOOST_CHECK_EQUAL(re[1], 0.0); Chris@366: BOOST_CHECK_EQUAL(re[2], 4.0); Chris@366: BOOST_CHECK_EQUAL(re[3], 0.0); Chris@366: COMPARE_CONST(im, 0.0); Chris@366: double back[4]; Chris@366: double backim[4]; Chris@366: FFT(4).process(true, re, im, back, backim); Chris@366: COMPARE_ARRAY(back, in); Chris@366: COMPARE_CONST(backim, 0.0); Chris@366: } Chris@366: Chris@366: BOOST_AUTO_TEST_CASE(r_nyquist) Chris@366: { Chris@366: double in[] = { 1, -1, 1, -1 }; Chris@366: double re[] = { 999, 999, 999, 999 }; Chris@366: double im[] = { 999, 999, 999, 999 }; Chris@366: FFTReal(4).forward(in, re, im); Chris@366: BOOST_CHECK_EQUAL(re[0], 0.0); Chris@366: BOOST_CHECK_EQUAL(re[1], 0.0); Chris@366: BOOST_CHECK_EQUAL(re[2], 4.0); Chris@366: BOOST_CHECK_EQUAL(re[3], 0.0); Chris@366: COMPARE_CONST(im, 0.0); Chris@366: double back[4]; Chris@366: // check conjugates are reconstructed Chris@366: re[3] = 999; Chris@366: im[3] = 999; Chris@366: FFTReal(4).inverse(re, im, back); Chris@366: COMPARE_ARRAY(back, in); Chris@366: } Chris@366: Chris@366: BOOST_AUTO_TEST_CASE(dirac) Chris@366: { Chris@366: double in[] = { 1, 0, 0, 0 }; Chris@366: double re[] = { 999, 999, 999, 999 }; Chris@366: double im[] = { 999, 999, 999, 999 }; Chris@366: FFT(4).process(false, in, 0, re, im); Chris@366: BOOST_CHECK_EQUAL(re[0], 1.0); Chris@366: BOOST_CHECK_EQUAL(re[1], 1.0); Chris@366: BOOST_CHECK_EQUAL(re[2], 1.0); Chris@366: BOOST_CHECK_EQUAL(re[3], 1.0); Chris@366: COMPARE_CONST(im, 0.0); Chris@366: double back[4]; Chris@366: double backim[4]; Chris@366: FFT(4).process(true, re, im, back, backim); Chris@366: COMPARE_ARRAY(back, in); Chris@366: COMPARE_CONST(backim, 0.0); Chris@366: } Chris@366: Chris@366: BOOST_AUTO_TEST_CASE(r_dirac) Chris@366: { Chris@366: double in[] = { 1, 0, 0, 0 }; Chris@366: double re[] = { 999, 999, 999, 999 }; Chris@366: double im[] = { 999, 999, 999, 999 }; Chris@366: FFTReal(4).forward(in, re, im); Chris@366: BOOST_CHECK_EQUAL(re[0], 1.0); Chris@366: BOOST_CHECK_EQUAL(re[1], 1.0); Chris@366: BOOST_CHECK_EQUAL(re[2], 1.0); Chris@366: BOOST_CHECK_EQUAL(re[3], 1.0); Chris@366: COMPARE_CONST(im, 0.0); Chris@366: double back[4]; Chris@366: // check conjugates are reconstructed Chris@366: re[3] = 999; Chris@366: im[3] = 999; Chris@366: FFTReal(4).inverse(re, im, back); Chris@366: COMPARE_ARRAY(back, in); Chris@366: } Chris@366: Chris@366: BOOST_AUTO_TEST_CASE(sizes) Chris@366: { Chris@366: // Complex supports any size. A single test with an odd size Chris@366: // will do here, without getting too much into our expectations Chris@366: // about supported butterflies etc Chris@366: Chris@366: double in[] = { 1, 1, 1 }; Chris@366: double re[] = { 999, 999, 999 }; Chris@366: double im[] = { 999, 999, 999 }; Chris@366: FFT(3).process(false, in, 0, re, im); Chris@366: BOOST_CHECK_EQUAL(re[0], 3.0); Chris@366: BOOST_CHECK_EQUAL(re[1], 0.0); Chris@366: BOOST_CHECK_EQUAL(re[2], 0.0); Chris@366: COMPARE_CONST(im, 0.0); Chris@366: double back[3]; Chris@366: double backim[3]; Chris@366: FFT(3).process(true, re, im, back, backim); Chris@366: COMPARE_ARRAY(back, in); Chris@366: COMPARE_CONST(backim, 0.0); Chris@366: } Chris@366: Chris@366: BOOST_AUTO_TEST_CASE(r_sizes) Chris@366: { Chris@366: // Real supports any even size, but not odd ones Chris@366: Chris@366: BOOST_CHECK_THROW(FFTReal r(3), std::invalid_argument); Chris@366: Chris@366: double in[] = { 1, 1, 1, 1, 1, 1 }; Chris@366: double re[] = { 999, 999, 999, 999, 999, 999 }; Chris@366: double im[] = { 999, 999, 999, 999, 999, 999 }; Chris@366: FFTReal(6).forward(in, re, im); Chris@366: BOOST_CHECK_EQUAL(re[0], 6.0); Chris@366: BOOST_CHECK_EQUAL(re[1], 0.0); Chris@366: BOOST_CHECK_EQUAL(re[2], 0.0); Chris@366: BOOST_CHECK_EQUAL(re[3], 0.0); Chris@366: BOOST_CHECK_EQUAL(re[4], 0.0); Chris@366: BOOST_CHECK_EQUAL(re[5], 0.0); Chris@366: COMPARE_CONST(im, 0.0); Chris@366: double back[6]; Chris@366: FFTReal(6).inverse(re, im, back); Chris@366: COMPARE_ARRAY(back, in); Chris@366: } Chris@366: Chris@366: BOOST_AUTO_TEST_SUITE_END() Chris@366: