|
1 /***************************************************************************** |
|
2 |
|
3 FFTRealPassDirect.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 (FFTRealPassDirect_CURRENT_CODEHEADER) |
|
27 #error Recursive inclusion of FFTRealPassDirect code header. |
|
28 #endif |
|
29 #define FFTRealPassDirect_CURRENT_CODEHEADER |
|
30 |
|
31 #if ! defined (FFTRealPassDirect_CODEHEADER_INCLUDED) |
|
32 #define FFTRealPassDirect_CODEHEADER_INCLUDED |
|
33 |
|
34 |
|
35 |
|
36 /*\\\ INCLUDE FILES \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/ |
|
37 |
|
38 #include "FFTRealUseTrigo.h" |
|
39 |
|
40 |
|
41 |
|
42 /*\\\ PUBLIC \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/ |
|
43 |
|
44 |
|
45 |
|
46 template <> |
|
47 void FFTRealPassDirect <1>::process (long len, DataType dest_ptr [], DataType src_ptr [], const DataType x_ptr [], const DataType cos_ptr [], long cos_len, const long br_ptr [], OscType osc_list []) |
|
48 { |
|
49 // First and second pass at once |
|
50 const long qlen = len >> 2; |
|
51 |
|
52 long coef_index = 0; |
|
53 do |
|
54 { |
|
55 // To do: unroll the loop (2x). |
|
56 const long ri_0 = br_ptr [coef_index >> 2]; |
|
57 const long ri_1 = ri_0 + 2 * qlen; // bit_rev_lut_ptr [coef_index + 1]; |
|
58 const long ri_2 = ri_0 + 1 * qlen; // bit_rev_lut_ptr [coef_index + 2]; |
|
59 const long ri_3 = ri_0 + 3 * qlen; // bit_rev_lut_ptr [coef_index + 3]; |
|
60 |
|
61 DataType * const df2 = dest_ptr + coef_index; |
|
62 df2 [1] = x_ptr [ri_0] - x_ptr [ri_1]; |
|
63 df2 [3] = x_ptr [ri_2] - x_ptr [ri_3]; |
|
64 |
|
65 const DataType sf_0 = x_ptr [ri_0] + x_ptr [ri_1]; |
|
66 const DataType sf_2 = x_ptr [ri_2] + x_ptr [ri_3]; |
|
67 |
|
68 df2 [0] = sf_0 + sf_2; |
|
69 df2 [2] = sf_0 - sf_2; |
|
70 |
|
71 coef_index += 4; |
|
72 } |
|
73 while (coef_index < len); |
|
74 } |
|
75 |
|
76 template <> |
|
77 void FFTRealPassDirect <2>::process (long len, DataType dest_ptr [], DataType src_ptr [], const DataType x_ptr [], const DataType cos_ptr [], long cos_len, const long br_ptr [], OscType osc_list []) |
|
78 { |
|
79 // Executes "previous" passes first. Inverts source and destination buffers |
|
80 FFTRealPassDirect <1>::process ( |
|
81 len, |
|
82 src_ptr, |
|
83 dest_ptr, |
|
84 x_ptr, |
|
85 cos_ptr, |
|
86 cos_len, |
|
87 br_ptr, |
|
88 osc_list |
|
89 ); |
|
90 |
|
91 // Third pass |
|
92 const DataType sqrt2_2 = DataType (SQRT2 * 0.5); |
|
93 |
|
94 long coef_index = 0; |
|
95 do |
|
96 { |
|
97 dest_ptr [coef_index ] = src_ptr [coef_index] + src_ptr [coef_index + 4]; |
|
98 dest_ptr [coef_index + 4] = src_ptr [coef_index] - src_ptr [coef_index + 4]; |
|
99 dest_ptr [coef_index + 2] = src_ptr [coef_index + 2]; |
|
100 dest_ptr [coef_index + 6] = src_ptr [coef_index + 6]; |
|
101 |
|
102 DataType v; |
|
103 |
|
104 v = (src_ptr [coef_index + 5] - src_ptr [coef_index + 7]) * sqrt2_2; |
|
105 dest_ptr [coef_index + 1] = src_ptr [coef_index + 1] + v; |
|
106 dest_ptr [coef_index + 3] = src_ptr [coef_index + 1] - v; |
|
107 |
|
108 v = (src_ptr [coef_index + 5] + src_ptr [coef_index + 7]) * sqrt2_2; |
|
109 dest_ptr [coef_index + 5] = v + src_ptr [coef_index + 3]; |
|
110 dest_ptr [coef_index + 7] = v - src_ptr [coef_index + 3]; |
|
111 |
|
112 coef_index += 8; |
|
113 } |
|
114 while (coef_index < len); |
|
115 } |
|
116 |
|
117 template <int PASS> |
|
118 void FFTRealPassDirect <PASS>::process (long len, DataType dest_ptr [], DataType src_ptr [], const DataType x_ptr [], const DataType cos_ptr [], long cos_len, const long br_ptr [], OscType osc_list []) |
|
119 { |
|
120 // Executes "previous" passes first. Inverts source and destination buffers |
|
121 FFTRealPassDirect <PASS - 1>::process ( |
|
122 len, |
|
123 src_ptr, |
|
124 dest_ptr, |
|
125 x_ptr, |
|
126 cos_ptr, |
|
127 cos_len, |
|
128 br_ptr, |
|
129 osc_list |
|
130 ); |
|
131 |
|
132 const long dist = 1L << (PASS - 1); |
|
133 const long c1_r = 0; |
|
134 const long c1_i = dist; |
|
135 const long c2_r = dist * 2; |
|
136 const long c2_i = dist * 3; |
|
137 const long cend = dist * 4; |
|
138 const long table_step = cos_len >> (PASS - 1); |
|
139 |
|
140 enum { TRIGO_OSC = PASS - FFTRealFixLenParam::TRIGO_BD_LIMIT }; |
|
141 enum { TRIGO_DIRECT = (TRIGO_OSC >= 0) ? 1 : 0 }; |
|
142 |
|
143 long coef_index = 0; |
|
144 do |
|
145 { |
|
146 const DataType * const sf = src_ptr + coef_index; |
|
147 DataType * const df = dest_ptr + coef_index; |
|
148 |
|
149 // Extreme coefficients are always real |
|
150 df [c1_r] = sf [c1_r] + sf [c2_r]; |
|
151 df [c2_r] = sf [c1_r] - sf [c2_r]; |
|
152 df [c1_i] = sf [c1_i]; |
|
153 df [c2_i] = sf [c2_i]; |
|
154 |
|
155 FFTRealUseTrigo <TRIGO_DIRECT>::prepare (osc_list [TRIGO_OSC]); |
|
156 |
|
157 // Others are conjugate complex numbers |
|
158 for (long i = 1; i < dist; ++ i) |
|
159 { |
|
160 DataType c; |
|
161 DataType s; |
|
162 FFTRealUseTrigo <TRIGO_DIRECT>::iterate ( |
|
163 osc_list [TRIGO_OSC], |
|
164 c, |
|
165 s, |
|
166 cos_ptr, |
|
167 i * table_step, |
|
168 (dist - i) * table_step |
|
169 ); |
|
170 |
|
171 const DataType sf_r_i = sf [c1_r + i]; |
|
172 const DataType sf_i_i = sf [c1_i + i]; |
|
173 |
|
174 const DataType v1 = sf [c2_r + i] * c - sf [c2_i + i] * s; |
|
175 df [c1_r + i] = sf_r_i + v1; |
|
176 df [c2_r - i] = sf_r_i - v1; |
|
177 |
|
178 const DataType v2 = sf [c2_r + i] * s + sf [c2_i + i] * c; |
|
179 df [c2_r + i] = v2 + sf_i_i; |
|
180 df [cend - i] = v2 - sf_i_i; |
|
181 } |
|
182 |
|
183 coef_index += cend; |
|
184 } |
|
185 while (coef_index < len); |
|
186 } |
|
187 |
|
188 |
|
189 |
|
190 /*\\\ PROTECTED \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/ |
|
191 |
|
192 |
|
193 |
|
194 /*\\\ PRIVATE \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/ |
|
195 |
|
196 |
|
197 |
|
198 #endif // FFTRealPassDirect_CODEHEADER_INCLUDED |
|
199 |
|
200 #undef FFTRealPassDirect_CURRENT_CODEHEADER |
|
201 |
|
202 |
|
203 |
|
204 /*\\\ EOF \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/ |