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