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 |
30 | |
31 | #if ! defined (FFTRealPassDirect_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 <> |
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 | dest_ptr: src_ptr, |
83 | src_ptr: 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 \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/ |
205 | |