arduino-audio-tools
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends Modules Pages
FFTReal.h
1/*****************************************************************************
2
3 FFTReal.h
4 By Laurent de Soras
5
6--- Legal stuff ---
7
8This program is free software. It comes without any warranty, to
9the extent permitted by applicable law. You can redistribute it
10and/or modify it under the terms of the Do What The Fuck You Want
11To Public License, Version 2, as published by Sam Hocevar. See
12http://sam.zoy.org/wtfpl/COPYING for more details.
13
14*Tab=3***********************************************************************/
15// #pragma GCC diagnostic ignored "-Wdeprecated-enum-enum-conversion" // commented out for platformio
16
17// We use the DefaultAllocator which supports PSRAM
18#include "AudioTools/CoreAudio/AudioBasic/Collections/Allocator.h"
19#define FFT_CUSTOM_ALLOC audio_tools::DefaultAllocator
20
21
22#if ! defined (ffft_FFTReal_HEADER_INCLUDED)
23#define ffft_FFTReal_HEADER_INCLUDED
24
25#if defined (_MSC_VER)
26 #pragma once
27 #pragma warning (4 : 4250) // "Inherits via dominance."
28#endif
29
30
31
32/*\\ INCLUDE FILES \\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
33
34/*****************************************************************************
35
36 def.h
37 By Laurent de Soras
38
39--- Legal stuff ---
40
41This program is free software. It comes without any warranty, to
42the extent permitted by applicable law. You can redistribute it
43and/or modify it under the terms of the Do What The Fuck You Want
44To Public License, Version 2, as published by Sam Hocevar. See
45http://sam.zoy.org/wtfpl/COPYING for more details.
46
47*Tab=3***********************************************************************/
48
49
50
51#if ! defined (ffft_def_HEADER_INCLUDED)
52#define ffft_def_HEADER_INCLUDED
53
54#if defined (_MSC_VER)
55 #pragma once
56 #pragma warning (4 : 4250) // "Inherits via dominance."
57#endif
58
59
60
61/*\\ INCLUDE FILES \\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
62
63
65namespace ffft {
66
67
68#ifndef PI
69const double PI = 3.1415926535897932384626433832795;
70#endif
71const double SQRT2 = 1.41421356237309514547462185873883;
72
73#if defined (_MSC_VER)
74
75 #define ffft_FORCEINLINE __forceinline
76
77#else
78
79 #define ffft_FORCEINLINE inline
80
81#endif
82
83
84
85} // namespace ffft
86
87
88
89#endif // ffft_def_HEADER_INCLUDED
90
91
92
93/*\\ EOF \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
94
95/*****************************************************************************
96
97 DynArray.h
98 By Laurent de Soras
99
100--- Legal stuff ---
101
102This program is free software. It comes without any warranty, to
103the extent permitted by applicable law. You can redistribute it
104and/or modify it under the terms of the Do What The Fuck You Want
105To Public License, Version 2, as published by Sam Hocevar. See
106http://sam.zoy.org/wtfpl/COPYING for more details.
107
108*Tab=3***********************************************************************/
109
110
111
112#if ! defined (ffft_DynArray_HEADER_INCLUDED)
113#define ffft_DynArray_HEADER_INCLUDED
114
115#if defined (_MSC_VER)
116 #pragma once
117 #pragma warning (4 : 4250) // "Inherits via dominance."
118#endif
119
120
121
122/*\\ INCLUDE FILES \\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
123
124
126namespace ffft
127{
128
129
130
131template <class T>
132class DynArray
133{
134
135/*\\ PUBLIC \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
136
137public:
138
139 typedef T DataType;
140
141 DynArray ();
142 explicit DynArray (long size);
143 ~DynArray ();
144
145 inline long size () const;
146 inline void resize (long size);
147
148 inline const DataType &
149 operator [] (long pos) const;
150 inline DataType &
151 operator [] (long pos);
152
153
154
155/*\\ PROTECTED \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
156
157protected:
158
159
160
161/*\\ PRIVATE \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
162
163private:
164
165 DataType * _data_ptr;
166 long _len;
167
168
169
170/*\\ FORBIDDEN MEMBER FUNCTIONS \\\\\\\\\\\\\\\\\\\\\\*/
171
172private:
173
174 DynArray (const DynArray &other);
175 DynArray & operator = (const DynArray &other);
176 bool operator == (const DynArray &other);
177 bool operator != (const DynArray &other);
178
179}; // class DynArray
180
181
182
183} // namespace ffft
184
185
186
187/*****************************************************************************
188
189 DynArray.hpp
190 By Laurent de Soras
191
192--- Legal stuff ---
193
194This program is free software. It comes without any warranty, to
195the extent permitted by applicable law. You can redistribute it
196and/or modify it under the terms of the Do What The Fuck You Want
197To Public License, Version 2, as published by Sam Hocevar. See
198http://sam.zoy.org/wtfpl/COPYING for more details.
199
200*Tab=3***********************************************************************/
201
202
203
204#if defined (ffft_DynArray_CURRENT_CODEHEADER)
205 #error Recursive inclusion of DynArray code header.
206#endif
207#define ffft_DynArray_CURRENT_CODEHEADER
208
209#if ! defined (ffft_DynArray_CODEHEADER_INCLUDED)
210#define ffft_DynArray_CODEHEADER_INCLUDED
211
212
213
214/*\\ INCLUDE FILES \\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
215
216#include <assert.h>
217
218
219
220namespace ffft
221{
222
223
224
225/*\\ PUBLIC \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
226
227
228
229template <class T>
230DynArray <T>::DynArray ()
231: _data_ptr (0)
232, _len (0)
233{
234 // Nothing
235}
236
237
238
239template <class T>
240DynArray <T>::DynArray (long size)
241: _data_ptr (0)
242, _len (0)
243{
244 assert (size >= 0);
245 if (size > 0)
246 {
247#ifdef FFT_CUSTOM_ALLOC
248 _data_ptr = FFT_CUSTOM_ALLOC.createArray<DataType>(size);
249#else
250 _data_ptr = new DataType [size];
251#endif
252 _len = size;
253 }
254}
255
256
257
258template <class T>
259DynArray <T>::~DynArray ()
260{
261#ifdef FFT_CUSTOM_ALLOC
262 FFT_CUSTOM_ALLOC.removeArray<DataType>(_data_ptr, _len);
263#else
264 delete [] _data_ptr;
265#endif
266 _data_ptr = 0;
267 _len = 0;
268}
269
270
271
272template <class T>
273long DynArray <T>::size () const
274{
275 return (_len);
276}
277
278
279
280template <class T>
281void DynArray <T>::resize (long size)
282{
283 assert (size >= 0);
284 if (size > 0)
285 {
286 DataType * old_data_ptr = _data_ptr;
287#ifdef FFT_CUSTOM_ALLOC
288 DataType * tmp_data_ptr = FFT_CUSTOM_ALLOC.createArray<DataType>(size);
289#else
290 DataType * tmp_data_ptr = new DataType [size];
291#endif
292 _data_ptr = tmp_data_ptr;
293 _len = size;
294
295#ifdef FFT_CUSTOM_ALLOC
296 FFT_CUSTOM_ALLOC.removeArray<DataType>(old_data_ptr, _len);
297#else
298 delete [] old_data_ptr;
299#endif
300 }
301}
302
303
304
305template <class T>
306const typename DynArray <T>::DataType & DynArray <T>::operator [] (long pos) const
307{
308 assert (pos >= 0);
309 assert (pos < _len);
310
311 return (_data_ptr [pos]);
312}
313
314
315
316template <class T>
317typename DynArray <T>::DataType & DynArray <T>::operator [] (long pos)
318{
319 assert (pos >= 0);
320 assert (pos < _len);
321
322 return (_data_ptr [pos]);
323}
324
325
326
327/*\\ PROTECTED \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
328
329
330
331/*\\ PRIVATE \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
332
333
334
335} // namespace ffft
336
337
338
339#endif // ffft_DynArray_CODEHEADER_INCLUDED
340
341#undef ffft_DynArray_CURRENT_CODEHEADER
342
343
344
345/*\\ EOF \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
346
347
348
349
350#endif // ffft_DynArray_HEADER_INCLUDED
351
352
353
354/*\\ EOF \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
355
356/*****************************************************************************
357
358 OscSinCos.h
359 By Laurent de Soras
360
361--- Legal stuff ---
362
363This program is free software. It comes without any warranty, to
364the extent permitted by applicable law. You can redistribute it
365and/or modify it under the terms of the Do What The Fuck You Want
366To Public License, Version 2, as published by Sam Hocevar. See
367http://sam.zoy.org/wtfpl/COPYING for more details.
368
369*Tab=3***********************************************************************/
370
371
372
373#if ! defined (ffft_OscSinCos_HEADER_INCLUDED)
374#define ffft_OscSinCos_HEADER_INCLUDED
375
376#if defined (_MSC_VER)
377 #pragma once
378 #pragma warning (4 : 4250) // "Inherits via dominance."
379#endif
380
381
382
383/*\\ INCLUDE FILES \\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
384
385
386
387
388
389namespace ffft
390{
391
392
393
394template <class T>
395class OscSinCos
396{
397
398/*\\ PUBLIC \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
399
400public:
401
402 typedef T DataType;
403
404 OscSinCos ();
405
406 ffft_FORCEINLINE void
407 set_step (double angle_rad);
408
409 ffft_FORCEINLINE DataType
410 get_cos () const;
411 ffft_FORCEINLINE DataType
412 get_sin () const;
413 ffft_FORCEINLINE void
414 step ();
415 ffft_FORCEINLINE void
416 clear_buffers ();
417
418
419
420/*\\ PROTECTED \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
421
422protected:
423
424
425
426/*\\ PRIVATE \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
427
428private:
429
430 DataType _pos_cos; // Current phase expressed with sin and cos. [-1 ; 1]
431 DataType _pos_sin; // -
432 DataType _step_cos; // Phase increment per step, [-1 ; 1]
433 DataType _step_sin; // -
434
435
436
437/*\\ FORBIDDEN MEMBER FUNCTIONS \\\\\\\\\\\\\\\\\\\\\\*/
438
439private:
440
441 OscSinCos (const OscSinCos &other);
442 OscSinCos & operator = (const OscSinCos &other);
443 bool operator == (const OscSinCos &other);
444 bool operator != (const OscSinCos &other);
445
446}; // class OscSinCos
447
448
449
450} // namespace ffft
451
452
453
454/*****************************************************************************
455
456 OscSinCos.hpp
457 By Laurent de Soras
458
459--- Legal stuff ---
460
461This program is free software. It comes without any warranty, to
462the extent permitted by applicable law. You can redistribute it
463and/or modify it under the terms of the Do What The Fuck You Want
464To Public License, Version 2, as published by Sam Hocevar. See
465http://sam.zoy.org/wtfpl/COPYING for more details.
466
467*Tab=3***********************************************************************/
468
469
470
471#if defined (ffft_OscSinCos_CURRENT_CODEHEADER)
472 #error Recursive inclusion of OscSinCos code header.
473#endif
474#define ffft_OscSinCos_CURRENT_CODEHEADER
475
476#if ! defined (ffft_OscSinCos_CODEHEADER_INCLUDED)
477#define ffft_OscSinCos_CODEHEADER_INCLUDED
478
479
480
481/*\\ INCLUDE FILES \\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
482#ifdef AVR
483# include <math.h>
484#else
485# include <cmath>
486#endif
487
488namespace std { }
489
490
491
492namespace ffft
493{
494
495
496
497/*\\ PUBLIC \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
498
499
500
501template <class T>
502OscSinCos <T>::OscSinCos ()
503: _pos_cos (1)
504, _pos_sin (0)
505, _step_cos (1)
506, _step_sin (0)
507{
508 // Nothing
509}
510
511
512
513template <class T>
514void OscSinCos <T>::set_step (double angle_rad)
515{
516 using namespace std;
517
518 _step_cos = static_cast <DataType> (cos (angle_rad));
519 _step_sin = static_cast <DataType> (sin (angle_rad));
520}
521
522
523
524template <class T>
525typename OscSinCos <T>::DataType OscSinCos <T>::get_cos () const
526{
527 return (_pos_cos);
528}
529
530
531
532template <class T>
533typename OscSinCos <T>::DataType OscSinCos <T>::get_sin () const
534{
535 return (_pos_sin);
536}
537
538
539
540template <class T>
541void OscSinCos <T>::step ()
542{
543 const DataType old_cos = _pos_cos;
544 const DataType old_sin = _pos_sin;
545
546 _pos_cos = old_cos * _step_cos - old_sin * _step_sin;
547 _pos_sin = old_cos * _step_sin + old_sin * _step_cos;
548}
549
550
551
552template <class T>
553void OscSinCos <T>::clear_buffers ()
554{
555 _pos_cos = static_cast <DataType> (1);
556 _pos_sin = static_cast <DataType> (0);
557}
558
559
560
561/*\\ PROTECTED \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
562
563
564
565/*\\ PRIVATE \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
566
567
568
569} // namespace ffft
570
571
572
573#endif // ffft_OscSinCos_CODEHEADER_INCLUDED
574
575#undef ffft_OscSinCos_CURRENT_CODEHEADER
576
577
578
579/*\\ EOF \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
580
581
582
583
584#endif // ffft_OscSinCos_HEADER_INCLUDED
585
586
587
588/*\\ EOF \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
589
590
591
592
593namespace ffft
594{
595
596
597
598template <class DT>
599class FFTReal
600{
601
602/*\\ PUBLIC \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
603
604public:
605
606 enum { MAX_BIT_DEPTH = 30 }; // So length can be represented as long int
607
608 typedef DT DataType;
609
610 explicit FFTReal (long length);
611 virtual ~FFTReal () {}
612
613 long get_length () const;
614 void do_fft (DataType f [], const DataType x []) const;
615 void do_ifft (const DataType f [], DataType x []) const;
616 void rescale (DataType x []) const;
617 DataType * use_buffer () const;
618
619
620
621/*\\ PROTECTED \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
622
623protected:
624
625
626
627/*\\ PRIVATE \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
628
629private:
630
631 // Over this bit depth, we use direct calculation for sin/cos
632 enum { TRIGO_BD_LIMIT = 12 };
633
634 typedef OscSinCos <DataType> OscType;
635
636 void init_br_lut ();
637 void init_trigo_lut ();
638 void init_trigo_osc ();
639
640 ffft_FORCEINLINE const long *
641 get_br_ptr () const;
642 ffft_FORCEINLINE const DataType *
643 get_trigo_ptr (int level) const;
644 ffft_FORCEINLINE long
645 get_trigo_level_index (int level) const;
646
647 inline void compute_fft_general (DataType f [], const DataType x []) const;
648 inline void compute_direct_pass_1_2 (DataType df [], const DataType x []) const;
649 inline void compute_direct_pass_3 (DataType df [], const DataType sf []) const;
650 inline void compute_direct_pass_n (DataType df [], const DataType sf [], int pass) const;
651 inline void compute_direct_pass_n_lut (DataType df [], const DataType sf [], int pass) const;
652 inline void compute_direct_pass_n_osc (DataType df [], const DataType sf [], int pass) const;
653
654 inline void compute_ifft_general (const DataType f [], DataType x []) const;
655 inline void compute_inverse_pass_n (DataType df [], const DataType sf [], int pass) const;
656 inline void compute_inverse_pass_n_osc (DataType df [], const DataType sf [], int pass) const;
657 inline void compute_inverse_pass_n_lut (DataType df [], const DataType sf [], int pass) const;
658 inline void compute_inverse_pass_3 (DataType df [], const DataType sf []) const;
659 inline void compute_inverse_pass_1_2 (DataType x [], const DataType sf []) const;
660
661 const long _length;
662 const int _nbr_bits;
663 DynArray <long>
664 _br_lut;
665 DynArray <DataType>
666 _trigo_lut;
667 mutable DynArray <DataType>
668 _buffer;
669 mutable DynArray <OscType>
670 _trigo_osc;
671
672
673
674/*\\ FORBIDDEN MEMBER FUNCTIONS \\\\\\\\\\\\\\\\\\\\\\*/
675
676private:
677
678 FFTReal ();
679 FFTReal (const FFTReal &other);
680 FFTReal & operator = (const FFTReal &other);
681 bool operator == (const FFTReal &other);
682 bool operator != (const FFTReal &other);
683
684}; // class FFTReal
685
686
687
688} // namespace ffft
689
690
691
692/*****************************************************************************
693
694 FFTReal.hpp
695 By Laurent de Soras
696
697--- Legal stuff ---
698
699This program is free software. It comes without any warranty, to
700the extent permitted by applicable law. You can redistribute it
701and/or modify it under the terms of the Do What The Fuck You Want
702To Public License, Version 2, as published by Sam Hocevar. See
703http://sam.zoy.org/wtfpl/COPYING for more details.
704
705*Tab=3***********************************************************************/
706
707
708
709#if defined (ffft_FFTReal_CURRENT_CODEHEADER)
710 #error Recursive inclusion of FFTReal code header.
711#endif
712#define ffft_FFTReal_CURRENT_CODEHEADER
713
714#if ! defined (ffft_FFTReal_CODEHEADER_INCLUDED)
715#define ffft_FFTReal_CODEHEADER_INCLUDED
716
717
718
719namespace ffft
720{
721
722
723
724static inline bool FFTReal_is_pow2 (long x)
725{
726 assert (x > 0);
727
728 return ((x & -x) == x);
729}
730
731
732
733static inline int FFTReal_get_next_pow2 (long x)
734{
735 --x;
736
737 int p = 0;
738 while ((x & ~0xFFFFL) != 0)
739 {
740 p += 16;
741 x >>= 16;
742 }
743 while ((x & ~0xFL) != 0)
744 {
745 p += 4;
746 x >>= 4;
747 }
748 while (x > 0)
749 {
750 ++p;
751 x >>= 1;
752 }
753
754 return (p);
755}
756
757
758
759/*\\ PUBLIC \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
760
761
762
763/*
764==============================================================================
765Name: ctor
766Input parameters:
767 - length: length of the array on which we want to do a FFT. Range: power of
768 2 only, > 0.
769Throws: std::bad_alloc
770==============================================================================
771*/
772
773template <class DT>
774FFTReal <DT>::FFTReal (long length)
775: _length (length)
776, _nbr_bits (FFTReal_get_next_pow2 (length))
777, _br_lut ()
778, _trigo_lut ()
779, _buffer (length)
780, _trigo_osc ()
781{
782 assert (FFTReal_is_pow2 (length));
783 assert (_nbr_bits <= MAX_BIT_DEPTH);
784
785 init_br_lut ();
786 init_trigo_lut ();
787 init_trigo_osc ();
788}
789
790
791
792/*
793==============================================================================
794Name: get_length
795Description:
796 Returns the number of points processed by this FFT object.
797Returns: The number of points, power of 2, > 0.
798Throws: Nothing
799==============================================================================
800*/
801
802template <class DT>
803long FFTReal <DT>::get_length () const
804{
805 return (_length);
806}
807
808
809
810/*
811==============================================================================
812Name: do_fft
813Description:
814 Compute the FFT of the array.
815Input parameters:
816 - x: pointer on the source array (time).
817Output parameters:
818 - f: pointer on the destination array (frequencies).
819 f [0...length(x)/2] = real values,
820 f [length(x)/2+1...length(x)-1] = negative imaginary values of
821 coefficents 1...length(x)/2-1.
822Throws: Nothing
823==============================================================================
824*/
825
826template <class DT>
827void FFTReal <DT>::do_fft (DataType f [], const DataType x []) const
828{
829 assert (f != 0);
830 assert (f != use_buffer ());
831 assert (x != 0);
832 assert (x != use_buffer ());
833 assert (x != f);
834
835 // General case
836 if (_nbr_bits > 2)
837 {
838 compute_fft_general (f, x);
839 }
840
841 // 4-point FFT
842 else if (_nbr_bits == 2)
843 {
844 f [1] = x [0] - x [2];
845 f [3] = x [1] - x [3];
846
847 const DataType b_0 = x [0] + x [2];
848 const DataType b_2 = x [1] + x [3];
849
850 f [0] = b_0 + b_2;
851 f [2] = b_0 - b_2;
852 }
853
854 // 2-point FFT
855 else if (_nbr_bits == 1)
856 {
857 f [0] = x [0] + x [1];
858 f [1] = x [0] - x [1];
859 }
860
861 // 1-point FFT
862 else
863 {
864 f [0] = x [0];
865 }
866}
867
868
869
870/*
871==============================================================================
872Name: do_ifft
873Description:
874 Compute the inverse FFT of the array. Note that data must be post-scaled:
875 IFFT (FFT (x)) = x * length (x).
876Input parameters:
877 - f: pointer on the source array (frequencies).
878 f [0...length(x)/2] = real values
879 f [length(x)/2+1...length(x)-1] = negative imaginary values of
880 coefficents 1...length(x)/2-1.
881Output parameters:
882 - x: pointer on the destination array (time).
883Throws: Nothing
884==============================================================================
885*/
886
887template <class DT>
888void FFTReal <DT>::do_ifft (const DataType f [], DataType x []) const
889{
890 assert (f != 0);
891 assert (f != use_buffer ());
892 assert (x != 0);
893 assert (x != use_buffer ());
894 assert (x != f);
895
896 // General case
897 if (_nbr_bits > 2)
898 {
899 compute_ifft_general (f, x);
900 }
901
902 // 4-point IFFT
903 else if (_nbr_bits == 2)
904 {
905 const DataType b_0 = f [0] + f [2];
906 const DataType b_2 = f [0] - f [2];
907
908 x [0] = b_0 + f [1] * 2;
909 x [2] = b_0 - f [1] * 2;
910 x [1] = b_2 + f [3] * 2;
911 x [3] = b_2 - f [3] * 2;
912 }
913
914 // 2-point IFFT
915 else if (_nbr_bits == 1)
916 {
917 x [0] = f [0] + f [1];
918 x [1] = f [0] - f [1];
919 }
920
921 // 1-point IFFT
922 else
923 {
924 x [0] = f [0];
925 }
926}
927
928
929
930/*
931==============================================================================
932Name: rescale
933Description:
934 Scale an array by divide each element by its length. This function should
935 be called after FFT + IFFT.
936Input parameters:
937 - x: pointer on array to rescale (time or frequency).
938Throws: Nothing
939==============================================================================
940*/
941
942template <class DT>
943void FFTReal <DT>::rescale (DataType x []) const
944{
945 const DataType mul = DataType (1.0 / _length);
946
947 if (_length < 4)
948 {
949 long i = _length - 1;
950 do
951 {
952 x [i] *= mul;
953 --i;
954 }
955 while (i >= 0);
956 }
957
958 else
959 {
960 assert ((_length & 3) == 0);
961
962 // Could be optimized with SIMD instruction sets (needs alignment check)
963 long i = _length - 4;
964 do
965 {
966 x [i + 0] *= mul;
967 x [i + 1] *= mul;
968 x [i + 2] *= mul;
969 x [i + 3] *= mul;
970 i -= 4;
971 }
972 while (i >= 0);
973 }
974}
975
976
977
978/*
979==============================================================================
980Name: use_buffer
981Description:
982 Access the internal buffer, whose length is the FFT one.
983 Buffer content will be erased at each do_fft() / do_ifft() call!
984 This buffer cannot be used as:
985 - source for FFT or IFFT done with this object
986 - destination for FFT or IFFT done with this object
987Returns:
988 Buffer start address
989Throws: Nothing
990==============================================================================
991*/
992
993template <class DT>
994typename FFTReal <DT>::DataType * FFTReal <DT>::use_buffer () const
995{
996 return (&_buffer [0]);
997}
998
999
1000
1001/*\\ PROTECTED \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
1002
1003
1004
1005/*\\ PRIVATE \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
1006
1007
1008
1009template <class DT>
1010void FFTReal <DT>::init_br_lut ()
1011{
1012 const long length = 1L << _nbr_bits;
1013 _br_lut.resize (length);
1014
1015 _br_lut [0] = 0;
1016 long br_index = 0;
1017 for (long cnt = 1; cnt < length; ++cnt)
1018 {
1019 // ++br_index (bit reversed)
1020 long bit = length >> 1;
1021 while (((br_index ^= bit) & bit) == 0)
1022 {
1023 bit >>= 1;
1024 }
1025
1026 _br_lut [cnt] = br_index;
1027 }
1028}
1029
1030
1031
1032template <class DT>
1033void FFTReal <DT>::init_trigo_lut ()
1034{
1035 using namespace std;
1036
1037 if (_nbr_bits > 3)
1038 {
1039 const long total_len = (1L << (_nbr_bits - 1)) - 4;
1040 _trigo_lut.resize (total_len);
1041
1042 for (int level = 3; level < _nbr_bits; ++level)
1043 {
1044 const long level_len = 1L << (level - 1);
1045 DataType * const level_ptr =
1046 &_trigo_lut [get_trigo_level_index (level)];
1047 const double mul = PI / (level_len << 1);
1048
1049 for (long i = 0; i < level_len; ++ i)
1050 {
1051 level_ptr [i] = static_cast <DataType> (cos (i * mul));
1052 }
1053 }
1054 }
1055}
1056
1057
1058
1059template <class DT>
1060void FFTReal <DT>::init_trigo_osc ()
1061{
1062 const int nbr_osc = _nbr_bits - TRIGO_BD_LIMIT;
1063 if (nbr_osc > 0)
1064 {
1065 _trigo_osc.resize (nbr_osc);
1066
1067 for (int osc_cnt = 0; osc_cnt < nbr_osc; ++osc_cnt)
1068 {
1069 OscType & osc = _trigo_osc [osc_cnt];
1070
1071 const long len = 1L << (TRIGO_BD_LIMIT + osc_cnt);
1072 const double mul = (0.5 * PI) / len;
1073 osc.set_step (mul);
1074 }
1075 }
1076}
1077
1078
1079
1080template <class DT>
1081const long * FFTReal <DT>::get_br_ptr () const
1082{
1083 return (&_br_lut [0]);
1084}
1085
1086
1087
1088template <class DT>
1089const typename FFTReal <DT>::DataType * FFTReal <DT>::get_trigo_ptr (int level) const
1090{
1091 assert (level >= 3);
1092
1093 return (&_trigo_lut [get_trigo_level_index (level)]);
1094}
1095
1096
1097
1098template <class DT>
1099long FFTReal <DT>::get_trigo_level_index (int level) const
1100{
1101 assert (level >= 3);
1102
1103 return ((1L << (level - 1)) - 4);
1104}
1105
1106
1107
1108// Transform in several passes
1109template <class DT>
1110void FFTReal <DT>::compute_fft_general (DataType f [], const DataType x []) const
1111{
1112 assert (f != 0);
1113 assert (f != use_buffer ());
1114 assert (x != 0);
1115 assert (x != use_buffer ());
1116 assert (x != f);
1117
1118 DataType * sf;
1119 DataType * df;
1120
1121 if ((_nbr_bits & 1) != 0)
1122 {
1123 df = use_buffer ();
1124 sf = f;
1125 }
1126 else
1127 {
1128 df = f;
1129 sf = use_buffer ();
1130 }
1131
1132 compute_direct_pass_1_2 (df, x);
1133 compute_direct_pass_3 (sf, df);
1134
1135 for (int pass = 3; pass < _nbr_bits; ++ pass)
1136 {
1137 compute_direct_pass_n (df, sf, pass);
1138
1139 DataType * const temp_ptr = df;
1140 df = sf;
1141 sf = temp_ptr;
1142 }
1143}
1144
1145
1146
1147template <class DT>
1148void FFTReal <DT>::compute_direct_pass_1_2 (DataType df [], const DataType x []) const
1149{
1150 assert (df != 0);
1151 assert (x != 0);
1152 assert (df != x);
1153
1154 const long * const bit_rev_lut_ptr = get_br_ptr ();
1155 long coef_index = 0;
1156 do
1157 {
1158 const long rev_index_0 = bit_rev_lut_ptr [coef_index];
1159 const long rev_index_1 = bit_rev_lut_ptr [coef_index + 1];
1160 const long rev_index_2 = bit_rev_lut_ptr [coef_index + 2];
1161 const long rev_index_3 = bit_rev_lut_ptr [coef_index + 3];
1162
1163 DataType * const df2 = df + coef_index;
1164 df2 [1] = x [rev_index_0] - x [rev_index_1];
1165 df2 [3] = x [rev_index_2] - x [rev_index_3];
1166
1167 const DataType sf_0 = x [rev_index_0] + x [rev_index_1];
1168 const DataType sf_2 = x [rev_index_2] + x [rev_index_3];
1169
1170 df2 [0] = sf_0 + sf_2;
1171 df2 [2] = sf_0 - sf_2;
1172
1173 coef_index += 4;
1174 }
1175 while (coef_index < _length);
1176}
1177
1178
1179
1180template <class DT>
1181void FFTReal <DT>::compute_direct_pass_3 (DataType df [], const DataType sf []) const
1182{
1183 assert (df != 0);
1184 assert (sf != 0);
1185 assert (df != sf);
1186
1187 const DataType sqrt2_2 = DataType (SQRT2 * 0.5);
1188 long coef_index = 0;
1189 do
1190 {
1191 DataType v;
1192
1193 df [coef_index] = sf [coef_index] + sf [coef_index + 4];
1194 df [coef_index + 4] = sf [coef_index] - sf [coef_index + 4];
1195 df [coef_index + 2] = sf [coef_index + 2];
1196 df [coef_index + 6] = sf [coef_index + 6];
1197
1198 v = (sf [coef_index + 5] - sf [coef_index + 7]) * sqrt2_2;
1199 df [coef_index + 1] = sf [coef_index + 1] + v;
1200 df [coef_index + 3] = sf [coef_index + 1] - v;
1201
1202 v = (sf [coef_index + 5] + sf [coef_index + 7]) * sqrt2_2;
1203 df [coef_index + 5] = v + sf [coef_index + 3];
1204 df [coef_index + 7] = v - sf [coef_index + 3];
1205
1206 coef_index += 8;
1207 }
1208 while (coef_index < _length);
1209}
1210
1211
1212
1213template <class DT>
1214void FFTReal <DT>::compute_direct_pass_n (DataType df [], const DataType sf [], int pass) const
1215{
1216 assert (df != 0);
1217 assert (sf != 0);
1218 assert (df != sf);
1219 assert (pass >= 3);
1220 assert (pass < _nbr_bits);
1221
1222 if (pass <= TRIGO_BD_LIMIT)
1223 {
1224 compute_direct_pass_n_lut (df, sf, pass);
1225 }
1226 else
1227 {
1228 compute_direct_pass_n_osc (df, sf, pass);
1229 }
1230}
1231
1232
1233
1234template <class DT>
1235void FFTReal <DT>::compute_direct_pass_n_lut (DataType df [], const DataType sf [], int pass) const
1236{
1237 assert (df != 0);
1238 assert (sf != 0);
1239 assert (df != sf);
1240 assert (pass >= 3);
1241 assert (pass < _nbr_bits);
1242
1243 const long nbr_coef = 1 << pass;
1244 const long h_nbr_coef = nbr_coef >> 1;
1245 const long d_nbr_coef = nbr_coef << 1;
1246 long coef_index = 0;
1247 const DataType * const cos_ptr = get_trigo_ptr (pass);
1248 do
1249 {
1250 const DataType * const sf1r = sf + coef_index;
1251 const DataType * const sf2r = sf1r + nbr_coef;
1252 DataType * const dfr = df + coef_index;
1253 DataType * const dfi = dfr + nbr_coef;
1254
1255 // Extreme coefficients are always real
1256 dfr [0] = sf1r [0] + sf2r [0];
1257 dfi [0] = sf1r [0] - sf2r [0]; // dfr [nbr_coef] =
1258 dfr [h_nbr_coef] = sf1r [h_nbr_coef];
1259 dfi [h_nbr_coef] = sf2r [h_nbr_coef];
1260
1261 // Others are conjugate complex numbers
1262 const DataType * const sf1i = sf1r + h_nbr_coef;
1263 const DataType * const sf2i = sf1i + nbr_coef;
1264 for (long i = 1; i < h_nbr_coef; ++ i)
1265 {
1266 const DataType c = cos_ptr [i]; // cos (i*PI/nbr_coef);
1267 const DataType s = cos_ptr [h_nbr_coef - i]; // sin (i*PI/nbr_coef);
1268 DataType v;
1269
1270 v = sf2r [i] * c - sf2i [i] * s;
1271 dfr [i] = sf1r [i] + v;
1272 dfi [-i] = sf1r [i] - v; // dfr [nbr_coef - i] =
1273
1274 v = sf2r [i] * s + sf2i [i] * c;
1275 dfi [i] = v + sf1i [i];
1276 dfi [nbr_coef - i] = v - sf1i [i];
1277 }
1278
1279 coef_index += d_nbr_coef;
1280 }
1281 while (coef_index < _length);
1282}
1283
1284
1285
1286template <class DT>
1287void FFTReal <DT>::compute_direct_pass_n_osc (DataType df [], const DataType sf [], int pass) const
1288{
1289 assert (df != 0);
1290 assert (sf != 0);
1291 assert (df != sf);
1292 assert (pass > TRIGO_BD_LIMIT);
1293 assert (pass < _nbr_bits);
1294
1295 const long nbr_coef = 1 << pass;
1296 const long h_nbr_coef = nbr_coef >> 1;
1297 const long d_nbr_coef = nbr_coef << 1;
1298 long coef_index = 0;
1299 OscType & osc = _trigo_osc [pass - (TRIGO_BD_LIMIT + 1)];
1300 do
1301 {
1302 const DataType * const sf1r = sf + coef_index;
1303 const DataType * const sf2r = sf1r + nbr_coef;
1304 DataType * const dfr = df + coef_index;
1305 DataType * const dfi = dfr + nbr_coef;
1306
1307 osc.clear_buffers ();
1308
1309 // Extreme coefficients are always real
1310 dfr [0] = sf1r [0] + sf2r [0];
1311 dfi [0] = sf1r [0] - sf2r [0]; // dfr [nbr_coef] =
1312 dfr [h_nbr_coef] = sf1r [h_nbr_coef];
1313 dfi [h_nbr_coef] = sf2r [h_nbr_coef];
1314
1315 // Others are conjugate complex numbers
1316 const DataType * const sf1i = sf1r + h_nbr_coef;
1317 const DataType * const sf2i = sf1i + nbr_coef;
1318 for (long i = 1; i < h_nbr_coef; ++ i)
1319 {
1320 osc.step ();
1321 const DataType c = osc.get_cos ();
1322 const DataType s = osc.get_sin ();
1323 DataType v;
1324
1325 v = sf2r [i] * c - sf2i [i] * s;
1326 dfr [i] = sf1r [i] + v;
1327 dfi [-i] = sf1r [i] - v; // dfr [nbr_coef - i] =
1328
1329 v = sf2r [i] * s + sf2i [i] * c;
1330 dfi [i] = v + sf1i [i];
1331 dfi [nbr_coef - i] = v - sf1i [i];
1332 }
1333
1334 coef_index += d_nbr_coef;
1335 }
1336 while (coef_index < _length);
1337}
1338
1339
1340
1341// Transform in several pass
1342template <class DT>
1343void FFTReal <DT>::compute_ifft_general (const DataType f [], DataType x []) const
1344{
1345 assert (f != 0);
1346 assert (f != use_buffer ());
1347 assert (x != 0);
1348 assert (x != use_buffer ());
1349 assert (x != f);
1350
1351 DataType * sf = const_cast <DataType *> (f);
1352 DataType * df;
1353 DataType * df_temp;
1354
1355 if (_nbr_bits & 1)
1356 {
1357 df = use_buffer ();
1358 df_temp = x;
1359 }
1360 else
1361 {
1362 df = x;
1363 df_temp = use_buffer ();
1364 }
1365
1366 for (int pass = _nbr_bits - 1; pass >= 3; -- pass)
1367 {
1368 compute_inverse_pass_n (df, sf, pass);
1369
1370 if (pass < _nbr_bits - 1)
1371 {
1372 DataType * const temp_ptr = df;
1373 df = sf;
1374 sf = temp_ptr;
1375 }
1376 else
1377 {
1378 sf = df;
1379 df = df_temp;
1380 }
1381 }
1382
1383 compute_inverse_pass_3 (df, sf);
1384 compute_inverse_pass_1_2 (x, df);
1385}
1386
1387
1388
1389template <class DT>
1390void FFTReal <DT>::compute_inverse_pass_n (DataType df [], const DataType sf [], int pass) const
1391{
1392 assert (df != 0);
1393 assert (sf != 0);
1394 assert (df != sf);
1395 assert (pass >= 3);
1396 assert (pass < _nbr_bits);
1397
1398 if (pass <= TRIGO_BD_LIMIT)
1399 {
1400 compute_inverse_pass_n_lut (df, sf, pass);
1401 }
1402 else
1403 {
1404 compute_inverse_pass_n_osc (df, sf, pass);
1405 }
1406}
1407
1408
1409
1410template <class DT>
1411void FFTReal <DT>::compute_inverse_pass_n_lut (DataType df [], const DataType sf [], int pass) const
1412{
1413 assert (df != 0);
1414 assert (sf != 0);
1415 assert (df != sf);
1416 assert (pass >= 3);
1417 assert (pass < _nbr_bits);
1418
1419 const long nbr_coef = 1 << pass;
1420 const long h_nbr_coef = nbr_coef >> 1;
1421 const long d_nbr_coef = nbr_coef << 1;
1422 long coef_index = 0;
1423 const DataType * const cos_ptr = get_trigo_ptr (pass);
1424 do
1425 {
1426 const DataType * const sfr = sf + coef_index;
1427 const DataType * const sfi = sfr + nbr_coef;
1428 DataType * const df1r = df + coef_index;
1429 DataType * const df2r = df1r + nbr_coef;
1430
1431 // Extreme coefficients are always real
1432 df1r [0] = sfr [0] + sfi [0]; // + sfr [nbr_coef]
1433 df2r [0] = sfr [0] - sfi [0]; // - sfr [nbr_coef]
1434 df1r [h_nbr_coef] = sfr [h_nbr_coef] * 2;
1435 df2r [h_nbr_coef] = sfi [h_nbr_coef] * 2;
1436
1437 // Others are conjugate complex numbers
1438 DataType * const df1i = df1r + h_nbr_coef;
1439 DataType * const df2i = df1i + nbr_coef;
1440 for (long i = 1; i < h_nbr_coef; ++ i)
1441 {
1442 df1r [i] = sfr [i] + sfi [-i]; // + sfr [nbr_coef - i]
1443 df1i [i] = sfi [i] - sfi [nbr_coef - i];
1444
1445 const DataType c = cos_ptr [i]; // cos (i*PI/nbr_coef);
1446 const DataType s = cos_ptr [h_nbr_coef - i]; // sin (i*PI/nbr_coef);
1447 const DataType vr = sfr [i] - sfi [-i]; // - sfr [nbr_coef - i]
1448 const DataType vi = sfi [i] + sfi [nbr_coef - i];
1449
1450 df2r [i] = vr * c + vi * s;
1451 df2i [i] = vi * c - vr * s;
1452 }
1453
1454 coef_index += d_nbr_coef;
1455 }
1456 while (coef_index < _length);
1457}
1458
1459
1460
1461template <class DT>
1462void FFTReal <DT>::compute_inverse_pass_n_osc (DataType df [], const DataType sf [], int pass) const
1463{
1464 assert (df != 0);
1465 assert (sf != 0);
1466 assert (df != sf);
1467 assert (pass > TRIGO_BD_LIMIT);
1468 assert (pass < _nbr_bits);
1469
1470 const long nbr_coef = 1 << pass;
1471 const long h_nbr_coef = nbr_coef >> 1;
1472 const long d_nbr_coef = nbr_coef << 1;
1473 long coef_index = 0;
1474 OscType & osc = _trigo_osc [pass - (TRIGO_BD_LIMIT + 1)];
1475 do
1476 {
1477 const DataType * const sfr = sf + coef_index;
1478 const DataType * const sfi = sfr + nbr_coef;
1479 DataType * const df1r = df + coef_index;
1480 DataType * const df2r = df1r + nbr_coef;
1481
1482 osc.clear_buffers ();
1483
1484 // Extreme coefficients are always real
1485 df1r [0] = sfr [0] + sfi [0]; // + sfr [nbr_coef]
1486 df2r [0] = sfr [0] - sfi [0]; // - sfr [nbr_coef]
1487 df1r [h_nbr_coef] = sfr [h_nbr_coef] * 2;
1488 df2r [h_nbr_coef] = sfi [h_nbr_coef] * 2;
1489
1490 // Others are conjugate complex numbers
1491 DataType * const df1i = df1r + h_nbr_coef;
1492 DataType * const df2i = df1i + nbr_coef;
1493 for (long i = 1; i < h_nbr_coef; ++ i)
1494 {
1495 df1r [i] = sfr [i] + sfi [-i]; // + sfr [nbr_coef - i]
1496 df1i [i] = sfi [i] - sfi [nbr_coef - i];
1497
1498 osc.step ();
1499 const DataType c = osc.get_cos ();
1500 const DataType s = osc.get_sin ();
1501 const DataType vr = sfr [i] - sfi [-i]; // - sfr [nbr_coef - i]
1502 const DataType vi = sfi [i] + sfi [nbr_coef - i];
1503
1504 df2r [i] = vr * c + vi * s;
1505 df2i [i] = vi * c - vr * s;
1506 }
1507
1508 coef_index += d_nbr_coef;
1509 }
1510 while (coef_index < _length);
1511}
1512
1513
1514
1515template <class DT>
1516void FFTReal <DT>::compute_inverse_pass_3 (DataType df [], const DataType sf []) const
1517{
1518 assert (df != 0);
1519 assert (sf != 0);
1520 assert (df != sf);
1521
1522 const DataType sqrt2_2 = DataType (SQRT2 * 0.5);
1523 long coef_index = 0;
1524 do
1525 {
1526 df [coef_index] = sf [coef_index] + sf [coef_index + 4];
1527 df [coef_index + 4] = sf [coef_index] - sf [coef_index + 4];
1528 df [coef_index + 2] = sf [coef_index + 2] * 2;
1529 df [coef_index + 6] = sf [coef_index + 6] * 2;
1530
1531 df [coef_index + 1] = sf [coef_index + 1] + sf [coef_index + 3];
1532 df [coef_index + 3] = sf [coef_index + 5] - sf [coef_index + 7];
1533
1534 const DataType vr = sf [coef_index + 1] - sf [coef_index + 3];
1535 const DataType vi = sf [coef_index + 5] + sf [coef_index + 7];
1536
1537 df [coef_index + 5] = (vr + vi) * sqrt2_2;
1538 df [coef_index + 7] = (vi - vr) * sqrt2_2;
1539
1540 coef_index += 8;
1541 }
1542 while (coef_index < _length);
1543}
1544
1545
1546
1547template <class DT>
1548void FFTReal <DT>::compute_inverse_pass_1_2 (DataType x [], const DataType sf []) const
1549{
1550 assert (x != 0);
1551 assert (sf != 0);
1552 assert (x != sf);
1553
1554 const long * bit_rev_lut_ptr = get_br_ptr ();
1555 const DataType * sf2 = sf;
1556 long coef_index = 0;
1557 do
1558 {
1559 {
1560 const DataType b_0 = sf2 [0] + sf2 [2];
1561 const DataType b_2 = sf2 [0] - sf2 [2];
1562 const DataType b_1 = sf2 [1] * 2;
1563 const DataType b_3 = sf2 [3] * 2;
1564
1565 x [bit_rev_lut_ptr [0]] = b_0 + b_1;
1566 x [bit_rev_lut_ptr [1]] = b_0 - b_1;
1567 x [bit_rev_lut_ptr [2]] = b_2 + b_3;
1568 x [bit_rev_lut_ptr [3]] = b_2 - b_3;
1569 }
1570 {
1571 const DataType b_0 = sf2 [4] + sf2 [6];
1572 const DataType b_2 = sf2 [4] - sf2 [6];
1573 const DataType b_1 = sf2 [5] * 2;
1574 const DataType b_3 = sf2 [7] * 2;
1575
1576 x [bit_rev_lut_ptr [4]] = b_0 + b_1;
1577 x [bit_rev_lut_ptr [5]] = b_0 - b_1;
1578 x [bit_rev_lut_ptr [6]] = b_2 + b_3;
1579 x [bit_rev_lut_ptr [7]] = b_2 - b_3;
1580 }
1581
1582 sf2 += 8;
1583 coef_index += 8;
1584 bit_rev_lut_ptr += 8;
1585 }
1586 while (coef_index < _length);
1587}
1588
1589
1590
1591} // namespace ffft
1592
1593
1594
1595#endif // ffft_FFTReal_CODEHEADER_INCLUDED
1596
1597#undef ffft_FFTReal_CURRENT_CODEHEADER
1598
1599
1600
1601/*\\ EOF \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
1602
1603
1604
1605
1606#endif // ffft_FFTReal_HEADER_INCLUDED
1607
1608
1609
1610/*\\ EOF \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
1611
1612/*****************************************************************************
1613
1614 FFTRealFixLen.h
1615 By Laurent de Soras
1616
1617--- Legal stuff ---
1618
1619This program is free software. It comes without any warranty, to
1620the extent permitted by applicable law. You can redistribute it
1621and/or modify it under the terms of the Do What The Fuck You Want
1622To Public License, Version 2, as published by Sam Hocevar. See
1623http://sam.zoy.org/wtfpl/COPYING for more details.
1624
1625*Tab=3***********************************************************************/
1626
1627
1628
1629#if ! defined (ffft_FFTRealFixLen_HEADER_INCLUDED)
1630#define ffft_FFTRealFixLen_HEADER_INCLUDED
1631
1632#if defined (_MSC_VER)
1633 #pragma once
1634 #pragma warning (4 : 4250) // "Inherits via dominance."
1635#endif
1636
1637
1638
1639/*\\ INCLUDE FILES \\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
1640
1641/*****************************************************************************
1642
1643 Array.h
1644 By Laurent de Soras
1645
1646--- Legal stuff ---
1647
1648This program is free software. It comes without any warranty, to
1649the extent permitted by applicable law. You can redistribute it
1650and/or modify it under the terms of the Do What The Fuck You Want
1651To Public License, Version 2, as published by Sam Hocevar. See
1652http://sam.zoy.org/wtfpl/COPYING for more details.
1653
1654*Tab=3***********************************************************************/
1655
1656
1657
1658#if ! defined (ffft_Array_HEADER_INCLUDED)
1659#define ffft_Array_HEADER_INCLUDED
1660
1661#if defined (_MSC_VER)
1662 #pragma once
1663 #pragma warning (4 : 4250) // "Inherits via dominance."
1664#endif
1665
1666
1667
1668/*\\ INCLUDE FILES \\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
1669
1670
1671
1672namespace ffft
1673{
1674
1675
1676
1677template <class T, long LEN>
1678class Array
1679{
1680
1681/*\\ PUBLIC \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
1682
1683public:
1684
1685 typedef T DataType;
1686
1687 Array ();
1688
1689 inline const DataType &
1690 operator [] (long pos) const;
1691 inline DataType &
1692 operator [] (long pos);
1693
1694 static inline long
1695 size ();
1696
1697
1698
1699/*\\ PROTECTED \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
1700
1701protected:
1702
1703
1704
1705/*\\ PRIVATE \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
1706
1707private:
1708
1709 DataType _data_arr [LEN];
1710
1711
1712
1713/*\\ FORBIDDEN MEMBER FUNCTIONS \\\\\\\\\\\\\\\\\\\\\\*/
1714
1715private:
1716
1717 Array (const Array &other);
1718 Array & operator = (const Array &other);
1719 bool operator == (const Array &other);
1720 bool operator != (const Array &other);
1721
1722}; // class Array
1723
1724
1725
1726} // namespace ffft
1727
1728
1729
1730/*****************************************************************************
1731
1732 Array.hpp
1733 By Laurent de Soras
1734
1735--- Legal stuff ---
1736
1737This program is free software. It comes without any warranty, to
1738the extent permitted by applicable law. You can redistribute it
1739and/or modify it under the terms of the Do What The Fuck You Want
1740To Public License, Version 2, as published by Sam Hocevar. See
1741http://sam.zoy.org/wtfpl/COPYING for more details.
1742
1743*Tab=3***********************************************************************/
1744
1745
1746
1747#if defined (ffft_Array_CURRENT_CODEHEADER)
1748 #error Recursive inclusion of Array code header.
1749#endif
1750#define ffft_Array_CURRENT_CODEHEADER
1751
1752#if ! defined (ffft_Array_CODEHEADER_INCLUDED)
1753#define ffft_Array_CODEHEADER_INCLUDED
1754
1755
1756
1757/*\\ INCLUDE FILES \\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
1758
1759
1760namespace ffft
1761{
1762
1763
1764
1765/*\\ PUBLIC \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
1766
1767
1768
1769template <class T, long LEN>
1770Array <T, LEN>::Array ()
1771{
1772 // Nothing
1773}
1774
1775
1776
1777template <class T, long LEN>
1778const typename Array <T, LEN>::DataType & Array <T, LEN>::operator [] (long pos) const
1779{
1780 assert (pos >= 0);
1781 assert (pos < LEN);
1782
1783 return (_data_arr [pos]);
1784}
1785
1786
1787
1788template <class T, long LEN>
1789typename Array <T, LEN>::DataType & Array <T, LEN>::operator [] (long pos)
1790{
1791 assert (pos >= 0);
1792 assert (pos < LEN);
1793
1794 return (_data_arr [pos]);
1795}
1796
1797
1798
1799template <class T, long LEN>
1800long Array <T, LEN>::size ()
1801{
1802 return (LEN);
1803}
1804
1805
1806
1807/*\\ PROTECTED \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
1808
1809
1810
1811/*\\ PRIVATE \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
1812
1813
1814
1815} // namespace ffft
1816
1817
1818
1819#endif // ffft_Array_CODEHEADER_INCLUDED
1820
1821#undef ffft_Array_CURRENT_CODEHEADER
1822
1823
1824
1825/*\\ EOF \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
1826
1827
1828
1829
1830#endif // ffft_Array_HEADER_INCLUDED
1831
1832
1833
1834/*\\ EOF \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
1835
1836/*****************************************************************************
1837
1838 DynArray.h
1839 By Laurent de Soras
1840
1841--- Legal stuff ---
1842
1843This program is free software. It comes without any warranty, to
1844the extent permitted by applicable law. You can redistribute it
1845and/or modify it under the terms of the Do What The Fuck You Want
1846To Public License, Version 2, as published by Sam Hocevar. See
1847http://sam.zoy.org/wtfpl/COPYING for more details.
1848
1849*Tab=3***********************************************************************/
1850
1851
1852
1853#if ! defined (ffft_DynArray_HEADER_INCLUDED)
1854#define ffft_DynArray_HEADER_INCLUDED
1855
1856#if defined (_MSC_VER)
1857 #pragma once
1858 #pragma warning (4 : 4250) // "Inherits via dominance."
1859#endif
1860
1861
1862
1863/*\\ INCLUDE FILES \\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
1864
1865
1866
1867namespace ffft
1868{
1869
1870
1871
1872template <class T>
1873class DynArray
1874{
1875
1876/*\\ PUBLIC \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
1877
1878public:
1879
1880 typedef T DataType;
1881
1882 DynArray ();
1883 explicit DynArray (long size);
1884 ~DynArray ();
1885
1886 inline long size () const;
1887 inline void resize (long size);
1888
1889 inline const DataType &
1890 operator [] (long pos) const;
1891 inline DataType &
1892 operator [] (long pos);
1893
1894
1895
1896/*\\ PROTECTED \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
1897
1898protected:
1899
1900
1901
1902/*\\ PRIVATE \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
1903
1904private:
1905
1906 DataType * _data_ptr;
1907 long _len;
1908
1909
1910
1911/*\\ FORBIDDEN MEMBER FUNCTIONS \\\\\\\\\\\\\\\\\\\\\\*/
1912
1913private:
1914
1915 DynArray (const DynArray &other);
1916 DynArray & operator = (const DynArray &other);
1917 bool operator == (const DynArray &other);
1918 bool operator != (const DynArray &other);
1919
1920}; // class DynArray
1921
1922
1923
1924} // namespace ffft
1925
1926
1927
1928/*****************************************************************************
1929
1930 DynArray.hpp
1931 By Laurent de Soras
1932
1933--- Legal stuff ---
1934
1935This program is free software. It comes without any warranty, to
1936the extent permitted by applicable law. You can redistribute it
1937and/or modify it under the terms of the Do What The Fuck You Want
1938To Public License, Version 2, as published by Sam Hocevar. See
1939http://sam.zoy.org/wtfpl/COPYING for more details.
1940
1941*Tab=3***********************************************************************/
1942
1943
1944
1945#if defined (ffft_DynArray_CURRENT_CODEHEADER)
1946 #error Recursive inclusion of DynArray code header.
1947#endif
1948#define ffft_DynArray_CURRENT_CODEHEADER
1949
1950#if ! defined (ffft_DynArray_CODEHEADER_INCLUDED)
1951#define ffft_DynArray_CODEHEADER_INCLUDED
1952
1953
1954
1955/*\\ INCLUDE FILES \\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
1956
1957
1958namespace ffft
1959{
1960
1961
1962
1963/*\\ PUBLIC \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
1964
1965
1966
1967template <class T>
1968DynArray <T>::DynArray ()
1969: _data_ptr (0)
1970, _len (0)
1971{
1972 // Nothing
1973}
1974
1975
1976
1977template <class T>
1978DynArray <T>::DynArray (long size)
1979: _data_ptr (0)
1980, _len (0)
1981{
1982 assert (size >= 0);
1983 if (size > 0)
1984 {
1985#ifdef FFT_CUSTOM_ALLOC
1986 _data_ptr = FFT_CUSTOM_ALLOC.createArray<DataType>(size);
1987#else
1988 _data_ptr = new DataType [size];
1989#endif
1990 _len = size;
1991 }
1992}
1993
1994
1995
1996template <class T>
1997DynArray <T>::~DynArray ()
1998{
1999#ifdef FFT_CUSTOM_ALLOC
2000 FFT_CUSTOM_ALLOC.removeArray<T>(_data_ptr, _len);
2001#else
2002 delete [] _data_ptr;
2003#endif
2004 _data_ptr = 0;
2005 _len = 0;
2006}
2007
2008
2009
2010template <class T>
2011long DynArray <T>::size () const
2012{
2013 return (_len);
2014}
2015
2016
2017
2018template <class T>
2019void DynArray <T>::resize (long size)
2020{
2021 assert (size >= 0);
2022 if (size > 0)
2023 {
2024 DataType * old_data_ptr = _data_ptr;
2025#ifdef FFT_CUSTOM_ALLOC
2026#else
2027 DataType * tmp_data_ptr = new DataType [size];
2028#endif
2029
2030 _data_ptr = tmp_data_ptr;
2031 _len = size;
2032
2033#ifdef FFT_CUSTOM_ALLOC
2034#else
2035 delete [] old_data_ptr;
2036#endif
2037 }
2038}
2039
2040
2041
2042template <class T>
2043const typename DynArray <T>::DataType & DynArray <T>::operator [] (long pos) const
2044{
2045 assert (pos >= 0);
2046 assert (pos < _len);
2047
2048 return (_data_ptr [pos]);
2049}
2050
2051
2052
2053template <class T>
2054typename DynArray <T>::DataType & DynArray <T>::operator [] (long pos)
2055{
2056 assert (pos >= 0);
2057 assert (pos < _len);
2058
2059 return (_data_ptr [pos]);
2060}
2061
2062
2063
2064/*\\ PROTECTED \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
2065
2066
2067
2068/*\\ PRIVATE \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
2069
2070
2071
2072} // namespace ffft
2073
2074
2075
2076#endif // ffft_DynArray_CODEHEADER_INCLUDED
2077
2078#undef ffft_DynArray_CURRENT_CODEHEADER
2079
2080
2081
2082/*\\ EOF \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
2083
2084
2085
2086
2087#endif // ffft_DynArray_HEADER_INCLUDED
2088
2089
2090
2091/*\\ EOF \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
2092
2093/*****************************************************************************
2094
2095 FFTRealFixLenParam.h
2096 By Laurent de Soras
2097
2098--- Legal stuff ---
2099
2100This program is free software. It comes without any warranty, to
2101the extent permitted by applicable law. You can redistribute it
2102and/or modify it under the terms of the Do What The Fuck You Want
2103To Public License, Version 2, as published by Sam Hocevar. See
2104http://sam.zoy.org/wtfpl/COPYING for more details.
2105
2106*Tab=3***********************************************************************/
2107
2108
2109
2110#if ! defined (ffft_FFTRealFixLenParam_HEADER_INCLUDED)
2111#define ffft_FFTRealFixLenParam_HEADER_INCLUDED
2112
2113#if defined (_MSC_VER)
2114 #pragma once
2115 #pragma warning (4 : 4250) // "Inherits via dominance."
2116#endif
2117
2118
2119
2120/*\\ INCLUDE FILES \\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
2121
2122
2123
2124namespace ffft
2125{
2126
2127
2128
2129class FFTRealFixLenParam
2130{
2131
2132/*\\ PUBLIC \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
2133
2134public:
2135
2136 // Over this bit depth, we use direct calculation for sin/cos
2137 enum { TRIGO_BD_LIMIT = 12 };
2138
2139 typedef float DataType;
2140
2141
2142
2143/*\\ PROTECTED \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
2144
2145protected:
2146
2147
2148
2149/*\\ PRIVATE \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
2150
2151private:
2152
2153
2154
2155/*\\ FORBIDDEN MEMBER FUNCTIONS \\\\\\\\\\\\\\\\\\\\\\*/
2156
2157private:
2158
2159 FFTRealFixLenParam ();
2160 FFTRealFixLenParam (const FFTRealFixLenParam &other);
2161 FFTRealFixLenParam &
2162 operator = (const FFTRealFixLenParam &other);
2163 bool operator == (const FFTRealFixLenParam &other);
2164 bool operator != (const FFTRealFixLenParam &other);
2165
2166}; // class FFTRealFixLenParam
2167
2168
2169
2170} // namespace ffft
2171
2172
2173
2174//#include "ffft/FFTRealFixLenParam.hpp"
2175
2176
2177
2178#endif // ffft_FFTRealFixLenParam_HEADER_INCLUDED
2179
2180
2181
2182/*\\ EOF \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
2183
2184/*****************************************************************************
2185
2186 OscSinCos.h
2187 By Laurent de Soras
2188
2189--- Legal stuff ---
2190
2191This program is free software. It comes without any warranty, to
2192the extent permitted by applicable law. You can redistribute it
2193and/or modify it under the terms of the Do What The Fuck You Want
2194To Public License, Version 2, as published by Sam Hocevar. See
2195http://sam.zoy.org/wtfpl/COPYING for more details.
2196
2197*Tab=3***********************************************************************/
2198
2199
2200
2201#if ! defined (ffft_OscSinCos_HEADER_INCLUDED)
2202#define ffft_OscSinCos_HEADER_INCLUDED
2203
2204#if defined (_MSC_VER)
2205 #pragma once
2206 #pragma warning (4 : 4250) // "Inherits via dominance."
2207#endif
2208
2209
2210
2211/*\\ INCLUDE FILES \\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
2212
2213
2214
2215
2216
2217namespace ffft
2218{
2219
2220
2221
2222template <class T>
2223class OscSinCos
2224{
2225
2226/*\\ PUBLIC \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
2227
2228public:
2229
2230 typedef T DataType;
2231
2232 OscSinCos ();
2233
2234 ffft_FORCEINLINE void
2235 set_step (double angle_rad);
2236
2237 ffft_FORCEINLINE DataType
2238 get_cos () const;
2239 ffft_FORCEINLINE DataType
2240 get_sin () const;
2241 ffft_FORCEINLINE void
2242 step ();
2243 ffft_FORCEINLINE void
2244 clear_buffers ();
2245
2246
2247
2248/*\\ PROTECTED \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
2249
2250protected:
2251
2252
2253
2254/*\\ PRIVATE \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
2255
2256private:
2257
2258 DataType _pos_cos; // Current phase expressed with sin and cos. [-1 ; 1]
2259 DataType _pos_sin; // -
2260 DataType _step_cos; // Phase increment per step, [-1 ; 1]
2261 DataType _step_sin; // -
2262
2263
2264
2265/*\\ FORBIDDEN MEMBER FUNCTIONS \\\\\\\\\\\\\\\\\\\\\\*/
2266
2267private:
2268
2269 OscSinCos (const OscSinCos &other);
2270 OscSinCos & operator = (const OscSinCos &other);
2271 bool operator == (const OscSinCos &other);
2272 bool operator != (const OscSinCos &other);
2273
2274}; // class OscSinCos
2275
2276
2277
2278} // namespace ffft
2279
2280
2281
2282/*****************************************************************************
2283
2284 OscSinCos.hpp
2285 By Laurent de Soras
2286
2287--- Legal stuff ---
2288
2289This program is free software. It comes without any warranty, to
2290the extent permitted by applicable law. You can redistribute it
2291and/or modify it under the terms of the Do What The Fuck You Want
2292To Public License, Version 2, as published by Sam Hocevar. See
2293http://sam.zoy.org/wtfpl/COPYING for more details.
2294
2295*Tab=3***********************************************************************/
2296
2297
2298
2299#if defined (ffft_OscSinCos_CURRENT_CODEHEADER)
2300 #error Recursive inclusion of OscSinCos code header.
2301#endif
2302#define ffft_OscSinCos_CURRENT_CODEHEADER
2303
2304#if ! defined (ffft_OscSinCos_CODEHEADER_INCLUDED)
2305#define ffft_OscSinCos_CODEHEADER_INCLUDED
2306
2307
2308
2309/*\\ INCLUDE FILES \\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
2310
2311#include <cmath>
2312
2313namespace std { }
2314
2315
2316
2317namespace ffft
2318{
2319
2320
2321
2322/*\\ PUBLIC \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
2323
2324
2325
2326template <class T>
2327OscSinCos <T>::OscSinCos ()
2328: _pos_cos (1)
2329, _pos_sin (0)
2330, _step_cos (1)
2331, _step_sin (0)
2332{
2333 // Nothing
2334}
2335
2336
2337
2338template <class T>
2339void OscSinCos <T>::set_step (double angle_rad)
2340{
2341 using namespace std;
2342
2343 _step_cos = static_cast <DataType> (cos (angle_rad));
2344 _step_sin = static_cast <DataType> (sin (angle_rad));
2345}
2346
2347
2348
2349template <class T>
2350typename OscSinCos <T>::DataType OscSinCos <T>::get_cos () const
2351{
2352 return (_pos_cos);
2353}
2354
2355
2356
2357template <class T>
2358typename OscSinCos <T>::DataType OscSinCos <T>::get_sin () const
2359{
2360 return (_pos_sin);
2361}
2362
2363
2364
2365template <class T>
2366void OscSinCos <T>::step ()
2367{
2368 const DataType old_cos = _pos_cos;
2369 const DataType old_sin = _pos_sin;
2370
2371 _pos_cos = old_cos * _step_cos - old_sin * _step_sin;
2372 _pos_sin = old_cos * _step_sin + old_sin * _step_cos;
2373}
2374
2375
2376
2377template <class T>
2378void OscSinCos <T>::clear_buffers ()
2379{
2380 _pos_cos = static_cast <DataType> (1);
2381 _pos_sin = static_cast <DataType> (0);
2382}
2383
2384
2385
2386/*\\ PROTECTED \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
2387
2388
2389
2390/*\\ PRIVATE \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
2391
2392
2393
2394} // namespace ffft
2395
2396
2397
2398#endif // ffft_OscSinCos_CODEHEADER_INCLUDED
2399
2400#undef ffft_OscSinCos_CURRENT_CODEHEADER
2401
2402
2403
2404/*\\ EOF \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
2405
2406
2407
2408
2409#endif // ffft_OscSinCos_HEADER_INCLUDED
2410
2411
2412
2413/*\\ EOF \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
2414
2415
2416
2417
2418namespace ffft
2419{
2420
2421
2422
2423template <int LL2>
2424class FFTRealFixLen
2425{
2426 typedef int CompileTimeCheck1 [(LL2 >= 0) ? 1 : -1];
2427 typedef int CompileTimeCheck2 [(LL2 <= 30) ? 1 : -1];
2428
2429/*\\ PUBLIC \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
2430
2431public:
2432
2433 typedef FFTRealFixLenParam::DataType DataType;
2434 typedef OscSinCos <DataType> OscType;
2435
2436 enum { FFT_LEN_L2 = LL2 };
2437 enum { FFT_LEN = 1 << FFT_LEN_L2 };
2438
2439 FFTRealFixLen ();
2440
2441 inline long get_length () const;
2442 void do_fft (DataType f [], const DataType x []);
2443 void do_ifft (const DataType f [], DataType x []);
2444 void rescale (DataType x []) const;
2445
2446
2447
2448/*\\ PROTECTED \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
2449
2450protected:
2451
2452
2453
2454/*\\ PRIVATE \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
2455
2456private:
2457
2458 enum { TRIGO_BD_LIMIT = FFTRealFixLenParam::TRIGO_BD_LIMIT };
2459
2460 enum { BR_ARR_SIZE_L2 = ((FFT_LEN_L2 - 3) < 0) ? 0 : (FFT_LEN_L2 - 2) };
2461 enum { BR_ARR_SIZE = 1 << BR_ARR_SIZE_L2 };
2462
2463 enum { TRIGO_BD = ((FFT_LEN_L2 - TRIGO_BD_LIMIT) < 0)
2464 ? (int)FFT_LEN_L2
2465 : (int)TRIGO_BD_LIMIT };
2466 enum { TRIGO_TABLE_ARR_SIZE_L2 = (LL2 < 4) ? 0 : (TRIGO_BD - 2) };
2467 enum { TRIGO_TABLE_ARR_SIZE = 1 << TRIGO_TABLE_ARR_SIZE_L2 };
2468
2469 enum { NBR_TRIGO_OSC = FFT_LEN_L2 - TRIGO_BD };
2470 enum { TRIGO_OSC_ARR_SIZE = (NBR_TRIGO_OSC > 0) ? NBR_TRIGO_OSC : 1 };
2471
2472 void build_br_lut ();
2473 void build_trigo_lut ();
2474 void build_trigo_osc ();
2475
2476 DynArray <DataType>
2477 _buffer;
2478 DynArray <long>
2479 _br_data;
2480 DynArray <DataType>
2481 _trigo_data;
2482 Array <OscType, TRIGO_OSC_ARR_SIZE>
2483 _trigo_osc;
2484
2485
2486
2487/*\\ FORBIDDEN MEMBER FUNCTIONS \\\\\\\\\\\\\\\\\\\\\\*/
2488
2489private:
2490
2491 FFTRealFixLen (const FFTRealFixLen &other);
2492 FFTRealFixLen& operator = (const FFTRealFixLen &other);
2493 bool operator == (const FFTRealFixLen &other);
2494 bool operator != (const FFTRealFixLen &other);
2495
2496}; // class FFTRealFixLen
2497
2498
2499
2500} // namespace ffft
2501
2502
2503
2504/*****************************************************************************
2505
2506 FFTRealFixLen.hpp
2507 By Laurent de Soras
2508
2509--- Legal stuff ---
2510
2511This program is free software. It comes without any warranty, to
2512the extent permitted by applicable law. You can redistribute it
2513and/or modify it under the terms of the Do What The Fuck You Want
2514To Public License, Version 2, as published by Sam Hocevar. See
2515http://sam.zoy.org/wtfpl/COPYING for more details.
2516
2517*Tab=3***********************************************************************/
2518
2519
2520
2521#if defined (ffft_FFTRealFixLen_CURRENT_CODEHEADER)
2522 #error Recursive inclusion of FFTRealFixLen code header.
2523#endif
2524#define ffft_FFTRealFixLen_CURRENT_CODEHEADER
2525
2526#if ! defined (ffft_FFTRealFixLen_CODEHEADER_INCLUDED)
2527#define ffft_FFTRealFixLen_CODEHEADER_INCLUDED
2528
2529
2530
2531/*\\ INCLUDE FILES \\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
2532
2533
2534/*****************************************************************************
2535
2536 FFTRealPassDirect.h
2537 By Laurent de Soras
2538
2539--- Legal stuff ---
2540
2541This program is free software. It comes without any warranty, to
2542the extent permitted by applicable law. You can redistribute it
2543and/or modify it under the terms of the Do What The Fuck You Want
2544To Public License, Version 2, as published by Sam Hocevar. See
2545http://sam.zoy.org/wtfpl/COPYING for more details.
2546
2547*Tab=3***********************************************************************/
2548
2549
2550
2551#if ! defined (ffft_FFTRealPassDirect_HEADER_INCLUDED)
2552#define ffft_FFTRealPassDirect_HEADER_INCLUDED
2553
2554#if defined (_MSC_VER)
2555 #pragma once
2556 #pragma warning (4 : 4250) // "Inherits via dominance."
2557#endif
2558
2559
2560
2561/*\\ INCLUDE FILES \\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
2562
2563
2564
2565
2566
2567
2568
2569namespace ffft
2570{
2571
2572
2573
2574template <int PASS>
2575class FFTRealPassDirect
2576{
2577
2578/*\\ PUBLIC \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
2579
2580public:
2581
2582 typedef FFTRealFixLenParam::DataType DataType;
2583 typedef OscSinCos <DataType> OscType;
2584
2585 ffft_FORCEINLINE static void
2586 process (long len, DataType dest_ptr [], DataType src_ptr [], const DataType x_ptr [], const DataType cos_ptr [], long cos_len, const long br_ptr [], OscType osc_list []);
2587
2588
2589
2590/*\\ PROTECTED \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
2591
2592protected:
2593
2594
2595
2596/*\\ PRIVATE \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
2597
2598private:
2599
2600
2601
2602/*\\ FORBIDDEN MEMBER FUNCTIONS \\\\\\\\\\\\\\\\\\\\\\*/
2603
2604private:
2605
2606 FFTRealPassDirect ();
2607 FFTRealPassDirect (const FFTRealPassDirect &other);
2608 FFTRealPassDirect &
2609 operator = (const FFTRealPassDirect &other);
2610 bool operator == (const FFTRealPassDirect &other);
2611 bool operator != (const FFTRealPassDirect &other);
2612
2613}; // class FFTRealPassDirect
2614
2615
2616
2617} // namespace ffft
2618
2619
2620
2621/*****************************************************************************
2622
2623 FFTRealPassDirect.hpp
2624 By Laurent de Soras
2625
2626--- Legal stuff ---
2627
2628This program is free software. It comes without any warranty, to
2629the extent permitted by applicable law. You can redistribute it
2630and/or modify it under the terms of the Do What The Fuck You Want
2631To Public License, Version 2, as published by Sam Hocevar. See
2632http://sam.zoy.org/wtfpl/COPYING for more details.
2633
2634*Tab=3***********************************************************************/
2635
2636
2637
2638#if defined (ffft_FFTRealPassDirect_CURRENT_CODEHEADER)
2639 #error Recursive inclusion of FFTRealPassDirect code header.
2640#endif
2641#define ffft_FFTRealPassDirect_CURRENT_CODEHEADER
2642
2643#if ! defined (ffft_FFTRealPassDirect_CODEHEADER_INCLUDED)
2644#define ffft_FFTRealPassDirect_CODEHEADER_INCLUDED
2645
2646
2647
2648/*\\ INCLUDE FILES \\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
2649
2650/*****************************************************************************
2651
2652 FFTRealUseTrigo.h
2653 By Laurent de Soras
2654
2655--- Legal stuff ---
2656
2657This program is free software. It comes without any warranty, to
2658the extent permitted by applicable law. You can redistribute it
2659and/or modify it under the terms of the Do What The Fuck You Want
2660To Public License, Version 2, as published by Sam Hocevar. See
2661http://sam.zoy.org/wtfpl/COPYING for more details.
2662
2663*Tab=3***********************************************************************/
2664
2665
2666
2667#if ! defined (ffft_FFTRealUseTrigo_HEADER_INCLUDED)
2668#define ffft_FFTRealUseTrigo_HEADER_INCLUDED
2669
2670#if defined (_MSC_VER)
2671 #pragma once
2672 #pragma warning (4 : 4250) // "Inherits via dominance."
2673#endif
2674
2675
2676
2677/*\\ INCLUDE FILES \\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
2678
2679
2680
2681
2682
2683
2684
2685namespace ffft
2686{
2687
2688
2689
2690template <int ALGO>
2691class FFTRealUseTrigo
2692{
2693
2694/*\\ PUBLIC \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
2695
2696public:
2697
2698 typedef FFTRealFixLenParam::DataType DataType;
2699 typedef OscSinCos <DataType> OscType;
2700
2701 ffft_FORCEINLINE static void
2702 prepare (OscType &osc);
2703 ffft_FORCEINLINE static void
2704 iterate (OscType &osc, DataType &c, DataType &s, const DataType cos_ptr [], long index_c, long index_s);
2705
2706
2707
2708/*\\ PROTECTED \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
2709
2710protected:
2711
2712
2713
2714/*\\ PRIVATE \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
2715
2716private:
2717
2718
2719
2720/*\\ FORBIDDEN MEMBER FUNCTIONS \\\\\\\\\\\\\\\\\\\\\\*/
2721
2722private:
2723
2724 FFTRealUseTrigo ();
2725 ~FFTRealUseTrigo ();
2726 FFTRealUseTrigo (const FFTRealUseTrigo &other);
2727 FFTRealUseTrigo &
2728 operator = (const FFTRealUseTrigo &other);
2729 bool operator == (const FFTRealUseTrigo &other);
2730 bool operator != (const FFTRealUseTrigo &other);
2731
2732}; // class FFTRealUseTrigo
2733
2734
2735
2736} // namespace ffft
2737
2738
2739
2740/*****************************************************************************
2741
2742 FFTRealUseTrigo.hpp
2743 By Laurent de Soras
2744
2745--- Legal stuff ---
2746
2747This program is free software. It comes without any warranty, to
2748the extent permitted by applicable law. You can redistribute it
2749and/or modify it under the terms of the Do What The Fuck You Want
2750To Public License, Version 2, as published by Sam Hocevar. See
2751http://sam.zoy.org/wtfpl/COPYING for more details.
2752
2753*Tab=3***********************************************************************/
2754
2755
2756
2757#if defined (ffft_FFTRealUseTrigo_CURRENT_CODEHEADER)
2758 #error Recursive inclusion of FFTRealUseTrigo code header.
2759#endif
2760#define ffft_FFTRealUseTrigo_CURRENT_CODEHEADER
2761
2762#if ! defined (ffft_FFTRealUseTrigo_CODEHEADER_INCLUDED)
2763#define ffft_FFTRealUseTrigo_CODEHEADER_INCLUDED
2764
2765
2766
2767/*\\ INCLUDE FILES \\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
2768
2769
2770
2771
2772
2773namespace ffft
2774{
2775
2776
2777
2778/*\\ PUBLIC \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
2779
2780
2781
2782template <int ALGO>
2783void FFTRealUseTrigo <ALGO>::prepare (OscType &osc)
2784{
2785 osc.clear_buffers ();
2786}
2787
2788template <>
2789inline void FFTRealUseTrigo <0>::prepare (OscType &osc)
2790{
2791 // Nothing
2792}
2793
2794
2795
2796template <int ALGO>
2797void FFTRealUseTrigo <ALGO>::iterate (OscType &osc, DataType &c, DataType &s, const DataType cos_ptr [], long index_c, long index_s)
2798{
2799 osc.step ();
2800 c = osc.get_cos ();
2801 s = osc.get_sin ();
2802}
2803
2804template <>
2805inline void FFTRealUseTrigo <0>::iterate (OscType &osc, DataType &c, DataType &s, const DataType cos_ptr [], long index_c, long index_s)
2806{
2807 c = cos_ptr [index_c];
2808 s = cos_ptr [index_s];
2809}
2810
2811
2812
2813/*\\ PROTECTED \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
2814
2815
2816
2817/*\\ PRIVATE \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
2818
2819
2820
2821} // namespace ffft
2822
2823
2824
2825#endif // ffft_FFTRealUseTrigo_CODEHEADER_INCLUDED
2826
2827#undef ffft_FFTRealUseTrigo_CURRENT_CODEHEADER
2828
2829
2830
2831/*\\ EOF \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
2832
2833
2834
2835
2836#endif // ffft_FFTRealUseTrigo_HEADER_INCLUDED
2837
2838
2839
2840/*\\ EOF \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
2841
2842
2843
2844
2845namespace ffft
2846{
2847
2848
2849
2850/*\\ PUBLIC \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
2851
2852
2853
2854template <>
2855inline void FFTRealPassDirect <1>::process (long len, DataType dest_ptr [], DataType src_ptr [], const DataType x_ptr [], const DataType cos_ptr [], long cos_len, const long br_ptr [], OscType osc_list [])
2856{
2857 // First and second pass at once
2858 const long qlen = len >> 2;
2859
2860 long coef_index = 0;
2861 do
2862 {
2863 // To do: unroll the loop (2x).
2864 const long ri_0 = br_ptr [coef_index >> 2];
2865 const long ri_1 = ri_0 + 2 * qlen; // bit_rev_lut_ptr [coef_index + 1];
2866 const long ri_2 = ri_0 + 1 * qlen; // bit_rev_lut_ptr [coef_index + 2];
2867 const long ri_3 = ri_0 + 3 * qlen; // bit_rev_lut_ptr [coef_index + 3];
2868
2869 DataType * const df2 = dest_ptr + coef_index;
2870 df2 [1] = x_ptr [ri_0] - x_ptr [ri_1];
2871 df2 [3] = x_ptr [ri_2] - x_ptr [ri_3];
2872
2873 const DataType sf_0 = x_ptr [ri_0] + x_ptr [ri_1];
2874 const DataType sf_2 = x_ptr [ri_2] + x_ptr [ri_3];
2875
2876 df2 [0] = sf_0 + sf_2;
2877 df2 [2] = sf_0 - sf_2;
2878
2879 coef_index += 4;
2880 }
2881 while (coef_index < len);
2882}
2883
2884template <>
2885inline void FFTRealPassDirect <2>::process (long len, DataType dest_ptr [], DataType src_ptr [], const DataType x_ptr [], const DataType cos_ptr [], long cos_len, const long br_ptr [], OscType osc_list [])
2886{
2887 // Executes "previous" passes first. Inverts source and destination buffers
2888 FFTRealPassDirect <1>::process (
2889 len,
2890 src_ptr,
2891 dest_ptr,
2892 x_ptr,
2893 cos_ptr,
2894 cos_len,
2895 br_ptr,
2896 osc_list
2897 );
2898
2899 // Third pass
2900 const DataType sqrt2_2 = DataType (SQRT2 * 0.5);
2901
2902 long coef_index = 0;
2903 do
2904 {
2905 dest_ptr [coef_index ] = src_ptr [coef_index] + src_ptr [coef_index + 4];
2906 dest_ptr [coef_index + 4] = src_ptr [coef_index] - src_ptr [coef_index + 4];
2907 dest_ptr [coef_index + 2] = src_ptr [coef_index + 2];
2908 dest_ptr [coef_index + 6] = src_ptr [coef_index + 6];
2909
2910 DataType v;
2911
2912 v = (src_ptr [coef_index + 5] - src_ptr [coef_index + 7]) * sqrt2_2;
2913 dest_ptr [coef_index + 1] = src_ptr [coef_index + 1] + v;
2914 dest_ptr [coef_index + 3] = src_ptr [coef_index + 1] - v;
2915
2916 v = (src_ptr [coef_index + 5] + src_ptr [coef_index + 7]) * sqrt2_2;
2917 dest_ptr [coef_index + 5] = v + src_ptr [coef_index + 3];
2918 dest_ptr [coef_index + 7] = v - src_ptr [coef_index + 3];
2919
2920 coef_index += 8;
2921 }
2922 while (coef_index < len);
2923}
2924
2925template <int PASS>
2926void FFTRealPassDirect <PASS>::process (long len, DataType dest_ptr [], DataType src_ptr [], const DataType x_ptr [], const DataType cos_ptr [], long cos_len, const long br_ptr [], OscType osc_list [])
2927{
2928 // Executes "previous" passes first. Inverts source and destination buffers
2929 FFTRealPassDirect <PASS - 1>::process (
2930 len,
2931 src_ptr,
2932 dest_ptr,
2933 x_ptr,
2934 cos_ptr,
2935 cos_len,
2936 br_ptr,
2937 osc_list
2938 );
2939
2940 const long dist = 1L << (PASS - 1);
2941 const long c1_r = 0;
2942 const long c1_i = dist;
2943 const long c2_r = dist * 2;
2944 const long c2_i = dist * 3;
2945 const long cend = dist * 4;
2946 const long table_step = cos_len >> (PASS - 1);
2947
2948 enum { TRIGO_OSC = PASS - FFTRealFixLenParam::TRIGO_BD_LIMIT };
2949 enum { TRIGO_DIRECT = (TRIGO_OSC >= 0) ? 1 : 0 };
2950
2951 long coef_index = 0;
2952 do
2953 {
2954 const DataType * const sf = src_ptr + coef_index;
2955 DataType * const df = dest_ptr + coef_index;
2956
2957 // Extreme coefficients are always real
2958 df [c1_r] = sf [c1_r] + sf [c2_r];
2959 df [c2_r] = sf [c1_r] - sf [c2_r];
2960 df [c1_i] = sf [c1_i];
2961 df [c2_i] = sf [c2_i];
2962
2963 FFTRealUseTrigo <TRIGO_DIRECT>::prepare (osc_list [TRIGO_OSC]);
2964
2965 // Others are conjugate complex numbers
2966 for (long i = 1; i < dist; ++ i)
2967 {
2968 DataType c;
2969 DataType s;
2970 FFTRealUseTrigo <TRIGO_DIRECT>::iterate (
2971 osc_list [TRIGO_OSC],
2972 c,
2973 s,
2974 cos_ptr,
2975 i * table_step,
2976 (dist - i) * table_step
2977 );
2978
2979 const DataType sf_r_i = sf [c1_r + i];
2980 const DataType sf_i_i = sf [c1_i + i];
2981
2982 const DataType v1 = sf [c2_r + i] * c - sf [c2_i + i] * s;
2983 df [c1_r + i] = sf_r_i + v1;
2984 df [c2_r - i] = sf_r_i - v1;
2985
2986 const DataType v2 = sf [c2_r + i] * s + sf [c2_i + i] * c;
2987 df [c2_r + i] = v2 + sf_i_i;
2988 df [cend - i] = v2 - sf_i_i;
2989 }
2990
2991 coef_index += cend;
2992 }
2993 while (coef_index < len);
2994}
2995
2996
2997
2998/*\\ PROTECTED \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
2999
3000
3001
3002/*\\ PRIVATE \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
3003
3004
3005
3006} // namespace ffft
3007
3008
3009
3010#endif // ffft_FFTRealPassDirect_CODEHEADER_INCLUDED
3011
3012#undef ffft_FFTRealPassDirect_CURRENT_CODEHEADER
3013
3014
3015
3016/*\\ EOF \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
3017
3018
3019
3020
3021#endif // ffft_FFTRealPassDirect_HEADER_INCLUDED
3022
3023
3024
3025/*\\ EOF \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
3026
3027/*****************************************************************************
3028
3029 FFTRealPassInverse.h
3030 By Laurent de Soras
3031
3032--- Legal stuff ---
3033
3034This program is free software. It comes without any warranty, to
3035the extent permitted by applicable law. You can redistribute it
3036and/or modify it under the terms of the Do What The Fuck You Want
3037To Public License, Version 2, as published by Sam Hocevar. See
3038http://sam.zoy.org/wtfpl/COPYING for more details.
3039
3040*Tab=3***********************************************************************/
3041
3042
3043
3044#if ! defined (ffft_FFTRealPassInverse_HEADER_INCLUDED)
3045#define ffft_FFTRealPassInverse_HEADER_INCLUDED
3046
3047#if defined (_MSC_VER)
3048 #pragma once
3049 #pragma warning (4 : 4250) // "Inherits via dominance."
3050#endif
3051
3052
3053
3054/*\\ INCLUDE FILES \\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
3055
3056
3057
3058
3059
3060
3061
3062
3063namespace ffft
3064{
3065
3066
3067
3068template <int PASS>
3069class FFTRealPassInverse
3070{
3071
3072/*\\ PUBLIC \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
3073
3074public:
3075
3076 typedef FFTRealFixLenParam::DataType DataType;
3077 typedef OscSinCos <DataType> OscType;
3078
3079 ffft_FORCEINLINE static void
3080 process (long len, DataType dest_ptr [], DataType src_ptr [], const DataType f_ptr [], const DataType cos_ptr [], long cos_len, const long br_ptr [], OscType osc_list []);
3081 ffft_FORCEINLINE static void
3082 process_rec (long len, DataType dest_ptr [], DataType src_ptr [], const DataType cos_ptr [], long cos_len, const long br_ptr [], OscType osc_list []);
3083 ffft_FORCEINLINE static void
3084 process_internal (long len, DataType dest_ptr [], const DataType src_ptr [], const DataType cos_ptr [], long cos_len, const long br_ptr [], OscType osc_list []);
3085
3086
3087
3088/*\\ PROTECTED \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
3089
3090protected:
3091
3092
3093
3094/*\\ PRIVATE \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
3095
3096private:
3097
3098
3099
3100/*\\ FORBIDDEN MEMBER FUNCTIONS \\\\\\\\\\\\\\\\\\\\\\*/
3101
3102private:
3103
3104 FFTRealPassInverse ();
3105 FFTRealPassInverse (const FFTRealPassInverse &other);
3106 FFTRealPassInverse &
3107 operator = (const FFTRealPassInverse &other);
3108 bool operator == (const FFTRealPassInverse &other);
3109 bool operator != (const FFTRealPassInverse &other);
3110
3111}; // class FFTRealPassInverse
3112
3113
3114
3115} // namespace ffft
3116
3117
3118
3119/*****************************************************************************
3120
3121 FFTRealPassInverse.hpp
3122 By Laurent de Soras
3123
3124--- Legal stuff ---
3125
3126This program is free software. It comes without any warranty, to
3127the extent permitted by applicable law. You can redistribute it
3128and/or modify it under the terms of the Do What The Fuck You Want
3129To Public License, Version 2, as published by Sam Hocevar. See
3130http://sam.zoy.org/wtfpl/COPYING for more details.
3131
3132*Tab=3***********************************************************************/
3133
3134
3135
3136#if defined (ffft_FFTRealPassInverse_CURRENT_CODEHEADER)
3137 #error Recursive inclusion of FFTRealPassInverse code header.
3138#endif
3139#define ffft_FFTRealPassInverse_CURRENT_CODEHEADER
3140
3141#if ! defined (ffft_FFTRealPassInverse_CODEHEADER_INCLUDED)
3142#define ffft_FFTRealPassInverse_CODEHEADER_INCLUDED
3143
3144
3145
3146/*\\ INCLUDE FILES \\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
3147
3148/*****************************************************************************
3149
3150 FFTRealUseTrigo.h
3151 By Laurent de Soras
3152
3153--- Legal stuff ---
3154
3155This program is free software. It comes without any warranty, to
3156the extent permitted by applicable law. You can redistribute it
3157and/or modify it under the terms of the Do What The Fuck You Want
3158To Public License, Version 2, as published by Sam Hocevar. See
3159http://sam.zoy.org/wtfpl/COPYING for more details.
3160
3161*Tab=3***********************************************************************/
3162
3163
3164
3165#if ! defined (ffft_FFTRealUseTrigo_HEADER_INCLUDED)
3166#define ffft_FFTRealUseTrigo_HEADER_INCLUDED
3167
3168#if defined (_MSC_VER)
3169 #pragma once
3170 #pragma warning (4 : 4250) // "Inherits via dominance."
3171#endif
3172
3173
3174
3175/*\\ INCLUDE FILES \\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
3176
3177
3178
3179
3180
3181
3182
3183namespace ffft
3184{
3185
3186
3187
3188template <int ALGO>
3189class FFTRealUseTrigo
3190{
3191
3192/*\\ PUBLIC \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
3193
3194public:
3195
3196 typedef FFTRealFixLenParam::DataType DataType;
3197 typedef OscSinCos <DataType> OscType;
3198
3199 ffft_FORCEINLINE static void
3200 prepare (OscType &osc);
3201 ffft_FORCEINLINE static void
3202 iterate (OscType &osc, DataType &c, DataType &s, const DataType cos_ptr [], long index_c, long index_s);
3203
3204
3205
3206/*\\ PROTECTED \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
3207
3208protected:
3209
3210
3211
3212/*\\ PRIVATE \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
3213
3214private:
3215
3216
3217
3218/*\\ FORBIDDEN MEMBER FUNCTIONS \\\\\\\\\\\\\\\\\\\\\\*/
3219
3220private:
3221
3222 FFTRealUseTrigo ();
3223 ~FFTRealUseTrigo ();
3224 FFTRealUseTrigo (const FFTRealUseTrigo &other);
3225 FFTRealUseTrigo &
3226 operator = (const FFTRealUseTrigo &other);
3227 bool operator == (const FFTRealUseTrigo &other);
3228 bool operator != (const FFTRealUseTrigo &other);
3229
3230}; // class FFTRealUseTrigo
3231
3232
3233
3234} // namespace ffft
3235
3236
3237
3238/*****************************************************************************
3239
3240 FFTRealUseTrigo.hpp
3241 By Laurent de Soras
3242
3243--- Legal stuff ---
3244
3245This program is free software. It comes without any warranty, to
3246the extent permitted by applicable law. You can redistribute it
3247and/or modify it under the terms of the Do What The Fuck You Want
3248To Public License, Version 2, as published by Sam Hocevar. See
3249http://sam.zoy.org/wtfpl/COPYING for more details.
3250
3251*Tab=3***********************************************************************/
3252
3253
3254
3255#if defined (ffft_FFTRealUseTrigo_CURRENT_CODEHEADER)
3256 #error Recursive inclusion of FFTRealUseTrigo code header.
3257#endif
3258#define ffft_FFTRealUseTrigo_CURRENT_CODEHEADER
3259
3260#if ! defined (ffft_FFTRealUseTrigo_CODEHEADER_INCLUDED)
3261#define ffft_FFTRealUseTrigo_CODEHEADER_INCLUDED
3262
3263
3264
3265/*\\ INCLUDE FILES \\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
3266
3267
3268
3269
3270
3271namespace ffft
3272{
3273
3274
3275
3276/*\\ PUBLIC \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
3277
3278
3279
3280template <int ALGO>
3281void FFTRealUseTrigo <ALGO>::prepare (OscType &osc)
3282{
3283 osc.clear_buffers ();
3284}
3285
3286template <>
3287inline void FFTRealUseTrigo <0>::prepare (OscType &osc)
3288{
3289 // Nothing
3290}
3291
3292
3293
3294template <int ALGO>
3295void FFTRealUseTrigo <ALGO>::iterate (OscType &osc, DataType &c, DataType &s, const DataType cos_ptr [], long index_c, long index_s)
3296{
3297 osc.step ();
3298 c = osc.get_cos ();
3299 s = osc.get_sin ();
3300}
3301
3302template <>
3303inline void FFTRealUseTrigo <0>::iterate (OscType &osc, DataType &c, DataType &s, const DataType cos_ptr [], long index_c, long index_s)
3304{
3305 c = cos_ptr [index_c];
3306 s = cos_ptr [index_s];
3307}
3308
3309
3310
3311/*\\ PROTECTED \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
3312
3313
3314
3315/*\\ PRIVATE \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
3316
3317
3318
3319} // namespace ffft
3320
3321
3322
3323#endif // ffft_FFTRealUseTrigo_CODEHEADER_INCLUDED
3324
3325#undef ffft_FFTRealUseTrigo_CURRENT_CODEHEADER
3326
3327
3328
3329/*\\ EOF \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
3330
3331
3332
3333
3334#endif // ffft_FFTRealUseTrigo_HEADER_INCLUDED
3335
3336
3337
3338/*\\ EOF \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
3339
3340
3341
3342
3343namespace ffft
3344{
3345
3346
3347
3348/*\\ PUBLIC \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
3349
3350
3351
3352template <int PASS>
3353void FFTRealPassInverse <PASS>::process (long len, DataType dest_ptr [], DataType src_ptr [], const DataType f_ptr [], const DataType cos_ptr [], long cos_len, const long br_ptr [], OscType osc_list [])
3354{
3355 process_internal (
3356 len,
3357 dest_ptr,
3358 f_ptr,
3359 cos_ptr,
3360 cos_len,
3361 br_ptr,
3362 osc_list
3363 );
3364 FFTRealPassInverse <PASS - 1>::process_rec (
3365 len,
3366 src_ptr,
3367 dest_ptr,
3368 cos_ptr,
3369 cos_len,
3370 br_ptr,
3371 osc_list
3372 );
3373}
3374
3375
3376
3377template <int PASS>
3378void FFTRealPassInverse <PASS>::process_rec (long len, DataType dest_ptr [], DataType src_ptr [], const DataType cos_ptr [], long cos_len, const long br_ptr [], OscType osc_list [])
3379{
3380 process_internal (
3381 len,
3382 dest_ptr,
3383 src_ptr,
3384 cos_ptr,
3385 cos_len,
3386 br_ptr,
3387 osc_list
3388 );
3389 FFTRealPassInverse <PASS - 1>::process_rec (
3390 len,
3391 src_ptr,
3392 dest_ptr,
3393 cos_ptr,
3394 cos_len,
3395 br_ptr,
3396 osc_list
3397 );
3398}
3399
3400template <>
3401inline void FFTRealPassInverse <0>::process_rec (long len, DataType dest_ptr [], DataType src_ptr [], const DataType cos_ptr [], long cos_len, const long br_ptr [], OscType osc_list [])
3402{
3403 // Stops recursion
3404}
3405
3406
3407
3408template <int PASS>
3409void FFTRealPassInverse <PASS>::process_internal (long len, DataType dest_ptr [], const DataType src_ptr [], const DataType cos_ptr [], long cos_len, const long br_ptr [], OscType osc_list [])
3410{
3411 const long dist = 1L << (PASS - 1);
3412 const long c1_r = 0;
3413 const long c1_i = dist;
3414 const long c2_r = dist * 2;
3415 const long c2_i = dist * 3;
3416 const long cend = dist * 4;
3417 const long table_step = cos_len >> (PASS - 1);
3418
3419 enum { TRIGO_OSC = PASS - FFTRealFixLenParam::TRIGO_BD_LIMIT };
3420 enum { TRIGO_DIRECT = (TRIGO_OSC >= 0) ? 1 : 0 };
3421
3422 long coef_index = 0;
3423 do
3424 {
3425 const DataType * const sf = src_ptr + coef_index;
3426 DataType * const df = dest_ptr + coef_index;
3427
3428 // Extreme coefficients are always real
3429 df [c1_r] = sf [c1_r] + sf [c2_r];
3430 df [c2_r] = sf [c1_r] - sf [c2_r];
3431 df [c1_i] = sf [c1_i] * 2;
3432 df [c2_i] = sf [c2_i] * 2;
3433
3434 FFTRealUseTrigo <TRIGO_DIRECT>::prepare (osc_list [TRIGO_OSC]);
3435
3436 // Others are conjugate complex numbers
3437 for (long i = 1; i < dist; ++ i)
3438 {
3439 df [c1_r + i] = sf [c1_r + i] + sf [c2_r - i];
3440 df [c1_i + i] = sf [c2_r + i] - sf [cend - i];
3441
3442 DataType c;
3443 DataType s;
3444 FFTRealUseTrigo <TRIGO_DIRECT>::iterate (
3445 osc_list [TRIGO_OSC],
3446 c,
3447 s,
3448 cos_ptr,
3449 i * table_step,
3450 (dist - i) * table_step
3451 );
3452
3453 const DataType vr = sf [c1_r + i] - sf [c2_r - i];
3454 const DataType vi = sf [c2_r + i] + sf [cend - i];
3455
3456 df [c2_r + i] = vr * c + vi * s;
3457 df [c2_i + i] = vi * c - vr * s;
3458 }
3459
3460 coef_index += cend;
3461 }
3462 while (coef_index < len);
3463}
3464
3465template <>
3466inline void FFTRealPassInverse <2>::process_internal (long len, DataType dest_ptr [], const DataType src_ptr [], const DataType cos_ptr [], long cos_len, const long br_ptr [], OscType osc_list [])
3467{
3468 // Antepenultimate pass
3469 const DataType sqrt2_2 = DataType (SQRT2 * 0.5);
3470
3471 long coef_index = 0;
3472 do
3473 {
3474 dest_ptr [coef_index ] = src_ptr [coef_index] + src_ptr [coef_index + 4];
3475 dest_ptr [coef_index + 4] = src_ptr [coef_index] - src_ptr [coef_index + 4];
3476 dest_ptr [coef_index + 2] = src_ptr [coef_index + 2] * 2;
3477 dest_ptr [coef_index + 6] = src_ptr [coef_index + 6] * 2;
3478
3479 dest_ptr [coef_index + 1] = src_ptr [coef_index + 1] + src_ptr [coef_index + 3];
3480 dest_ptr [coef_index + 3] = src_ptr [coef_index + 5] - src_ptr [coef_index + 7];
3481
3482 const DataType vr = src_ptr [coef_index + 1] - src_ptr [coef_index + 3];
3483 const DataType vi = src_ptr [coef_index + 5] + src_ptr [coef_index + 7];
3484
3485 dest_ptr [coef_index + 5] = (vr + vi) * sqrt2_2;
3486 dest_ptr [coef_index + 7] = (vi - vr) * sqrt2_2;
3487
3488 coef_index += 8;
3489 }
3490 while (coef_index < len);
3491}
3492
3493template <>
3494inline void FFTRealPassInverse <1>::process_internal (long len, DataType dest_ptr [], const DataType src_ptr [], const DataType cos_ptr [], long cos_len, const long br_ptr [], OscType osc_list [])
3495{
3496 // Penultimate and last pass at once
3497 const long qlen = len >> 2;
3498
3499 long coef_index = 0;
3500 do
3501 {
3502 const long ri_0 = br_ptr [coef_index >> 2];
3503
3504 const DataType b_0 = src_ptr [coef_index ] + src_ptr [coef_index + 2];
3505 const DataType b_2 = src_ptr [coef_index ] - src_ptr [coef_index + 2];
3506 const DataType b_1 = src_ptr [coef_index + 1] * 2;
3507 const DataType b_3 = src_ptr [coef_index + 3] * 2;
3508
3509 dest_ptr [ri_0 ] = b_0 + b_1;
3510 dest_ptr [ri_0 + 2 * qlen] = b_0 - b_1;
3511 dest_ptr [ri_0 + 1 * qlen] = b_2 + b_3;
3512 dest_ptr [ri_0 + 3 * qlen] = b_2 - b_3;
3513
3514 coef_index += 4;
3515 }
3516 while (coef_index < len);
3517}
3518
3519
3520
3521/*\\ PROTECTED \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
3522
3523
3524
3525/*\\ PRIVATE \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
3526
3527
3528
3529} // namespace ffft
3530
3531
3532
3533#endif // ffft_FFTRealPassInverse_CODEHEADER_INCLUDED
3534
3535#undef ffft_FFTRealPassInverse_CURRENT_CODEHEADER
3536
3537
3538
3539/*\\ EOF \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
3540
3541
3542
3543
3544#endif // ffft_FFTRealPassInverse_HEADER_INCLUDED
3545
3546
3547
3548/*\\ EOF \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
3549
3550/*****************************************************************************
3551
3552 FFTRealSelect.h
3553 By Laurent de Soras
3554
3555--- Legal stuff ---
3556
3557This program is free software. It comes without any warranty, to
3558the extent permitted by applicable law. You can redistribute it
3559and/or modify it under the terms of the Do What The Fuck You Want
3560To Public License, Version 2, as published by Sam Hocevar. See
3561http://sam.zoy.org/wtfpl/COPYING for more details.
3562
3563*Tab=3***********************************************************************/
3564
3565
3566
3567#if ! defined (ffft_FFTRealSelect_HEADER_INCLUDED)
3568#define ffft_FFTRealSelect_HEADER_INCLUDED
3569
3570#if defined (_MSC_VER)
3571 #pragma once
3572#endif
3573
3574
3575
3576/*\\ INCLUDE FILES \\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
3577
3578
3579
3580
3581
3582namespace ffft
3583{
3584
3585
3586
3587template <int P>
3588class FFTRealSelect
3589{
3590
3591/*\\ PUBLIC \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
3592
3593public:
3594
3595 ffft_FORCEINLINE static float *
3596 sel_bin (float *e_ptr, float *o_ptr);
3597
3598
3599
3600/*\\ FORBIDDEN MEMBER FUNCTIONS \\\\\\\\\\\\\\\\\\\\\\*/
3601
3602private:
3603
3604 FFTRealSelect ();
3605 ~FFTRealSelect ();
3606 FFTRealSelect (const FFTRealSelect &other);
3607 FFTRealSelect& operator = (const FFTRealSelect &other);
3608 bool operator == (const FFTRealSelect &other);
3609 bool operator != (const FFTRealSelect &other);
3610
3611}; // class FFTRealSelect
3612
3613
3614
3615} // namespace ffft
3616
3617
3618
3619/*****************************************************************************
3620
3621 FFTRealSelect.hpp
3622 By Laurent de Soras
3623
3624--- Legal stuff ---
3625
3626This program is free software. It comes without any warranty, to
3627the extent permitted by applicable law. You can redistribute it
3628and/or modify it under the terms of the Do What The Fuck You Want
3629To Public License, Version 2, as published by Sam Hocevar. See
3630http://sam.zoy.org/wtfpl/COPYING for more details.
3631
3632*Tab=3***********************************************************************/
3633
3634
3635
3636#if defined (ffft_FFTRealSelect_CURRENT_CODEHEADER)
3637 #error Recursive inclusion of FFTRealSelect code header.
3638#endif
3639#define ffft_FFTRealSelect_CURRENT_CODEHEADER
3640
3641#if ! defined (ffft_FFTRealSelect_CODEHEADER_INCLUDED)
3642#define ffft_FFTRealSelect_CODEHEADER_INCLUDED
3643
3644
3645
3646namespace ffft
3647{
3648
3649
3650
3651/*\\ PUBLIC \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
3652
3653
3654
3655template <int P>
3656float * FFTRealSelect <P>::sel_bin (float *e_ptr, float *o_ptr)
3657{
3658 return (o_ptr);
3659}
3660
3661
3662
3663template <>
3664inline float * FFTRealSelect <0>::sel_bin (float *e_ptr, float *o_ptr)
3665{
3666 return (e_ptr);
3667}
3668
3669
3670
3671} // namespace ffft
3672
3673
3674
3675#endif // ffft_FFTRealSelect_CODEHEADER_INCLUDED
3676
3677#undef ffft_FFTRealSelect_CURRENT_CODEHEADER
3678
3679
3680
3681/*\\ EOF \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
3682
3683
3684
3685
3686#endif // ffft_FFTRealSelect_HEADER_INCLUDED
3687
3688
3689
3690/*\\ EOF \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
3691
3692
3693namespace std { }
3694
3695
3696
3697namespace ffft
3698{
3699
3700
3701
3702/*\\ PUBLIC \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
3703
3704
3705
3706template <int LL2>
3707FFTRealFixLen <LL2>::FFTRealFixLen ()
3708: _buffer (FFT_LEN)
3709, _br_data (BR_ARR_SIZE)
3710, _trigo_data (TRIGO_TABLE_ARR_SIZE)
3711, _trigo_osc ()
3712{
3713 build_br_lut ();
3714 build_trigo_lut ();
3715 build_trigo_osc ();
3716}
3717
3718
3719
3720template <int LL2>
3721long FFTRealFixLen <LL2>::get_length () const
3722{
3723 return (FFT_LEN);
3724}
3725
3726
3727
3728// General case
3729template <int LL2>
3730void FFTRealFixLen <LL2>::do_fft (DataType f [], const DataType x [])
3731{
3732 assert (f != 0);
3733 assert (x != 0);
3734 assert (x != f);
3735 assert (FFT_LEN_L2 >= 3);
3736
3737 // Do the transform in several passes
3738 const DataType * cos_ptr = &_trigo_data [0];
3739 const long * br_ptr = &_br_data [0];
3740
3741 FFTRealPassDirect <FFT_LEN_L2 - 1>::process (
3742 FFT_LEN,
3743 f,
3744 &_buffer [0],
3745 x,
3746 cos_ptr,
3747 TRIGO_TABLE_ARR_SIZE,
3748 br_ptr,
3749 &_trigo_osc [0]
3750 );
3751}
3752
3753// 4-point FFT
3754template <>
3755inline void FFTRealFixLen <2>::do_fft (DataType f [], const DataType x [])
3756{
3757 assert (f != 0);
3758 assert (x != 0);
3759 assert (x != f);
3760
3761 f [1] = x [0] - x [2];
3762 f [3] = x [1] - x [3];
3763
3764 const DataType b_0 = x [0] + x [2];
3765 const DataType b_2 = x [1] + x [3];
3766
3767 f [0] = b_0 + b_2;
3768 f [2] = b_0 - b_2;
3769}
3770
3771// 2-point FFT
3772template <>
3773inline void FFTRealFixLen <1>::do_fft (DataType f [], const DataType x [])
3774{
3775 assert (f != 0);
3776 assert (x != 0);
3777 assert (x != f);
3778
3779 f [0] = x [0] + x [1];
3780 f [1] = x [0] - x [1];
3781}
3782
3783// 1-point FFT
3784template <>
3785inline void FFTRealFixLen <0>::do_fft (DataType f [], const DataType x [])
3786{
3787 assert (f != 0);
3788 assert (x != 0);
3789
3790 f [0] = x [0];
3791}
3792
3793
3794
3795// General case
3796template <int LL2>
3797void FFTRealFixLen <LL2>::do_ifft (const DataType f [], DataType x [])
3798{
3799 assert (f != 0);
3800 assert (x != 0);
3801 assert (x != f);
3802 assert (FFT_LEN_L2 >= 3);
3803
3804 // Do the transform in several passes
3805 DataType * s_ptr =
3806 FFTRealSelect <FFT_LEN_L2 & 1>::sel_bin (&_buffer [0], x);
3807 DataType * d_ptr =
3808 FFTRealSelect <FFT_LEN_L2 & 1>::sel_bin (x, &_buffer [0]);
3809 const DataType * cos_ptr = &_trigo_data [0];
3810 const long * br_ptr = &_br_data [0];
3811
3812 FFTRealPassInverse <FFT_LEN_L2 - 1>::process (
3813 FFT_LEN,
3814 d_ptr,
3815 s_ptr,
3816 f,
3817 cos_ptr,
3818 TRIGO_TABLE_ARR_SIZE,
3819 br_ptr,
3820 &_trigo_osc [0]
3821 );
3822}
3823
3824// 4-point IFFT
3825template <>
3826inline void FFTRealFixLen <2>::do_ifft (const DataType f [], DataType x [])
3827{
3828 assert (f != 0);
3829 assert (x != 0);
3830 assert (x != f);
3831
3832 const DataType b_0 = f [0] + f [2];
3833 const DataType b_2 = f [0] - f [2];
3834
3835 x [0] = b_0 + f [1] * 2;
3836 x [2] = b_0 - f [1] * 2;
3837 x [1] = b_2 + f [3] * 2;
3838 x [3] = b_2 - f [3] * 2;
3839}
3840
3841// 2-point IFFT
3842template <>
3843inline void FFTRealFixLen <1>::do_ifft (const DataType f [], DataType x [])
3844{
3845 assert (f != 0);
3846 assert (x != 0);
3847 assert (x != f);
3848
3849 x [0] = f [0] + f [1];
3850 x [1] = f [0] - f [1];
3851}
3852
3853// 1-point IFFT
3854template <>
3855inline void FFTRealFixLen <0>::do_ifft (const DataType f [], DataType x [])
3856{
3857 assert (f != 0);
3858 assert (x != 0);
3859 assert (x != f);
3860
3861 x [0] = f [0];
3862}
3863
3864
3865
3866
3867template <int LL2>
3868void FFTRealFixLen <LL2>::rescale (DataType x []) const
3869{
3870 assert (x != 0);
3871
3872 const DataType mul = DataType (1.0 / FFT_LEN);
3873
3874 if (FFT_LEN < 4)
3875 {
3876 long i = FFT_LEN - 1;
3877 do
3878 {
3879 x [i] *= mul;
3880 --i;
3881 }
3882 while (i >= 0);
3883 }
3884
3885 else
3886 {
3887 assert ((FFT_LEN & 3) == 0);
3888
3889 // Could be optimized with SIMD instruction sets (needs alignment check)
3890 long i = FFT_LEN - 4;
3891 do
3892 {
3893 x [i + 0] *= mul;
3894 x [i + 1] *= mul;
3895 x [i + 2] *= mul;
3896 x [i + 3] *= mul;
3897 i -= 4;
3898 }
3899 while (i >= 0);
3900 }
3901}
3902
3903
3904
3905/*\\ PROTECTED \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
3906
3907
3908
3909/*\\ PRIVATE \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
3910
3911
3912
3913template <int LL2>
3914void FFTRealFixLen <LL2>::build_br_lut ()
3915{
3916 _br_data [0] = 0;
3917 for (long cnt = 1; cnt < BR_ARR_SIZE; ++cnt)
3918 {
3919 long index = cnt << 2;
3920 long br_index = 0;
3921
3922 int bit_cnt = FFT_LEN_L2;
3923 do
3924 {
3925 br_index <<= 1;
3926 br_index += (index & 1);
3927 index >>= 1;
3928
3929 -- bit_cnt;
3930 }
3931 while (bit_cnt > 0);
3932
3933 _br_data [cnt] = br_index;
3934 }
3935}
3936
3937
3938
3939template <int LL2>
3940void FFTRealFixLen <LL2>::build_trigo_lut ()
3941{
3942 const double mul = (0.5 * PI) / TRIGO_TABLE_ARR_SIZE;
3943 for (long i = 0; i < TRIGO_TABLE_ARR_SIZE; ++ i)
3944 {
3945 using namespace std;
3946
3947 _trigo_data [i] = DataType (cos (i * mul));
3948 }
3949}
3950
3951
3952
3953template <int LL2>
3954void FFTRealFixLen <LL2>::build_trigo_osc ()
3955{
3956 for (int i = 0; i < NBR_TRIGO_OSC; ++i)
3957 {
3958 OscType & osc = _trigo_osc [i];
3959
3960 const long len = static_cast <long> (TRIGO_TABLE_ARR_SIZE) << (i + 1);
3961 const double mul = (0.5 * PI) / len;
3962 osc.set_step (mul);
3963 }
3964}
3965
3966
3967
3968} // namespace ffft
3969
3970
3971
3972#endif // ffft_FFTRealFixLen_CODEHEADER_INCLUDED
3973
3974#undef ffft_FFTRealFixLen_CURRENT_CODEHEADER
3975
3976
3977
3978/*\\ EOF \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
3979
3980
3981
3982
3983#endif // ffft_FFTRealFixLen_HEADER_INCLUDED
3984
3985
3986
3987/*\\ EOF \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
ffft::Array
Definition FFTReal.h:1679
ffft::DynArray
Definition FFTReal.h:133
ffft::FFTRealFixLen
Definition FFTReal.h:2425
ffft::FFTRealFixLenParam
Definition FFTReal.h:2130
ffft::FFTReal
Definition FFTReal.h:600
ffft::FFTRealPassDirect
Definition FFTReal.h:2576
ffft::FFTRealPassInverse
Definition FFTReal.h:3070
ffft::FFTRealSelect
Definition FFTReal.h:3589
ffft::FFTRealUseTrigo
Definition FFTReal.h:2692
ffft::OscSinCos
Definition FFTReal.h:396
ffft
AudioTools internal: FFT.
Definition FFTReal.h:65
Generated by doxygen 1.9.8