1 /*****************************************************************************
2
3 FFTRealPassDirect.hpp
4 By Laurent de Soras
5
6 --- Legal stuff ---
7
8 This program is free software. It comes without any warranty, to
9 the extent permitted by applicable law. You can redistribute it
10 and/or modify it under the terms of the Do What The Fuck You Want
11 To Public License, Version 2, as published by Sam Hocevar. See
12 http://sam.zoy.org/wtfpl/COPYING for more details.
13
14 *Tab=3***********************************************************************/
15
16
17
18 #if defined (ffft_FFTRealPassDirect_CURRENT_CODEHEADER)
19 #error Recursive inclusion of FFTRealPassDirect code header.
20 #endif
21 #define ffft_FFTRealPassDirect_CURRENT_CODEHEADER
22
23 #if ! defined (ffft_FFTRealPassDirect_CODEHEADER_INCLUDED)
24 #define ffft_FFTRealPassDirect_CODEHEADER_INCLUDED
25
26
27
28 /*\\\ INCLUDE FILES \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
29
30 #include "FFTRealUseTrigo.h"
31
32
33
34 namespace ffft
35 {
36
37
38
39 /*\\\ PUBLIC \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
40
41
42
43 template <>
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[])44 inline 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 [])
45 {
46 // First and second pass at once
47 const long qlen = len >> 2;
48
49 long coef_index = 0;
50 do
51 {
52 // To do: unroll the loop (2x).
53 const long ri_0 = br_ptr [coef_index >> 2];
54 const long ri_1 = ri_0 + 2 * qlen; // bit_rev_lut_ptr [coef_index + 1];
55 const long ri_2 = ri_0 + 1 * qlen; // bit_rev_lut_ptr [coef_index + 2];
56 const long ri_3 = ri_0 + 3 * qlen; // bit_rev_lut_ptr [coef_index + 3];
57
58 DataType * const df2 = dest_ptr + coef_index;
59 df2 [1] = x_ptr [ri_0] - x_ptr [ri_1];
60 df2 [3] = x_ptr [ri_2] - x_ptr [ri_3];
61
62 const DataType sf_0 = x_ptr [ri_0] + x_ptr [ri_1];
63 const DataType sf_2 = x_ptr [ri_2] + x_ptr [ri_3];
64
65 df2 [0] = sf_0 + sf_2;
66 df2 [2] = sf_0 - sf_2;
67
68 coef_index += 4;
69 }
70 while (coef_index < len);
71 }
72
73 template <>
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[])74 inline 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 [])
75 {
76 // Executes "previous" passes first. Inverts source and destination buffers
77 FFTRealPassDirect <1>::process (
78 len,
79 src_ptr,
80 dest_ptr,
81 x_ptr,
82 cos_ptr,
83 cos_len,
84 br_ptr,
85 osc_list
86 );
87
88 // Third pass
89 const DataType sqrt2_2 = DataType (SQRT2 * 0.5);
90
91 long coef_index = 0;
92 do
93 {
94 dest_ptr [coef_index ] = src_ptr [coef_index] + src_ptr [coef_index + 4];
95 dest_ptr [coef_index + 4] = src_ptr [coef_index] - src_ptr [coef_index + 4];
96 dest_ptr [coef_index + 2] = src_ptr [coef_index + 2];
97 dest_ptr [coef_index + 6] = src_ptr [coef_index + 6];
98
99 DataType v;
100
101 v = (src_ptr [coef_index + 5] - src_ptr [coef_index + 7]) * sqrt2_2;
102 dest_ptr [coef_index + 1] = src_ptr [coef_index + 1] + v;
103 dest_ptr [coef_index + 3] = src_ptr [coef_index + 1] - v;
104
105 v = (src_ptr [coef_index + 5] + src_ptr [coef_index + 7]) * sqrt2_2;
106 dest_ptr [coef_index + 5] = v + src_ptr [coef_index + 3];
107 dest_ptr [coef_index + 7] = v - src_ptr [coef_index + 3];
108
109 coef_index += 8;
110 }
111 while (coef_index < len);
112 }
113
114 template <int 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[])115 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 [])
116 {
117 // Executes "previous" passes first. Inverts source and destination buffers
118 FFTRealPassDirect <PASS - 1>::process (
119 len,
120 src_ptr,
121 dest_ptr,
122 x_ptr,
123 cos_ptr,
124 cos_len,
125 br_ptr,
126 osc_list
127 );
128
129 const long dist = 1L << (PASS - 1);
130 const long c1_r = 0;
131 const long c1_i = dist;
132 const long c2_r = dist * 2;
133 const long c2_i = dist * 3;
134 const long cend = dist * 4;
135 const long table_step = cos_len >> (PASS - 1);
136
137 enum { TRIGO_OSC = PASS - FFTRealFixLenParam::TRIGO_BD_LIMIT };
138 enum { TRIGO_DIRECT = (TRIGO_OSC >= 0) ? 1 : 0 };
139
140 long coef_index = 0;
141 do
142 {
143 const DataType * const sf = src_ptr + coef_index;
144 DataType * const df = dest_ptr + coef_index;
145
146 // Extreme coefficients are always real
147 df [c1_r] = sf [c1_r] + sf [c2_r];
148 df [c2_r] = sf [c1_r] - sf [c2_r];
149 df [c1_i] = sf [c1_i];
150 df [c2_i] = sf [c2_i];
151
152 FFTRealUseTrigo <TRIGO_DIRECT>::prepare (osc_list [TRIGO_OSC]);
153
154 // Others are conjugate complex numbers
155 for (long i = 1; i < dist; ++ i)
156 {
157 DataType c;
158 DataType s;
159 FFTRealUseTrigo <TRIGO_DIRECT>::iterate (
160 osc_list [TRIGO_OSC],
161 c,
162 s,
163 cos_ptr,
164 i * table_step,
165 (dist - i) * table_step
166 );
167
168 const DataType sf_r_i = sf [c1_r + i];
169 const DataType sf_i_i = sf [c1_i + i];
170
171 const DataType v1 = sf [c2_r + i] * c - sf [c2_i + i] * s;
172 df [c1_r + i] = sf_r_i + v1;
173 df [c2_r - i] = sf_r_i - v1;
174
175 const DataType v2 = sf [c2_r + i] * s + sf [c2_i + i] * c;
176 df [c2_r + i] = v2 + sf_i_i;
177 df [cend - i] = v2 - sf_i_i;
178 }
179
180 coef_index += cend;
181 }
182 while (coef_index < len);
183 }
184
185
186
187 /*\\\ PROTECTED \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
188
189
190
191 /*\\\ PRIVATE \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
192
193
194
195 } // namespace ffft
196
197
198
199 #endif // ffft_FFTRealPassDirect_CODEHEADER_INCLUDED
200
201 #undef ffft_FFTRealPassDirect_CURRENT_CODEHEADER
202
203
204
205 /*\\\ EOF \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
206