|
1 /***************************************************************************** |
|
2 |
|
3 FFTRealPassInverse.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 (FFTRealPassInverse_CURRENT_CODEHEADER) |
|
27 #error Recursive inclusion of FFTRealPassInverse code header. |
|
28 #endif |
|
29 #define FFTRealPassInverse_CURRENT_CODEHEADER |
|
30 |
|
31 #if ! defined (FFTRealPassInverse_CODEHEADER_INCLUDED) |
|
32 #define FFTRealPassInverse_CODEHEADER_INCLUDED |
|
33 |
|
34 |
|
35 |
|
36 /*\\\ INCLUDE FILES \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/ |
|
37 |
|
38 #include "FFTRealUseTrigo.h" |
|
39 |
|
40 |
|
41 |
|
42 /*\\\ PUBLIC \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/ |
|
43 |
|
44 |
|
45 |
|
46 template <int PASS> |
|
47 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 []) |
|
48 { |
|
49 process_internal ( |
|
50 len, |
|
51 dest_ptr, |
|
52 f_ptr, |
|
53 cos_ptr, |
|
54 cos_len, |
|
55 br_ptr, |
|
56 osc_list |
|
57 ); |
|
58 FFTRealPassInverse <PASS - 1>::process_rec ( |
|
59 len, |
|
60 src_ptr, |
|
61 dest_ptr, |
|
62 cos_ptr, |
|
63 cos_len, |
|
64 br_ptr, |
|
65 osc_list |
|
66 ); |
|
67 } |
|
68 |
|
69 |
|
70 |
|
71 template <int PASS> |
|
72 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 []) |
|
73 { |
|
74 process_internal ( |
|
75 len, |
|
76 dest_ptr, |
|
77 src_ptr, |
|
78 cos_ptr, |
|
79 cos_len, |
|
80 br_ptr, |
|
81 osc_list |
|
82 ); |
|
83 FFTRealPassInverse <PASS - 1>::process_rec ( |
|
84 len, |
|
85 src_ptr, |
|
86 dest_ptr, |
|
87 cos_ptr, |
|
88 cos_len, |
|
89 br_ptr, |
|
90 osc_list |
|
91 ); |
|
92 } |
|
93 |
|
94 template <> |
|
95 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 []) |
|
96 { |
|
97 // Stops recursion |
|
98 } |
|
99 |
|
100 |
|
101 |
|
102 template <int PASS> |
|
103 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 []) |
|
104 { |
|
105 const long dist = 1L << (PASS - 1); |
|
106 const long c1_r = 0; |
|
107 const long c1_i = dist; |
|
108 const long c2_r = dist * 2; |
|
109 const long c2_i = dist * 3; |
|
110 const long cend = dist * 4; |
|
111 const long table_step = cos_len >> (PASS - 1); |
|
112 |
|
113 enum { TRIGO_OSC = PASS - FFTRealFixLenParam::TRIGO_BD_LIMIT }; |
|
114 enum { TRIGO_DIRECT = (TRIGO_OSC >= 0) ? 1 : 0 }; |
|
115 |
|
116 long coef_index = 0; |
|
117 do |
|
118 { |
|
119 const DataType * const sf = src_ptr + coef_index; |
|
120 DataType * const df = dest_ptr + coef_index; |
|
121 |
|
122 // Extreme coefficients are always real |
|
123 df [c1_r] = sf [c1_r] + sf [c2_r]; |
|
124 df [c2_r] = sf [c1_r] - sf [c2_r]; |
|
125 df [c1_i] = sf [c1_i] * 2; |
|
126 df [c2_i] = sf [c2_i] * 2; |
|
127 |
|
128 FFTRealUseTrigo <TRIGO_DIRECT>::prepare (osc_list [TRIGO_OSC]); |
|
129 |
|
130 // Others are conjugate complex numbers |
|
131 for (long i = 1; i < dist; ++ i) |
|
132 { |
|
133 df [c1_r + i] = sf [c1_r + i] + sf [c2_r - i]; |
|
134 df [c1_i + i] = sf [c2_r + i] - sf [cend - i]; |
|
135 |
|
136 DataType c; |
|
137 DataType s; |
|
138 FFTRealUseTrigo <TRIGO_DIRECT>::iterate ( |
|
139 osc_list [TRIGO_OSC], |
|
140 c, |
|
141 s, |
|
142 cos_ptr, |
|
143 i * table_step, |
|
144 (dist - i) * table_step |
|
145 ); |
|
146 |
|
147 const DataType vr = sf [c1_r + i] - sf [c2_r - i]; |
|
148 const DataType vi = sf [c2_r + i] + sf [cend - i]; |
|
149 |
|
150 df [c2_r + i] = vr * c + vi * s; |
|
151 df [c2_i + i] = vi * c - vr * s; |
|
152 } |
|
153 |
|
154 coef_index += cend; |
|
155 } |
|
156 while (coef_index < len); |
|
157 } |
|
158 |
|
159 template <> |
|
160 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 []) |
|
161 { |
|
162 // Antepenultimate pass |
|
163 const DataType sqrt2_2 = DataType (SQRT2 * 0.5); |
|
164 |
|
165 long coef_index = 0; |
|
166 do |
|
167 { |
|
168 dest_ptr [coef_index ] = src_ptr [coef_index] + src_ptr [coef_index + 4]; |
|
169 dest_ptr [coef_index + 4] = src_ptr [coef_index] - src_ptr [coef_index + 4]; |
|
170 dest_ptr [coef_index + 2] = src_ptr [coef_index + 2] * 2; |
|
171 dest_ptr [coef_index + 6] = src_ptr [coef_index + 6] * 2; |
|
172 |
|
173 dest_ptr [coef_index + 1] = src_ptr [coef_index + 1] + src_ptr [coef_index + 3]; |
|
174 dest_ptr [coef_index + 3] = src_ptr [coef_index + 5] - src_ptr [coef_index + 7]; |
|
175 |
|
176 const DataType vr = src_ptr [coef_index + 1] - src_ptr [coef_index + 3]; |
|
177 const DataType vi = src_ptr [coef_index + 5] + src_ptr [coef_index + 7]; |
|
178 |
|
179 dest_ptr [coef_index + 5] = (vr + vi) * sqrt2_2; |
|
180 dest_ptr [coef_index + 7] = (vi - vr) * sqrt2_2; |
|
181 |
|
182 coef_index += 8; |
|
183 } |
|
184 while (coef_index < len); |
|
185 } |
|
186 |
|
187 template <> |
|
188 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 []) |
|
189 { |
|
190 // Penultimate and last pass at once |
|
191 const long qlen = len >> 2; |
|
192 |
|
193 long coef_index = 0; |
|
194 do |
|
195 { |
|
196 const long ri_0 = br_ptr [coef_index >> 2]; |
|
197 |
|
198 const DataType b_0 = src_ptr [coef_index ] + src_ptr [coef_index + 2]; |
|
199 const DataType b_2 = src_ptr [coef_index ] - src_ptr [coef_index + 2]; |
|
200 const DataType b_1 = src_ptr [coef_index + 1] * 2; |
|
201 const DataType b_3 = src_ptr [coef_index + 3] * 2; |
|
202 |
|
203 dest_ptr [ri_0 ] = b_0 + b_1; |
|
204 dest_ptr [ri_0 + 2 * qlen] = b_0 - b_1; |
|
205 dest_ptr [ri_0 + 1 * qlen] = b_2 + b_3; |
|
206 dest_ptr [ri_0 + 3 * qlen] = b_2 - b_3; |
|
207 |
|
208 coef_index += 4; |
|
209 } |
|
210 while (coef_index < len); |
|
211 } |
|
212 |
|
213 |
|
214 |
|
215 /*\\\ PROTECTED \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/ |
|
216 |
|
217 |
|
218 |
|
219 /*\\\ PRIVATE \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/ |
|
220 |
|
221 |
|
222 |
|
223 #endif // FFTRealPassInverse_CODEHEADER_INCLUDED |
|
224 |
|
225 #undef FFTRealPassInverse_CURRENT_CODEHEADER |
|
226 |
|
227 |
|
228 |
|
229 /*\\\ EOF \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/ |