1 /*****************************************************************************
2
3 FFTRealPassInverse.hpp
4 By Laurent de Soras
5
6 --- Legal stuff ---
7
8 This program is free software. It comes without any warranty, to
9 the extent permitted by applicable law. You can redistribute it
10 and/or modify it under the terms of the Do What The Fuck You Want
11 To Public License, Version 2, as published by Sam Hocevar. See
12 http://sam.zoy.org/wtfpl/COPYING for more details.
13
14 *Tab=3***********************************************************************/
15
16
17
18 #if defined (ffft_FFTRealPassInverse_CURRENT_CODEHEADER)
19 #error Recursive inclusion of FFTRealPassInverse code header.
20 #endif
21 #define ffft_FFTRealPassInverse_CURRENT_CODEHEADER
22
23 #if ! defined (ffft_FFTRealPassInverse_CODEHEADER_INCLUDED)
24 #define ffft_FFTRealPassInverse_CODEHEADER_INCLUDED
25
26
27
28 /*\\\ INCLUDE FILES \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
29
30 #include "FFTRealUseTrigo.h"
31
32
33
34 namespace ffft
35 {
36
37
38
39 /*\\\ PUBLIC \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
40
41
42
43 template <int 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[])44 void 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 [])
45 {
46 process_internal (
47 len,
48 dest_ptr,
49 f_ptr,
50 cos_ptr,
51 cos_len,
52 br_ptr,
53 osc_list
54 );
55 FFTRealPassInverse <PASS - 1>::process_rec (
56 len,
57 src_ptr,
58 dest_ptr,
59 cos_ptr,
60 cos_len,
61 br_ptr,
62 osc_list
63 );
64 }
65
66
67
68 template <int 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[])69 void 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 [])
70 {
71 process_internal (
72 len,
73 dest_ptr,
74 src_ptr,
75 cos_ptr,
76 cos_len,
77 br_ptr,
78 osc_list
79 );
80 FFTRealPassInverse <PASS - 1>::process_rec (
81 len,
82 src_ptr,
83 dest_ptr,
84 cos_ptr,
85 cos_len,
86 br_ptr,
87 osc_list
88 );
89 }
90
91 template <>
process_rec(long len,DataType dest_ptr[],DataType src_ptr[],const DataType cos_ptr[],long cos_len,const long br_ptr[],OscType osc_list[])92 inline 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 [])
93 {
94 // Stops recursion
95 }
96
97
98
99 template <int 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[])100 void 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 [])
101 {
102 const long dist = 1L << (PASS - 1);
103 const long c1_r = 0;
104 const long c1_i = dist;
105 const long c2_r = dist * 2;
106 const long c2_i = dist * 3;
107 const long cend = dist * 4;
108 const long table_step = cos_len >> (PASS - 1);
109
110 enum { TRIGO_OSC = PASS - FFTRealFixLenParam::TRIGO_BD_LIMIT };
111 enum { TRIGO_DIRECT = (TRIGO_OSC >= 0) ? 1 : 0 };
112
113 long coef_index = 0;
114 do
115 {
116 const DataType * const sf = src_ptr + coef_index;
117 DataType * const df = dest_ptr + coef_index;
118
119 // Extreme coefficients are always real
120 df [c1_r] = sf [c1_r] + sf [c2_r];
121 df [c2_r] = sf [c1_r] - sf [c2_r];
122 df [c1_i] = sf [c1_i] * 2;
123 df [c2_i] = sf [c2_i] * 2;
124
125 FFTRealUseTrigo <TRIGO_DIRECT>::prepare (osc_list [TRIGO_OSC]);
126
127 // Others are conjugate complex numbers
128 for (long i = 1; i < dist; ++ i)
129 {
130 df [c1_r + i] = sf [c1_r + i] + sf [c2_r - i];
131 df [c1_i + i] = sf [c2_r + i] - sf [cend - i];
132
133 DataType c;
134 DataType s;
135 FFTRealUseTrigo <TRIGO_DIRECT>::iterate (
136 osc_list [TRIGO_OSC],
137 c,
138 s,
139 cos_ptr,
140 i * table_step,
141 (dist - i) * table_step
142 );
143
144 const DataType vr = sf [c1_r + i] - sf [c2_r - i];
145 const DataType vi = sf [c2_r + i] + sf [cend - i];
146
147 df [c2_r + i] = vr * c + vi * s;
148 df [c2_i + i] = vi * c - vr * s;
149 }
150
151 coef_index += cend;
152 }
153 while (coef_index < len);
154 }
155
156 template <>
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[])157 inline 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 [])
158 {
159 // Antepenultimate pass
160 const DataType sqrt2_2 = DataType (SQRT2 * 0.5);
161
162 long coef_index = 0;
163 do
164 {
165 dest_ptr [coef_index ] = src_ptr [coef_index] + src_ptr [coef_index + 4];
166 dest_ptr [coef_index + 4] = src_ptr [coef_index] - src_ptr [coef_index + 4];
167 dest_ptr [coef_index + 2] = src_ptr [coef_index + 2] * 2;
168 dest_ptr [coef_index + 6] = src_ptr [coef_index + 6] * 2;
169
170 dest_ptr [coef_index + 1] = src_ptr [coef_index + 1] + src_ptr [coef_index + 3];
171 dest_ptr [coef_index + 3] = src_ptr [coef_index + 5] - src_ptr [coef_index + 7];
172
173 const DataType vr = src_ptr [coef_index + 1] - src_ptr [coef_index + 3];
174 const DataType vi = src_ptr [coef_index + 5] + src_ptr [coef_index + 7];
175
176 dest_ptr [coef_index + 5] = (vr + vi) * sqrt2_2;
177 dest_ptr [coef_index + 7] = (vi - vr) * sqrt2_2;
178
179 coef_index += 8;
180 }
181 while (coef_index < len);
182 }
183
184 template <>
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[])185 inline 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 [])
186 {
187 // Penultimate and last pass at once
188 const long qlen = len >> 2;
189
190 long coef_index = 0;
191 do
192 {
193 const long ri_0 = br_ptr [coef_index >> 2];
194
195 const DataType b_0 = src_ptr [coef_index ] + src_ptr [coef_index + 2];
196 const DataType b_2 = src_ptr [coef_index ] - src_ptr [coef_index + 2];
197 const DataType b_1 = src_ptr [coef_index + 1] * 2;
198 const DataType b_3 = src_ptr [coef_index + 3] * 2;
199
200 dest_ptr [ri_0 ] = b_0 + b_1;
201 dest_ptr [ri_0 + 2 * qlen] = b_0 - b_1;
202 dest_ptr [ri_0 + 1 * qlen] = b_2 + b_3;
203 dest_ptr [ri_0 + 3 * qlen] = b_2 - b_3;
204
205 coef_index += 4;
206 }
207 while (coef_index < len);
208 }
209
210
211
212 /*\\\ PROTECTED \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
213
214
215
216 /*\\\ PRIVATE \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
217
218
219
220 } // namespace ffft
221
222
223
224 #endif // ffft_FFTRealPassInverse_CODEHEADER_INCLUDED
225
226 #undef ffft_FFTRealPassInverse_CURRENT_CODEHEADER
227
228
229
230 /*\\\ EOF \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
231