Mercurial > hg > qm-dsp

annotate tests/TestFFT.cpp @ 343:24d8ea972643

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Don't use tabs
author Chris Cannam <c.cannam@qmul.ac.uk>
date Wed, 02 Oct 2013 18:22:06 +0100
parents 9c8ee77db9de
children 6ec45e85ed81
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@335 10 BOOST_AUTO_TEST_SUITE(TestFFT)
c@335 11
c@335 12 #define COMPARE_CONST(a, n) \
c@335 13 for (int cmp_i = 0; cmp_i < (int)(sizeof(a)/sizeof(a[0])); ++cmp_i) { \
c@335 14 BOOST_CHECK_SMALL(a[cmp_i] - n, 1e-14); \
c@335 15 }
c@335 16
c@335 17 #define COMPARE_ARRAY(a, b) \
c@335 18 for (int cmp_i = 0; cmp_i < (int)(sizeof(a)/sizeof(a[0])); ++cmp_i) { \
c@335 19 BOOST_CHECK_SMALL(a[cmp_i] - b[cmp_i], 1e-14); \
c@335 20 }
c@335 21
c@339 22 //!!! need at least one test with complex time-domain signal
c@335 23
c@335 24 BOOST_AUTO_TEST_CASE(forwardArrayBounds)
c@335 25 {
c@335 26 // initialise bins to something recognisable, so we can tell
c@335 27 // if they haven't been written
c@335 28 double in[] = { 1, 1, -1, -1 };
c@335 29 double re[] = { 999, 999, 999, 999, 999, 999 };
c@335 30 double im[] = { 999, 999, 999, 999, 999, 999 };
c@335 31 FFT(4).process(false, in, 0, re+1, im+1);
c@335 32 // And check we haven't overrun the arrays
c@335 33 BOOST_CHECK_EQUAL(re[0], 999.0);
c@335 34 BOOST_CHECK_EQUAL(im[0], 999.0);
c@335 35 BOOST_CHECK_EQUAL(re[5], 999.0);
c@335 36 BOOST_CHECK_EQUAL(im[5], 999.0);
c@335 37 }
c@335 38
c@339 39 BOOST_AUTO_TEST_CASE(r_forwardArrayBounds)
c@339 40 {
c@339 41 // initialise bins to something recognisable, so we can tell
c@339 42 // if they haven't been written
c@339 43 double in[] = { 1, 1, -1, -1 };
c@339 44 double re[] = { 999, 999, 999, 999, 999, 999 };
c@339 45 double im[] = { 999, 999, 999, 999, 999, 999 };
c@339 46 FFTReal(4).forward(in, re+1, im+1);
c@339 47 // And check we haven't overrun the arrays
c@339 48 BOOST_CHECK_EQUAL(re[0], 999.0);
c@339 49 BOOST_CHECK_EQUAL(im[0], 999.0);
c@339 50 BOOST_CHECK_EQUAL(re[5], 999.0);
c@339 51 BOOST_CHECK_EQUAL(im[5], 999.0);
c@339 52 }
c@339 53
c@335 54 BOOST_AUTO_TEST_CASE(inverseArrayBounds)
c@335 55 {
c@335 56 // initialise bins to something recognisable, so we can tell
c@335 57 // if they haven't been written
c@339 58 double re[] = { 0, 1, 0, 1 };
c@339 59 double im[] = { 0, -2, 0, 2 };
c@335 60 double outre[] = { 999, 999, 999, 999, 999, 999 };
c@335 61 double outim[] = { 999, 999, 999, 999, 999, 999 };
c@339 62 FFT(4).process(true, re, im, outre+1, outim+1);
c@335 63 // And check we haven't overrun the arrays
c@335 64 BOOST_CHECK_EQUAL(outre[0], 999.0);
c@335 65 BOOST_CHECK_EQUAL(outim[0], 999.0);
c@335 66 BOOST_CHECK_EQUAL(outre[5], 999.0);
c@335 67 BOOST_CHECK_EQUAL(outim[5], 999.0);
c@335 68 }
c@335 69
c@339 70 BOOST_AUTO_TEST_CASE(r_inverseArrayBounds)
c@339 71 {
c@339 72 // initialise bins to something recognisable, so we can tell
c@339 73 // if they haven't been written
c@339 74 double re[] = { 0, 1, 0 };
c@339 75 double im[] = { 0, -2, 0 };
c@339 76 double outre[] = { 999, 999, 999, 999, 999, 999 };
c@339 77 FFTReal(4).inverse(re, im, outre+1);
c@339 78 // And check we haven't overrun the arrays
c@339 79 BOOST_CHECK_EQUAL(outre[0], 999.0);
c@339 80 BOOST_CHECK_EQUAL(outre[5], 999.0);
c@339 81 }
c@339 82
c@339 83 BOOST_AUTO_TEST_CASE(dc)
c@339 84 {
c@339 85 // DC-only signal. The DC bin is purely real
c@339 86 double in[] = { 1, 1, 1, 1 };
c@339 87 double re[] = { 999, 999, 999, 999 };
c@339 88 double im[] = { 999, 999, 999, 999 };
c@339 89 FFT(4).process(false, in, 0, re, im);
c@339 90 BOOST_CHECK_EQUAL(re[0], 4.0);
c@339 91 BOOST_CHECK_EQUAL(re[1], 0.0);
c@339 92 BOOST_CHECK_EQUAL(re[2], 0.0);
c@339 93 BOOST_CHECK_EQUAL(re[3], 0.0);
c@339 94 COMPARE_CONST(im, 0.0);
c@339 95 double back[4];
c@339 96 double backim[4];
c@339 97 FFT(4).process(true, re, im, back, backim);
c@339 98 COMPARE_ARRAY(back, in);
c@339 99 COMPARE_CONST(backim, 0.0);
c@339 100 }
c@339 101
c@339 102 BOOST_AUTO_TEST_CASE(r_dc)
c@339 103 {
c@339 104 // DC-only signal. The DC bin is purely real
c@339 105 double in[] = { 1, 1, 1, 1 };
c@339 106 double re[] = { 999, 999, 999, 999 };
c@339 107 double im[] = { 999, 999, 999, 999 };
c@339 108 FFTReal(4).forward(in, re, im);
c@339 109 BOOST_CHECK_EQUAL(re[0], 4.0);
c@339 110 BOOST_CHECK_EQUAL(re[1], 0.0);
c@339 111 BOOST_CHECK_EQUAL(re[2], 0.0);
c@339 112 BOOST_CHECK_EQUAL(re[3], 0.0);
c@339 113 COMPARE_CONST(im, 0.0);
c@339 114 double back[4];
c@339 115 // check conjugates are reconstructed
c@339 116 re[3] = 999;
c@339 117 im[3] = 999;
c@339 118 FFTReal(4).inverse(re, im, back);
c@339 119 COMPARE_ARRAY(back, in);
c@339 120 }
c@339 121
c@339 122 BOOST_AUTO_TEST_CASE(sine)
c@339 123 {
c@339 124 // Sine. Output is purely imaginary
c@339 125 double in[] = { 0, 1, 0, -1 };
c@339 126 double re[] = { 999, 999, 999, 999 };
c@339 127 double im[] = { 999, 999, 999, 999 };
c@339 128 FFT(4).process(false, in, 0, re, im);
c@339 129 COMPARE_CONST(re, 0.0);
c@339 130 BOOST_CHECK_EQUAL(im[0], 0.0);
c@339 131 BOOST_CHECK_EQUAL(im[1], -2.0);
c@339 132 BOOST_CHECK_EQUAL(im[2], 0.0);
c@339 133 BOOST_CHECK_EQUAL(im[3], 2.0);
c@339 134 double back[4];
c@339 135 double backim[4];
c@339 136 FFT(4).process(true, re, im, back, backim);
c@339 137 COMPARE_ARRAY(back, in);
c@339 138 COMPARE_CONST(backim, 0.0);
c@339 139 }
c@339 140
c@339 141 BOOST_AUTO_TEST_CASE(r_sine)
c@339 142 {
c@339 143 // Sine. Output is purely imaginary
c@339 144 double in[] = { 0, 1, 0, -1 };
c@339 145 double re[] = { 999, 999, 999, 999 };
c@339 146 double im[] = { 999, 999, 999, 999 };
c@339 147 FFTReal(4).forward(in, re, im);
c@339 148 COMPARE_CONST(re, 0.0);
c@339 149 BOOST_CHECK_EQUAL(im[0], 0.0);
c@339 150 BOOST_CHECK_EQUAL(im[1], -2.0);
c@339 151 BOOST_CHECK_EQUAL(im[2], 0.0);
c@339 152 BOOST_CHECK_EQUAL(im[3], 2.0);
c@339 153 double back[4];
c@339 154 // check conjugates are reconstructed
c@339 155 re[3] = 999;
c@339 156 im[3] = 999;
c@339 157 FFTReal(4).inverse(re, im, back);
c@339 158 COMPARE_ARRAY(back, in);
c@339 159 }
c@339 160
c@339 161 BOOST_AUTO_TEST_CASE(cosine)
c@339 162 {
c@339 163 // Cosine. Output is purely real
c@339 164 double in[] = { 1, 0, -1, 0 };
c@339 165 double re[] = { 999, 999, 999, 999 };
c@339 166 double im[] = { 999, 999, 999, 999 };
c@339 167 FFT(4).process(false, in, 0, re, im);
c@339 168 BOOST_CHECK_EQUAL(re[0], 0.0);
c@339 169 BOOST_CHECK_EQUAL(re[1], 2.0);
c@339 170 BOOST_CHECK_EQUAL(re[2], 0.0);
c@339 171 BOOST_CHECK_EQUAL(re[3], 2.0);
c@339 172 COMPARE_CONST(im, 0.0);
c@339 173 double back[4];
c@339 174 double backim[4];
c@339 175 FFT(4).process(true, re, im, back, backim);
c@339 176 COMPARE_ARRAY(back, in);
c@339 177 COMPARE_CONST(backim, 0.0);
c@339 178 }
c@339 179
c@339 180 BOOST_AUTO_TEST_CASE(r_cosine)
c@339 181 {
c@339 182 // Cosine. Output is purely real
c@339 183 double in[] = { 1, 0, -1, 0 };
c@339 184 double re[] = { 999, 999, 999, 999 };
c@339 185 double im[] = { 999, 999, 999, 999 };
c@339 186 FFTReal(4).forward(in, re, im);
c@339 187 BOOST_CHECK_EQUAL(re[0], 0.0);
c@339 188 BOOST_CHECK_EQUAL(re[1], 2.0);
c@339 189 BOOST_CHECK_EQUAL(re[2], 0.0);
c@339 190 BOOST_CHECK_EQUAL(re[3], 2.0);
c@339 191 COMPARE_CONST(im, 0.0);
c@339 192 double back[4];
c@339 193 // check conjugates are reconstructed
c@339 194 re[3] = 999;
c@339 195 im[3] = 999;
c@339 196 FFTReal(4).inverse(re, im, back);
c@339 197 COMPARE_ARRAY(back, in);
c@339 198 }
c@339 199
c@339 200 BOOST_AUTO_TEST_CASE(sineCosine)
c@339 201 {
c@339 202 // Sine and cosine mixed
c@339 203 double in[] = { 0.5, 1, -0.5, -1 };
c@339 204 double re[] = { 999, 999, 999, 999 };
c@339 205 double im[] = { 999, 999, 999, 999 };
c@339 206 FFT(4).process(false, in, 0, re, im);
c@339 207 BOOST_CHECK_EQUAL(re[0], 0.0);
c@339 208 BOOST_CHECK_CLOSE(re[1], 1.0, 1e-12);
c@339 209 BOOST_CHECK_EQUAL(re[2], 0.0);
c@339 210 BOOST_CHECK_CLOSE(re[3], 1.0, 1e-12);
c@339 211 BOOST_CHECK_EQUAL(im[0], 0.0);
c@339 212 BOOST_CHECK_CLOSE(im[1], -2.0, 1e-12);
c@339 213 BOOST_CHECK_EQUAL(im[2], 0.0);
c@339 214 BOOST_CHECK_CLOSE(im[3], 2.0, 1e-12);
c@339 215 double back[4];
c@339 216 double backim[4];
c@339 217 FFT(4).process(true, re, im, back, backim);
c@339 218 COMPARE_ARRAY(back, in);
c@339 219 COMPARE_CONST(backim, 0.0);
c@339 220 }
c@339 221
c@339 222 BOOST_AUTO_TEST_CASE(r_sineCosine)
c@339 223 {
c@339 224 // Sine and cosine mixed
c@339 225 double in[] = { 0.5, 1, -0.5, -1 };
c@339 226 double re[] = { 999, 999, 999, 999 };
c@339 227 double im[] = { 999, 999, 999, 999 };
c@339 228 FFTReal(4).forward(in, re, im);
c@339 229 BOOST_CHECK_EQUAL(re[0], 0.0);
c@339 230 BOOST_CHECK_CLOSE(re[1], 1.0, 1e-12);
c@339 231 BOOST_CHECK_EQUAL(re[2], 0.0);
c@339 232 BOOST_CHECK_CLOSE(re[3], 1.0, 1e-12);
c@339 233 BOOST_CHECK_EQUAL(im[0], 0.0);
c@339 234 BOOST_CHECK_CLOSE(im[1], -2.0, 1e-12);
c@339 235 BOOST_CHECK_EQUAL(im[2], 0.0);
c@339 236 BOOST_CHECK_CLOSE(im[3], 2.0, 1e-12);
c@339 237 double back[4];
c@339 238 // check conjugates are reconstructed
c@339 239 re[3] = 999;
c@339 240 im[3] = 999;
c@339 241 FFTReal(4).inverse(re, im, back);
c@339 242 COMPARE_ARRAY(back, in);
c@339 243 }
c@339 244
c@339 245 BOOST_AUTO_TEST_CASE(nyquist)
c@339 246 {
c@339 247 double in[] = { 1, -1, 1, -1 };
c@339 248 double re[] = { 999, 999, 999, 999 };
c@339 249 double im[] = { 999, 999, 999, 999 };
c@339 250 FFT(4).process(false, in, 0, re, im);
c@339 251 BOOST_CHECK_EQUAL(re[0], 0.0);
c@339 252 BOOST_CHECK_EQUAL(re[1], 0.0);
c@339 253 BOOST_CHECK_EQUAL(re[2], 4.0);
c@339 254 BOOST_CHECK_EQUAL(re[3], 0.0);
c@339 255 COMPARE_CONST(im, 0.0);
c@339 256 double back[4];
c@339 257 double backim[4];
c@339 258 FFT(4).process(true, re, im, back, backim);
c@339 259 COMPARE_ARRAY(back, in);
c@339 260 COMPARE_CONST(backim, 0.0);
c@339 261 }
c@339 262
c@339 263 BOOST_AUTO_TEST_CASE(r_nyquist)
c@339 264 {
c@339 265 double in[] = { 1, -1, 1, -1 };
c@339 266 double re[] = { 999, 999, 999, 999 };
c@339 267 double im[] = { 999, 999, 999, 999 };
c@339 268 FFTReal(4).forward(in, re, im);
c@339 269 BOOST_CHECK_EQUAL(re[0], 0.0);
c@339 270 BOOST_CHECK_EQUAL(re[1], 0.0);
c@339 271 BOOST_CHECK_EQUAL(re[2], 4.0);
c@339 272 BOOST_CHECK_EQUAL(re[3], 0.0);
c@339 273 COMPARE_CONST(im, 0.0);
c@339 274 double back[4];
c@339 275 // check conjugates are reconstructed
c@339 276 re[3] = 999;
c@339 277 im[3] = 999;
c@339 278 FFTReal(4).inverse(re, im, back);
c@339 279 COMPARE_ARRAY(back, in);
c@339 280 }
c@339 281
c@339 282 BOOST_AUTO_TEST_CASE(dirac)
c@339 283 {
c@339 284 double in[] = { 1, 0, 0, 0 };
c@339 285 double re[] = { 999, 999, 999, 999 };
c@339 286 double im[] = { 999, 999, 999, 999 };
c@339 287 FFT(4).process(false, in, 0, re, im);
c@339 288 BOOST_CHECK_EQUAL(re[0], 1.0);
c@339 289 BOOST_CHECK_EQUAL(re[1], 1.0);
c@339 290 BOOST_CHECK_EQUAL(re[2], 1.0);
c@339 291 BOOST_CHECK_EQUAL(re[3], 1.0);
c@339 292 COMPARE_CONST(im, 0.0);
c@339 293 double back[4];
c@339 294 double backim[4];
c@339 295 FFT(4).process(true, re, im, back, backim);
c@339 296 COMPARE_ARRAY(back, in);
c@339 297 COMPARE_CONST(backim, 0.0);
c@339 298 }
c@339 299
c@339 300 BOOST_AUTO_TEST_CASE(r_dirac)
c@339 301 {
c@339 302 double in[] = { 1, 0, 0, 0 };
c@339 303 double re[] = { 999, 999, 999, 999 };
c@339 304 double im[] = { 999, 999, 999, 999 };
c@339 305 FFTReal(4).forward(in, re, im);
c@339 306 BOOST_CHECK_EQUAL(re[0], 1.0);
c@339 307 BOOST_CHECK_EQUAL(re[1], 1.0);
c@339 308 BOOST_CHECK_EQUAL(re[2], 1.0);
c@339 309 BOOST_CHECK_EQUAL(re[3], 1.0);
c@339 310 COMPARE_CONST(im, 0.0);
c@339 311 double back[4];
c@339 312 // check conjugates are reconstructed
c@339 313 re[3] = 999;
c@339 314 im[3] = 999;
c@339 315 FFTReal(4).inverse(re, im, back);
c@339 316 COMPARE_ARRAY(back, in);
c@339 317 }
c@339 318
c@335 319 BOOST_AUTO_TEST_SUITE_END()
c@335 320