|
1 /***************************************************************************** |
|
2 |
|
3 TestAccuracy.hpp |
|
4 Copyright (c) 2005 Laurent de Soras |
|
5 |
|
6 --- Legal stuff --- |
|
7 |
|
8 This library is free software; you can redistribute it and/or |
|
9 modify it under the terms of the GNU Lesser General Public |
|
10 License as published by the Free Software Foundation; either |
|
11 version 2.1 of the License, or (at your option) any later version. |
|
12 |
|
13 This library is distributed in the hope that it will be useful, |
|
14 but WITHOUT ANY WARRANTY; without even the implied warranty of |
|
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
|
16 Lesser General Public License for more details. |
|
17 |
|
18 You should have received a copy of the GNU Lesser General Public |
|
19 License along with this library; if not, write to the Free Software |
|
20 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA |
|
21 |
|
22 *Tab=3***********************************************************************/ |
|
23 |
|
24 |
|
25 |
|
26 #if defined (TestAccuracy_CURRENT_CODEHEADER) |
|
27 #error Recursive inclusion of TestAccuracy code header. |
|
28 #endif |
|
29 #define TestAccuracy_CURRENT_CODEHEADER |
|
30 |
|
31 #if ! defined (TestAccuracy_CODEHEADER_INCLUDED) |
|
32 #define TestAccuracy_CODEHEADER_INCLUDED |
|
33 |
|
34 |
|
35 |
|
36 /*\\\ INCLUDE FILES \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/ |
|
37 |
|
38 #include "def.h" |
|
39 #include "test_fnc.h" |
|
40 #include "TestWhiteNoiseGen.h" |
|
41 |
|
42 #include <typeinfo> |
|
43 #include <vector> |
|
44 |
|
45 #include <cmath> |
|
46 #include <cstdio> |
|
47 |
|
48 |
|
49 |
|
50 static const double TestAccuracy_LN10 = 2.3025850929940456840179914546844; |
|
51 |
|
52 |
|
53 |
|
54 /*\\\ PUBLIC \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/ |
|
55 |
|
56 |
|
57 |
|
58 template <class FO> |
|
59 int TestAccuracy <FO>::perform_test_single_object (FO &fft) |
|
60 { |
|
61 assert (&fft != 0); |
|
62 |
|
63 using namespace std; |
|
64 |
|
65 int ret_val = 0; |
|
66 |
|
67 const std::type_info & ti = typeid (fft); |
|
68 const char * class_name_0 = ti.name (); |
|
69 |
|
70 if (ret_val == 0) |
|
71 { |
|
72 ret_val = perform_test_d (fft, class_name_0); |
|
73 } |
|
74 if (ret_val == 0) |
|
75 { |
|
76 ret_val = perform_test_i (fft, class_name_0); |
|
77 } |
|
78 if (ret_val == 0) |
|
79 { |
|
80 ret_val = perform_test_di (fft, class_name_0); |
|
81 } |
|
82 |
|
83 if (ret_val == 0) |
|
84 { |
|
85 printf ("\n"); |
|
86 } |
|
87 |
|
88 return (ret_val); |
|
89 } |
|
90 |
|
91 |
|
92 |
|
93 template <class FO> |
|
94 int TestAccuracy <FO>::perform_test_d (FO &fft, const char *class_name_0) |
|
95 { |
|
96 assert (&fft != 0); |
|
97 assert (class_name_0 != 0); |
|
98 |
|
99 using namespace std; |
|
100 |
|
101 int ret_val = 0; |
|
102 const long len = fft.get_length (); |
|
103 const long nbr_tests = limit ( |
|
104 NBR_ACC_TESTS / len / len, |
|
105 1L, |
|
106 static_cast <long> (MAX_NBR_TESTS) |
|
107 ); |
|
108 |
|
109 printf ("Testing %s::do_fft () [%ld samples]... ", class_name_0, len); |
|
110 fflush (stdout); |
|
111 TestWhiteNoiseGen <DataType> noise; |
|
112 std::vector <DataType> x (len); |
|
113 std::vector <DataType> s1 (len); |
|
114 std::vector <DataType> s2 (len); |
|
115 BigFloat err_avg = 0; |
|
116 |
|
117 for (long test = 0; test < nbr_tests && ret_val == 0; ++ test) |
|
118 { |
|
119 noise.generate (&x [0], len); |
|
120 fft.do_fft (&s1 [0], &x [0]); |
|
121 compute_tf (&s2 [0], &x [0], len); |
|
122 |
|
123 BigFloat max_err; |
|
124 compare_vect_display (&s1 [0], &s2 [0], len, max_err); |
|
125 err_avg += max_err; |
|
126 } |
|
127 err_avg /= NBR_ACC_TESTS; |
|
128 |
|
129 printf ("done.\n"); |
|
130 printf ( |
|
131 "Average maximum error: %.6f %% (%f dB)\n", |
|
132 static_cast <double> (err_avg * 100), |
|
133 static_cast <double> ((20 / TestAccuracy_LN10) * log (err_avg)) |
|
134 ); |
|
135 |
|
136 return (ret_val); |
|
137 } |
|
138 |
|
139 |
|
140 |
|
141 template <class FO> |
|
142 int TestAccuracy <FO>::perform_test_i (FO &fft, const char *class_name_0) |
|
143 { |
|
144 assert (&fft != 0); |
|
145 assert (class_name_0 != 0); |
|
146 |
|
147 using namespace std; |
|
148 |
|
149 int ret_val = 0; |
|
150 const long len = fft.get_length (); |
|
151 const long nbr_tests = limit ( |
|
152 NBR_ACC_TESTS / len / len, |
|
153 10L, |
|
154 static_cast <long> (MAX_NBR_TESTS) |
|
155 ); |
|
156 |
|
157 printf ("Testing %s::do_ifft () [%ld samples]... ", class_name_0, len); |
|
158 fflush (stdout); |
|
159 TestWhiteNoiseGen <DataType> noise; |
|
160 std::vector <DataType> s (len); |
|
161 std::vector <DataType> x1 (len); |
|
162 std::vector <DataType> x2 (len); |
|
163 BigFloat err_avg = 0; |
|
164 |
|
165 for (long test = 0; test < nbr_tests && ret_val == 0; ++ test) |
|
166 { |
|
167 noise.generate (&s [0], len); |
|
168 fft.do_ifft (&s [0], &x1 [0]); |
|
169 compute_itf (&x2 [0], &s [0], len); |
|
170 |
|
171 BigFloat max_err; |
|
172 compare_vect_display (&x1 [0], &x2 [0], len, max_err); |
|
173 err_avg += max_err; |
|
174 } |
|
175 err_avg /= NBR_ACC_TESTS; |
|
176 |
|
177 printf ("done.\n"); |
|
178 printf ( |
|
179 "Average maximum error: %.6f %% (%f dB)\n", |
|
180 static_cast <double> (err_avg * 100), |
|
181 static_cast <double> ((20 / TestAccuracy_LN10) * log (err_avg)) |
|
182 ); |
|
183 |
|
184 return (ret_val); |
|
185 } |
|
186 |
|
187 |
|
188 |
|
189 template <class FO> |
|
190 int TestAccuracy <FO>::perform_test_di (FO &fft, const char *class_name_0) |
|
191 { |
|
192 assert (&fft != 0); |
|
193 assert (class_name_0 != 0); |
|
194 |
|
195 using namespace std; |
|
196 |
|
197 int ret_val = 0; |
|
198 const long len = fft.get_length (); |
|
199 const long nbr_tests = limit ( |
|
200 NBR_ACC_TESTS / len / len, |
|
201 1L, |
|
202 static_cast <long> (MAX_NBR_TESTS) |
|
203 ); |
|
204 |
|
205 printf ( |
|
206 "Testing %s::do_fft () / do_ifft () / rescale () [%ld samples]... ", |
|
207 class_name_0, |
|
208 len |
|
209 ); |
|
210 fflush (stdout); |
|
211 TestWhiteNoiseGen <DataType> noise; |
|
212 std::vector <DataType> x (len); |
|
213 std::vector <DataType> s (len); |
|
214 std::vector <DataType> y (len); |
|
215 BigFloat err_avg = 0; |
|
216 |
|
217 for (long test = 0; test < nbr_tests && ret_val == 0; ++ test) |
|
218 { |
|
219 noise.generate (&x [0], len); |
|
220 fft.do_fft (&s [0], &x [0]); |
|
221 fft.do_ifft (&s [0], &y [0]); |
|
222 fft.rescale (&y [0]); |
|
223 |
|
224 BigFloat max_err; |
|
225 compare_vect_display (&x [0], &y [0], len, max_err); |
|
226 err_avg += max_err; |
|
227 } |
|
228 err_avg /= NBR_ACC_TESTS; |
|
229 |
|
230 printf ("done.\n"); |
|
231 printf ( |
|
232 "Average maximum error: %.6f %% (%f dB)\n", |
|
233 static_cast <double> (err_avg * 100), |
|
234 static_cast <double> ((20 / TestAccuracy_LN10) * log (err_avg)) |
|
235 ); |
|
236 |
|
237 return (ret_val); |
|
238 } |
|
239 |
|
240 |
|
241 |
|
242 /*\\\ PROTECTED \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/ |
|
243 |
|
244 |
|
245 |
|
246 /*\\\ PRIVATE \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/ |
|
247 |
|
248 |
|
249 |
|
250 // Positive transform |
|
251 template <class FO> |
|
252 void TestAccuracy <FO>::compute_tf (DataType s [], const DataType x [], long length) |
|
253 { |
|
254 assert (s != 0); |
|
255 assert (x != 0); |
|
256 assert (length >= 2); |
|
257 assert ((length & 1) == 0); |
|
258 |
|
259 const long nbr_bins = length >> 1; |
|
260 |
|
261 // DC and Nyquist |
|
262 BigFloat dc = 0; |
|
263 BigFloat ny = 0; |
|
264 for (long pos = 0; pos < length; pos += 2) |
|
265 { |
|
266 const BigFloat even = x [pos ]; |
|
267 const BigFloat odd = x [pos + 1]; |
|
268 dc += even + odd; |
|
269 ny += even - odd; |
|
270 } |
|
271 s [0 ] = static_cast <DataType> (dc); |
|
272 s [nbr_bins] = static_cast <DataType> (ny); |
|
273 |
|
274 // Regular bins |
|
275 for (long bin = 1; bin < nbr_bins; ++ bin) |
|
276 { |
|
277 BigFloat sum_r = 0; |
|
278 BigFloat sum_i = 0; |
|
279 |
|
280 const BigFloat m = bin * static_cast <BigFloat> (2 * PI) / length; |
|
281 |
|
282 for (long pos = 0; pos < length; ++pos) |
|
283 { |
|
284 using namespace std; |
|
285 |
|
286 const BigFloat phase = pos * m; |
|
287 const BigFloat e_r = cos (phase); |
|
288 const BigFloat e_i = sin (phase); |
|
289 |
|
290 sum_r += x [pos] * e_r; |
|
291 sum_i += x [pos] * e_i; |
|
292 } |
|
293 |
|
294 s [ bin] = static_cast <DataType> (sum_r); |
|
295 s [nbr_bins + bin] = static_cast <DataType> (sum_i); |
|
296 } |
|
297 } |
|
298 |
|
299 |
|
300 |
|
301 // Negative transform |
|
302 template <class FO> |
|
303 void TestAccuracy <FO>::compute_itf (DataType x [], const DataType s [], long length) |
|
304 { |
|
305 assert (s != 0); |
|
306 assert (x != 0); |
|
307 assert (length >= 2); |
|
308 assert ((length & 1) == 0); |
|
309 |
|
310 const long nbr_bins = length >> 1; |
|
311 |
|
312 // DC and Nyquist |
|
313 BigFloat dc = s [0 ]; |
|
314 BigFloat ny = s [nbr_bins]; |
|
315 |
|
316 // Regular bins |
|
317 for (long pos = 0; pos < length; ++pos) |
|
318 { |
|
319 BigFloat sum = dc + ny * (1 - 2 * (pos & 1)); |
|
320 |
|
321 const BigFloat m = pos * static_cast <BigFloat> (-2 * PI) / length; |
|
322 |
|
323 for (long bin = 1; bin < nbr_bins; ++ bin) |
|
324 { |
|
325 using namespace std; |
|
326 |
|
327 const BigFloat phase = bin * m; |
|
328 const BigFloat e_r = cos (phase); |
|
329 const BigFloat e_i = sin (phase); |
|
330 |
|
331 sum += 2 * ( e_r * s [bin ] |
|
332 - e_i * s [bin + nbr_bins]); |
|
333 } |
|
334 |
|
335 x [pos] = static_cast <DataType> (sum); |
|
336 } |
|
337 } |
|
338 |
|
339 |
|
340 |
|
341 template <class FO> |
|
342 int TestAccuracy <FO>::compare_vect_display (const DataType x_ptr [], const DataType y_ptr [], long len, BigFloat &max_err_rel) |
|
343 { |
|
344 assert (x_ptr != 0); |
|
345 assert (y_ptr != 0); |
|
346 assert (len > 0); |
|
347 assert (&max_err_rel != 0); |
|
348 |
|
349 using namespace std; |
|
350 |
|
351 int ret_val = 0; |
|
352 |
|
353 BigFloat power = compute_power (&x_ptr [0], &y_ptr [0], len); |
|
354 BigFloat power_dif; |
|
355 long max_err_pos; |
|
356 compare_vect (&x_ptr [0], &y_ptr [0], power_dif, max_err_pos, len); |
|
357 |
|
358 if (power == 0) |
|
359 { |
|
360 power = power_dif; |
|
361 } |
|
362 const BigFloat power_err_rel = power_dif / power; |
|
363 |
|
364 BigFloat max_err = 0; |
|
365 max_err_rel = 0; |
|
366 if (max_err_pos >= 0) |
|
367 { |
|
368 max_err = y_ptr [max_err_pos] - x_ptr [max_err_pos]; |
|
369 max_err_rel = 2 * fabs (max_err) / ( fabs (y_ptr [max_err_pos]) |
|
370 + fabs (x_ptr [max_err_pos])); |
|
371 } |
|
372 |
|
373 if (power_err_rel > 0.001) |
|
374 { |
|
375 printf ("Power error : %f (%.6f %%)\n", |
|
376 static_cast <double> (power_err_rel), |
|
377 static_cast <double> (power_err_rel * 100) |
|
378 ); |
|
379 if (max_err_pos >= 0) |
|
380 { |
|
381 printf ( |
|
382 "Maximum error: %f - %f = %f (%f)\n", |
|
383 static_cast <double> (y_ptr [max_err_pos]), |
|
384 static_cast <double> (x_ptr [max_err_pos]), |
|
385 static_cast <double> (max_err), |
|
386 static_cast <double> (max_err_pos) |
|
387 ); |
|
388 } |
|
389 } |
|
390 |
|
391 return (ret_val); |
|
392 } |
|
393 |
|
394 |
|
395 |
|
396 template <class FO> |
|
397 typename TestAccuracy <FO>::BigFloat TestAccuracy <FO>::compute_power (const DataType x_ptr [], long len) |
|
398 { |
|
399 assert (x_ptr != 0); |
|
400 assert (len > 0); |
|
401 |
|
402 BigFloat power = 0; |
|
403 for (long pos = 0; pos < len; ++pos) |
|
404 { |
|
405 const BigFloat val = x_ptr [pos]; |
|
406 |
|
407 power += val * val; |
|
408 } |
|
409 |
|
410 using namespace std; |
|
411 |
|
412 power = sqrt (power) / len; |
|
413 |
|
414 return (power); |
|
415 } |
|
416 |
|
417 |
|
418 |
|
419 template <class FO> |
|
420 typename TestAccuracy <FO>::BigFloat TestAccuracy <FO>::compute_power (const DataType x_ptr [], const DataType y_ptr [], long len) |
|
421 { |
|
422 assert (x_ptr != 0); |
|
423 assert (y_ptr != 0); |
|
424 assert (len > 0); |
|
425 |
|
426 return ((compute_power (x_ptr, len) + compute_power (y_ptr, len)) * 0.5); |
|
427 } |
|
428 |
|
429 |
|
430 |
|
431 template <class FO> |
|
432 void TestAccuracy <FO>::compare_vect (const DataType x_ptr [], const DataType y_ptr [], BigFloat &power, long &max_err_pos, long len) |
|
433 { |
|
434 assert (x_ptr != 0); |
|
435 assert (y_ptr != 0); |
|
436 assert (len > 0); |
|
437 assert (&power != 0); |
|
438 assert (&max_err_pos != 0); |
|
439 |
|
440 power = 0; |
|
441 BigFloat max_dif2 = 0; |
|
442 max_err_pos = -1; |
|
443 |
|
444 for (long pos = 0; pos < len; ++pos) |
|
445 { |
|
446 const BigFloat x = x_ptr [pos]; |
|
447 const BigFloat y = y_ptr [pos]; |
|
448 const BigFloat dif = y - x; |
|
449 const BigFloat dif2 = dif * dif; |
|
450 |
|
451 power += dif2; |
|
452 if (dif2 > max_dif2) |
|
453 { |
|
454 max_err_pos = pos; |
|
455 max_dif2 = dif2; |
|
456 } |
|
457 } |
|
458 |
|
459 using namespace std; |
|
460 |
|
461 power = sqrt (power) / len; |
|
462 } |
|
463 |
|
464 |
|
465 |
|
466 #endif // TestAccuracy_CODEHEADER_INCLUDED |
|
467 |
|
468 #undef TestAccuracy_CURRENT_CODEHEADER |
|
469 |
|
470 |
|
471 |
|
472 /*\\\ EOF \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/ |