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 \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
|
---|