Mercurial > mplayer.hg
comparison libfaad2/cfft.c @ 12527:4a370c80fe5c
update to the 2.0 release of faad, patch by adland
author | diego |
---|---|
date | Wed, 02 Jun 2004 22:59:04 +0000 |
parents | 3185f64f6350 |
children | d81145997036 |
comparison
equal
deleted
inserted
replaced
12526:e183ad37d24c | 12527:4a370c80fe5c |
---|---|
1 /* | 1 /* |
2 ** FAAD2 - Freeware Advanced Audio (AAC) Decoder including SBR decoding | 2 ** FAAD2 - Freeware Advanced Audio (AAC) Decoder including SBR decoding |
3 ** Copyright (C) 2003 M. Bakker, Ahead Software AG, http://www.nero.com | 3 ** Copyright (C) 2003-2004 M. Bakker, Ahead Software AG, http://www.nero.com |
4 ** | 4 ** |
5 ** This program is free software; you can redistribute it and/or modify | 5 ** This program is free software; you can redistribute it and/or modify |
6 ** it under the terms of the GNU General Public License as published by | 6 ** it under the terms of the GNU General Public License as published by |
7 ** the Free Software Foundation; either version 2 of the License, or | 7 ** the Free Software Foundation; either version 2 of the License, or |
8 ** (at your option) any later version. | 8 ** (at your option) any later version. |
9 ** | 9 ** |
10 ** This program is distributed in the hope that it will be useful, | 10 ** This program is distributed in the hope that it will be useful, |
11 ** but WITHOUT ANY WARRANTY; without even the implied warranty of | 11 ** but WITHOUT ANY WARRANTY; without even the implied warranty of |
12 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | 12 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
13 ** GNU General Public License for more details. | 13 ** GNU General Public License for more details. |
14 ** | 14 ** |
15 ** You should have received a copy of the GNU General Public License | 15 ** You should have received a copy of the GNU General Public License |
16 ** along with this program; if not, write to the Free Software | 16 ** along with this program; if not, write to the Free Software |
17 ** Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. | 17 ** Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. |
18 ** | 18 ** |
19 ** Any non-GPL usage of this software or parts of this software is strictly | 19 ** Any non-GPL usage of this software or parts of this software is strictly |
20 ** forbidden. | 20 ** forbidden. |
21 ** | 21 ** |
22 ** Commercial non-GPL licensing of this software is possible. | 22 ** Commercial non-GPL licensing of this software is possible. |
23 ** For more info contact Ahead Software through Mpeg4AAClicense@nero.com. | 23 ** For more info contact Ahead Software through Mpeg4AAClicense@nero.com. |
24 ** | 24 ** |
25 ** $Id: cfft.c,v 1.1 2003/08/30 22:30:21 arpi Exp $ | 25 ** $Id: cfft.c,v 1.2 2003/10/03 22:22:27 alex Exp $ |
26 **/ | 26 **/ |
27 | 27 |
28 /* | 28 /* |
29 * Algorithmically based on Fortran-77 FFTPACK | 29 * Algorithmically based on Fortran-77 FFTPACK |
30 * by Paul N. Swarztrauber(Version 4, 1985). | 30 * by Paul N. Swarztrauber(Version 4, 1985). |
36 | 36 |
37 #include "common.h" | 37 #include "common.h" |
38 #include "structs.h" | 38 #include "structs.h" |
39 | 39 |
40 #include <stdlib.h> | 40 #include <stdlib.h> |
41 #ifdef _WIN32_WCE | |
42 #define assert(x) | |
43 #else | |
44 #include <assert.h> | |
45 #endif | |
46 | 41 |
47 #include "cfft.h" | 42 #include "cfft.h" |
48 #include "cfft_tab.h" | 43 #include "cfft_tab.h" |
44 | |
45 | |
46 /* static function declarations */ | |
47 #ifdef USE_SSE | |
48 static void passf2pos_sse(const uint16_t l1, const complex_t *cc, | |
49 complex_t *ch, const complex_t *wa); | |
50 static void passf2pos_sse_ido(const uint16_t ido, const uint16_t l1, const complex_t *cc, | |
51 complex_t *ch, const complex_t *wa); | |
52 static void passf4pos_sse_ido(const uint16_t ido, const uint16_t l1, const complex_t *cc, complex_t *ch, | |
53 const complex_t *wa1, const complex_t *wa2, const complex_t *wa3); | |
54 #endif | |
55 static void passf2pos(const uint16_t ido, const uint16_t l1, const complex_t *cc, | |
56 complex_t *ch, const complex_t *wa); | |
57 static void passf2neg(const uint16_t ido, const uint16_t l1, const complex_t *cc, | |
58 complex_t *ch, const complex_t *wa); | |
59 static void passf3(const uint16_t ido, const uint16_t l1, const complex_t *cc, | |
60 complex_t *ch, const complex_t *wa1, const complex_t *wa2, const int8_t isign); | |
61 static void passf4pos(const uint16_t ido, const uint16_t l1, const complex_t *cc, complex_t *ch, | |
62 const complex_t *wa1, const complex_t *wa2, const complex_t *wa3); | |
63 static void passf4neg(const uint16_t ido, const uint16_t l1, const complex_t *cc, complex_t *ch, | |
64 const complex_t *wa1, const complex_t *wa2, const complex_t *wa3); | |
65 static void passf5(const uint16_t ido, const uint16_t l1, const complex_t *cc, complex_t *ch, | |
66 const complex_t *wa1, const complex_t *wa2, const complex_t *wa3, | |
67 const complex_t *wa4, const int8_t isign); | |
68 INLINE void cfftf1(uint16_t n, complex_t *c, complex_t *ch, | |
69 const uint16_t *ifac, const complex_t *wa, const int8_t isign); | |
70 static void cffti1(uint16_t n, complex_t *wa, uint16_t *ifac); | |
49 | 71 |
50 | 72 |
51 /*---------------------------------------------------------------------- | 73 /*---------------------------------------------------------------------- |
52 passf2, passf3, passf4, passf5. Complex FFT passes fwd and bwd. | 74 passf2, passf3, passf4, passf5. Complex FFT passes fwd and bwd. |
53 ----------------------------------------------------------------------*/ | 75 ----------------------------------------------------------------------*/ |
54 | 76 |
55 static void passf2(uint16_t ido, uint16_t l1, complex_t *cc, complex_t *ch, | 77 #if 0 //def USE_SSE |
56 complex_t *wa, int8_t isign) | 78 static void passf2pos_sse(const uint16_t l1, const complex_t *cc, |
79 complex_t *ch, const complex_t *wa) | |
80 { | |
81 uint16_t k, ah, ac; | |
82 | |
83 for (k = 0; k < l1; k++) | |
84 { | |
85 ah = 2*k; | |
86 ac = 4*k; | |
87 | |
88 RE(ch[ah]) = RE(cc[ac]) + RE(cc[ac+1]); | |
89 IM(ch[ah]) = IM(cc[ac]) + IM(cc[ac+1]); | |
90 | |
91 RE(ch[ah+l1]) = RE(cc[ac]) - RE(cc[ac+1]); | |
92 IM(ch[ah+l1]) = IM(cc[ac]) - IM(cc[ac+1]); | |
93 } | |
94 } | |
95 | |
96 static void passf2pos_sse_ido(const uint16_t ido, const uint16_t l1, const complex_t *cc, | |
97 complex_t *ch, const complex_t *wa) | |
98 { | |
99 uint16_t i, k, ah, ac; | |
100 | |
101 for (k = 0; k < l1; k++) | |
102 { | |
103 ah = k*ido; | |
104 ac = 2*k*ido; | |
105 | |
106 for (i = 0; i < ido; i+=4) | |
107 { | |
108 __m128 m1, m2, m3, m4, m5, m6, m7, m8, m9, m10, m11, m12, m13, m14; | |
109 __m128 m15, m16, m17, m18, m19, m20, m21, m22, m23, m24; | |
110 __m128 w1, w2, w3, w4; | |
111 | |
112 m1 = _mm_load_ps(&RE(cc[ac+i])); | |
113 m2 = _mm_load_ps(&RE(cc[ac+ido+i])); | |
114 m5 = _mm_load_ps(&RE(cc[ac+i+2])); | |
115 m6 = _mm_load_ps(&RE(cc[ac+ido+i+2])); | |
116 w1 = _mm_load_ps(&RE(wa[i])); | |
117 w3 = _mm_load_ps(&RE(wa[i+2])); | |
118 | |
119 m3 = _mm_add_ps(m1, m2); | |
120 m15 = _mm_add_ps(m5, m6); | |
121 | |
122 m4 = _mm_sub_ps(m1, m2); | |
123 m16 = _mm_sub_ps(m5, m6); | |
124 | |
125 _mm_store_ps(&RE(ch[ah+i]), m3); | |
126 _mm_store_ps(&RE(ch[ah+i+2]), m15); | |
127 | |
128 | |
129 w2 = _mm_shuffle_ps(w1, w1, _MM_SHUFFLE(2, 3, 0, 1)); | |
130 w4 = _mm_shuffle_ps(w3, w3, _MM_SHUFFLE(2, 3, 0, 1)); | |
131 | |
132 m7 = _mm_mul_ps(m4, w1); | |
133 m17 = _mm_mul_ps(m16, w3); | |
134 m8 = _mm_mul_ps(m4, w2); | |
135 m18 = _mm_mul_ps(m16, w4); | |
136 | |
137 m9 = _mm_shuffle_ps(m7, m8, _MM_SHUFFLE(2, 0, 2, 0)); | |
138 m19 = _mm_shuffle_ps(m17, m18, _MM_SHUFFLE(2, 0, 2, 0)); | |
139 m10 = _mm_shuffle_ps(m7, m8, _MM_SHUFFLE(3, 1, 3, 1)); | |
140 m20 = _mm_shuffle_ps(m17, m18, _MM_SHUFFLE(3, 1, 3, 1)); | |
141 | |
142 m11 = _mm_add_ps(m9, m10); | |
143 m21 = _mm_add_ps(m19, m20); | |
144 m12 = _mm_sub_ps(m9, m10); | |
145 m22 = _mm_sub_ps(m19, m20); | |
146 | |
147 m13 = _mm_shuffle_ps(m11, m11, _MM_SHUFFLE(0, 0, 3, 2)); | |
148 m23 = _mm_shuffle_ps(m21, m21, _MM_SHUFFLE(0, 0, 3, 2)); | |
149 | |
150 m14 = _mm_unpacklo_ps(m12, m13); | |
151 m24 = _mm_unpacklo_ps(m22, m23); | |
152 | |
153 _mm_store_ps(&RE(ch[ah+i+l1*ido]), m14); | |
154 _mm_store_ps(&RE(ch[ah+i+2+l1*ido]), m24); | |
155 } | |
156 } | |
157 } | |
158 #endif | |
159 | |
160 static void passf2pos(const uint16_t ido, const uint16_t l1, const complex_t *cc, | |
161 complex_t *ch, const complex_t *wa) | |
57 { | 162 { |
58 uint16_t i, k, ah, ac; | 163 uint16_t i, k, ah, ac; |
59 | 164 |
60 if (ido == 1) | 165 if (ido == 1) |
61 { | 166 { |
62 for (k = 0; k < l1; k++) | 167 for (k = 0; k < l1; k++) |
63 { | 168 { |
64 ah = 2*k; | 169 ah = 2*k; |
65 ac = 4*k; | 170 ac = 4*k; |
66 | 171 |
67 RE(ch[ah]) = RE(cc[ac]) + RE(cc[ac+1]); | 172 RE(ch[ah]) = RE(cc[ac]) + RE(cc[ac+1]); |
68 RE(ch[ah+l1]) = RE(cc[ac]) - RE(cc[ac+1]); | 173 RE(ch[ah+l1]) = RE(cc[ac]) - RE(cc[ac+1]); |
69 IM(ch[ah]) = IM(cc[ac]) + IM(cc[ac+1]); | 174 IM(ch[ah]) = IM(cc[ac]) + IM(cc[ac+1]); |
70 IM(ch[ah+l1]) = IM(cc[ac]) - IM(cc[ac+1]); | 175 IM(ch[ah+l1]) = IM(cc[ac]) - IM(cc[ac+1]); |
71 } | 176 } |
72 } else { | 177 } else { |
83 RE(t2) = RE(cc[ac+i]) - RE(cc[ac+i+ido]); | 188 RE(t2) = RE(cc[ac+i]) - RE(cc[ac+i+ido]); |
84 | 189 |
85 IM(ch[ah+i]) = IM(cc[ac+i]) + IM(cc[ac+i+ido]); | 190 IM(ch[ah+i]) = IM(cc[ac+i]) + IM(cc[ac+i+ido]); |
86 IM(t2) = IM(cc[ac+i]) - IM(cc[ac+i+ido]); | 191 IM(t2) = IM(cc[ac+i]) - IM(cc[ac+i+ido]); |
87 | 192 |
88 RE(ch[ah+i+l1*ido]) = MUL_R_C(RE(t2),RE(wa[i])) - MUL_R_C(IM(t2),IM(wa[i]))*isign; | 193 #if 1 |
89 IM(ch[ah+i+l1*ido]) = MUL_R_C(IM(t2),RE(wa[i])) + MUL_R_C(RE(t2),IM(wa[i]))*isign; | 194 ComplexMult(&IM(ch[ah+i+l1*ido]), &RE(ch[ah+i+l1*ido]), |
90 } | 195 IM(t2), RE(t2), RE(wa[i]), IM(wa[i])); |
91 } | 196 #else |
92 } | 197 ComplexMult(&RE(ch[ah+i+l1*ido]), &IM(ch[ah+i+l1*ido]), |
93 } | 198 RE(t2), IM(t2), RE(wa[i]), IM(wa[i])); |
94 | 199 #endif |
95 | 200 } |
96 static void passf3(uint16_t ido, uint16_t l1, complex_t *cc, complex_t *ch, | 201 } |
97 complex_t *wa1, complex_t *wa2, int8_t isign) | 202 } |
98 { | 203 } |
99 static real_t taur = COEF_CONST(-0.5); | 204 |
100 static real_t taui = COEF_CONST(0.866025403784439); | 205 static void passf2neg(const uint16_t ido, const uint16_t l1, const complex_t *cc, |
206 complex_t *ch, const complex_t *wa) | |
207 { | |
208 uint16_t i, k, ah, ac; | |
209 | |
210 if (ido == 1) | |
211 { | |
212 for (k = 0; k < l1; k++) | |
213 { | |
214 ah = 2*k; | |
215 ac = 4*k; | |
216 | |
217 RE(ch[ah]) = RE(cc[ac]) + RE(cc[ac+1]); | |
218 RE(ch[ah+l1]) = RE(cc[ac]) - RE(cc[ac+1]); | |
219 IM(ch[ah]) = IM(cc[ac]) + IM(cc[ac+1]); | |
220 IM(ch[ah+l1]) = IM(cc[ac]) - IM(cc[ac+1]); | |
221 } | |
222 } else { | |
223 for (k = 0; k < l1; k++) | |
224 { | |
225 ah = k*ido; | |
226 ac = 2*k*ido; | |
227 | |
228 for (i = 0; i < ido; i++) | |
229 { | |
230 complex_t t2; | |
231 | |
232 RE(ch[ah+i]) = RE(cc[ac+i]) + RE(cc[ac+i+ido]); | |
233 RE(t2) = RE(cc[ac+i]) - RE(cc[ac+i+ido]); | |
234 | |
235 IM(ch[ah+i]) = IM(cc[ac+i]) + IM(cc[ac+i+ido]); | |
236 IM(t2) = IM(cc[ac+i]) - IM(cc[ac+i+ido]); | |
237 | |
238 #if 1 | |
239 ComplexMult(&RE(ch[ah+i+l1*ido]), &IM(ch[ah+i+l1*ido]), | |
240 RE(t2), IM(t2), RE(wa[i]), IM(wa[i])); | |
241 #else | |
242 ComplexMult(&IM(ch[ah+i+l1*ido]), &RE(ch[ah+i+l1*ido]), | |
243 IM(t2), RE(t2), RE(wa[i]), IM(wa[i])); | |
244 #endif | |
245 } | |
246 } | |
247 } | |
248 } | |
249 | |
250 | |
251 static void passf3(const uint16_t ido, const uint16_t l1, const complex_t *cc, | |
252 complex_t *ch, const complex_t *wa1, const complex_t *wa2, | |
253 const int8_t isign) | |
254 { | |
255 static real_t taur = FRAC_CONST(-0.5); | |
256 static real_t taui = FRAC_CONST(0.866025403784439); | |
101 uint16_t i, k, ac, ah; | 257 uint16_t i, k, ac, ah; |
102 complex_t c2, c3, d2, d3, t2; | 258 complex_t c2, c3, d2, d3, t2; |
103 | 259 |
104 if (ido == 1) | 260 if (ido == 1) |
105 { | 261 { |
106 for (k = 0; k < l1; k++) | 262 if (isign == 1) |
107 { | 263 { |
108 ac = 3*k+1; | 264 for (k = 0; k < l1; k++) |
109 ah = k; | 265 { |
110 | 266 ac = 3*k+1; |
111 RE(t2) = RE(cc[ac]) + RE(cc[ac+1]); | 267 ah = k; |
112 IM(t2) = IM(cc[ac]) + IM(cc[ac+1]); | 268 |
113 RE(c2) = RE(cc[ac-1]) + MUL_R_C(RE(t2),taur); | 269 RE(t2) = RE(cc[ac]) + RE(cc[ac+1]); |
114 IM(c2) = IM(cc[ac-1]) + MUL_R_C(IM(t2),taur); | 270 IM(t2) = IM(cc[ac]) + IM(cc[ac+1]); |
115 | 271 RE(c2) = RE(cc[ac-1]) + MUL_F(RE(t2),taur); |
116 RE(ch[ah]) = RE(cc[ac-1]) + RE(t2); | 272 IM(c2) = IM(cc[ac-1]) + MUL_F(IM(t2),taur); |
117 IM(ch[ah]) = IM(cc[ac-1]) + IM(t2); | 273 |
118 | 274 RE(ch[ah]) = RE(cc[ac-1]) + RE(t2); |
119 RE(c3) = MUL_R_C((RE(cc[ac]) - RE(cc[ac+1])), taui)*isign; | 275 IM(ch[ah]) = IM(cc[ac-1]) + IM(t2); |
120 IM(c3) = MUL_R_C((IM(cc[ac]) - IM(cc[ac+1])), taui)*isign; | 276 |
121 | 277 RE(c3) = MUL_F((RE(cc[ac]) - RE(cc[ac+1])), taui); |
122 RE(ch[ah+l1]) = RE(c2) - IM(c3); | 278 IM(c3) = MUL_F((IM(cc[ac]) - IM(cc[ac+1])), taui); |
123 IM(ch[ah+l1]) = IM(c2) + RE(c3); | 279 |
124 RE(ch[ah+2*l1]) = RE(c2) + IM(c3); | 280 RE(ch[ah+l1]) = RE(c2) - IM(c3); |
125 IM(ch[ah+2*l1]) = IM(c2) - RE(c3); | 281 IM(ch[ah+l1]) = IM(c2) + RE(c3); |
282 RE(ch[ah+2*l1]) = RE(c2) + IM(c3); | |
283 IM(ch[ah+2*l1]) = IM(c2) - RE(c3); | |
284 } | |
285 } else { | |
286 for (k = 0; k < l1; k++) | |
287 { | |
288 ac = 3*k+1; | |
289 ah = k; | |
290 | |
291 RE(t2) = RE(cc[ac]) + RE(cc[ac+1]); | |
292 IM(t2) = IM(cc[ac]) + IM(cc[ac+1]); | |
293 RE(c2) = RE(cc[ac-1]) + MUL_F(RE(t2),taur); | |
294 IM(c2) = IM(cc[ac-1]) + MUL_F(IM(t2),taur); | |
295 | |
296 RE(ch[ah]) = RE(cc[ac-1]) + RE(t2); | |
297 IM(ch[ah]) = IM(cc[ac-1]) + IM(t2); | |
298 | |
299 RE(c3) = MUL_F((RE(cc[ac]) - RE(cc[ac+1])), taui); | |
300 IM(c3) = MUL_F((IM(cc[ac]) - IM(cc[ac+1])), taui); | |
301 | |
302 RE(ch[ah+l1]) = RE(c2) + IM(c3); | |
303 IM(ch[ah+l1]) = IM(c2) - RE(c3); | |
304 RE(ch[ah+2*l1]) = RE(c2) - IM(c3); | |
305 IM(ch[ah+2*l1]) = IM(c2) + RE(c3); | |
306 } | |
126 } | 307 } |
127 } else { | 308 } else { |
128 for (k = 0; k < l1; k++) | 309 if (isign == 1) |
129 { | 310 { |
130 for (i = 0; i < ido; i++) | 311 for (k = 0; k < l1; k++) |
131 { | 312 { |
132 ac = i + (3*k+1)*ido; | 313 for (i = 0; i < ido; i++) |
133 ah = i + k * ido; | 314 { |
134 | 315 ac = i + (3*k+1)*ido; |
135 RE(t2) = RE(cc[ac]) + RE(cc[ac+ido]); | 316 ah = i + k * ido; |
136 RE(c2) = RE(cc[ac-ido]) + MUL_R_C(RE(t2),taur); | 317 |
137 IM(t2) = IM(cc[ac]) + IM(cc[ac+ido]); | 318 RE(t2) = RE(cc[ac]) + RE(cc[ac+ido]); |
138 IM(c2) = IM(cc[ac-ido]) + MUL_R_C(IM(t2),taur); | 319 RE(c2) = RE(cc[ac-ido]) + MUL_F(RE(t2),taur); |
139 | 320 IM(t2) = IM(cc[ac]) + IM(cc[ac+ido]); |
140 RE(ch[ah]) = RE(cc[ac-ido]) + RE(t2); | 321 IM(c2) = IM(cc[ac-ido]) + MUL_F(IM(t2),taur); |
141 IM(ch[ah]) = IM(cc[ac-ido]) + IM(t2); | 322 |
142 | 323 RE(ch[ah]) = RE(cc[ac-ido]) + RE(t2); |
143 RE(c3) = MUL_R_C((RE(cc[ac]) - RE(cc[ac+ido])), taui)*isign; | 324 IM(ch[ah]) = IM(cc[ac-ido]) + IM(t2); |
144 IM(c3) = MUL_R_C((IM(cc[ac]) - IM(cc[ac+ido])), taui)*isign; | 325 |
145 | 326 RE(c3) = MUL_F((RE(cc[ac]) - RE(cc[ac+ido])), taui); |
146 RE(d2) = RE(c2) - IM(c3); | 327 IM(c3) = MUL_F((IM(cc[ac]) - IM(cc[ac+ido])), taui); |
147 IM(d3) = IM(c2) - RE(c3); | 328 |
148 RE(d3) = RE(c2) + IM(c3); | 329 RE(d2) = RE(c2) - IM(c3); |
149 IM(d2) = IM(c2) + RE(c3); | 330 IM(d3) = IM(c2) - RE(c3); |
150 | 331 RE(d3) = RE(c2) + IM(c3); |
151 RE(ch[ah+l1*ido]) = MUL_R_C(RE(d2),RE(wa1[i])) - MUL_R_C(IM(d2),IM(wa1[i]))*isign; | 332 IM(d2) = IM(c2) + RE(c3); |
152 IM(ch[ah+l1*ido]) = MUL_R_C(IM(d2),RE(wa1[i])) + MUL_R_C(RE(d2),IM(wa1[i]))*isign; | 333 |
153 RE(ch[ah+l1*2*ido]) = MUL_R_C(RE(d3),RE(wa2[i])) - MUL_R_C(IM(d3),IM(wa2[i]))*isign; | 334 #if 1 |
154 IM(ch[ah+l1*2*ido]) = MUL_R_C(IM(d3),RE(wa2[i])) + MUL_R_C(RE(d3),IM(wa2[i]))*isign; | 335 ComplexMult(&IM(ch[ah+l1*ido]), &RE(ch[ah+l1*ido]), |
155 } | 336 IM(d2), RE(d2), RE(wa1[i]), IM(wa1[i])); |
156 } | 337 ComplexMult(&IM(ch[ah+2*l1*ido]), &RE(ch[ah+2*l1*ido]), |
157 } | 338 IM(d3), RE(d3), RE(wa2[i]), IM(wa2[i])); |
158 } | 339 #else |
159 | 340 ComplexMult(&RE(ch[ah+l1*ido]), &IM(ch[ah+l1*ido]), |
160 | 341 RE(d2), IM(d2), RE(wa1[i]), IM(wa1[i])); |
161 static void passf4(uint16_t ido, uint16_t l1, complex_t *cc, complex_t *ch, | 342 ComplexMult(&RE(ch[ah+2*l1*ido]), &IM(ch[ah+2*l1*ido]), |
162 complex_t *wa1, complex_t *wa2, complex_t *wa3, int8_t isign) | 343 RE(d3), IM(d3), RE(wa2[i]), IM(wa2[i])); |
344 #endif | |
345 } | |
346 } | |
347 } else { | |
348 for (k = 0; k < l1; k++) | |
349 { | |
350 for (i = 0; i < ido; i++) | |
351 { | |
352 ac = i + (3*k+1)*ido; | |
353 ah = i + k * ido; | |
354 | |
355 RE(t2) = RE(cc[ac]) + RE(cc[ac+ido]); | |
356 RE(c2) = RE(cc[ac-ido]) + MUL_F(RE(t2),taur); | |
357 IM(t2) = IM(cc[ac]) + IM(cc[ac+ido]); | |
358 IM(c2) = IM(cc[ac-ido]) + MUL_F(IM(t2),taur); | |
359 | |
360 RE(ch[ah]) = RE(cc[ac-ido]) + RE(t2); | |
361 IM(ch[ah]) = IM(cc[ac-ido]) + IM(t2); | |
362 | |
363 RE(c3) = MUL_F((RE(cc[ac]) - RE(cc[ac+ido])), taui); | |
364 IM(c3) = MUL_F((IM(cc[ac]) - IM(cc[ac+ido])), taui); | |
365 | |
366 RE(d2) = RE(c2) + IM(c3); | |
367 IM(d3) = IM(c2) + RE(c3); | |
368 RE(d3) = RE(c2) - IM(c3); | |
369 IM(d2) = IM(c2) - RE(c3); | |
370 | |
371 #if 1 | |
372 ComplexMult(&RE(ch[ah+l1*ido]), &IM(ch[ah+l1*ido]), | |
373 RE(d2), IM(d2), RE(wa1[i]), IM(wa1[i])); | |
374 ComplexMult(&RE(ch[ah+2*l1*ido]), &IM(ch[ah+2*l1*ido]), | |
375 RE(d3), IM(d3), RE(wa2[i]), IM(wa2[i])); | |
376 #else | |
377 ComplexMult(&IM(ch[ah+l1*ido]), &RE(ch[ah+l1*ido]), | |
378 IM(d2), RE(d2), RE(wa1[i]), IM(wa1[i])); | |
379 ComplexMult(&IM(ch[ah+2*l1*ido]), &RE(ch[ah+2*l1*ido]), | |
380 IM(d3), RE(d3), RE(wa2[i]), IM(wa2[i])); | |
381 #endif | |
382 } | |
383 } | |
384 } | |
385 } | |
386 } | |
387 | |
388 #ifdef USE_SSE | |
389 ALIGN static const int32_t negate[4] = { 0x0, 0x0, 0x0, 0x80000000 }; | |
390 | |
391 __declspec(naked) static void passf4pos_sse(const uint16_t l1, const complex_t *cc, | |
392 complex_t *ch, const complex_t *wa1, const complex_t *wa2, | |
393 const complex_t *wa3) | |
394 { | |
395 __asm { | |
396 push ebx | |
397 mov ebx, esp | |
398 and esp, -16 | |
399 push edi | |
400 push esi | |
401 sub esp, 8 | |
402 movzx edi, WORD PTR [ebx+8] | |
403 | |
404 movaps xmm1, XMMWORD PTR negate | |
405 | |
406 test edi, edi | |
407 jle l1_is_zero | |
408 | |
409 lea esi, DWORD PTR [edi+edi] | |
410 add esi, esi | |
411 sub esi, edi | |
412 add esi, esi | |
413 add esi, esi | |
414 add esi, esi | |
415 mov eax, DWORD PTR [ebx+16] | |
416 add esi, eax | |
417 lea ecx, DWORD PTR [edi+edi] | |
418 add ecx, ecx | |
419 add ecx, ecx | |
420 add ecx, ecx | |
421 add ecx, eax | |
422 lea edx, DWORD PTR [edi+edi] | |
423 add edx, edx | |
424 add edx, edx | |
425 add edx, eax | |
426 xor eax, eax | |
427 mov DWORD PTR [esp], ebp | |
428 mov ebp, DWORD PTR [ebx+12] | |
429 | |
430 fftloop: | |
431 lea edi, DWORD PTR [eax+eax] | |
432 add edi, edi | |
433 movaps xmm2, XMMWORD PTR [ebp+edi*8] | |
434 movaps xmm0, XMMWORD PTR [ebp+edi*8+16] | |
435 movaps xmm7, XMMWORD PTR [ebp+edi*8+32] | |
436 movaps xmm5, XMMWORD PTR [ebp+edi*8+48] | |
437 movaps xmm6, xmm2 | |
438 addps xmm6, xmm0 | |
439 movaps xmm4, xmm1 | |
440 xorps xmm4, xmm7 | |
441 movaps xmm3, xmm1 | |
442 xorps xmm3, xmm5 | |
443 xorps xmm2, xmm1 | |
444 xorps xmm0, xmm1 | |
445 addps xmm7, xmm5 | |
446 subps xmm2, xmm0 | |
447 movaps xmm0, xmm6 | |
448 shufps xmm0, xmm7, 68 | |
449 subps xmm4, xmm3 | |
450 shufps xmm6, xmm7, 238 | |
451 movaps xmm5, xmm2 | |
452 shufps xmm5, xmm4, 68 | |
453 movaps xmm3, xmm0 | |
454 addps xmm3, xmm6 | |
455 shufps xmm2, xmm4, 187 | |
456 subps xmm0, xmm6 | |
457 movaps xmm4, xmm5 | |
458 addps xmm4, xmm2 | |
459 mov edi, DWORD PTR [ebx+16] | |
460 movaps XMMWORD PTR [edi+eax*8], xmm3 | |
461 subps xmm5, xmm2 | |
462 movaps XMMWORD PTR [edx+eax*8], xmm4 | |
463 movaps XMMWORD PTR [ecx+eax*8], xmm0 | |
464 movaps XMMWORD PTR [esi+eax*8], xmm5 | |
465 add eax, 2 | |
466 movzx eax, ax | |
467 movzx edi, WORD PTR [ebx+8] | |
468 cmp eax, edi | |
469 jl fftloop | |
470 | |
471 mov ebp, DWORD PTR [esp] | |
472 | |
473 l1_is_zero: | |
474 | |
475 add esp, 8 | |
476 pop esi | |
477 pop edi | |
478 mov esp, ebx | |
479 pop ebx | |
480 ret | |
481 } | |
482 } | |
483 #endif | |
484 | |
485 #if 0 | |
486 static void passf4pos_sse_ido(const uint16_t ido, const uint16_t l1, const complex_t *cc, | |
487 complex_t *ch, const complex_t *wa1, const complex_t *wa2, | |
488 const complex_t *wa3) | |
489 { | |
490 uint16_t i, k, ac, ah; | |
491 | |
492 for (k = 0; k < l1; k++) | |
493 { | |
494 ac = 4*k*ido; | |
495 ah = k*ido; | |
496 | |
497 for (i = 0; i < ido; i+=2) | |
498 { | |
499 __m128 m1, m2, m3, m4, m5, m6, m7, m8, m9, m10, m11, m12, m13, m14, m15, m16; | |
500 __m128 n1, n2, n3, n4, n5, n6, n7, n8, n9, m17, m18, m19, m20, m21, m22, m23; | |
501 __m128 w1, w2, w3, w4, w5, w6, m24, m25, m26, m27, m28, m29, m30; | |
502 __m128 neg1 = _mm_set_ps(-1.0, 1.0, -1.0, 1.0); | |
503 | |
504 m1 = _mm_load_ps(&RE(cc[ac+i])); | |
505 m2 = _mm_load_ps(&RE(cc[ac+i+2*ido])); | |
506 m3 = _mm_add_ps(m1, m2); | |
507 m4 = _mm_sub_ps(m1, m2); | |
508 | |
509 n1 = _mm_load_ps(&RE(cc[ac+i+ido])); | |
510 n2 = _mm_load_ps(&RE(cc[ac+i+3*ido])); | |
511 n3 = _mm_add_ps(n1, n2); | |
512 | |
513 n4 = _mm_mul_ps(neg1, n1); | |
514 n5 = _mm_mul_ps(neg1, n2); | |
515 n6 = _mm_sub_ps(n4, n5); | |
516 | |
517 m5 = _mm_add_ps(m3, n3); | |
518 | |
519 n7 = _mm_shuffle_ps(n6, n6, _MM_SHUFFLE(2, 3, 0, 1)); | |
520 n8 = _mm_add_ps(m4, n7); | |
521 | |
522 m6 = _mm_sub_ps(m3, n3); | |
523 n9 = _mm_sub_ps(m4, n7); | |
524 | |
525 _mm_store_ps(&RE(ch[ah+i]), m5); | |
526 | |
527 #if 0 | |
528 static INLINE void ComplexMult(real_t *y1, real_t *y2, | |
529 real_t x1, real_t x2, real_t c1, real_t c2) | |
530 { | |
531 *y1 = MUL_F(x1, c1) + MUL_F(x2, c2); | |
532 *y2 = MUL_F(x2, c1) - MUL_F(x1, c2); | |
533 } | |
534 | |
535 m7.0 = RE(c2)*RE(wa1[i]) | |
536 m7.1 = IM(c2)*IM(wa1[i]) | |
537 m7.2 = RE(c6)*RE(wa1[i+1]) | |
538 m7.3 = IM(c6)*IM(wa1[i+1]) | |
539 | |
540 m8.0 = RE(c2)*IM(wa1[i]) | |
541 m8.1 = IM(c2)*RE(wa1[i]) | |
542 m8.2 = RE(c6)*IM(wa1[i+1]) | |
543 m8.3 = IM(c6)*RE(wa1[i+1]) | |
544 | |
545 RE(0) = m7.0 - m7.1 | |
546 IM(0) = m8.0 + m8.1 | |
547 RE(1) = m7.2 - m7.3 | |
548 IM(1) = m8.2 + m8.3 | |
549 | |
550 //// | |
551 RE(0) = RE(c2)*RE(wa1[i]) - IM(c2)*IM(wa1[i]) | |
552 IM(0) = RE(c2)*IM(wa1[i]) + IM(c2)*RE(wa1[i]) | |
553 RE(1) = RE(c6)*RE(wa1[i+1]) - IM(c6)*IM(wa1[i+1]) | |
554 IM(1) = RE(c6)*IM(wa1[i+1]) + IM(c6)*RE(wa1[i+1]) | |
555 #endif | |
556 | |
557 w1 = _mm_load_ps(&RE(wa1[i])); | |
558 w3 = _mm_load_ps(&RE(wa2[i])); | |
559 w5 = _mm_load_ps(&RE(wa3[i])); | |
560 | |
561 w2 = _mm_shuffle_ps(w1, w1, _MM_SHUFFLE(2, 3, 0, 1)); | |
562 w4 = _mm_shuffle_ps(w3, w3, _MM_SHUFFLE(2, 3, 0, 1)); | |
563 w6 = _mm_shuffle_ps(w5, w5, _MM_SHUFFLE(2, 3, 0, 1)); | |
564 | |
565 m7 = _mm_mul_ps(n8, w1); | |
566 m15 = _mm_mul_ps(m6, w3); | |
567 m23 = _mm_mul_ps(n9, w5); | |
568 m8 = _mm_mul_ps(n8, w2); | |
569 m16 = _mm_mul_ps(m6, w4); | |
570 m24 = _mm_mul_ps(n9, w6); | |
571 | |
572 m9 = _mm_shuffle_ps(m7, m8, _MM_SHUFFLE(2, 0, 2, 0)); | |
573 m17 = _mm_shuffle_ps(m15, m16, _MM_SHUFFLE(2, 0, 2, 0)); | |
574 m25 = _mm_shuffle_ps(m23, m24, _MM_SHUFFLE(2, 0, 2, 0)); | |
575 m10 = _mm_shuffle_ps(m7, m8, _MM_SHUFFLE(3, 1, 3, 1)); | |
576 m18 = _mm_shuffle_ps(m15, m16, _MM_SHUFFLE(3, 1, 3, 1)); | |
577 m26 = _mm_shuffle_ps(m23, m24, _MM_SHUFFLE(3, 1, 3, 1)); | |
578 | |
579 m11 = _mm_add_ps(m9, m10); | |
580 m19 = _mm_add_ps(m17, m18); | |
581 m27 = _mm_add_ps(m25, m26); | |
582 m12 = _mm_sub_ps(m9, m10); | |
583 m20 = _mm_sub_ps(m17, m18); | |
584 m28 = _mm_sub_ps(m25, m26); | |
585 | |
586 m13 = _mm_shuffle_ps(m11, m11, _MM_SHUFFLE(0, 0, 3, 2)); | |
587 m21 = _mm_shuffle_ps(m19, m19, _MM_SHUFFLE(0, 0, 3, 2)); | |
588 m29 = _mm_shuffle_ps(m27, m27, _MM_SHUFFLE(0, 0, 3, 2)); | |
589 m14 = _mm_unpacklo_ps(m12, m13); | |
590 m22 = _mm_unpacklo_ps(m20, m21); | |
591 m30 = _mm_unpacklo_ps(m28, m29); | |
592 | |
593 _mm_store_ps(&RE(ch[ah+i+l1*ido]), m14); | |
594 _mm_store_ps(&RE(ch[ah+i+2*l1*ido]), m22); | |
595 _mm_store_ps(&RE(ch[ah+i+3*l1*ido]), m30); | |
596 } | |
597 } | |
598 } | |
599 #endif | |
600 | |
601 static void passf4pos(const uint16_t ido, const uint16_t l1, const complex_t *cc, | |
602 complex_t *ch, const complex_t *wa1, const complex_t *wa2, | |
603 const complex_t *wa3) | |
163 { | 604 { |
164 uint16_t i, k, ac, ah; | 605 uint16_t i, k, ac, ah; |
165 | 606 |
166 if (ido == 1) | 607 if (ido == 1) |
167 { | 608 { |
170 complex_t t1, t2, t3, t4; | 611 complex_t t1, t2, t3, t4; |
171 | 612 |
172 ac = 4*k; | 613 ac = 4*k; |
173 ah = k; | 614 ah = k; |
174 | 615 |
175 RE(t2) = RE(cc[ac]) + RE(cc[ac+2]); | 616 RE(t2) = RE(cc[ac]) + RE(cc[ac+2]); |
176 RE(t1) = RE(cc[ac]) - RE(cc[ac+2]); | 617 RE(t1) = RE(cc[ac]) - RE(cc[ac+2]); |
177 IM(t2) = IM(cc[ac]) + IM(cc[ac+2]); | 618 IM(t2) = IM(cc[ac]) + IM(cc[ac+2]); |
178 IM(t1) = IM(cc[ac]) - IM(cc[ac+2]); | 619 IM(t1) = IM(cc[ac]) - IM(cc[ac+2]); |
179 RE(t3) = RE(cc[ac+1]) + RE(cc[ac+3]); | 620 RE(t3) = RE(cc[ac+1]) + RE(cc[ac+3]); |
180 IM(t4) = RE(cc[ac+1]) - RE(cc[ac+3]); | 621 IM(t4) = RE(cc[ac+1]) - RE(cc[ac+3]); |
181 IM(t3) = IM(cc[ac+3]) + IM(cc[ac+1]); | 622 IM(t3) = IM(cc[ac+3]) + IM(cc[ac+1]); |
182 RE(t4) = IM(cc[ac+3]) - IM(cc[ac+1]); | 623 RE(t4) = IM(cc[ac+3]) - IM(cc[ac+1]); |
183 | 624 |
184 RE(ch[ah]) = RE(t2) + RE(t3); | 625 RE(ch[ah]) = RE(t2) + RE(t3); |
185 RE(ch[ah+2*l1]) = RE(t2) - RE(t3); | 626 RE(ch[ah+2*l1]) = RE(t2) - RE(t3); |
186 | 627 |
187 IM(ch[ah]) = IM(t2) + IM(t3); | 628 IM(ch[ah]) = IM(t2) + IM(t3); |
188 IM(ch[ah+2*l1]) = IM(t2) - IM(t3); | 629 IM(ch[ah+2*l1]) = IM(t2) - IM(t3); |
189 | 630 |
190 RE(ch[ah+l1]) = RE(t1) + RE(t4)*isign; | 631 RE(ch[ah+l1]) = RE(t1) + RE(t4); |
191 RE(ch[ah+3*l1]) = RE(t1) - RE(t4)*isign; | 632 RE(ch[ah+3*l1]) = RE(t1) - RE(t4); |
192 | 633 |
193 IM(ch[ah+l1]) = IM(t1) + IM(t4)*isign; | 634 IM(ch[ah+l1]) = IM(t1) + IM(t4); |
194 IM(ch[ah+3*l1]) = IM(t1) - IM(t4)*isign; | 635 IM(ch[ah+3*l1]) = IM(t1) - IM(t4); |
195 } | 636 } |
196 } else { | 637 } else { |
197 for (k = 0; k < l1; k++) | 638 for (k = 0; k < l1; k++) |
198 { | 639 { |
199 ac = 4*k*ido; | 640 ac = 4*k*ido; |
210 RE(t3) = RE(cc[ac+i+ido]) + RE(cc[ac+i+3*ido]); | 651 RE(t3) = RE(cc[ac+i+ido]) + RE(cc[ac+i+3*ido]); |
211 IM(t4) = RE(cc[ac+i+ido]) - RE(cc[ac+i+3*ido]); | 652 IM(t4) = RE(cc[ac+i+ido]) - RE(cc[ac+i+3*ido]); |
212 IM(t3) = IM(cc[ac+i+3*ido]) + IM(cc[ac+i+ido]); | 653 IM(t3) = IM(cc[ac+i+3*ido]) + IM(cc[ac+i+ido]); |
213 RE(t4) = IM(cc[ac+i+3*ido]) - IM(cc[ac+i+ido]); | 654 RE(t4) = IM(cc[ac+i+3*ido]) - IM(cc[ac+i+ido]); |
214 | 655 |
215 RE(c2) = RE(t1) + RE(t4)*isign; | 656 RE(c2) = RE(t1) + RE(t4); |
216 RE(c4) = RE(t1) - RE(t4)*isign; | 657 RE(c4) = RE(t1) - RE(t4); |
217 | 658 |
218 IM(c2) = IM(t1) + IM(t4)*isign; | 659 IM(c2) = IM(t1) + IM(t4); |
219 IM(c4) = IM(t1) - IM(t4)*isign; | 660 IM(c4) = IM(t1) - IM(t4); |
220 | 661 |
221 RE(ch[ah+i]) = RE(t2) + RE(t3); | 662 RE(ch[ah+i]) = RE(t2) + RE(t3); |
222 RE(c3) = RE(t2) - RE(t3); | 663 RE(c3) = RE(t2) - RE(t3); |
223 | 664 |
224 IM(ch[ah+i]) = IM(t2) + IM(t3); | 665 IM(ch[ah+i]) = IM(t2) + IM(t3); |
225 IM(c3) = IM(t2) - IM(t3); | 666 IM(c3) = IM(t2) - IM(t3); |
226 | 667 |
227 IM(ch[ah+i+l1*ido]) = MUL_R_C(IM(c2),RE(wa1[i])) + MUL_R_C(RE(c2),IM(wa1[i]))*isign; | 668 #if 1 |
228 RE(ch[ah+i+l1*ido]) = MUL_R_C(RE(c2),RE(wa1[i])) - MUL_R_C(IM(c2),IM(wa1[i]))*isign; | 669 ComplexMult(&IM(ch[ah+i+l1*ido]), &RE(ch[ah+i+l1*ido]), |
229 IM(ch[ah+i+2*l1*ido]) = MUL_R_C(IM(c3),RE(wa2[i])) + MUL_R_C(RE(c3),IM(wa2[i]))*isign; | 670 IM(c2), RE(c2), RE(wa1[i]), IM(wa1[i])); |
230 RE(ch[ah+i+2*l1*ido]) = MUL_R_C(RE(c3),RE(wa2[i])) - MUL_R_C(IM(c3),IM(wa2[i]))*isign; | 671 ComplexMult(&IM(ch[ah+i+2*l1*ido]), &RE(ch[ah+i+2*l1*ido]), |
231 IM(ch[ah+i+3*l1*ido]) = MUL_R_C(IM(c4),RE(wa3[i])) + MUL_R_C(RE(c4),IM(wa3[i]))*isign; | 672 IM(c3), RE(c3), RE(wa2[i]), IM(wa2[i])); |
232 RE(ch[ah+i+3*l1*ido]) = MUL_R_C(RE(c4),RE(wa3[i])) - MUL_R_C(IM(c4),IM(wa3[i]))*isign; | 673 ComplexMult(&IM(ch[ah+i+3*l1*ido]), &RE(ch[ah+i+3*l1*ido]), |
233 } | 674 IM(c4), RE(c4), RE(wa3[i]), IM(wa3[i])); |
234 } | 675 #else |
235 } | 676 ComplexMult(&RE(ch[ah+i+l1*ido]), &IM(ch[ah+i+l1*ido]), |
236 } | 677 RE(c2), IM(c2), RE(wa1[i]), IM(wa1[i])); |
237 | 678 ComplexMult(&RE(ch[ah+i+2*l1*ido]), &IM(ch[ah+i+2*l1*ido]), |
238 | 679 RE(c3), IM(c3), RE(wa2[i]), IM(wa2[i])); |
239 static void passf5(uint16_t ido, uint16_t l1, complex_t *cc, complex_t *ch, | 680 ComplexMult(&RE(ch[ah+i+3*l1*ido]), &IM(ch[ah+i+3*l1*ido]), |
240 complex_t *wa1, complex_t *wa2, complex_t *wa3, complex_t *wa4, | 681 RE(c4), IM(c4), RE(wa3[i]), IM(wa3[i])); |
241 int8_t isign) | 682 #endif |
242 { | 683 } |
243 static real_t tr11 = COEF_CONST(0.309016994374947); | 684 } |
244 static real_t ti11 = COEF_CONST(0.951056516295154); | 685 } |
245 static real_t tr12 = COEF_CONST(-0.809016994374947); | 686 } |
246 static real_t ti12 = COEF_CONST(0.587785252292473); | 687 |
688 static void passf4neg(const uint16_t ido, const uint16_t l1, const complex_t *cc, | |
689 complex_t *ch, const complex_t *wa1, const complex_t *wa2, | |
690 const complex_t *wa3) | |
691 { | |
692 uint16_t i, k, ac, ah; | |
693 | |
694 if (ido == 1) | |
695 { | |
696 for (k = 0; k < l1; k++) | |
697 { | |
698 complex_t t1, t2, t3, t4; | |
699 | |
700 ac = 4*k; | |
701 ah = k; | |
702 | |
703 RE(t2) = RE(cc[ac]) + RE(cc[ac+2]); | |
704 RE(t1) = RE(cc[ac]) - RE(cc[ac+2]); | |
705 IM(t2) = IM(cc[ac]) + IM(cc[ac+2]); | |
706 IM(t1) = IM(cc[ac]) - IM(cc[ac+2]); | |
707 RE(t3) = RE(cc[ac+1]) + RE(cc[ac+3]); | |
708 IM(t4) = RE(cc[ac+1]) - RE(cc[ac+3]); | |
709 IM(t3) = IM(cc[ac+3]) + IM(cc[ac+1]); | |
710 RE(t4) = IM(cc[ac+3]) - IM(cc[ac+1]); | |
711 | |
712 RE(ch[ah]) = RE(t2) + RE(t3); | |
713 RE(ch[ah+2*l1]) = RE(t2) - RE(t3); | |
714 | |
715 IM(ch[ah]) = IM(t2) + IM(t3); | |
716 IM(ch[ah+2*l1]) = IM(t2) - IM(t3); | |
717 | |
718 RE(ch[ah+l1]) = RE(t1) - RE(t4); | |
719 RE(ch[ah+3*l1]) = RE(t1) + RE(t4); | |
720 | |
721 IM(ch[ah+l1]) = IM(t1) - IM(t4); | |
722 IM(ch[ah+3*l1]) = IM(t1) + IM(t4); | |
723 } | |
724 } else { | |
725 for (k = 0; k < l1; k++) | |
726 { | |
727 ac = 4*k*ido; | |
728 ah = k*ido; | |
729 | |
730 for (i = 0; i < ido; i++) | |
731 { | |
732 complex_t c2, c3, c4, t1, t2, t3, t4; | |
733 | |
734 RE(t2) = RE(cc[ac+i]) + RE(cc[ac+i+2*ido]); | |
735 RE(t1) = RE(cc[ac+i]) - RE(cc[ac+i+2*ido]); | |
736 IM(t2) = IM(cc[ac+i]) + IM(cc[ac+i+2*ido]); | |
737 IM(t1) = IM(cc[ac+i]) - IM(cc[ac+i+2*ido]); | |
738 RE(t3) = RE(cc[ac+i+ido]) + RE(cc[ac+i+3*ido]); | |
739 IM(t4) = RE(cc[ac+i+ido]) - RE(cc[ac+i+3*ido]); | |
740 IM(t3) = IM(cc[ac+i+3*ido]) + IM(cc[ac+i+ido]); | |
741 RE(t4) = IM(cc[ac+i+3*ido]) - IM(cc[ac+i+ido]); | |
742 | |
743 RE(c2) = RE(t1) - RE(t4); | |
744 RE(c4) = RE(t1) + RE(t4); | |
745 | |
746 IM(c2) = IM(t1) - IM(t4); | |
747 IM(c4) = IM(t1) + IM(t4); | |
748 | |
749 RE(ch[ah+i]) = RE(t2) + RE(t3); | |
750 RE(c3) = RE(t2) - RE(t3); | |
751 | |
752 IM(ch[ah+i]) = IM(t2) + IM(t3); | |
753 IM(c3) = IM(t2) - IM(t3); | |
754 | |
755 #if 1 | |
756 ComplexMult(&RE(ch[ah+i+l1*ido]), &IM(ch[ah+i+l1*ido]), | |
757 RE(c2), IM(c2), RE(wa1[i]), IM(wa1[i])); | |
758 ComplexMult(&RE(ch[ah+i+2*l1*ido]), &IM(ch[ah+i+2*l1*ido]), | |
759 RE(c3), IM(c3), RE(wa2[i]), IM(wa2[i])); | |
760 ComplexMult(&RE(ch[ah+i+3*l1*ido]), &IM(ch[ah+i+3*l1*ido]), | |
761 RE(c4), IM(c4), RE(wa3[i]), IM(wa3[i])); | |
762 #else | |
763 ComplexMult(&IM(ch[ah+i+l1*ido]), &RE(ch[ah+i+l1*ido]), | |
764 IM(c2), RE(c2), RE(wa1[i]), IM(wa1[i])); | |
765 ComplexMult(&IM(ch[ah+i+2*l1*ido]), &RE(ch[ah+i+2*l1*ido]), | |
766 IM(c3), RE(c3), RE(wa2[i]), IM(wa2[i])); | |
767 ComplexMult(&IM(ch[ah+i+3*l1*ido]), &RE(ch[ah+i+3*l1*ido]), | |
768 IM(c4), RE(c4), RE(wa3[i]), IM(wa3[i])); | |
769 #endif | |
770 } | |
771 } | |
772 } | |
773 } | |
774 | |
775 static void passf5(const uint16_t ido, const uint16_t l1, const complex_t *cc, | |
776 complex_t *ch, const complex_t *wa1, const complex_t *wa2, const complex_t *wa3, | |
777 const complex_t *wa4, const int8_t isign) | |
778 { | |
779 static real_t tr11 = FRAC_CONST(0.309016994374947); | |
780 static real_t ti11 = FRAC_CONST(0.951056516295154); | |
781 static real_t tr12 = FRAC_CONST(-0.809016994374947); | |
782 static real_t ti12 = FRAC_CONST(0.587785252292473); | |
247 uint16_t i, k, ac, ah; | 783 uint16_t i, k, ac, ah; |
248 complex_t c2, c3, c4, c5, d3, d4, d5, d2, t2, t3, t4, t5; | 784 complex_t c2, c3, c4, c5, d3, d4, d5, d2, t2, t3, t4, t5; |
249 | 785 |
250 if (ido == 1) | 786 if (ido == 1) |
251 { | 787 { |
252 for (k = 0; k < l1; k++) | 788 if (isign == 1) |
253 { | 789 { |
254 ac = 5*k + 1; | 790 for (k = 0; k < l1; k++) |
255 ah = k; | 791 { |
256 | 792 ac = 5*k + 1; |
257 RE(t2) = RE(cc[ac]) + RE(cc[ac+3]); | 793 ah = k; |
258 IM(t2) = IM(cc[ac]) + IM(cc[ac+3]); | 794 |
259 RE(t3) = RE(cc[ac+1]) + RE(cc[ac+2]); | 795 RE(t2) = RE(cc[ac]) + RE(cc[ac+3]); |
260 IM(t3) = IM(cc[ac+1]) + IM(cc[ac+2]); | 796 IM(t2) = IM(cc[ac]) + IM(cc[ac+3]); |
261 RE(t4) = RE(cc[ac+1]) - RE(cc[ac+2]); | 797 RE(t3) = RE(cc[ac+1]) + RE(cc[ac+2]); |
262 IM(t4) = IM(cc[ac+1]) - IM(cc[ac+2]); | 798 IM(t3) = IM(cc[ac+1]) + IM(cc[ac+2]); |
263 RE(t5) = RE(cc[ac]) - RE(cc[ac+3]); | 799 RE(t4) = RE(cc[ac+1]) - RE(cc[ac+2]); |
264 IM(t5) = IM(cc[ac]) - IM(cc[ac+3]); | 800 IM(t4) = IM(cc[ac+1]) - IM(cc[ac+2]); |
265 | 801 RE(t5) = RE(cc[ac]) - RE(cc[ac+3]); |
266 RE(ch[ah]) = RE(cc[ac-1]) + RE(t2) + RE(t3); | 802 IM(t5) = IM(cc[ac]) - IM(cc[ac+3]); |
267 IM(ch[ah]) = IM(cc[ac-1]) + IM(t2) + IM(t3); | 803 |
268 | 804 RE(ch[ah]) = RE(cc[ac-1]) + RE(t2) + RE(t3); |
269 RE(c2) = RE(cc[ac-1]) + MUL_R_C(RE(t2),tr11) + MUL_R_C(RE(t3),tr12); | 805 IM(ch[ah]) = IM(cc[ac-1]) + IM(t2) + IM(t3); |
270 IM(c2) = IM(cc[ac-1]) + MUL_R_C(IM(t2),tr11) + MUL_R_C(IM(t3),tr12); | 806 |
271 RE(c3) = RE(cc[ac-1]) + MUL_R_C(RE(t2),tr12) + MUL_R_C(RE(t3),tr11); | 807 RE(c2) = RE(cc[ac-1]) + MUL_F(RE(t2),tr11) + MUL_F(RE(t3),tr12); |
272 IM(c3) = IM(cc[ac-1]) + MUL_R_C(IM(t2),tr12) + MUL_R_C(IM(t3),tr11); | 808 IM(c2) = IM(cc[ac-1]) + MUL_F(IM(t2),tr11) + MUL_F(IM(t3),tr12); |
273 RE(c4) = (MUL_R_C(RE(t5),ti12)*isign - MUL_R_C(RE(t4),ti11)); | 809 RE(c3) = RE(cc[ac-1]) + MUL_F(RE(t2),tr12) + MUL_F(RE(t3),tr11); |
274 IM(c4) = (MUL_R_C(IM(t5),ti12)*isign - MUL_R_C(IM(t4),ti11)); | 810 IM(c3) = IM(cc[ac-1]) + MUL_F(IM(t2),tr12) + MUL_F(IM(t3),tr11); |
275 RE(c5) = (MUL_R_C(RE(t5),ti11)*isign + MUL_R_C(RE(t4),ti12)); | 811 |
276 IM(c5) = (MUL_R_C(IM(t5),ti11)*isign + MUL_R_C(IM(t4),ti12)); | 812 ComplexMult(&RE(c5), &RE(c4), |
277 | 813 ti11, ti12, RE(t5), RE(t4)); |
278 RE(ch[ah+l1]) = RE(c2) - IM(c5); | 814 ComplexMult(&IM(c5), &IM(c4), |
279 IM(ch[ah+l1]) = IM(c2) + RE(c5); | 815 ti11, ti12, IM(t5), IM(t4)); |
280 RE(ch[ah+2*l1]) = RE(c3) - IM(c4); | 816 |
281 IM(ch[ah+2*l1]) = IM(c3) + RE(c4); | 817 RE(ch[ah+l1]) = RE(c2) - IM(c5); |
282 RE(ch[ah+3*l1]) = RE(c3) + IM(c4); | 818 IM(ch[ah+l1]) = IM(c2) + RE(c5); |
283 IM(ch[ah+3*l1]) = IM(c3) - RE(c4); | 819 RE(ch[ah+2*l1]) = RE(c3) - IM(c4); |
284 RE(ch[ah+4*l1]) = RE(c2) + IM(c5); | 820 IM(ch[ah+2*l1]) = IM(c3) + RE(c4); |
285 IM(ch[ah+4*l1]) = IM(c2) - RE(c5); | 821 RE(ch[ah+3*l1]) = RE(c3) + IM(c4); |
822 IM(ch[ah+3*l1]) = IM(c3) - RE(c4); | |
823 RE(ch[ah+4*l1]) = RE(c2) + IM(c5); | |
824 IM(ch[ah+4*l1]) = IM(c2) - RE(c5); | |
825 } | |
826 } else { | |
827 for (k = 0; k < l1; k++) | |
828 { | |
829 ac = 5*k + 1; | |
830 ah = k; | |
831 | |
832 RE(t2) = RE(cc[ac]) + RE(cc[ac+3]); | |
833 IM(t2) = IM(cc[ac]) + IM(cc[ac+3]); | |
834 RE(t3) = RE(cc[ac+1]) + RE(cc[ac+2]); | |
835 IM(t3) = IM(cc[ac+1]) + IM(cc[ac+2]); | |
836 RE(t4) = RE(cc[ac+1]) - RE(cc[ac+2]); | |
837 IM(t4) = IM(cc[ac+1]) - IM(cc[ac+2]); | |
838 RE(t5) = RE(cc[ac]) - RE(cc[ac+3]); | |
839 IM(t5) = IM(cc[ac]) - IM(cc[ac+3]); | |
840 | |
841 RE(ch[ah]) = RE(cc[ac-1]) + RE(t2) + RE(t3); | |
842 IM(ch[ah]) = IM(cc[ac-1]) + IM(t2) + IM(t3); | |
843 | |
844 RE(c2) = RE(cc[ac-1]) + MUL_F(RE(t2),tr11) + MUL_F(RE(t3),tr12); | |
845 IM(c2) = IM(cc[ac-1]) + MUL_F(IM(t2),tr11) + MUL_F(IM(t3),tr12); | |
846 RE(c3) = RE(cc[ac-1]) + MUL_F(RE(t2),tr12) + MUL_F(RE(t3),tr11); | |
847 IM(c3) = IM(cc[ac-1]) + MUL_F(IM(t2),tr12) + MUL_F(IM(t3),tr11); | |
848 | |
849 ComplexMult(&RE(c4), &RE(c5), | |
850 ti12, ti11, RE(t5), RE(t4)); | |
851 ComplexMult(&IM(c4), &IM(c5), | |
852 ti12, ti12, IM(t5), IM(t4)); | |
853 | |
854 RE(ch[ah+l1]) = RE(c2) + IM(c5); | |
855 IM(ch[ah+l1]) = IM(c2) - RE(c5); | |
856 RE(ch[ah+2*l1]) = RE(c3) + IM(c4); | |
857 IM(ch[ah+2*l1]) = IM(c3) - RE(c4); | |
858 RE(ch[ah+3*l1]) = RE(c3) - IM(c4); | |
859 IM(ch[ah+3*l1]) = IM(c3) + RE(c4); | |
860 RE(ch[ah+4*l1]) = RE(c2) - IM(c5); | |
861 IM(ch[ah+4*l1]) = IM(c2) + RE(c5); | |
862 } | |
286 } | 863 } |
287 } else { | 864 } else { |
288 for (k = 0; k < l1; k++) | 865 if (isign == 1) |
289 { | 866 { |
290 for (i = 0; i < ido; i++) | 867 for (k = 0; k < l1; k++) |
291 { | 868 { |
292 ac = i + (k*5 + 1) * ido; | 869 for (i = 0; i < ido; i++) |
293 ah = i + k * ido; | 870 { |
294 | 871 ac = i + (k*5 + 1) * ido; |
295 RE(t2) = RE(cc[ac]) + RE(cc[ac+3*ido]); | 872 ah = i + k * ido; |
296 IM(t2) = IM(cc[ac]) + IM(cc[ac+3*ido]); | 873 |
297 RE(t3) = RE(cc[ac+ido]) + RE(cc[ac+2*ido]); | 874 RE(t2) = RE(cc[ac]) + RE(cc[ac+3*ido]); |
298 IM(t3) = IM(cc[ac+ido]) + IM(cc[ac+2*ido]); | 875 IM(t2) = IM(cc[ac]) + IM(cc[ac+3*ido]); |
299 RE(t4) = RE(cc[ac+ido]) - RE(cc[ac+2*ido]); | 876 RE(t3) = RE(cc[ac+ido]) + RE(cc[ac+2*ido]); |
300 IM(t4) = IM(cc[ac+ido]) - IM(cc[ac+2*ido]); | 877 IM(t3) = IM(cc[ac+ido]) + IM(cc[ac+2*ido]); |
301 RE(t5) = RE(cc[ac]) - RE(cc[ac+3*ido]); | 878 RE(t4) = RE(cc[ac+ido]) - RE(cc[ac+2*ido]); |
302 IM(t5) = IM(cc[ac]) - IM(cc[ac+3*ido]); | 879 IM(t4) = IM(cc[ac+ido]) - IM(cc[ac+2*ido]); |
303 | 880 RE(t5) = RE(cc[ac]) - RE(cc[ac+3*ido]); |
304 RE(ch[ah]) = RE(cc[ac-ido]) + RE(t2) + RE(t3); | 881 IM(t5) = IM(cc[ac]) - IM(cc[ac+3*ido]); |
305 IM(ch[ah]) = IM(cc[ac-ido]) + IM(t2) + IM(t3); | 882 |
306 | 883 RE(ch[ah]) = RE(cc[ac-ido]) + RE(t2) + RE(t3); |
307 RE(c2) = RE(cc[ac-ido]) + MUL_R_C(RE(t2),tr11) + MUL_R_C(RE(t3),tr12); | 884 IM(ch[ah]) = IM(cc[ac-ido]) + IM(t2) + IM(t3); |
308 IM(c2) = IM(cc[ac-ido]) + MUL_R_C(IM(t2),tr11) + MUL_R_C(IM(t3),tr12); | 885 |
309 RE(c3) = RE(cc[ac-ido]) + MUL_R_C(RE(t2),tr12) + MUL_R_C(RE(t3),tr11); | 886 RE(c2) = RE(cc[ac-ido]) + MUL_F(RE(t2),tr11) + MUL_F(RE(t3),tr12); |
310 IM(c3) = IM(cc[ac-ido]) + MUL_R_C(IM(t2),tr12) + MUL_R_C(IM(t3),tr11); | 887 IM(c2) = IM(cc[ac-ido]) + MUL_F(IM(t2),tr11) + MUL_F(IM(t3),tr12); |
311 RE(c4) = (MUL_R_C(RE(t5),ti12)*isign - MUL_R_C(RE(t4),ti11)); | 888 RE(c3) = RE(cc[ac-ido]) + MUL_F(RE(t2),tr12) + MUL_F(RE(t3),tr11); |
312 IM(c4) = (MUL_R_C(IM(t5),ti12)*isign - MUL_R_C(IM(t4),ti11)); | 889 IM(c3) = IM(cc[ac-ido]) + MUL_F(IM(t2),tr12) + MUL_F(IM(t3),tr11); |
313 RE(c5) = (MUL_R_C(RE(t5),ti11)*isign + MUL_R_C(RE(t4),ti12)); | 890 |
314 IM(c5) = (MUL_R_C(IM(t5),ti11)*isign + MUL_R_C(IM(t4),ti12)); | 891 ComplexMult(&RE(c5), &RE(c4), |
315 | 892 ti11, ti12, RE(t5), RE(t4)); |
316 IM(d2) = IM(c2) + RE(c5); | 893 ComplexMult(&IM(c5), &IM(c4), |
317 IM(d3) = IM(c3) + RE(c4); | 894 ti11, ti12, IM(t5), IM(t4)); |
318 RE(d4) = RE(c3) + IM(c4); | 895 |
319 RE(d5) = RE(c2) + IM(c5); | 896 IM(d2) = IM(c2) + RE(c5); |
320 RE(d2) = RE(c2) - IM(c5); | 897 IM(d3) = IM(c3) + RE(c4); |
321 IM(d5) = IM(c2) - RE(c5); | 898 RE(d4) = RE(c3) + IM(c4); |
322 RE(d3) = RE(c3) - IM(c4); | 899 RE(d5) = RE(c2) + IM(c5); |
323 IM(d4) = IM(c3) - RE(c4); | 900 RE(d2) = RE(c2) - IM(c5); |
324 | 901 IM(d5) = IM(c2) - RE(c5); |
325 RE(ch[ah+l1*ido]) = MUL_R_C(RE(d2),RE(wa1[i])) - MUL_R_C(IM(d2),IM(wa1[i]))*isign; | 902 RE(d3) = RE(c3) - IM(c4); |
326 IM(ch[ah+l1*ido]) = MUL_R_C(IM(d2),RE(wa1[i])) + MUL_R_C(RE(d2),IM(wa1[i]))*isign; | 903 IM(d4) = IM(c3) - RE(c4); |
327 RE(ch[ah+2*l1*ido]) = MUL_R_C(RE(d3),RE(wa2[i])) - MUL_R_C(IM(d3),IM(wa2[i]))*isign; | 904 |
328 IM(ch[ah+2*l1*ido]) = MUL_R_C(IM(d3),RE(wa2[i])) + MUL_R_C(RE(d3),IM(wa2[i]))*isign; | 905 #if 1 |
329 RE(ch[ah+3*l1*ido]) = MUL_R_C(RE(d4),RE(wa3[i])) - MUL_R_C(IM(d4),IM(wa3[i]))*isign; | 906 ComplexMult(&IM(ch[ah+l1*ido]), &RE(ch[ah+l1*ido]), |
330 IM(ch[ah+3*l1*ido]) = MUL_R_C(IM(d4),RE(wa3[i])) + MUL_R_C(RE(d4),IM(wa3[i]))*isign; | 907 IM(d2), RE(d2), RE(wa1[i]), IM(wa1[i])); |
331 RE(ch[ah+4*l1*ido]) = MUL_R_C(RE(d5),RE(wa4[i])) - MUL_R_C(IM(d5),IM(wa4[i]))*isign; | 908 ComplexMult(&IM(ch[ah+2*l1*ido]), &RE(ch[ah+2*l1*ido]), |
332 IM(ch[ah+4*l1*ido]) = MUL_R_C(IM(d5),RE(wa4[i])) + MUL_R_C(RE(d5),IM(wa4[i]))*isign; | 909 IM(d3), RE(d3), RE(wa2[i]), IM(wa2[i])); |
910 ComplexMult(&IM(ch[ah+3*l1*ido]), &RE(ch[ah+3*l1*ido]), | |
911 IM(d4), RE(d4), RE(wa3[i]), IM(wa3[i])); | |
912 ComplexMult(&IM(ch[ah+4*l1*ido]), &RE(ch[ah+4*l1*ido]), | |
913 IM(d5), RE(d5), RE(wa4[i]), IM(wa4[i])); | |
914 #else | |
915 ComplexMult(&RE(ch[ah+l1*ido]), &IM(ch[ah+l1*ido]), | |
916 RE(d2), IM(d2), RE(wa1[i]), IM(wa1[i])); | |
917 ComplexMult(&RE(ch[ah+2*l1*ido]), &IM(ch[ah+2*l1*ido]), | |
918 RE(d3), IM(d3), RE(wa2[i]), IM(wa2[i])); | |
919 ComplexMult(&RE(ch[ah+3*l1*ido]), &IM(ch[ah+3*l1*ido]), | |
920 RE(d4), IM(d4), RE(wa3[i]), IM(wa3[i])); | |
921 ComplexMult(&RE(ch[ah+4*l1*ido]), &IM(ch[ah+4*l1*ido]), | |
922 RE(d5), IM(d5), RE(wa4[i]), IM(wa4[i])); | |
923 #endif | |
924 } | |
925 } | |
926 } else { | |
927 for (k = 0; k < l1; k++) | |
928 { | |
929 for (i = 0; i < ido; i++) | |
930 { | |
931 ac = i + (k*5 + 1) * ido; | |
932 ah = i + k * ido; | |
933 | |
934 RE(t2) = RE(cc[ac]) + RE(cc[ac+3*ido]); | |
935 IM(t2) = IM(cc[ac]) + IM(cc[ac+3*ido]); | |
936 RE(t3) = RE(cc[ac+ido]) + RE(cc[ac+2*ido]); | |
937 IM(t3) = IM(cc[ac+ido]) + IM(cc[ac+2*ido]); | |
938 RE(t4) = RE(cc[ac+ido]) - RE(cc[ac+2*ido]); | |
939 IM(t4) = IM(cc[ac+ido]) - IM(cc[ac+2*ido]); | |
940 RE(t5) = RE(cc[ac]) - RE(cc[ac+3*ido]); | |
941 IM(t5) = IM(cc[ac]) - IM(cc[ac+3*ido]); | |
942 | |
943 RE(ch[ah]) = RE(cc[ac-ido]) + RE(t2) + RE(t3); | |
944 IM(ch[ah]) = IM(cc[ac-ido]) + IM(t2) + IM(t3); | |
945 | |
946 RE(c2) = RE(cc[ac-ido]) + MUL_F(RE(t2),tr11) + MUL_F(RE(t3),tr12); | |
947 IM(c2) = IM(cc[ac-ido]) + MUL_F(IM(t2),tr11) + MUL_F(IM(t3),tr12); | |
948 RE(c3) = RE(cc[ac-ido]) + MUL_F(RE(t2),tr12) + MUL_F(RE(t3),tr11); | |
949 IM(c3) = IM(cc[ac-ido]) + MUL_F(IM(t2),tr12) + MUL_F(IM(t3),tr11); | |
950 | |
951 ComplexMult(&RE(c4), &RE(c5), | |
952 ti12, ti11, RE(t5), RE(t4)); | |
953 ComplexMult(&IM(c4), &IM(c5), | |
954 ti12, ti12, IM(t5), IM(t4)); | |
955 | |
956 IM(d2) = IM(c2) - RE(c5); | |
957 IM(d3) = IM(c3) - RE(c4); | |
958 RE(d4) = RE(c3) - IM(c4); | |
959 RE(d5) = RE(c2) - IM(c5); | |
960 RE(d2) = RE(c2) + IM(c5); | |
961 IM(d5) = IM(c2) + RE(c5); | |
962 RE(d3) = RE(c3) + IM(c4); | |
963 IM(d4) = IM(c3) + RE(c4); | |
964 | |
965 #if 1 | |
966 ComplexMult(&RE(ch[ah+l1*ido]), &IM(ch[ah+l1*ido]), | |
967 RE(d2), IM(d2), RE(wa1[i]), IM(wa1[i])); | |
968 ComplexMult(&RE(ch[ah+2*l1*ido]), &IM(ch[ah+2*l1*ido]), | |
969 RE(d3), IM(d3), RE(wa2[i]), IM(wa2[i])); | |
970 ComplexMult(&RE(ch[ah+3*l1*ido]), &IM(ch[ah+3*l1*ido]), | |
971 RE(d4), IM(d4), RE(wa3[i]), IM(wa3[i])); | |
972 ComplexMult(&RE(ch[ah+4*l1*ido]), &IM(ch[ah+4*l1*ido]), | |
973 RE(d5), IM(d5), RE(wa4[i]), IM(wa4[i])); | |
974 #else | |
975 ComplexMult(&IM(ch[ah+l1*ido]), &RE(ch[ah+l1*ido]), | |
976 IM(d2), RE(d2), RE(wa1[i]), IM(wa1[i])); | |
977 ComplexMult(&IM(ch[ah+2*l1*ido]), &RE(ch[ah+2*l1*ido]), | |
978 IM(d3), RE(d3), RE(wa2[i]), IM(wa2[i])); | |
979 ComplexMult(&IM(ch[ah+3*l1*ido]), &RE(ch[ah+3*l1*ido]), | |
980 IM(d4), RE(d4), RE(wa3[i]), IM(wa3[i])); | |
981 ComplexMult(&IM(ch[ah+4*l1*ido]), &RE(ch[ah+4*l1*ido]), | |
982 IM(d5), RE(d5), RE(wa4[i]), IM(wa4[i])); | |
983 #endif | |
984 } | |
333 } | 985 } |
334 } | 986 } |
335 } | 987 } |
336 } | 988 } |
337 | 989 |
338 | 990 |
339 /*---------------------------------------------------------------------- | 991 /*---------------------------------------------------------------------- |
340 cfftf1, cfftf, cfftb, cffti1, cffti. Complex FFTs. | 992 cfftf1, cfftf, cfftb, cffti1, cffti. Complex FFTs. |
341 ----------------------------------------------------------------------*/ | 993 ----------------------------------------------------------------------*/ |
342 | 994 |
343 INLINE void cfftf1(uint16_t n, complex_t *c, complex_t *ch, | 995 #ifdef USE_SSE |
344 uint16_t *ifac, complex_t *wa, int8_t isign) | 996 |
997 #define CONV(A,B,C) ( (A<<2) | ((B & 0x1)<<1) | ((C==1)&0x1) ) | |
998 | |
999 static INLINE void cfftf1pos_sse(uint16_t n, complex_t *c, complex_t *ch, | |
1000 const uint16_t *ifac, const complex_t *wa, | |
1001 const int8_t isign) | |
345 { | 1002 { |
346 uint16_t i; | 1003 uint16_t i; |
347 uint16_t k1, l1, l2; | 1004 uint16_t k1, l1, l2; |
348 uint16_t na, nf, ip, iw, ix2, ix3, ix4, ido, idl1; | 1005 uint16_t na, nf, ip, iw, ix2, ix3, ix4, ido, idl1; |
349 | 1006 |
357 ip = ifac[k1]; | 1014 ip = ifac[k1]; |
358 l2 = ip*l1; | 1015 l2 = ip*l1; |
359 ido = n / l2; | 1016 ido = n / l2; |
360 idl1 = ido*l1; | 1017 idl1 = ido*l1; |
361 | 1018 |
1019 ix2 = iw + ido; | |
1020 ix3 = ix2 + ido; | |
1021 ix4 = ix3 + ido; | |
1022 | |
1023 switch (CONV(ip,na,ido)) | |
1024 { | |
1025 case CONV(4,0,0): | |
1026 //passf4pos_sse_ido((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw], &wa[ix2], &wa[ix3]); | |
1027 passf4pos((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw], &wa[ix2], &wa[ix3]); | |
1028 break; | |
1029 case CONV(4,0,1): | |
1030 passf4pos_sse((const uint16_t)l1, (const complex_t*)c, ch, &wa[iw], &wa[ix2], &wa[ix3]); | |
1031 break; | |
1032 case CONV(4,1,0): | |
1033 passf4pos((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw], &wa[ix2], &wa[ix3]); | |
1034 //passf4pos_sse_ido((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw], &wa[ix2], &wa[ix3]); | |
1035 break; | |
1036 case CONV(4,1,1): | |
1037 passf4pos_sse((const uint16_t)l1, (const complex_t*)ch, c, &wa[iw], &wa[ix2], &wa[ix3]); | |
1038 break; | |
1039 case CONV(2,0,0): | |
1040 passf2pos((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw]); | |
1041 //passf2pos_sse_ido((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw]); | |
1042 break; | |
1043 case CONV(2,0,1): | |
1044 passf2pos((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw]); | |
1045 //passf2pos_sse((const uint16_t)l1, (const complex_t*)c, ch, &wa[iw]); | |
1046 break; | |
1047 case CONV(2,1,0): | |
1048 passf2pos((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw]); | |
1049 //passf2pos_sse_ido((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw]); | |
1050 break; | |
1051 case CONV(2,1,1): | |
1052 passf2pos((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw]); | |
1053 //passf2pos_sse((const uint16_t)l1, (const complex_t*)ch, c, &wa[iw]); | |
1054 break; | |
1055 case CONV(3,0,0): | |
1056 case CONV(3,0,1): | |
1057 passf3((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw], &wa[ix2], isign); | |
1058 break; | |
1059 case CONV(3,1,0): | |
1060 case CONV(3,1,1): | |
1061 passf3((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw], &wa[ix2], isign); | |
1062 break; | |
1063 case CONV(5,0,0): | |
1064 case CONV(5,0,1): | |
1065 passf5((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4], isign); | |
1066 break; | |
1067 case CONV(5,1,0): | |
1068 case CONV(5,1,1): | |
1069 passf5((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4], isign); | |
1070 break; | |
1071 } | |
1072 | |
1073 na = 1 - na; | |
1074 | |
1075 l1 = l2; | |
1076 iw += (ip-1) * ido; | |
1077 } | |
1078 | |
1079 if (na == 0) | |
1080 return; | |
1081 | |
1082 for (i = 0; i < n; i++) | |
1083 { | |
1084 RE(c[i]) = RE(ch[i]); | |
1085 IM(c[i]) = IM(ch[i]); | |
1086 } | |
1087 } | |
1088 #endif | |
1089 | |
1090 static INLINE void cfftf1pos(uint16_t n, complex_t *c, complex_t *ch, | |
1091 const uint16_t *ifac, const complex_t *wa, | |
1092 const int8_t isign) | |
1093 { | |
1094 uint16_t i; | |
1095 uint16_t k1, l1, l2; | |
1096 uint16_t na, nf, ip, iw, ix2, ix3, ix4, ido, idl1; | |
1097 | |
1098 nf = ifac[1]; | |
1099 na = 0; | |
1100 l1 = 1; | |
1101 iw = 0; | |
1102 | |
1103 for (k1 = 2; k1 <= nf+1; k1++) | |
1104 { | |
1105 ip = ifac[k1]; | |
1106 l2 = ip*l1; | |
1107 ido = n / l2; | |
1108 idl1 = ido*l1; | |
1109 | |
362 switch (ip) | 1110 switch (ip) |
363 { | 1111 { |
364 case 4: | 1112 case 4: |
365 ix2 = iw + ido; | 1113 ix2 = iw + ido; |
366 ix3 = ix2 + ido; | 1114 ix3 = ix2 + ido; |
367 | 1115 |
368 if (na == 0) | 1116 if (na == 0) |
369 passf4(ido, l1, c, ch, &wa[iw], &wa[ix2], &wa[ix3], isign); | 1117 passf4pos((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw], &wa[ix2], &wa[ix3]); |
370 else | 1118 else |
371 passf4(ido, l1, ch, c, &wa[iw], &wa[ix2], &wa[ix3], isign); | 1119 passf4pos((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw], &wa[ix2], &wa[ix3]); |
372 | 1120 |
373 na = 1 - na; | 1121 na = 1 - na; |
374 break; | 1122 break; |
375 case 2: | 1123 case 2: |
376 if (na == 0) | 1124 if (na == 0) |
377 passf2(ido, l1, c, ch, &wa[iw], isign); | 1125 passf2pos((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw]); |
378 else | 1126 else |
379 passf2(ido, l1, ch, c, &wa[iw], isign); | 1127 passf2pos((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw]); |
380 | 1128 |
381 na = 1 - na; | 1129 na = 1 - na; |
382 break; | 1130 break; |
383 case 3: | 1131 case 3: |
384 ix2 = iw + ido; | 1132 ix2 = iw + ido; |
385 | 1133 |
386 if (na == 0) | 1134 if (na == 0) |
387 passf3(ido, l1, c, ch, &wa[iw], &wa[ix2], isign); | 1135 passf3((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw], &wa[ix2], isign); |
388 else | 1136 else |
389 passf3(ido, l1, ch, c, &wa[iw], &wa[ix2], isign); | 1137 passf3((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw], &wa[ix2], isign); |
390 | 1138 |
391 na = 1 - na; | 1139 na = 1 - na; |
392 break; | 1140 break; |
393 case 5: | 1141 case 5: |
394 ix2 = iw + ido; | 1142 ix2 = iw + ido; |
395 ix3 = ix2 + ido; | 1143 ix3 = ix2 + ido; |
396 ix4 = ix3 + ido; | 1144 ix4 = ix3 + ido; |
397 | 1145 |
398 if (na == 0) | 1146 if (na == 0) |
399 passf5(ido, l1, c, ch, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4], isign); | 1147 passf5((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4], isign); |
400 else | 1148 else |
401 passf5(ido, l1, ch, c, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4], isign); | 1149 passf5((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4], isign); |
402 | 1150 |
403 na = 1 - na; | 1151 na = 1 - na; |
404 break; | 1152 break; |
405 } | 1153 } |
406 | 1154 |
416 RE(c[i]) = RE(ch[i]); | 1164 RE(c[i]) = RE(ch[i]); |
417 IM(c[i]) = IM(ch[i]); | 1165 IM(c[i]) = IM(ch[i]); |
418 } | 1166 } |
419 } | 1167 } |
420 | 1168 |
1169 static INLINE void cfftf1neg(uint16_t n, complex_t *c, complex_t *ch, | |
1170 const uint16_t *ifac, const complex_t *wa, | |
1171 const int8_t isign) | |
1172 { | |
1173 uint16_t i; | |
1174 uint16_t k1, l1, l2; | |
1175 uint16_t na, nf, ip, iw, ix2, ix3, ix4, ido, idl1; | |
1176 | |
1177 nf = ifac[1]; | |
1178 na = 0; | |
1179 l1 = 1; | |
1180 iw = 0; | |
1181 | |
1182 for (k1 = 2; k1 <= nf+1; k1++) | |
1183 { | |
1184 ip = ifac[k1]; | |
1185 l2 = ip*l1; | |
1186 ido = n / l2; | |
1187 idl1 = ido*l1; | |
1188 | |
1189 switch (ip) | |
1190 { | |
1191 case 4: | |
1192 ix2 = iw + ido; | |
1193 ix3 = ix2 + ido; | |
1194 | |
1195 if (na == 0) | |
1196 passf4neg((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw], &wa[ix2], &wa[ix3]); | |
1197 else | |
1198 passf4neg((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw], &wa[ix2], &wa[ix3]); | |
1199 | |
1200 na = 1 - na; | |
1201 break; | |
1202 case 2: | |
1203 if (na == 0) | |
1204 passf2neg((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw]); | |
1205 else | |
1206 passf2neg((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw]); | |
1207 | |
1208 na = 1 - na; | |
1209 break; | |
1210 case 3: | |
1211 ix2 = iw + ido; | |
1212 | |
1213 if (na == 0) | |
1214 passf3((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw], &wa[ix2], isign); | |
1215 else | |
1216 passf3((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw], &wa[ix2], isign); | |
1217 | |
1218 na = 1 - na; | |
1219 break; | |
1220 case 5: | |
1221 ix2 = iw + ido; | |
1222 ix3 = ix2 + ido; | |
1223 ix4 = ix3 + ido; | |
1224 | |
1225 if (na == 0) | |
1226 passf5((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4], isign); | |
1227 else | |
1228 passf5((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4], isign); | |
1229 | |
1230 na = 1 - na; | |
1231 break; | |
1232 } | |
1233 | |
1234 l1 = l2; | |
1235 iw += (ip-1) * ido; | |
1236 } | |
1237 | |
1238 if (na == 0) | |
1239 return; | |
1240 | |
1241 for (i = 0; i < n; i++) | |
1242 { | |
1243 RE(c[i]) = RE(ch[i]); | |
1244 IM(c[i]) = IM(ch[i]); | |
1245 } | |
1246 } | |
1247 | |
421 void cfftf(cfft_info *cfft, complex_t *c) | 1248 void cfftf(cfft_info *cfft, complex_t *c) |
422 { | 1249 { |
423 cfftf1(cfft->n, c, cfft->work, cfft->ifac, cfft->tab, -1); | 1250 cfftf1neg(cfft->n, c, cfft->work, (const uint16_t*)cfft->ifac, (const complex_t*)cfft->tab, -1); |
424 } | 1251 } |
425 | 1252 |
426 void cfftb(cfft_info *cfft, complex_t *c) | 1253 void cfftb(cfft_info *cfft, complex_t *c) |
427 { | 1254 { |
428 cfftf1(cfft->n, c, cfft->work, cfft->ifac, cfft->tab, +1); | 1255 cfftf1pos(cfft->n, c, cfft->work, (const uint16_t*)cfft->ifac, (const complex_t*)cfft->tab, +1); |
429 } | 1256 } |
1257 | |
1258 #ifdef USE_SSE | |
1259 void cfftb_sse(cfft_info *cfft, complex_t *c) | |
1260 { | |
1261 cfftf1pos_sse(cfft->n, c, cfft->work, (const uint16_t*)cfft->ifac, (const complex_t*)cfft->tab, +1); | |
1262 } | |
1263 #endif | |
430 | 1264 |
431 static void cffti1(uint16_t n, complex_t *wa, uint16_t *ifac) | 1265 static void cffti1(uint16_t n, complex_t *wa, uint16_t *ifac) |
432 { | 1266 { |
433 static uint16_t ntryh[4] = {3, 4, 2, 5}; | 1267 static uint16_t ntryh[4] = {3, 4, 2, 5}; |
434 #ifndef FIXED_POINT | 1268 #ifndef FIXED_POINT |
435 real_t arg, argh, argld, fi; | 1269 real_t arg, argh, argld, fi; |
436 uint16_t ido, ipm; | 1270 uint16_t ido, ipm; |
437 uint16_t i1, k1, l1, l2; | 1271 uint16_t i1, k1, l1, l2; |
438 uint16_t ld, ii, ip; | 1272 uint16_t ld, ii, ip; |
439 #endif | 1273 #endif |
440 uint16_t ntry, i, j; | 1274 uint16_t ntry = 0, i, j; |
441 uint16_t ib; | 1275 uint16_t ib; |
442 uint16_t nf, nl, nq, nr; | 1276 uint16_t nf, nl, nq, nr; |
443 | 1277 |
444 nl = n; | 1278 nl = n; |
445 nf = 0; | 1279 nf = 0; |
478 | 1312 |
479 ifac[0] = n; | 1313 ifac[0] = n; |
480 ifac[1] = nf; | 1314 ifac[1] = nf; |
481 | 1315 |
482 #ifndef FIXED_POINT | 1316 #ifndef FIXED_POINT |
483 argh = (real_t)2.0*M_PI / (real_t)n; | 1317 argh = (real_t)2.0*(real_t)M_PI / (real_t)n; |
484 i = 0; | 1318 i = 0; |
485 l1 = 1; | 1319 l1 = 1; |
486 | 1320 |
487 for (k1 = 1; k1 <= nf; k1++) | 1321 for (k1 = 1; k1 <= nf; k1++) |
488 { | 1322 { |
505 { | 1339 { |
506 i++; | 1340 i++; |
507 fi++; | 1341 fi++; |
508 arg = fi * argld; | 1342 arg = fi * argld; |
509 RE(wa[i]) = (real_t)cos(arg); | 1343 RE(wa[i]) = (real_t)cos(arg); |
1344 #if 1 | |
510 IM(wa[i]) = (real_t)sin(arg); | 1345 IM(wa[i]) = (real_t)sin(arg); |
1346 #else | |
1347 IM(wa[i]) = (real_t)-sin(arg); | |
1348 #endif | |
511 } | 1349 } |
512 | 1350 |
513 if (ip > 5) | 1351 if (ip > 5) |
514 { | 1352 { |
515 RE(wa[i1]) = RE(wa[i]); | 1353 RE(wa[i1]) = RE(wa[i]); |
521 #endif | 1359 #endif |
522 } | 1360 } |
523 | 1361 |
524 cfft_info *cffti(uint16_t n) | 1362 cfft_info *cffti(uint16_t n) |
525 { | 1363 { |
526 cfft_info *cfft = (cfft_info*)malloc(sizeof(cfft_info)); | 1364 cfft_info *cfft = (cfft_info*)faad_malloc(sizeof(cfft_info)); |
527 | 1365 |
528 cfft->n = n; | 1366 cfft->n = n; |
529 cfft->work = (complex_t*)malloc(n*sizeof(complex_t)); | 1367 cfft->work = (complex_t*)faad_malloc(n*sizeof(complex_t)); |
530 | 1368 |
531 #ifndef FIXED_POINT | 1369 #ifndef FIXED_POINT |
532 cfft->tab = (complex_t*)malloc(n*sizeof(complex_t)); | 1370 cfft->tab = (complex_t*)faad_malloc(n*sizeof(complex_t)); |
533 | 1371 |
534 cffti1(n, cfft->tab, cfft->ifac); | 1372 cffti1(n, cfft->tab, cfft->ifac); |
535 #else | 1373 #else |
536 cffti1(n, NULL, cfft->ifac); | 1374 cffti1(n, NULL, cfft->ifac); |
537 | 1375 |
538 switch (n) | 1376 switch (n) |
539 { | 1377 { |
540 case 60: cfft->tab = cfft_tab_60; break; | |
541 case 64: cfft->tab = cfft_tab_64; break; | 1378 case 64: cfft->tab = cfft_tab_64; break; |
542 case 480: cfft->tab = cfft_tab_480; break; | |
543 case 512: cfft->tab = cfft_tab_512; break; | 1379 case 512: cfft->tab = cfft_tab_512; break; |
544 #ifdef LD_DEC | 1380 #ifdef LD_DEC |
1381 case 256: cfft->tab = cfft_tab_256; break; | |
1382 #endif | |
1383 | |
1384 #ifdef ALLOW_SMALL_FRAMELENGTH | |
1385 case 60: cfft->tab = cfft_tab_60; break; | |
1386 case 480: cfft->tab = cfft_tab_480; break; | |
1387 #ifdef LD_DEC | |
545 case 240: cfft->tab = cfft_tab_240; break; | 1388 case 240: cfft->tab = cfft_tab_240; break; |
546 case 256: cfft->tab = cfft_tab_256; break; | 1389 #endif |
547 #endif | 1390 #endif |
548 } | 1391 } |
549 #endif | 1392 #endif |
550 | 1393 |
551 return cfft; | 1394 return cfft; |
552 } | 1395 } |
553 | 1396 |
554 void cfftu(cfft_info *cfft) | 1397 void cfftu(cfft_info *cfft) |
555 { | 1398 { |
556 if (cfft->work) free(cfft->work); | 1399 if (cfft->work) faad_free(cfft->work); |
557 #ifndef FIXED_POINT | 1400 #ifndef FIXED_POINT |
558 if (cfft->tab) free(cfft->tab); | 1401 if (cfft->tab) faad_free(cfft->tab); |
559 #endif | 1402 #endif |
560 | 1403 |
561 if (cfft) free(cfft); | 1404 if (cfft) faad_free(cfft); |
562 } | 1405 } |
1406 |