source: trunk/demos/spectrum/3rdparty/fftreal/TestAccuracy.hpp

Last change on this file was 769, checked in by Dmitry A. Kuminov, 15 years ago

trunk: Merged in qt 4.6.3 sources from branches/vendor/nokia/qt.
File size: 10.2 KB
Line 
1/*****************************************************************************
2
3 TestAccuracy.hpp
4 Copyright (c) 2005 Laurent de Soras
5
6--- Legal stuff ---
7
8This library is free software; you can redistribute it and/or
9modify it under the terms of the GNU Lesser General Public
10License as published by the Free Software Foundation; either
11version 2.1 of the License, or (at your option) any later version.
12
13This library is distributed in the hope that it will be useful,
14but WITHOUT ANY WARRANTY; without even the implied warranty of
15MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16Lesser General Public License for more details.
17
18You should have received a copy of the GNU Lesser General Public
19License along with this library; if not, write to the Free Software
20Foundation, 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
50static const double TestAccuracy_LN10 = 2.3025850929940456840179914546844;
51
52
53
54/*\\\ PUBLIC \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
55
56
57
58template <class FO>
59int 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
93template <class FO>
94int 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
141template <class FO>
142int 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
189template <class FO>
190int 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
251template <class FO>
252void 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
302template <class FO>
303void 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
341template <class FO>
342int 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
396template <class FO>
397typename 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
419template <class FO>
420typename 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
431template <class FO>
432void 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 \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
Note: See TracBrowser for help on using the repository browser.