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@355
|
10 #include <stdexcept>
|
c@355
|
11
|
c@335
|
12 BOOST_AUTO_TEST_SUITE(TestFFT)
|
c@335
|
13
|
c@335
|
14 #define COMPARE_CONST(a, n) \
|
c@335
|
15 for (int cmp_i = 0; cmp_i < (int)(sizeof(a)/sizeof(a[0])); ++cmp_i) { \
|
c@335
|
16 BOOST_CHECK_SMALL(a[cmp_i] - n, 1e-14); \
|
c@335
|
17 }
|
c@335
|
18
|
c@335
|
19 #define COMPARE_ARRAY(a, b) \
|
c@335
|
20 for (int cmp_i = 0; cmp_i < (int)(sizeof(a)/sizeof(a[0])); ++cmp_i) { \
|
c@335
|
21 BOOST_CHECK_SMALL(a[cmp_i] - b[cmp_i], 1e-14); \
|
c@335
|
22 }
|
c@335
|
23
|
c@339
|
24 //!!! need at least one test with complex time-domain signal
|
c@335
|
25
|
c@335
|
26 BOOST_AUTO_TEST_CASE(forwardArrayBounds)
|
c@335
|
27 {
|
c@335
|
28 // initialise bins to something recognisable, so we can tell
|
c@335
|
29 // if they haven't been written
|
c@335
|
30 double in[] = { 1, 1, -1, -1 };
|
c@335
|
31 double re[] = { 999, 999, 999, 999, 999, 999 };
|
c@335
|
32 double im[] = { 999, 999, 999, 999, 999, 999 };
|
c@335
|
33 FFT(4).process(false, in, 0, re+1, im+1);
|
c@335
|
34 // And check we haven't overrun the arrays
|
c@335
|
35 BOOST_CHECK_EQUAL(re[0], 999.0);
|
c@335
|
36 BOOST_CHECK_EQUAL(im[0], 999.0);
|
c@335
|
37 BOOST_CHECK_EQUAL(re[5], 999.0);
|
c@335
|
38 BOOST_CHECK_EQUAL(im[5], 999.0);
|
c@335
|
39 }
|
c@335
|
40
|
c@339
|
41 BOOST_AUTO_TEST_CASE(r_forwardArrayBounds)
|
c@339
|
42 {
|
c@339
|
43 // initialise bins to something recognisable, so we can tell
|
c@339
|
44 // if they haven't been written
|
c@339
|
45 double in[] = { 1, 1, -1, -1 };
|
c@339
|
46 double re[] = { 999, 999, 999, 999, 999, 999 };
|
c@339
|
47 double im[] = { 999, 999, 999, 999, 999, 999 };
|
c@339
|
48 FFTReal(4).forward(in, re+1, im+1);
|
c@339
|
49 // And check we haven't overrun the arrays
|
c@339
|
50 BOOST_CHECK_EQUAL(re[0], 999.0);
|
c@339
|
51 BOOST_CHECK_EQUAL(im[0], 999.0);
|
c@339
|
52 BOOST_CHECK_EQUAL(re[5], 999.0);
|
c@339
|
53 BOOST_CHECK_EQUAL(im[5], 999.0);
|
c@339
|
54 }
|
c@339
|
55
|
c@335
|
56 BOOST_AUTO_TEST_CASE(inverseArrayBounds)
|
c@335
|
57 {
|
c@335
|
58 // initialise bins to something recognisable, so we can tell
|
c@335
|
59 // if they haven't been written
|
c@339
|
60 double re[] = { 0, 1, 0, 1 };
|
c@339
|
61 double im[] = { 0, -2, 0, 2 };
|
c@335
|
62 double outre[] = { 999, 999, 999, 999, 999, 999 };
|
c@335
|
63 double outim[] = { 999, 999, 999, 999, 999, 999 };
|
c@339
|
64 FFT(4).process(true, re, im, outre+1, outim+1);
|
c@335
|
65 // And check we haven't overrun the arrays
|
c@335
|
66 BOOST_CHECK_EQUAL(outre[0], 999.0);
|
c@335
|
67 BOOST_CHECK_EQUAL(outim[0], 999.0);
|
c@335
|
68 BOOST_CHECK_EQUAL(outre[5], 999.0);
|
c@335
|
69 BOOST_CHECK_EQUAL(outim[5], 999.0);
|
c@335
|
70 }
|
c@335
|
71
|
c@339
|
72 BOOST_AUTO_TEST_CASE(r_inverseArrayBounds)
|
c@339
|
73 {
|
c@339
|
74 // initialise bins to something recognisable, so we can tell
|
c@339
|
75 // if they haven't been written
|
c@339
|
76 double re[] = { 0, 1, 0 };
|
c@339
|
77 double im[] = { 0, -2, 0 };
|
c@339
|
78 double outre[] = { 999, 999, 999, 999, 999, 999 };
|
c@339
|
79 FFTReal(4).inverse(re, im, outre+1);
|
c@339
|
80 // And check we haven't overrun the arrays
|
c@339
|
81 BOOST_CHECK_EQUAL(outre[0], 999.0);
|
c@339
|
82 BOOST_CHECK_EQUAL(outre[5], 999.0);
|
c@339
|
83 }
|
c@339
|
84
|
c@339
|
85 BOOST_AUTO_TEST_CASE(dc)
|
c@339
|
86 {
|
c@339
|
87 // DC-only signal. The DC bin is purely real
|
c@339
|
88 double in[] = { 1, 1, 1, 1 };
|
c@339
|
89 double re[] = { 999, 999, 999, 999 };
|
c@339
|
90 double im[] = { 999, 999, 999, 999 };
|
c@339
|
91 FFT(4).process(false, in, 0, re, im);
|
c@339
|
92 BOOST_CHECK_EQUAL(re[0], 4.0);
|
c@339
|
93 BOOST_CHECK_EQUAL(re[1], 0.0);
|
c@339
|
94 BOOST_CHECK_EQUAL(re[2], 0.0);
|
c@339
|
95 BOOST_CHECK_EQUAL(re[3], 0.0);
|
c@339
|
96 COMPARE_CONST(im, 0.0);
|
c@339
|
97 double back[4];
|
c@339
|
98 double backim[4];
|
c@339
|
99 FFT(4).process(true, re, im, back, backim);
|
c@339
|
100 COMPARE_ARRAY(back, in);
|
c@339
|
101 COMPARE_CONST(backim, 0.0);
|
c@339
|
102 }
|
c@339
|
103
|
c@339
|
104 BOOST_AUTO_TEST_CASE(r_dc)
|
c@339
|
105 {
|
c@339
|
106 // DC-only signal. The DC bin is purely real
|
c@339
|
107 double in[] = { 1, 1, 1, 1 };
|
c@339
|
108 double re[] = { 999, 999, 999, 999 };
|
c@339
|
109 double im[] = { 999, 999, 999, 999 };
|
c@339
|
110 FFTReal(4).forward(in, re, im);
|
c@339
|
111 BOOST_CHECK_EQUAL(re[0], 4.0);
|
c@339
|
112 BOOST_CHECK_EQUAL(re[1], 0.0);
|
c@339
|
113 BOOST_CHECK_EQUAL(re[2], 0.0);
|
c@339
|
114 BOOST_CHECK_EQUAL(re[3], 0.0);
|
c@339
|
115 COMPARE_CONST(im, 0.0);
|
c@339
|
116 double back[4];
|
c@339
|
117 // check conjugates are reconstructed
|
c@339
|
118 re[3] = 999;
|
c@339
|
119 im[3] = 999;
|
c@339
|
120 FFTReal(4).inverse(re, im, back);
|
c@339
|
121 COMPARE_ARRAY(back, in);
|
c@339
|
122 }
|
c@339
|
123
|
c@339
|
124 BOOST_AUTO_TEST_CASE(sine)
|
c@339
|
125 {
|
c@339
|
126 // Sine. Output is purely imaginary
|
c@339
|
127 double in[] = { 0, 1, 0, -1 };
|
c@339
|
128 double re[] = { 999, 999, 999, 999 };
|
c@339
|
129 double im[] = { 999, 999, 999, 999 };
|
c@339
|
130 FFT(4).process(false, in, 0, re, im);
|
c@339
|
131 COMPARE_CONST(re, 0.0);
|
c@339
|
132 BOOST_CHECK_EQUAL(im[0], 0.0);
|
c@339
|
133 BOOST_CHECK_EQUAL(im[1], -2.0);
|
c@339
|
134 BOOST_CHECK_EQUAL(im[2], 0.0);
|
c@339
|
135 BOOST_CHECK_EQUAL(im[3], 2.0);
|
c@339
|
136 double back[4];
|
c@339
|
137 double backim[4];
|
c@339
|
138 FFT(4).process(true, re, im, back, backim);
|
c@339
|
139 COMPARE_ARRAY(back, in);
|
c@339
|
140 COMPARE_CONST(backim, 0.0);
|
c@339
|
141 }
|
c@339
|
142
|
c@339
|
143 BOOST_AUTO_TEST_CASE(r_sine)
|
c@339
|
144 {
|
c@339
|
145 // Sine. Output is purely imaginary
|
c@339
|
146 double in[] = { 0, 1, 0, -1 };
|
c@339
|
147 double re[] = { 999, 999, 999, 999 };
|
c@339
|
148 double im[] = { 999, 999, 999, 999 };
|
c@339
|
149 FFTReal(4).forward(in, re, im);
|
c@339
|
150 COMPARE_CONST(re, 0.0);
|
c@339
|
151 BOOST_CHECK_EQUAL(im[0], 0.0);
|
c@339
|
152 BOOST_CHECK_EQUAL(im[1], -2.0);
|
c@339
|
153 BOOST_CHECK_EQUAL(im[2], 0.0);
|
c@339
|
154 BOOST_CHECK_EQUAL(im[3], 2.0);
|
c@339
|
155 double back[4];
|
c@339
|
156 // check conjugates are reconstructed
|
c@339
|
157 re[3] = 999;
|
c@339
|
158 im[3] = 999;
|
c@339
|
159 FFTReal(4).inverse(re, im, back);
|
c@339
|
160 COMPARE_ARRAY(back, in);
|
c@339
|
161 }
|
c@339
|
162
|
c@339
|
163 BOOST_AUTO_TEST_CASE(cosine)
|
c@339
|
164 {
|
c@339
|
165 // Cosine. Output is purely real
|
c@339
|
166 double in[] = { 1, 0, -1, 0 };
|
c@339
|
167 double re[] = { 999, 999, 999, 999 };
|
c@339
|
168 double im[] = { 999, 999, 999, 999 };
|
c@339
|
169 FFT(4).process(false, in, 0, re, im);
|
c@339
|
170 BOOST_CHECK_EQUAL(re[0], 0.0);
|
c@339
|
171 BOOST_CHECK_EQUAL(re[1], 2.0);
|
c@339
|
172 BOOST_CHECK_EQUAL(re[2], 0.0);
|
c@339
|
173 BOOST_CHECK_EQUAL(re[3], 2.0);
|
c@339
|
174 COMPARE_CONST(im, 0.0);
|
c@339
|
175 double back[4];
|
c@339
|
176 double backim[4];
|
c@339
|
177 FFT(4).process(true, re, im, back, backim);
|
c@339
|
178 COMPARE_ARRAY(back, in);
|
c@339
|
179 COMPARE_CONST(backim, 0.0);
|
c@339
|
180 }
|
c@339
|
181
|
c@339
|
182 BOOST_AUTO_TEST_CASE(r_cosine)
|
c@339
|
183 {
|
c@339
|
184 // Cosine. Output is purely real
|
c@339
|
185 double in[] = { 1, 0, -1, 0 };
|
c@339
|
186 double re[] = { 999, 999, 999, 999 };
|
c@339
|
187 double im[] = { 999, 999, 999, 999 };
|
c@339
|
188 FFTReal(4).forward(in, re, im);
|
c@339
|
189 BOOST_CHECK_EQUAL(re[0], 0.0);
|
c@339
|
190 BOOST_CHECK_EQUAL(re[1], 2.0);
|
c@339
|
191 BOOST_CHECK_EQUAL(re[2], 0.0);
|
c@339
|
192 BOOST_CHECK_EQUAL(re[3], 2.0);
|
c@339
|
193 COMPARE_CONST(im, 0.0);
|
c@339
|
194 double back[4];
|
c@339
|
195 // check conjugates are reconstructed
|
c@339
|
196 re[3] = 999;
|
c@339
|
197 im[3] = 999;
|
c@339
|
198 FFTReal(4).inverse(re, im, back);
|
c@339
|
199 COMPARE_ARRAY(back, in);
|
c@339
|
200 }
|
c@339
|
201
|
c@339
|
202 BOOST_AUTO_TEST_CASE(sineCosine)
|
c@339
|
203 {
|
c@339
|
204 // Sine and cosine mixed
|
c@339
|
205 double in[] = { 0.5, 1, -0.5, -1 };
|
c@339
|
206 double re[] = { 999, 999, 999, 999 };
|
c@339
|
207 double im[] = { 999, 999, 999, 999 };
|
c@339
|
208 FFT(4).process(false, in, 0, re, im);
|
c@339
|
209 BOOST_CHECK_EQUAL(re[0], 0.0);
|
c@339
|
210 BOOST_CHECK_CLOSE(re[1], 1.0, 1e-12);
|
c@339
|
211 BOOST_CHECK_EQUAL(re[2], 0.0);
|
c@339
|
212 BOOST_CHECK_CLOSE(re[3], 1.0, 1e-12);
|
c@339
|
213 BOOST_CHECK_EQUAL(im[0], 0.0);
|
c@339
|
214 BOOST_CHECK_CLOSE(im[1], -2.0, 1e-12);
|
c@339
|
215 BOOST_CHECK_EQUAL(im[2], 0.0);
|
c@339
|
216 BOOST_CHECK_CLOSE(im[3], 2.0, 1e-12);
|
c@339
|
217 double back[4];
|
c@339
|
218 double backim[4];
|
c@339
|
219 FFT(4).process(true, re, im, back, backim);
|
c@339
|
220 COMPARE_ARRAY(back, in);
|
c@339
|
221 COMPARE_CONST(backim, 0.0);
|
c@339
|
222 }
|
c@339
|
223
|
c@339
|
224 BOOST_AUTO_TEST_CASE(r_sineCosine)
|
c@339
|
225 {
|
c@339
|
226 // Sine and cosine mixed
|
c@339
|
227 double in[] = { 0.5, 1, -0.5, -1 };
|
c@339
|
228 double re[] = { 999, 999, 999, 999 };
|
c@339
|
229 double im[] = { 999, 999, 999, 999 };
|
c@339
|
230 FFTReal(4).forward(in, re, im);
|
c@339
|
231 BOOST_CHECK_EQUAL(re[0], 0.0);
|
c@339
|
232 BOOST_CHECK_CLOSE(re[1], 1.0, 1e-12);
|
c@339
|
233 BOOST_CHECK_EQUAL(re[2], 0.0);
|
c@339
|
234 BOOST_CHECK_CLOSE(re[3], 1.0, 1e-12);
|
c@339
|
235 BOOST_CHECK_EQUAL(im[0], 0.0);
|
c@339
|
236 BOOST_CHECK_CLOSE(im[1], -2.0, 1e-12);
|
c@339
|
237 BOOST_CHECK_EQUAL(im[2], 0.0);
|
c@339
|
238 BOOST_CHECK_CLOSE(im[3], 2.0, 1e-12);
|
c@339
|
239 double back[4];
|
c@339
|
240 // check conjugates are reconstructed
|
c@339
|
241 re[3] = 999;
|
c@339
|
242 im[3] = 999;
|
c@339
|
243 FFTReal(4).inverse(re, im, back);
|
c@339
|
244 COMPARE_ARRAY(back, in);
|
c@339
|
245 }
|
c@339
|
246
|
c@339
|
247 BOOST_AUTO_TEST_CASE(nyquist)
|
c@339
|
248 {
|
c@339
|
249 double in[] = { 1, -1, 1, -1 };
|
c@339
|
250 double re[] = { 999, 999, 999, 999 };
|
c@339
|
251 double im[] = { 999, 999, 999, 999 };
|
c@339
|
252 FFT(4).process(false, in, 0, re, im);
|
c@339
|
253 BOOST_CHECK_EQUAL(re[0], 0.0);
|
c@339
|
254 BOOST_CHECK_EQUAL(re[1], 0.0);
|
c@339
|
255 BOOST_CHECK_EQUAL(re[2], 4.0);
|
c@339
|
256 BOOST_CHECK_EQUAL(re[3], 0.0);
|
c@339
|
257 COMPARE_CONST(im, 0.0);
|
c@339
|
258 double back[4];
|
c@339
|
259 double backim[4];
|
c@339
|
260 FFT(4).process(true, re, im, back, backim);
|
c@339
|
261 COMPARE_ARRAY(back, in);
|
c@339
|
262 COMPARE_CONST(backim, 0.0);
|
c@339
|
263 }
|
c@339
|
264
|
c@339
|
265 BOOST_AUTO_TEST_CASE(r_nyquist)
|
c@339
|
266 {
|
c@339
|
267 double in[] = { 1, -1, 1, -1 };
|
c@339
|
268 double re[] = { 999, 999, 999, 999 };
|
c@339
|
269 double im[] = { 999, 999, 999, 999 };
|
c@339
|
270 FFTReal(4).forward(in, re, im);
|
c@339
|
271 BOOST_CHECK_EQUAL(re[0], 0.0);
|
c@339
|
272 BOOST_CHECK_EQUAL(re[1], 0.0);
|
c@339
|
273 BOOST_CHECK_EQUAL(re[2], 4.0);
|
c@339
|
274 BOOST_CHECK_EQUAL(re[3], 0.0);
|
c@339
|
275 COMPARE_CONST(im, 0.0);
|
c@339
|
276 double back[4];
|
c@339
|
277 // check conjugates are reconstructed
|
c@339
|
278 re[3] = 999;
|
c@339
|
279 im[3] = 999;
|
c@339
|
280 FFTReal(4).inverse(re, im, back);
|
c@339
|
281 COMPARE_ARRAY(back, in);
|
c@339
|
282 }
|
c@339
|
283
|
c@339
|
284 BOOST_AUTO_TEST_CASE(dirac)
|
c@339
|
285 {
|
c@339
|
286 double in[] = { 1, 0, 0, 0 };
|
c@339
|
287 double re[] = { 999, 999, 999, 999 };
|
c@339
|
288 double im[] = { 999, 999, 999, 999 };
|
c@339
|
289 FFT(4).process(false, in, 0, re, im);
|
c@339
|
290 BOOST_CHECK_EQUAL(re[0], 1.0);
|
c@339
|
291 BOOST_CHECK_EQUAL(re[1], 1.0);
|
c@339
|
292 BOOST_CHECK_EQUAL(re[2], 1.0);
|
c@339
|
293 BOOST_CHECK_EQUAL(re[3], 1.0);
|
c@339
|
294 COMPARE_CONST(im, 0.0);
|
c@339
|
295 double back[4];
|
c@339
|
296 double backim[4];
|
c@339
|
297 FFT(4).process(true, re, im, back, backim);
|
c@339
|
298 COMPARE_ARRAY(back, in);
|
c@339
|
299 COMPARE_CONST(backim, 0.0);
|
c@339
|
300 }
|
c@339
|
301
|
c@339
|
302 BOOST_AUTO_TEST_CASE(r_dirac)
|
c@339
|
303 {
|
c@339
|
304 double in[] = { 1, 0, 0, 0 };
|
c@339
|
305 double re[] = { 999, 999, 999, 999 };
|
c@339
|
306 double im[] = { 999, 999, 999, 999 };
|
c@339
|
307 FFTReal(4).forward(in, re, im);
|
c@339
|
308 BOOST_CHECK_EQUAL(re[0], 1.0);
|
c@339
|
309 BOOST_CHECK_EQUAL(re[1], 1.0);
|
c@339
|
310 BOOST_CHECK_EQUAL(re[2], 1.0);
|
c@339
|
311 BOOST_CHECK_EQUAL(re[3], 1.0);
|
c@339
|
312 COMPARE_CONST(im, 0.0);
|
c@339
|
313 double back[4];
|
c@339
|
314 // check conjugates are reconstructed
|
c@339
|
315 re[3] = 999;
|
c@339
|
316 im[3] = 999;
|
c@339
|
317 FFTReal(4).inverse(re, im, back);
|
c@339
|
318 COMPARE_ARRAY(back, in);
|
c@339
|
319 }
|
c@339
|
320
|
c@355
|
321 BOOST_AUTO_TEST_CASE(sizes)
|
c@355
|
322 {
|
c@355
|
323 // Complex supports any size. A single test with an odd size
|
c@355
|
324 // will do here, without getting too much into our expectations
|
c@355
|
325 // about supported butterflies etc
|
c@355
|
326
|
c@355
|
327 double in[] = { 1, 1, 1 };
|
c@355
|
328 double re[] = { 999, 999, 999 };
|
c@355
|
329 double im[] = { 999, 999, 999 };
|
c@355
|
330 FFT(3).process(false, in, 0, re, im);
|
c@355
|
331 BOOST_CHECK_EQUAL(re[0], 3.0);
|
c@355
|
332 BOOST_CHECK_EQUAL(re[1], 0.0);
|
c@355
|
333 BOOST_CHECK_EQUAL(re[2], 0.0);
|
c@355
|
334 COMPARE_CONST(im, 0.0);
|
c@355
|
335 double back[3];
|
c@355
|
336 double backim[3];
|
c@355
|
337 FFT(3).process(true, re, im, back, backim);
|
c@355
|
338 COMPARE_ARRAY(back, in);
|
c@355
|
339 COMPARE_CONST(backim, 0.0);
|
c@355
|
340 }
|
c@355
|
341
|
c@355
|
342 BOOST_AUTO_TEST_CASE(r_sizes)
|
c@355
|
343 {
|
c@355
|
344 // Real supports any even size, but not odd ones
|
c@355
|
345
|
c@355
|
346 BOOST_CHECK_THROW(FFTReal r(3), std::invalid_argument);
|
c@355
|
347
|
c@355
|
348 double in[] = { 1, 1, 1, 1, 1, 1 };
|
c@355
|
349 double re[] = { 999, 999, 999, 999, 999, 999 };
|
c@355
|
350 double im[] = { 999, 999, 999, 999, 999, 999 };
|
c@355
|
351 FFTReal(6).forward(in, re, im);
|
c@355
|
352 BOOST_CHECK_EQUAL(re[0], 6.0);
|
c@355
|
353 BOOST_CHECK_EQUAL(re[1], 0.0);
|
c@355
|
354 BOOST_CHECK_EQUAL(re[2], 0.0);
|
c@355
|
355 BOOST_CHECK_EQUAL(re[3], 0.0);
|
c@355
|
356 BOOST_CHECK_EQUAL(re[4], 0.0);
|
c@355
|
357 BOOST_CHECK_EQUAL(re[5], 0.0);
|
c@355
|
358 COMPARE_CONST(im, 0.0);
|
c@355
|
359 double back[6];
|
c@355
|
360 FFTReal(6).inverse(re, im, back);
|
c@355
|
361 COMPARE_ARRAY(back, in);
|
c@355
|
362 }
|
c@355
|
363
|
c@335
|
364 BOOST_AUTO_TEST_SUITE_END()
|
c@335
|
365
|