Mercurial > hg > qm-dsp

comparison tests/TestFFT.cpp @ 114:f6ccde089491 pvoc

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