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 |
30 | |
31 | #if ! defined (FFTRealPassInverse_CODEHEADER_INCLUDED) |
32 | #define |
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 | src_ptr: 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 \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/ |
230 | |