To check out this repository please hg clone the following URL, or open the URL using EasyMercurial or your preferred Mercurial client.

The primary repository for this project is hosted at https://github.com/cannam/constant-q-cpp/ .
This repository is a read-only copy which is updated automatically every hour.

Statistics Download as Zip
| Branch: | Revision:

root / test / TestFFT.cpp

History | View | Annotate | Download (13.4 KB)

1
/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
2

    
3
#include "dsp/FFT.h"
4

    
5
#define BOOST_TEST_DYN_LINK
6
#define BOOST_TEST_MAIN
7

    
8
#include <boost/test/unit_test.hpp>
9

    
10
#include <stdexcept>
11

    
12
BOOST_AUTO_TEST_SUITE(TestFFT)
13

    
14
#define COMPARE_CONST(a, n) \
15
    for (int cmp_i = 0; cmp_i < (int)(sizeof(a)/sizeof(a[0])); ++cmp_i) { \
16
        BOOST_CHECK_SMALL(a[cmp_i] - n, 1e-14);                                \
17
    }
18

    
19
#define COMPARE_ARRAY(a, b)                                                \
20
    for (int cmp_i = 0; cmp_i < (int)(sizeof(a)/sizeof(a[0])); ++cmp_i) { \
21
        BOOST_CHECK_SMALL(a[cmp_i] - b[cmp_i], 1e-14);                        \
22
    }
