Mercurial > hg > qm-dsp

annotate tests/TestFFT.cpp @ 339:9c8ee77db9de

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