|
1 /***************************************************************************** |
|
2 |
|
3 FFTRealFixLen.hpp |
|
4 Copyright (c) 2005 Laurent de Soras |
|
5 |
|
6 --- Legal stuff --- |
|
7 |
|
8 This library is free software; you can redistribute it and/or |
|
9 modify it under the terms of the GNU Lesser General Public |
|
10 License as published by the Free Software Foundation; either |
|
11 version 2.1 of the License, or (at your option) any later version. |
|
12 |
|
13 This library is distributed in the hope that it will be useful, |
|
14 but WITHOUT ANY WARRANTY; without even the implied warranty of |
|
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
|
16 Lesser General Public License for more details. |
|
17 |
|
18 You should have received a copy of the GNU Lesser General Public |
|
19 License along with this library; if not, write to the Free Software |
|
20 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA |
|
21 |
|
22 *Tab=3***********************************************************************/ |
|
23 |
|
24 |
|
25 |
|
26 #if defined (FFTRealFixLen_CURRENT_CODEHEADER) |
|
27 #error Recursive inclusion of FFTRealFixLen code header. |
|
28 #endif |
|
29 #define FFTRealFixLen_CURRENT_CODEHEADER |
|
30 |
|
31 #if ! defined (FFTRealFixLen_CODEHEADER_INCLUDED) |
|
32 #define FFTRealFixLen_CODEHEADER_INCLUDED |
|
33 |
|
34 |
|
35 |
|
36 /*\\\ INCLUDE FILES \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/ |
|
37 |
|
38 #include "def.h" |
|
39 #include "FFTRealPassDirect.h" |
|
40 #include "FFTRealPassInverse.h" |
|
41 #include "FFTRealSelect.h" |
|
42 |
|
43 #include <cassert> |
|
44 #include <cmath> |
|
45 |
|
46 namespace std { } |
|
47 |
|
48 |
|
49 |
|
50 /*\\\ PUBLIC \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/ |
|
51 |
|
52 |
|
53 |
|
54 template <int LL2> |
|
55 FFTRealFixLen <LL2>::FFTRealFixLen () |
|
56 : _buffer (FFT_LEN) |
|
57 , _br_data (BR_ARR_SIZE) |
|
58 , _trigo_data (TRIGO_TABLE_ARR_SIZE) |
|
59 , _trigo_osc () |
|
60 { |
|
61 build_br_lut (); |
|
62 build_trigo_lut (); |
|
63 build_trigo_osc (); |
|
64 } |
|
65 |
|
66 |
|
67 |
|
68 template <int LL2> |
|
69 long FFTRealFixLen <LL2>::get_length () const |
|
70 { |
|
71 return (FFT_LEN); |
|
72 } |
|
73 |
|
74 |
|
75 |
|
76 // General case |
|
77 template <int LL2> |
|
78 void FFTRealFixLen <LL2>::do_fft (DataType f [], const DataType x []) |
|
79 { |
|
80 assert (f != 0); |
|
81 assert (x != 0); |
|
82 assert (x != f); |
|
83 assert (FFT_LEN_L2 >= 3); |
|
84 |
|
85 // Do the transform in several passes |
|
86 const DataType * cos_ptr = &_trigo_data [0]; |
|
87 const long * br_ptr = &_br_data [0]; |
|
88 |
|
89 FFTRealPassDirect <FFT_LEN_L2 - 1>::process ( |
|
90 FFT_LEN, |
|
91 f, |
|
92 &_buffer [0], |
|
93 x, |
|
94 cos_ptr, |
|
95 TRIGO_TABLE_ARR_SIZE, |
|
96 br_ptr, |
|
97 &_trigo_osc [0] |
|
98 ); |
|
99 } |
|
100 |
|
101 // 4-point FFT |
|
102 template <> |
|
103 void FFTRealFixLen <2>::do_fft (DataType f [], const DataType x []) |
|
104 { |
|
105 assert (f != 0); |
|
106 assert (x != 0); |
|
107 assert (x != f); |
|
108 |
|
109 f [1] = x [0] - x [2]; |
|
110 f [3] = x [1] - x [3]; |
|
111 |
|
112 const DataType b_0 = x [0] + x [2]; |
|
113 const DataType b_2 = x [1] + x [3]; |
|
114 |
|
115 f [0] = b_0 + b_2; |
|
116 f [2] = b_0 - b_2; |
|
117 } |
|
118 |
|
119 // 2-point FFT |
|
120 template <> |
|
121 void FFTRealFixLen <1>::do_fft (DataType f [], const DataType x []) |
|
122 { |
|
123 assert (f != 0); |
|
124 assert (x != 0); |
|
125 assert (x != f); |
|
126 |
|
127 f [0] = x [0] + x [1]; |
|
128 f [1] = x [0] - x [1]; |
|
129 } |
|
130 |
|
131 // 1-point FFT |
|
132 template <> |
|
133 void FFTRealFixLen <0>::do_fft (DataType f [], const DataType x []) |
|
134 { |
|
135 assert (f != 0); |
|
136 assert (x != 0); |
|
137 |
|
138 f [0] = x [0]; |
|
139 } |
|
140 |
|
141 |
|
142 |
|
143 // General case |
|
144 template <int LL2> |
|
145 void FFTRealFixLen <LL2>::do_ifft (const DataType f [], DataType x []) |
|
146 { |
|
147 assert (f != 0); |
|
148 assert (x != 0); |
|
149 assert (x != f); |
|
150 assert (FFT_LEN_L2 >= 3); |
|
151 |
|
152 // Do the transform in several passes |
|
153 DataType * s_ptr = |
|
154 FFTRealSelect <FFT_LEN_L2 & 1>::sel_bin (&_buffer [0], x); |
|
155 DataType * d_ptr = |
|
156 FFTRealSelect <FFT_LEN_L2 & 1>::sel_bin (x, &_buffer [0]); |
|
157 const DataType * cos_ptr = &_trigo_data [0]; |
|
158 const long * br_ptr = &_br_data [0]; |
|
159 |
|
160 FFTRealPassInverse <FFT_LEN_L2 - 1>::process ( |
|
161 FFT_LEN, |
|
162 d_ptr, |
|
163 s_ptr, |
|
164 f, |
|
165 cos_ptr, |
|
166 TRIGO_TABLE_ARR_SIZE, |
|
167 br_ptr, |
|
168 &_trigo_osc [0] |
|
169 ); |
|
170 } |
|
171 |
|
172 // 4-point IFFT |
|
173 template <> |
|
174 void FFTRealFixLen <2>::do_ifft (const DataType f [], DataType x []) |
|
175 { |
|
176 assert (f != 0); |
|
177 assert (x != 0); |
|
178 assert (x != f); |
|
179 |
|
180 const DataType b_0 = f [0] + f [2]; |
|
181 const DataType b_2 = f [0] - f [2]; |
|
182 |
|
183 x [0] = b_0 + f [1] * 2; |
|
184 x [2] = b_0 - f [1] * 2; |
|
185 x [1] = b_2 + f [3] * 2; |
|
186 x [3] = b_2 - f [3] * 2; |
|
187 } |
|
188 |
|
189 // 2-point IFFT |
|
190 template <> |
|
191 void FFTRealFixLen <1>::do_ifft (const DataType f [], DataType x []) |
|
192 { |
|
193 assert (f != 0); |
|
194 assert (x != 0); |
|
195 assert (x != f); |
|
196 |
|
197 x [0] = f [0] + f [1]; |
|
198 x [1] = f [0] - f [1]; |
|
199 } |
|
200 |
|
201 // 1-point IFFT |
|
202 template <> |
|
203 void FFTRealFixLen <0>::do_ifft (const DataType f [], DataType x []) |
|
204 { |
|
205 assert (f != 0); |
|
206 assert (x != 0); |
|
207 assert (x != f); |
|
208 |
|
209 x [0] = f [0]; |
|
210 } |
|
211 |
|
212 |
|
213 |
|
214 |
|
215 template <int LL2> |
|
216 void FFTRealFixLen <LL2>::rescale (DataType x []) const |
|
217 { |
|
218 assert (x != 0); |
|
219 |
|
220 const DataType mul = DataType (1.0 / FFT_LEN); |
|
221 |
|
222 if (FFT_LEN < 4) |
|
223 { |
|
224 long i = FFT_LEN - 1; |
|
225 do |
|
226 { |
|
227 x [i] *= mul; |
|
228 --i; |
|
229 } |
|
230 while (i >= 0); |
|
231 } |
|
232 |
|
233 else |
|
234 { |
|
235 assert ((FFT_LEN & 3) == 0); |
|
236 |
|
237 // Could be optimized with SIMD instruction sets (needs alignment check) |
|
238 long i = FFT_LEN - 4; |
|
239 do |
|
240 { |
|
241 x [i + 0] *= mul; |
|
242 x [i + 1] *= mul; |
|
243 x [i + 2] *= mul; |
|
244 x [i + 3] *= mul; |
|
245 i -= 4; |
|
246 } |
|
247 while (i >= 0); |
|
248 } |
|
249 } |
|
250 |
|
251 |
|
252 |
|
253 /*\\\ PROTECTED \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/ |
|
254 |
|
255 |
|
256 |
|
257 /*\\\ PRIVATE \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/ |
|
258 |
|
259 |
|
260 |
|
261 template <int LL2> |
|
262 void FFTRealFixLen <LL2>::build_br_lut () |
|
263 { |
|
264 _br_data [0] = 0; |
|
265 for (long cnt = 1; cnt < BR_ARR_SIZE; ++cnt) |
|
266 { |
|
267 long index = cnt << 2; |
|
268 long br_index = 0; |
|
269 |
|
270 int bit_cnt = FFT_LEN_L2; |
|
271 do |
|
272 { |
|
273 br_index <<= 1; |
|
274 br_index += (index & 1); |
|
275 index >>= 1; |
|
276 |
|
277 -- bit_cnt; |
|
278 } |
|
279 while (bit_cnt > 0); |
|
280 |
|
281 _br_data [cnt] = br_index; |
|
282 } |
|
283 } |
|
284 |
|
285 |
|
286 |
|
287 template <int LL2> |
|
288 void FFTRealFixLen <LL2>::build_trigo_lut () |
|
289 { |
|
290 const double mul = (0.5 * PI) / TRIGO_TABLE_ARR_SIZE; |
|
291 for (long i = 0; i < TRIGO_TABLE_ARR_SIZE; ++ i) |
|
292 { |
|
293 using namespace std; |
|
294 |
|
295 _trigo_data [i] = DataType (cos (i * mul)); |
|
296 } |
|
297 } |
|
298 |
|
299 |
|
300 |
|
301 template <int LL2> |
|
302 void FFTRealFixLen <LL2>::build_trigo_osc () |
|
303 { |
|
304 for (int i = 0; i < NBR_TRIGO_OSC; ++i) |
|
305 { |
|
306 OscType & osc = _trigo_osc [i]; |
|
307 |
|
308 const long len = static_cast <long> (TRIGO_TABLE_ARR_SIZE) << (i + 1); |
|
309 const double mul = (0.5 * PI) / len; |
|
310 osc.set_step (mul); |
|
311 } |
|
312 } |
|
313 |
|
314 |
|
315 |
|
316 #endif // FFTRealFixLen_CODEHEADER_INCLUDED |
|
317 |
|
318 #undef FFTRealFixLen_CURRENT_CODEHEADER |
|
319 |
|
320 |
|
321 |
|
322 /*\\\ EOF \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/ |