23

    
24
//!!! need at least one test with complex time-domain signal
25

    
26
BOOST_AUTO_TEST_CASE(forwardArrayBounds)
27
{
28
    // initialise bins to something recognisable, so we can tell if
29
    // they haven't been written; and allocate the inputs on the heap
30
    // so that, if running under valgrind, we get warnings about
31
    // overruns
32
    double *in = new double[4];
33
    in[0] = 1;
34
    in[1] = 1;
35
    in[2] = -1;
36
    in[3] = -1;
37
    double re[] = { 999, 999, 999, 999, 999, 999 };
38
    double im[] = { 999, 999, 999, 999, 999, 999 };
39
    FFT(4).process(false, in, 0, re+1, im+1);
40
    // And check we haven't overrun the arrays
41
    BOOST_CHECK_EQUAL(re[0], 999.0);
42
    BOOST_CHECK_EQUAL(im[0], 999.0);
43
    BOOST_CHECK_EQUAL(re[5], 999.0);
44
    BOOST_CHECK_EQUAL(im[5], 999.0);
45
    delete[] in;
46
}
47

    
48
BOOST_AUTO_TEST_CASE(r_forwardArrayBounds)
49
{
50
    // initialise bins to something recognisable, so we can tell if
51
    // they haven't been written; and allocate the inputs on the heap
52
    // so that, if running under valgrind, we get warnings about
53
    // overruns
54
    double *in = new double[4];
55
    in[0] = 1;
56
    in[1] = 1;
57
    in[2] = -1;
58
    in[3] = -1;
59
    double re[] = { 999, 999, 999, 999, 999, 999 };
60
    double im[] = { 999, 999, 999, 999, 999, 999 };
61
    FFTReal(4).forward(in, re+1, im+1);
62
    // And check we haven't overrun the arrays
63
    BOOST_CHECK_EQUAL(re[0], 999.0);
64
    BOOST_CHECK_EQUAL(im[0], 999.0);
65
    BOOST_CHECK_EQUAL(re[5], 999.0);
66
    BOOST_CHECK_EQUAL(im[5], 999.0);
67
    delete[] in;
68
}
69

    
70
BOOST_AUTO_TEST_CASE(inverseArrayBounds)
71
{
72
    // initialise bins to something recognisable, so we can tell if
73
    // they haven't been written; and allocate the inputs on the heap
74
    // so that, if running under valgrind, we get warnings about
75
    // overruns
76
    double *re = new double[4];
77
    double *im = new double[4];
78
    re[0] = 0;
79
    re[1] = 1;
80
    re[2] = 0;
81
    re[3] = 1;
82
    im[0] = 0;
83
    im[1] = -2;
84
    im[2] = 0;
85
    im[3] = 2;
86
    double outre[] = { 999, 999, 999, 999, 999, 999 };
87
    double outim[] = { 999, 999, 999, 999, 999, 999 };
88
    FFT(4).process(true, re, im, outre+1, outim+1);
89
    // And check we haven't overrun the arrays
90
    BOOST_CHECK_EQUAL(outre[0], 999.0);
91
    BOOST_CHECK_EQUAL(outim[0], 999.0);
92
    BOOST_CHECK_EQUAL(outre[5], 999.0);
93
    BOOST_CHECK_EQUAL(outim[5], 999.0);
94
    delete[] re;
95
    delete[] im;
96
}
97

    
98
BOOST_AUTO_TEST_CASE(r_inverseArrayBounds)
99
{
100
    // initialise bins to something recognisable, so we can tell if
101
    // they haven't been written; and allocate the inputs on the heap
102
    // so that, if running under valgrind, we get warnings about
103
    // overruns
104
    double *re = new double[3];
105
    double *im = new double[3];
106
    re[0] = 0;
107
    re[1] = 1;
108
    re[2] = 0;
109
    im[0] = 0;
110
    im[1] = -2;
111
    im[2] = 0;
112
    double outre[] = { 999, 999, 999, 999, 999, 999 };
113
    FFTReal(4).inverse(re, im, outre+1);
114
    // And check we haven't overrun the arrays
115
    BOOST_CHECK_EQUAL(outre[0], 999.0);
116
    BOOST_CHECK_EQUAL(outre[5], 999.0);
117
    delete[] re;
118
    delete[] im;
119
}
120

    
121
BOOST_AUTO_TEST_CASE(dc)
122
{
123
    // DC-only signal. The DC bin is purely real
124
    double in[] = { 1, 1, 1, 1 };
125
    double re[] = { 999, 999, 999, 999 };
126
    double im[] = { 999, 999, 999, 999 };
127
    FFT(4).process(false, in, 0, re, im);
128
    BOOST_CHECK_EQUAL(re[0], 4.0);
129
    BOOST_CHECK_EQUAL(re[1], 0.0);
130
    BOOST_CHECK_EQUAL(re[2], 0.0);
131
    BOOST_CHECK_EQUAL(re[3], 0.0);
132
    COMPARE_CONST(im, 0.0);
133
    double back[4];
134
    double backim[4];
135
    FFT(4).process(true, re, im, back, backim);
136
    COMPARE_ARRAY(back, in);
137
    COMPARE_CONST(backim, 0.0);
138
}
139

    
140
BOOST_AUTO_TEST_CASE(r_dc)
141
{
142
    // DC-only signal. The DC bin is purely real
143
    double in[] = { 1, 1, 1, 1 };
144
    double re[] = { 999, 999, 999, 999 };
145
    double im[] = { 999, 999, 999, 999 };
146
    FFTReal(4).forward(in, re, im);
147
    BOOST_CHECK_EQUAL(re[0], 4.0);
148
    BOOST_CHECK_EQUAL(re[1], 0.0);
149
    BOOST_CHECK_EQUAL(re[2], 0.0);
150
    BOOST_CHECK_EQUAL(re[3], 0.0);
151
    COMPARE_CONST(im, 0.0);
152
    double back[4];
153
    // check conjugates are reconstructed
154
    re[3] = 999;
155
    im[3] = 999;
156
    FFTReal(4).inverse(re, im, back);
157
    COMPARE_ARRAY(back, in);
158
}
159

    
160
BOOST_AUTO_TEST_CASE(c_dc)
161
{
162
    // DC-only signal. The DC bin is purely real
163
    double rin[] = { 1, 1, 1, 1 };
164
    double iin[] = { 1, 1, 1, 1 };
165
    double re[] = { 999, 999, 999, 999 };
166
    double im[] = { 999, 999, 999, 999 };
167
    FFT(4).process(false, rin, iin, re, im);
168
    BOOST_CHECK_EQUAL(re[0], 4.0);
169
    BOOST_CHECK_EQUAL(re[1], 0.0);
170
    BOOST_CHECK_EQUAL(re[2], 0.0);
171
    BOOST_CHECK_EQUAL(re[3], 0.0);
172
    BOOST_CHECK_EQUAL(im[0], 4.0);
173
    BOOST_CHECK_EQUAL(im[1], 0.0);
174
    BOOST_CHECK_EQUAL(im[2], 0.0);
175
    BOOST_CHECK_EQUAL(im[3], 0.0);
176
    double back[4];
177
    double backim[4];
178
    FFT(4).process(true, re, im, back, backim);
179
    COMPARE_ARRAY(back, rin);
180
    COMPARE_ARRAY(backim, iin);
181
}
182

    
183
BOOST_AUTO_TEST_CASE(sine)
184
{
185
    // Sine. Output is purely imaginary
186
    double in[] = { 0, 1, 0, -1 };
187
    double re[] = { 999, 999, 999, 999 };
188
    double im[] = { 999, 999, 999, 999 };
189
    FFT(4).process(false, in, 0, re, im);
190
    COMPARE_CONST(re, 0.0);
191
    BOOST_CHECK_EQUAL(im[0], 0.0);
192
    BOOST_CHECK_EQUAL(im[1], -2.0);
193
    BOOST_CHECK_EQUAL(im[2], 0.0);
194
    BOOST_CHECK_EQUAL(im[3], 2.0);
195
    double back[4];
196
    double backim[4];
197
    FFT(4).process(true, re, im, back, backim);
198
    COMPARE_ARRAY(back, in);
199
    COMPARE_CONST(backim, 0.0);
200
}
201

    
202
BOOST_AUTO_TEST_CASE(r_sine)
203
{
204
    // Sine. Output is purely imaginary
205
    double in[] = { 0, 1, 0, -1 };
206
    double re[] = { 999, 999, 999, 999 };
207
    double im[] = { 999, 999, 999, 999 };
208
    FFTReal(4).forward(in, re, im);
209
    COMPARE_CONST(re, 0.0);
210
    BOOST_CHECK_EQUAL(im[0], 0.0);
211
    BOOST_CHECK_EQUAL(im[1], -2.0);
212
    BOOST_CHECK_EQUAL(im[2], 0.0);
213
    BOOST_CHECK_EQUAL(im[3], 2.0);
214
    double back[4];
215
    // check conjugates are reconstructed
216
    re[3] = 999;
217
    im[3] = 999;
218
    FFTReal(4).inverse(re, im, back);
219
    COMPARE_ARRAY(back, in);
220
}
221

    
222
BOOST_AUTO_TEST_CASE(cosine)
223
{
224
    // Cosine. Output is purely real
225
    double in[] = { 1, 0, -1, 0 };
226
    double re[] = { 999, 999, 999, 999 };
227
    double im[] = { 999, 999, 999, 999 };
228
    FFT(4).process(false, in, 0, re, im);
229
    BOOST_CHECK_EQUAL(re[0], 0.0);
230
    BOOST_CHECK_EQUAL(re[1], 2.0);
231
    BOOST_CHECK_EQUAL(re[2], 0.0);
232
    BOOST_CHECK_EQUAL(re[3], 2.0);
233
    COMPARE_CONST(im, 0.0);
234
    double back[4];
235
    double backim[4];
236
    FFT(4).process(true, re, im, back, backim);
237
    COMPARE_ARRAY(back, in);
238
    COMPARE_CONST(backim, 0.0);
239
}
240

    
241
BOOST_AUTO_TEST_CASE(r_cosine)
242
{
243
    // Cosine. Output is purely real
244
    double in[] = { 1, 0, -1, 0 };
245
    double re[] = { 999, 999, 999, 999 };
246
    double im[] = { 999, 999, 999, 999 };
247
    FFTReal(4).forward(in, re, im);
248
    BOOST_CHECK_EQUAL(re[0], 0.0);
249
    BOOST_CHECK_EQUAL(re[1], 2.0);
250
    BOOST_CHECK_EQUAL(re[2], 0.0);
251
    BOOST_CHECK_EQUAL(re[3], 2.0);
252
    COMPARE_CONST(im, 0.0);
253
    double back[4];
254
    // check conjugates are reconstructed
255
    re[3] = 999;
256
    im[3] = 999;
257
    FFTReal(4).inverse(re, im, back);
258
    COMPARE_ARRAY(back, in);
259
}
260

    
261
BOOST_AUTO_TEST_CASE(c_cosine)
262
{
263
    // Cosine. Output is purely real
264
    double rin[] = { 1, 0, -1, 0 };
265
    double iin[] = { 1, 0, -1, 0 };
266
    double re[] = { 999, 999, 999, 999 };
267
    double im[] = { 999, 999, 999, 999 };
268
    FFT(4).process(false, rin, iin, re, im);
269
    BOOST_CHECK_EQUAL(re[0], 0.0);
270
    BOOST_CHECK_EQUAL(re[1], 2.0);
271
    BOOST_CHECK_EQUAL(re[2], 0.0);
272
    BOOST_CHECK_EQUAL(re[3], 2.0);
273
    BOOST_CHECK_EQUAL(im[0], 0.0);
274
    BOOST_CHECK_EQUAL(im[1], 2.0);
275
    BOOST_CHECK_EQUAL(im[2], 0.0);
276
    BOOST_CHECK_EQUAL(im[3], 2.0);
277
    double back[4];
278
    double backim[4];
279
    FFT(4).process(true, re, im, back, backim);
280
    COMPARE_ARRAY(back, rin);
281
    COMPARE_ARRAY(backim, iin);
282
}
283
        
