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>
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[])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>
process_rec(long len,DataType dest_ptr[],DataType src_ptr[],const DataType cos_ptr[],long cos_len,const long br_ptr[],OscType osc_list[])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 <>
process_rec(long len,DataType dest_ptr[],DataType src_ptr[],const DataType cos_ptr[],long cos_len,const long br_ptr[],OscType osc_list[])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>
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[])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 <>
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[])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 <>
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[])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