source: trunk/demos/spectrum/3rdparty/fftreal/FFTReal.hpp@ 1004

Last change on this file since 1004 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: 20.0 KB
Line 
1/*****************************************************************************
2
3 FFTReal.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 (FFTReal_CURRENT_CODEHEADER)
27 #error Recursive inclusion of FFTReal code header.
28#endif
29#define FFTReal_CURRENT_CODEHEADER
30
31#if ! defined (FFTReal_CODEHEADER_INCLUDED)
32#define FFTReal_CODEHEADER_INCLUDED
33
34
35
36/*\\\ INCLUDE FILES \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
37
38#include <cassert>
39#include <cmath>
40
41
42
43static inline bool FFTReal_is_pow2 (long x)
44{
45 assert (x > 0);
46
47 return ((x & -x) == x);
48}
49
50
51
52static inline int FFTReal_get_next_pow2 (long x)
53{
54 --x;
55
56 int p = 0;
57 while ((x & ~0xFFFFL) != 0)
58 {
59 p += 16;
60 x >>= 16;
61 }
62 while ((x & ~0xFL) != 0)
63 {
64 p += 4;
65 x >>= 4;
66 }
67 while (x > 0)
68 {
69 ++p;
70 x >>= 1;
71 }
72
73 return (p);
74}
75
76
77
78/*\\\ PUBLIC \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
79
80
81
82/*
83==============================================================================
84Name: ctor
85Input parameters:
86 - length: length of the array on which we want to do a FFT. Range: power of
87 2 only, > 0.
88Throws: std::bad_alloc
89==============================================================================
90*/
91
92template <class DT>
93FFTReal <DT>::FFTReal (long length)
94: _length (length)
95, _nbr_bits (FFTReal_get_next_pow2 (length))
96, _br_lut ()
97, _trigo_lut ()
98, _buffer (length)
99, _trigo_osc ()
100{
101 assert (FFTReal_is_pow2 (length));
102 assert (_nbr_bits <= MAX_BIT_DEPTH);
103
104 init_br_lut ();
105 init_trigo_lut ();
106 init_trigo_osc ();
107}
108
109
110
111/*
112==============================================================================
113Name: get_length
114Description:
115 Returns the number of points processed by this FFT object.
116Returns: The number of points, power of 2, > 0.
117Throws: Nothing
118==============================================================================
119*/
120
121template <class DT>
122long FFTReal <DT>::get_length () const
123{
124 return (_length);
125}
126
127
128
129/*
130==============================================================================
131Name: do_fft
132Description:
133 Compute the FFT of the array.
134Input parameters:
135 - x: pointer on the source array (time).
136Output parameters:
137 - f: pointer on the destination array (frequencies).
138 f [0...length(x)/2] = real values,
139 f [length(x)/2+1...length(x)-1] = negative imaginary values of
140 coefficents 1...length(x)/2-1.
141Throws: Nothing
142==============================================================================
143*/
144
145template <class DT>
146void FFTReal <DT>::do_fft (DataType f [], const DataType x []) const
147{
148 assert (f != 0);
149 assert (f != use_buffer ());
150 assert (x != 0);
151 assert (x != use_buffer ());
152 assert (x != f);
153
154 // General case
155 if (_nbr_bits > 2)
156 {
157 compute_fft_general (f, x);
158 }
159
160 // 4-point FFT
161 else if (_nbr_bits == 2)
162 {
163 f [1] = x [0] - x [2];
164 f [3] = x [1] - x [3];
165
166 const DataType b_0 = x [0] + x [2];
167 const DataType b_2 = x [1] + x [3];
168
169 f [0] = b_0 + b_2;
170 f [2] = b_0 - b_2;
171 }
172
173 // 2-point FFT
174 else if (_nbr_bits == 1)
175 {
176 f [0] = x [0] + x [1];
177 f [1] = x [0] - x [1];
178 }
179
180 // 1-point FFT
181 else
182 {
183 f [0] = x [0];
184 }
185}
186
187
188
189/*
190==============================================================================
191Name: do_ifft
192Description:
193 Compute the inverse FFT of the array. Note that data must be post-scaled:
194 IFFT (FFT (x)) = x * length (x).
195Input parameters:
196 - f: pointer on the source array (frequencies).
197 f [0...length(x)/2] = real values
198 f [length(x)/2+1...length(x)-1] = negative imaginary values of
199 coefficents 1...length(x)/2-1.
200Output parameters:
201 - x: pointer on the destination array (time).
202Throws: Nothing
203==============================================================================
204*/
205
206template <class DT>
207void FFTReal <DT>::do_ifft (const DataType f [], DataType x []) const
208{
209 assert (f != 0);
210 assert (f != use_buffer ());
211 assert (x != 0);
212 assert (x != use_buffer ());
213 assert (x != f);
214
215 // General case
216 if (_nbr_bits > 2)
217 {
218 compute_ifft_general (f, x);
219 }
220
221 // 4-point IFFT
222 else if (_nbr_bits == 2)
223 {
224 const DataType b_0 = f [0] + f [2];
225 const DataType b_2 = f [0] - f [2];
226
227 x [0] = b_0 + f [1] * 2;
228 x [2] = b_0 - f [1] * 2;
229 x [1] = b_2 + f [3] * 2;
230 x [3] = b_2 - f [3] * 2;
231 }
232
233 // 2-point IFFT
234 else if (_nbr_bits == 1)
235 {
236 x [0] = f [0] + f [1];
237 x [1] = f [0] - f [1];
238 }
239
240 // 1-point IFFT
241 else
242 {
243 x [0] = f [0];
244 }
245}
246
247
248
249/*
250==============================================================================
251Name: rescale
252Description:
253 Scale an array by divide each element by its length. This function should
254 be called after FFT + IFFT.
255Input parameters:
256 - x: pointer on array to rescale (time or frequency).
257Throws: Nothing
258==============================================================================
259*/
260
261template <class DT>
262void FFTReal <DT>::rescale (DataType x []) const
263{
264 const DataType mul = DataType (1.0 / _length);
265
266 if (_length < 4)
267 {
268 long i = _length - 1;
269 do
270 {
271 x [i] *= mul;
272 --i;
273 }
274 while (i >= 0);
275 }
276
277 else
278 {
279 assert ((_length & 3) == 0);
280
281 // Could be optimized with SIMD instruction sets (needs alignment check)
282 long i = _length - 4;
283 do
284 {
285 x [i + 0] *= mul;
286 x [i + 1] *= mul;
287 x [i + 2] *= mul;
288 x [i + 3] *= mul;
289 i -= 4;
290 }
291 while (i >= 0);
292 }
293}
294
295
296
297/*
298==============================================================================
299Name: use_buffer
300Description:
301 Access the internal buffer, whose length is the FFT one.
302 Buffer content will be erased at each do_fft() / do_ifft() call!
303 This buffer cannot be used as:
304 - source for FFT or IFFT done with this object
305 - destination for FFT or IFFT done with this object
306Returns:
307 Buffer start address
308Throws: Nothing
309==============================================================================
310*/
311
312template <class DT>
313typename FFTReal <DT>::DataType * FFTReal <DT>::use_buffer () const
314{
315 return (&_buffer [0]);
316}
317
318
319
320/*\\\ PROTECTED \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
321
322
323
324/*\\\ PRIVATE \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
325
326
327
328template <class DT>
329void FFTReal <DT>::init_br_lut ()
330{
331 const long length = 1L << _nbr_bits;
332 _br_lut.resize (length);
333
334 _br_lut [0] = 0;
335 long br_index = 0;
336 for (long cnt = 1; cnt < length; ++cnt)
337 {
338 // ++br_index (bit reversed)
339 long bit = length >> 1;
340 while (((br_index ^= bit) & bit) == 0)
341 {
342 bit >>= 1;
343 }
344
345 _br_lut [cnt] = br_index;
346 }
347}
348
349
350
351template <class DT>
352void FFTReal <DT>::init_trigo_lut ()
353{
354 using namespace std;
355
356 if (_nbr_bits > 3)
357 {
358 const long total_len = (1L << (_nbr_bits - 1)) - 4;
359 _trigo_lut.resize (total_len);
360
361 for (int level = 3; level < _nbr_bits; ++level)
362 {
363 const long level_len = 1L << (level - 1);
364 DataType * const level_ptr =
365 &_trigo_lut [get_trigo_level_index (level)];
366 const double mul = PI / (level_len << 1);
367
368 for (long i = 0; i < level_len; ++ i)
369 {
370 level_ptr [i] = static_cast <DataType> (cos (i * mul));
371 }
372 }
373 }
374}
375
376
377
378template <class DT>
379void FFTReal <DT>::init_trigo_osc ()
380{
381 const int nbr_osc = _nbr_bits - TRIGO_BD_LIMIT;
382 if (nbr_osc > 0)
383 {
384 _trigo_osc.resize (nbr_osc);
385
386 for (int osc_cnt = 0; osc_cnt < nbr_osc; ++osc_cnt)
387 {
388 OscType & osc = _trigo_osc [osc_cnt];
389
390 const long len = 1L << (TRIGO_BD_LIMIT + osc_cnt);
391 const double mul = (0.5 * PI) / len;
392 osc.set_step (mul);
393 }
394 }
395}
396
397
398
399template <class DT>
400const long * FFTReal <DT>::get_br_ptr () const
401{
402 return (&_br_lut [0]);
403}
404
405
406
407template <class DT>
408const typename FFTReal <DT>::DataType * FFTReal <DT>::get_trigo_ptr (int level) const
409{
410 assert (level >= 3);
411
412 return (&_trigo_lut [get_trigo_level_index (level)]);
413}
414
415
416
417template <class DT>
418long FFTReal <DT>::get_trigo_level_index (int level) const
419{
420 assert (level >= 3);
421
422 return ((1L << (level - 1)) - 4);
423}
424
425
426
427// Transform in several passes
428template <class DT>
429void FFTReal <DT>::compute_fft_general (DataType f [], const DataType x []) const
430{
431 assert (f != 0);
432 assert (f != use_buffer ());
433 assert (x != 0);
434 assert (x != use_buffer ());
435 assert (x != f);
436
437 DataType * sf;
438 DataType * df;
439
440 if ((_nbr_bits & 1) != 0)
441 {
442 df = use_buffer ();
443 sf = f;
444 }
445 else
446 {
447 df = f;
448 sf = use_buffer ();
449 }
450
451 compute_direct_pass_1_2 (df, x);
452 compute_direct_pass_3 (sf, df);
453
454 for (int pass = 3; pass < _nbr_bits; ++ pass)
455 {
456 compute_direct_pass_n (df, sf, pass);
457
458 DataType * const temp_ptr = df;
459 df = sf;
460 sf = temp_ptr;
461 }
462}
463
464
465
466template <class DT>
467void FFTReal <DT>::compute_direct_pass_1_2 (DataType df [], const DataType x []) const
468{
469 assert (df != 0);
470 assert (x != 0);
471 assert (df != x);
472
473 const long * const bit_rev_lut_ptr = get_br_ptr ();
474 long coef_index = 0;
475 do
476 {
477 const long rev_index_0 = bit_rev_lut_ptr [coef_index];
478 const long rev_index_1 = bit_rev_lut_ptr [coef_index + 1];
479 const long rev_index_2 = bit_rev_lut_ptr [coef_index + 2];
480 const long rev_index_3 = bit_rev_lut_ptr [coef_index + 3];
481
482 DataType * const df2 = df + coef_index;
483 df2 [1] = x [rev_index_0] - x [rev_index_1];
484 df2 [3] = x [rev_index_2] - x [rev_index_3];
485
486 const DataType sf_0 = x [rev_index_0] + x [rev_index_1];
487 const DataType sf_2 = x [rev_index_2] + x [rev_index_3];
488
489 df2 [0] = sf_0 + sf_2;
490 df2 [2] = sf_0 - sf_2;
491
492 coef_index += 4;
493 }
494 while (coef_index < _length);
495}
496
497
498
499template <class DT>
500void FFTReal <DT>::compute_direct_pass_3 (DataType df [], const DataType sf []) const
501{
502 assert (df != 0);
503 assert (sf != 0);
504 assert (df != sf);
505
506 const DataType sqrt2_2 = DataType (SQRT2 * 0.5);
507 long coef_index = 0;
508 do
509 {
510 DataType v;
511
512 df [coef_index] = sf [coef_index] + sf [coef_index + 4];
513 df [coef_index + 4] = sf [coef_index] - sf [coef_index + 4];
514 df [coef_index + 2] = sf [coef_index + 2];
515 df [coef_index + 6] = sf [coef_index + 6];
516
517 v = (sf [coef_index + 5] - sf [coef_index + 7]) * sqrt2_2;
518 df [coef_index + 1] = sf [coef_index + 1] + v;
519 df [coef_index + 3] = sf [coef_index + 1] - v;
520
521 v = (sf [coef_index + 5] + sf [coef_index + 7]) * sqrt2_2;
522 df [coef_index + 5] = v + sf [coef_index + 3];
523 df [coef_index + 7] = v - sf [coef_index + 3];
524
525 coef_index += 8;
526 }
527 while (coef_index < _length);
528}
529
530
531
532template <class DT>
533void FFTReal <DT>::compute_direct_pass_n (DataType df [], const DataType sf [], int pass) const
534{
535 assert (df != 0);
536 assert (sf != 0);
537 assert (df != sf);
538 assert (pass >= 3);
539 assert (pass < _nbr_bits);
540
541 if (pass <= TRIGO_BD_LIMIT)
542 {
543 compute_direct_pass_n_lut (df, sf, pass);
544 }
545 else
546 {
547 compute_direct_pass_n_osc (df, sf, pass);
548 }
549}
550
551
552
553template <class DT>
554void FFTReal <DT>::compute_direct_pass_n_lut (DataType df [], const DataType sf [], int pass) const
555{
556 assert (df != 0);
557 assert (sf != 0);
558 assert (df != sf);
559 assert (pass >= 3);
560 assert (pass < _nbr_bits);
561
562 const long nbr_coef = 1 << pass;
563 const long h_nbr_coef = nbr_coef >> 1;
564 const long d_nbr_coef = nbr_coef << 1;
565 long coef_index = 0;
566 const DataType * const cos_ptr = get_trigo_ptr (pass);
567 do
568 {
569 const DataType * const sf1r = sf + coef_index;
570 const DataType * const sf2r = sf1r + nbr_coef;
571 DataType * const dfr = df + coef_index;
572 DataType * const dfi = dfr + nbr_coef;
573
574 // Extreme coefficients are always real
575 dfr [0] = sf1r [0] + sf2r [0];
576 dfi [0] = sf1r [0] - sf2r [0]; // dfr [nbr_coef] =
577 dfr [h_nbr_coef] = sf1r [h_nbr_coef];
578 dfi [h_nbr_coef] = sf2r [h_nbr_coef];
579
580 // Others are conjugate complex numbers
581 const DataType * const sf1i = sf1r + h_nbr_coef;
582 const DataType * const sf2i = sf1i + nbr_coef;
583 for (long i = 1; i < h_nbr_coef; ++ i)
584 {
585 const DataType c = cos_ptr [i]; // cos (i*PI/nbr_coef);
586 const DataType s = cos_ptr [h_nbr_coef - i]; // sin (i*PI/nbr_coef);
587 DataType v;
588
589 v = sf2r [i] * c - sf2i [i] * s;
590 dfr [i] = sf1r [i] + v;
591 dfi [-i] = sf1r [i] - v; // dfr [nbr_coef - i] =
592
593 v = sf2r [i] * s + sf2i [i] * c;
594 dfi [i] = v + sf1i [i];
595 dfi [nbr_coef - i] = v - sf1i [i];
596 }
597
598 coef_index += d_nbr_coef;
599 }
600 while (coef_index < _length);
601}
602
603
604
605template <class DT>
606void FFTReal <DT>::compute_direct_pass_n_osc (DataType df [], const DataType sf [], int pass) const
607{
608 assert (df != 0);
609 assert (sf != 0);
610 assert (df != sf);
611 assert (pass > TRIGO_BD_LIMIT);
612 assert (pass < _nbr_bits);
613
614 const long nbr_coef = 1 << pass;
615 const long h_nbr_coef = nbr_coef >> 1;
616 const long d_nbr_coef = nbr_coef << 1;
617 long coef_index = 0;
618 OscType & osc = _trigo_osc [pass - (TRIGO_BD_LIMIT + 1)];
619 do
620 {
621 const DataType * const sf1r = sf + coef_index;
622 const DataType * const sf2r = sf1r + nbr_coef;
623 DataType * const dfr = df + coef_index;
624 DataType * const dfi = dfr + nbr_coef;
625
626 osc.clear_buffers ();
627
628 // Extreme coefficients are always real
629 dfr [0] = sf1r [0] + sf2r [0];
630 dfi [0] = sf1r [0] - sf2r [0]; // dfr [nbr_coef] =
631 dfr [h_nbr_coef] = sf1r [h_nbr_coef];
632 dfi [h_nbr_coef] = sf2r [h_nbr_coef];
633
634 // Others are conjugate complex numbers
635 const DataType * const sf1i = sf1r + h_nbr_coef;
636 const DataType * const sf2i = sf1i + nbr_coef;
637 for (long i = 1; i < h_nbr_coef; ++ i)
638 {
639 osc.step ();
640 const DataType c = osc.get_cos ();
641 const DataType s = osc.get_sin ();
642 DataType v;
643
644 v = sf2r [i] * c - sf2i [i] * s;
645 dfr [i] = sf1r [i] + v;
646 dfi [-i] = sf1r [i] - v; // dfr [nbr_coef - i] =
647
648 v = sf2r [i] * s + sf2i [i] * c;
649 dfi [i] = v + sf1i [i];
650 dfi [nbr_coef - i] = v - sf1i [i];
651 }
652
653 coef_index += d_nbr_coef;
654 }
655 while (coef_index < _length);
656}
657
658
659
660// Transform in several pass
661template <class DT>
662void FFTReal <DT>::compute_ifft_general (const DataType f [], DataType x []) const
663{
664 assert (f != 0);
665 assert (f != use_buffer ());
666 assert (x != 0);
667 assert (x != use_buffer ());
668 assert (x != f);
669
670 DataType * sf = const_cast <DataType *> (f);
671 DataType * df;
672 DataType * df_temp;
673
674 if (_nbr_bits & 1)
675 {
676 df = use_buffer ();
677 df_temp = x;
678 }
679 else
680 {
681 df = x;
682 df_temp = use_buffer ();
683 }
684
685 for (int pass = _nbr_bits - 1; pass >= 3; -- pass)
686 {
687 compute_inverse_pass_n (df, sf, pass);
688
689 if (pass < _nbr_bits - 1)
690 {
691 DataType * const temp_ptr = df;
692 df = sf;
693 sf = temp_ptr;
694 }
695 else
696 {
697 sf = df;
698 df = df_temp;
699 }
700 }
701
702 compute_inverse_pass_3 (df, sf);
703 compute_inverse_pass_1_2 (x, df);
704}
705
706
707
708template <class DT>
709void FFTReal <DT>::compute_inverse_pass_n (DataType df [], const DataType sf [], int pass) const
710{
711 assert (df != 0);
712 assert (sf != 0);
713 assert (df != sf);
714 assert (pass >= 3);
715 assert (pass < _nbr_bits);
716
717 if (pass <= TRIGO_BD_LIMIT)
718 {
719 compute_inverse_pass_n_lut (df, sf, pass);
720 }
721 else
722 {
723 compute_inverse_pass_n_osc (df, sf, pass);
724 }
725}
726
727
728
729template <class DT>
730void FFTReal <DT>::compute_inverse_pass_n_lut (DataType df [], const DataType sf [], int pass) const
731{
732 assert (df != 0);
733 assert (sf != 0);
734 assert (df != sf);
735 assert (pass >= 3);
736 assert (pass < _nbr_bits);
737
738 const long nbr_coef = 1 << pass;
739 const long h_nbr_coef = nbr_coef >> 1;
740 const long d_nbr_coef = nbr_coef << 1;
741 long coef_index = 0;
742 const DataType * const cos_ptr = get_trigo_ptr (pass);
743 do
744 {
745 const DataType * const sfr = sf + coef_index;
746 const DataType * const sfi = sfr + nbr_coef;
747 DataType * const df1r = df + coef_index;
748 DataType * const df2r = df1r + nbr_coef;
749
750 // Extreme coefficients are always real
751 df1r [0] = sfr [0] + sfi [0]; // + sfr [nbr_coef]
752 df2r [0] = sfr [0] - sfi [0]; // - sfr [nbr_coef]
753 df1r [h_nbr_coef] = sfr [h_nbr_coef] * 2;
754 df2r [h_nbr_coef] = sfi [h_nbr_coef] * 2;
755
756 // Others are conjugate complex numbers
757 DataType * const df1i = df1r + h_nbr_coef;
758 DataType * const df2i = df1i + nbr_coef;
759 for (long i = 1; i < h_nbr_coef; ++ i)
760 {
761 df1r [i] = sfr [i] + sfi [-i]; // + sfr [nbr_coef - i]
762 df1i [i] = sfi [i] - sfi [nbr_coef - i];
763
764 const DataType c = cos_ptr [i]; // cos (i*PI/nbr_coef);
765 const DataType s = cos_ptr [h_nbr_coef - i]; // sin (i*PI/nbr_coef);
766 const DataType vr = sfr [i] - sfi [-i]; // - sfr [nbr_coef - i]
767 const DataType vi = sfi [i] + sfi [nbr_coef - i];
768
769 df2r [i] = vr * c + vi * s;
770 df2i [i] = vi * c - vr * s;
771 }
772
773 coef_index += d_nbr_coef;
774 }
775 while (coef_index < _length);
776}
777
778
779
780template <class DT>
781void FFTReal <DT>::compute_inverse_pass_n_osc (DataType df [], const DataType sf [], int pass) const
782{
783 assert (df != 0);
784 assert (sf != 0);
785 assert (df != sf);
786 assert (pass > TRIGO_BD_LIMIT);
787 assert (pass < _nbr_bits);
788
789 const long nbr_coef = 1 << pass;
790 const long h_nbr_coef = nbr_coef >> 1;
791 const long d_nbr_coef = nbr_coef << 1;
792 long coef_index = 0;
793 OscType & osc = _trigo_osc [pass - (TRIGO_BD_LIMIT + 1)];
794 do
795 {
796 const DataType * const sfr = sf + coef_index;
797 const DataType * const sfi = sfr + nbr_coef;
798 DataType * const df1r = df + coef_index;
799 DataType * const df2r = df1r + nbr_coef;
800
801 osc.clear_buffers ();
802
803 // Extreme coefficients are always real
804 df1r [0] = sfr [0] + sfi [0]; // + sfr [nbr_coef]
805 df2r [0] = sfr [0] - sfi [0]; // - sfr [nbr_coef]
806 df1r [h_nbr_coef] = sfr [h_nbr_coef] * 2;
807 df2r [h_nbr_coef] = sfi [h_nbr_coef] * 2;
808
809 // Others are conjugate complex numbers
810 DataType * const df1i = df1r + h_nbr_coef;
811 DataType * const df2i = df1i + nbr_coef;
812 for (long i = 1; i < h_nbr_coef; ++ i)
813 {
814 df1r [i] = sfr [i] + sfi [-i]; // + sfr [nbr_coef - i]
815 df1i [i] = sfi [i] - sfi [nbr_coef - i];
816
817 osc.step ();
818 const DataType c = osc.get_cos ();
819 const DataType s = osc.get_sin ();
820 const DataType vr = sfr [i] - sfi [-i]; // - sfr [nbr_coef - i]
821 const DataType vi = sfi [i] + sfi [nbr_coef - i];
822
823 df2r [i] = vr * c + vi * s;
824 df2i [i] = vi * c - vr * s;
825 }
826
827 coef_index += d_nbr_coef;
828 }
829 while (coef_index < _length);
830}
831
832
833
834template <class DT>
835void FFTReal <DT>::compute_inverse_pass_3 (DataType df [], const DataType sf []) const
836{
837 assert (df != 0);
838 assert (sf != 0);
839 assert (df != sf);
840
841 const DataType sqrt2_2 = DataType (SQRT2 * 0.5);
842 long coef_index = 0;
843 do
844 {
845 df [coef_index] = sf [coef_index] + sf [coef_index + 4];
846 df [coef_index + 4] = sf [coef_index] - sf [coef_index + 4];
847 df [coef_index + 2] = sf [coef_index + 2] * 2;
848 df [coef_index + 6] = sf [coef_index + 6] * 2;
849
850 df [coef_index + 1] = sf [coef_index + 1] + sf [coef_index + 3];
851 df [coef_index + 3] = sf [coef_index + 5] - sf [coef_index + 7];
852
853 const DataType vr = sf [coef_index + 1] - sf [coef_index + 3];
854 const DataType vi = sf [coef_index + 5] + sf [coef_index + 7];
855
856 df [coef_index + 5] = (vr + vi) * sqrt2_2;
857 df [coef_index + 7] = (vi - vr) * sqrt2_2;
858
859 coef_index += 8;
860 }
861 while (coef_index < _length);
862}
863
864
865
866template <class DT>
867void FFTReal <DT>::compute_inverse_pass_1_2 (DataType x [], const DataType sf []) const
868{
869 assert (x != 0);
870 assert (sf != 0);
871 assert (x != sf);
872
873 const long * bit_rev_lut_ptr = get_br_ptr ();
874 const DataType * sf2 = sf;
875 long coef_index = 0;
876 do
877 {
878 {
879 const DataType b_0 = sf2 [0] + sf2 [2];
880 const DataType b_2 = sf2 [0] - sf2 [2];
881 const DataType b_1 = sf2 [1] * 2;
882 const DataType b_3 = sf2 [3] * 2;
883
884 x [bit_rev_lut_ptr [0]] = b_0 + b_1;
885 x [bit_rev_lut_ptr [1]] = b_0 - b_1;
886 x [bit_rev_lut_ptr [2]] = b_2 + b_3;
887 x [bit_rev_lut_ptr [3]] = b_2 - b_3;
888 }
889 {
890 const DataType b_0 = sf2 [4] + sf2 [6];
891 const DataType b_2 = sf2 [4] - sf2 [6];
892 const DataType b_1 = sf2 [5] * 2;
893 const DataType b_3 = sf2 [7] * 2;
894
895 x [bit_rev_lut_ptr [4]] = b_0 + b_1;
896 x [bit_rev_lut_ptr [5]] = b_0 - b_1;
897 x [bit_rev_lut_ptr [6]] = b_2 + b_3;
898 x [bit_rev_lut_ptr [7]] = b_2 - b_3;
899 }
900
901 sf2 += 8;
902 coef_index += 8;
903 bit_rev_lut_ptr += 8;
904 }
905 while (coef_index < _length);
906}
907
908
909
910#endif // FFTReal_CODEHEADER_INCLUDED
911
912#undef FFTReal_CURRENT_CODEHEADER
913
914
915
916/*\\\ EOF \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
Note: See TracBrowser for help on using the repository browser.