284
BOOST_AUTO_TEST_CASE(sineCosine)
285
{
286
    // Sine and cosine mixed
287
    double in[] = { 0.5, 1, -0.5, -1 };
288
    double re[] = { 999, 999, 999, 999 };
289
    double im[] = { 999, 999, 999, 999 };
290
    FFT(4).process(false, in, 0, re, im);
291
    BOOST_CHECK_EQUAL(re[0], 0.0);
292
    BOOST_CHECK_CLOSE(re[1], 1.0, 1e-12);
293
    BOOST_CHECK_EQUAL(re[2], 0.0);
294
    BOOST_CHECK_CLOSE(re[3], 1.0, 1e-12);
295
    BOOST_CHECK_EQUAL(im[0], 0.0);
296
    BOOST_CHECK_CLOSE(im[1], -2.0, 1e-12);
297
    BOOST_CHECK_EQUAL(im[2], 0.0);
298
    BOOST_CHECK_CLOSE(im[3], 2.0, 1e-12);
299
    double back[4];
300
    double backim[4];
301
    FFT(4).process(true, re, im, back, backim);
302
    COMPARE_ARRAY(back, in);
303
    COMPARE_CONST(backim, 0.0);
304
}
305
        
306
BOOST_AUTO_TEST_CASE(r_sineCosine)
307
{
308
    // Sine and cosine mixed
309
    double in[] = { 0.5, 1, -0.5, -1 };
310
    double re[] = { 999, 999, 999, 999 };
311
    double im[] = { 999, 999, 999, 999 };
312
    FFTReal(4).forward(in, re, im);
313
    BOOST_CHECK_EQUAL(re[0], 0.0);
314
    BOOST_CHECK_CLOSE(re[1], 1.0, 1e-12);
315
    BOOST_CHECK_EQUAL(re[2], 0.0);
316
    BOOST_CHECK_CLOSE(re[3], 1.0, 1e-12);
317
    BOOST_CHECK_EQUAL(im[0], 0.0);
318
    BOOST_CHECK_CLOSE(im[1], -2.0, 1e-12);
319
    BOOST_CHECK_EQUAL(im[2], 0.0);
320
    BOOST_CHECK_CLOSE(im[3], 2.0, 1e-12);
321
    double back[4];
322
    // check conjugates are reconstructed
323
    re[3] = 999;
324
    im[3] = 999;
325
    FFTReal(4).inverse(re, im, back);
326
    COMPARE_ARRAY(back, in);
327
}
328
        
329
BOOST_AUTO_TEST_CASE(c_sineCosine)
330
{
331
    double rin[] = { 1, 0, -1, 0 };
332
    double iin[] = { 0, 1, 0, -1 };
333
    double re[] = { 999, 999, 999, 999 };
334
    double im[] = { 999, 999, 999, 999 };
335
    FFT(4).process(false, rin, iin, re, im);
336
    BOOST_CHECK_EQUAL(re[0], 0.0);
337
    BOOST_CHECK_EQUAL(re[1], 4.0);
338
    BOOST_CHECK_EQUAL(re[2], 0.0);
339
    BOOST_CHECK_EQUAL(re[3], 0.0);
340
    COMPARE_CONST(im, 0.0);
341
    double back[4];
342
    double backim[4];
343
    FFT(4).process(true, re, im, back, backim);
344
    COMPARE_ARRAY(back, rin);
345
    COMPARE_ARRAY(backim, iin);
346
}
347

    
348
BOOST_AUTO_TEST_CASE(nyquist)
349
{
350
    double in[] = { 1, -1, 1, -1 };
351
    double re[] = { 999, 999, 999, 999 };
352
    double im[] = { 999, 999, 999, 999 };
353
    FFT(4).process(false, in, 0, re, im);
354
    BOOST_CHECK_EQUAL(re[0], 0.0);
355
    BOOST_CHECK_EQUAL(re[1], 0.0);
356
    BOOST_CHECK_EQUAL(re[2], 4.0);
357
    BOOST_CHECK_EQUAL(re[3], 0.0);
358
    COMPARE_CONST(im, 0.0);
359
    double back[4];
360
    double backim[4];
361
    FFT(4).process(true, re, im, back, backim);
362
    COMPARE_ARRAY(back, in);
363
    COMPARE_CONST(backim, 0.0);
364
}
365

    
366
BOOST_AUTO_TEST_CASE(r_nyquist)
367
{
368
    double in[] = { 1, -1, 1, -1 };
369
    double re[] = { 999, 999, 999, 999 };
370
    double im[] = { 999, 999, 999, 999 };
371
    FFTReal(4).forward(in, re, im);
372
    BOOST_CHECK_EQUAL(re[0], 0.0);
373
    BOOST_CHECK_EQUAL(re[1], 0.0);
374
    BOOST_CHECK_EQUAL(re[2], 4.0);
375
    BOOST_CHECK_EQUAL(re[3], 0.0);
376
    COMPARE_CONST(im, 0.0);
377
    double back[4];
378
    // check conjugates are reconstructed
379
    re[3] = 999;
380
    im[3] = 999;
381
    FFTReal(4).inverse(re, im, back);
382
    COMPARE_ARRAY(back, in);
383
}
384

    
385
BOOST_AUTO_TEST_CASE(dirac)
386
{
387
    double in[] = { 1, 0, 0, 0 };
388
    double re[] = { 999, 999, 999, 999 };
389
    double im[] = { 999, 999, 999, 999 };
390
    FFT(4).process(false, in, 0, re, im);
391
    BOOST_CHECK_EQUAL(re[0], 1.0);
392
    BOOST_CHECK_EQUAL(re[1], 1.0);
393
    BOOST_CHECK_EQUAL(re[2], 1.0);
394
    BOOST_CHECK_EQUAL(re[3], 1.0);
395
    COMPARE_CONST(im, 0.0);
396
    double back[4];
397
    double backim[4];
398
    FFT(4).process(true, re, im, back, backim);
399
    COMPARE_ARRAY(back, in);
400
    COMPARE_CONST(backim, 0.0);
401
}
402

    
403
BOOST_AUTO_TEST_CASE(r_dirac)
404
{
405
    double in[] = { 1, 0, 0, 0 };
406
    double re[] = { 999, 999, 999, 999 };
407
    double im[] = { 999, 999, 999, 999 };
408
    FFTReal(4).forward(in, re, im);
409
    BOOST_CHECK_EQUAL(re[0], 1.0);
410
    BOOST_CHECK_EQUAL(re[1], 1.0);
411
    BOOST_CHECK_EQUAL(re[2], 1.0);
412
    BOOST_CHECK_EQUAL(re[3], 1.0);
413
    COMPARE_CONST(im, 0.0);
414
    double back[4];
415
    // check conjugates are reconstructed
416
    re[3] = 999;
417
    im[3] = 999;
418
    FFTReal(4).inverse(re, im, back);
419
    COMPARE_ARRAY(back, in);
420
}
421

    
422
BOOST_AUTO_TEST_CASE(sizes) 
423
{
424
    // Complex supports any size. A single test with an odd size
425
    // will do here, without getting too much into our expectations
426
    // about supported butterflies etc
427

    
428
    double in[] = { 1, 1, 1 };
429
    double re[] = { 999, 999, 999 };
430
    double im[] = { 999, 999, 999 };
431
    FFT(3).process(false, in, 0, re, im);
432
    BOOST_CHECK_EQUAL(re[0], 3.0);
433
    BOOST_CHECK_EQUAL(re[1], 0.0);
434
    BOOST_CHECK_EQUAL(re[2], 0.0);
435
    COMPARE_CONST(im, 0.0);
436
    double back[3];
437
    double backim[3];
438
    FFT(3).process(true, re, im, back, backim);
439
    COMPARE_ARRAY(back, in);
440
    COMPARE_CONST(backim, 0.0);
441
}
442

    
443
BOOST_AUTO_TEST_CASE(r_sizes)
444
{
445
    // Real supports any even size, but not odd ones
446

    
447
    BOOST_CHECK_THROW(FFTReal r(3), std::invalid_argument);
448

    
449
    double in[] = { 1, 1, 1, 1, 1, 1 };
450
    double re[] = { 999, 999, 999, 999, 999, 999 };
451
    double im[] = { 999, 999, 999, 999, 999, 999 };
452
    FFTReal(6).forward(in, re, im);
453
    BOOST_CHECK_EQUAL(re[0], 6.0);
454
    BOOST_CHECK_EQUAL(re[1], 0.0);
455
    BOOST_CHECK_EQUAL(re[2], 0.0);
456
    BOOST_CHECK_EQUAL(re[3], 0.0);
457
    BOOST_CHECK_EQUAL(re[4], 0.0);
458
    BOOST_CHECK_EQUAL(re[5], 0.0);
459
    COMPARE_CONST(im, 0.0);
460
    double back[6];
461
    FFTReal(6).inverse(re, im, back);
462
    COMPARE_ARRAY(back, in);
463
}
464

    
465
BOOST_AUTO_TEST_SUITE_END()
466