Mercurial > mplayer.hg
annotate libfaad2/cfft.c @ 29269:4d9de809b174
Add a hack to detect when we are writing into a Windows pipe since the fseek
incorrectly does not fail like it should.
This ensures we will not incorrectly append the file header at the end.
Based on patch by Zhou Zongyi [zhouzongyi at pset.suntec.net]
author | reimar |
---|---|
date | Sat, 16 May 2009 13:59:53 +0000 |
parents | 59b6fa5b4201 |
children |
rev | line source |
---|---|
10725 | 1 /* |
2 ** FAAD2 - Freeware Advanced Audio (AAC) Decoder including SBR decoding | |
12527 | 3 ** Copyright (C) 2003-2004 M. Bakker, Ahead Software AG, http://www.nero.com |
4 ** | |
10725 | 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 | |
7 ** the Free Software Foundation; either version 2 of the License, or | |
8 ** (at your option) any later version. | |
12527 | 9 ** |
10725 | 10 ** This program is distributed in the hope that it will be useful, |
11 ** but WITHOUT ANY WARRANTY; without even the implied warranty of | |
12 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | |
13 ** GNU General Public License for more details. | |
12527 | 14 ** |
10725 | 15 ** You should have received a copy of the GNU General Public License |
12527 | 16 ** along with this program; if not, write to the Free Software |
10725 | 17 ** Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. |
18 ** | |
19 ** Any non-GPL usage of this software or parts of this software is strictly | |
20 ** forbidden. | |
21 ** | |
22 ** Commercial non-GPL licensing of this software is possible. | |
23 ** For more info contact Ahead Software through Mpeg4AAClicense@nero.com. | |
24 ** | |
18141 | 25 ** $Id: cfft.c,v 1.30 2004/09/08 09:43:11 gcp Exp $ |
10725 | 26 **/ |
27 | |
28 /* | |
29 * Algorithmically based on Fortran-77 FFTPACK | |
30 * by Paul N. Swarztrauber(Version 4, 1985). | |
31 * | |
32 * Does even sized fft only | |
33 */ | |
34 | |
35 /* isign is +1 for backward and -1 for forward transforms */ | |
36 | |
37 #include "common.h" | |
38 #include "structs.h" | |
39 | |
40 #include <stdlib.h> | |
41 | |
42 #include "cfft.h" | |
43 #include "cfft_tab.h" | |
44 | |
45 | |
12527 | 46 /* static function declarations */ |
47 static void passf2pos(const uint16_t ido, const uint16_t l1, const complex_t *cc, | |
48 complex_t *ch, const complex_t *wa); | |
49 static void passf2neg(const uint16_t ido, const uint16_t l1, const complex_t *cc, | |
50 complex_t *ch, const complex_t *wa); | |
51 static void passf3(const uint16_t ido, const uint16_t l1, const complex_t *cc, | |
52 complex_t *ch, const complex_t *wa1, const complex_t *wa2, const int8_t isign); | |
53 static void passf4pos(const uint16_t ido, const uint16_t l1, const complex_t *cc, complex_t *ch, | |
54 const complex_t *wa1, const complex_t *wa2, const complex_t *wa3); | |
55 static void passf4neg(const uint16_t ido, const uint16_t l1, const complex_t *cc, complex_t *ch, | |
56 const complex_t *wa1, const complex_t *wa2, const complex_t *wa3); | |
57 static void passf5(const uint16_t ido, const uint16_t l1, const complex_t *cc, complex_t *ch, | |
58 const complex_t *wa1, const complex_t *wa2, const complex_t *wa3, | |
59 const complex_t *wa4, const int8_t isign); | |
60 INLINE void cfftf1(uint16_t n, complex_t *c, complex_t *ch, | |
61 const uint16_t *ifac, const complex_t *wa, const int8_t isign); | |
62 static void cffti1(uint16_t n, complex_t *wa, uint16_t *ifac); | |
63 | |
64 | |
10725 | 65 /*---------------------------------------------------------------------- |
66 passf2, passf3, passf4, passf5. Complex FFT passes fwd and bwd. | |
67 ----------------------------------------------------------------------*/ | |
68 | |
12527 | 69 static void passf2pos(const uint16_t ido, const uint16_t l1, const complex_t *cc, |
70 complex_t *ch, const complex_t *wa) | |
10725 | 71 { |
72 uint16_t i, k, ah, ac; | |
73 | |
74 if (ido == 1) | |
75 { | |
76 for (k = 0; k < l1; k++) | |
77 { | |
78 ah = 2*k; | |
79 ac = 4*k; | |
80 | |
12527 | 81 RE(ch[ah]) = RE(cc[ac]) + RE(cc[ac+1]); |
82 RE(ch[ah+l1]) = RE(cc[ac]) - RE(cc[ac+1]); | |
83 IM(ch[ah]) = IM(cc[ac]) + IM(cc[ac+1]); | |
84 IM(ch[ah+l1]) = IM(cc[ac]) - IM(cc[ac+1]); | |
85 } | |
86 } else { | |
87 for (k = 0; k < l1; k++) | |
88 { | |
89 ah = k*ido; | |
90 ac = 2*k*ido; | |
91 | |
92 for (i = 0; i < ido; i++) | |
93 { | |
94 complex_t t2; | |
95 | |
96 RE(ch[ah+i]) = RE(cc[ac+i]) + RE(cc[ac+i+ido]); | |
97 RE(t2) = RE(cc[ac+i]) - RE(cc[ac+i+ido]); | |
98 | |
99 IM(ch[ah+i]) = IM(cc[ac+i]) + IM(cc[ac+i+ido]); | |
100 IM(t2) = IM(cc[ac+i]) - IM(cc[ac+i+ido]); | |
101 | |
102 #if 1 | |
103 ComplexMult(&IM(ch[ah+i+l1*ido]), &RE(ch[ah+i+l1*ido]), | |
104 IM(t2), RE(t2), RE(wa[i]), IM(wa[i])); | |
105 #else | |
106 ComplexMult(&RE(ch[ah+i+l1*ido]), &IM(ch[ah+i+l1*ido]), | |
107 RE(t2), IM(t2), RE(wa[i]), IM(wa[i])); | |
108 #endif | |
109 } | |
110 } | |
111 } | |
112 } | |
113 | |
114 static void passf2neg(const uint16_t ido, const uint16_t l1, const complex_t *cc, | |
115 complex_t *ch, const complex_t *wa) | |
116 { | |
117 uint16_t i, k, ah, ac; | |
118 | |
119 if (ido == 1) | |
120 { | |
121 for (k = 0; k < l1; k++) | |
122 { | |
123 ah = 2*k; | |
124 ac = 4*k; | |
125 | |
126 RE(ch[ah]) = RE(cc[ac]) + RE(cc[ac+1]); | |
10725 | 127 RE(ch[ah+l1]) = RE(cc[ac]) - RE(cc[ac+1]); |
10989 | 128 IM(ch[ah]) = IM(cc[ac]) + IM(cc[ac+1]); |
10725 | 129 IM(ch[ah+l1]) = IM(cc[ac]) - IM(cc[ac+1]); |
130 } | |
131 } else { | |
132 for (k = 0; k < l1; k++) | |
133 { | |
134 ah = k*ido; | |
135 ac = 2*k*ido; | |
136 | |
137 for (i = 0; i < ido; i++) | |
138 { | |
139 complex_t t2; | |
140 | |
10989 | 141 RE(ch[ah+i]) = RE(cc[ac+i]) + RE(cc[ac+i+ido]); |
142 RE(t2) = RE(cc[ac+i]) - RE(cc[ac+i+ido]); | |
10725 | 143 |
10989 | 144 IM(ch[ah+i]) = IM(cc[ac+i]) + IM(cc[ac+i+ido]); |
145 IM(t2) = IM(cc[ac+i]) - IM(cc[ac+i+ido]); | |
10725 | 146 |
12527 | 147 #if 1 |
148 ComplexMult(&RE(ch[ah+i+l1*ido]), &IM(ch[ah+i+l1*ido]), | |
149 RE(t2), IM(t2), RE(wa[i]), IM(wa[i])); | |
150 #else | |
151 ComplexMult(&IM(ch[ah+i+l1*ido]), &RE(ch[ah+i+l1*ido]), | |
152 IM(t2), RE(t2), RE(wa[i]), IM(wa[i])); | |
153 #endif | |
10725 | 154 } |
155 } | |
156 } | |
157 } | |
158 | |
159 | |
12527 | 160 static void passf3(const uint16_t ido, const uint16_t l1, const complex_t *cc, |
161 complex_t *ch, const complex_t *wa1, const complex_t *wa2, | |
162 const int8_t isign) | |
10725 | 163 { |
12527 | 164 static real_t taur = FRAC_CONST(-0.5); |
165 static real_t taui = FRAC_CONST(0.866025403784439); | |
10725 | 166 uint16_t i, k, ac, ah; |
167 complex_t c2, c3, d2, d3, t2; | |
168 | |
169 if (ido == 1) | |
170 { | |
12527 | 171 if (isign == 1) |
10725 | 172 { |
12527 | 173 for (k = 0; k < l1; k++) |
174 { | |
175 ac = 3*k+1; | |
176 ah = k; | |
10725 | 177 |
12527 | 178 RE(t2) = RE(cc[ac]) + RE(cc[ac+1]); |
179 IM(t2) = IM(cc[ac]) + IM(cc[ac+1]); | |
180 RE(c2) = RE(cc[ac-1]) + MUL_F(RE(t2),taur); | |
181 IM(c2) = IM(cc[ac-1]) + MUL_F(IM(t2),taur); | |
182 | |
183 RE(ch[ah]) = RE(cc[ac-1]) + RE(t2); | |
184 IM(ch[ah]) = IM(cc[ac-1]) + IM(t2); | |
185 | |
186 RE(c3) = MUL_F((RE(cc[ac]) - RE(cc[ac+1])), taui); | |
187 IM(c3) = MUL_F((IM(cc[ac]) - IM(cc[ac+1])), taui); | |
10725 | 188 |
12527 | 189 RE(ch[ah+l1]) = RE(c2) - IM(c3); |
190 IM(ch[ah+l1]) = IM(c2) + RE(c3); | |
191 RE(ch[ah+2*l1]) = RE(c2) + IM(c3); | |
192 IM(ch[ah+2*l1]) = IM(c2) - RE(c3); | |
193 } | |
194 } else { | |
195 for (k = 0; k < l1; k++) | |
196 { | |
197 ac = 3*k+1; | |
198 ah = k; | |
10725 | 199 |
12527 | 200 RE(t2) = RE(cc[ac]) + RE(cc[ac+1]); |
201 IM(t2) = IM(cc[ac]) + IM(cc[ac+1]); | |
202 RE(c2) = RE(cc[ac-1]) + MUL_F(RE(t2),taur); | |
203 IM(c2) = IM(cc[ac-1]) + MUL_F(IM(t2),taur); | |
204 | |
205 RE(ch[ah]) = RE(cc[ac-1]) + RE(t2); | |
206 IM(ch[ah]) = IM(cc[ac-1]) + IM(t2); | |
10725 | 207 |
12527 | 208 RE(c3) = MUL_F((RE(cc[ac]) - RE(cc[ac+1])), taui); |
209 IM(c3) = MUL_F((IM(cc[ac]) - IM(cc[ac+1])), taui); | |
210 | |
211 RE(ch[ah+l1]) = RE(c2) + IM(c3); | |
212 IM(ch[ah+l1]) = IM(c2) - RE(c3); | |
213 RE(ch[ah+2*l1]) = RE(c2) - IM(c3); | |
214 IM(ch[ah+2*l1]) = IM(c2) + RE(c3); | |
215 } | |
10725 | 216 } |
217 } else { | |
12527 | 218 if (isign == 1) |
10725 | 219 { |
12527 | 220 for (k = 0; k < l1; k++) |
10725 | 221 { |
12527 | 222 for (i = 0; i < ido; i++) |
223 { | |
224 ac = i + (3*k+1)*ido; | |
225 ah = i + k * ido; | |
226 | |
227 RE(t2) = RE(cc[ac]) + RE(cc[ac+ido]); | |
228 RE(c2) = RE(cc[ac-ido]) + MUL_F(RE(t2),taur); | |
229 IM(t2) = IM(cc[ac]) + IM(cc[ac+ido]); | |
230 IM(c2) = IM(cc[ac-ido]) + MUL_F(IM(t2),taur); | |
10725 | 231 |
12527 | 232 RE(ch[ah]) = RE(cc[ac-ido]) + RE(t2); |
233 IM(ch[ah]) = IM(cc[ac-ido]) + IM(t2); | |
234 | |
235 RE(c3) = MUL_F((RE(cc[ac]) - RE(cc[ac+ido])), taui); | |
236 IM(c3) = MUL_F((IM(cc[ac]) - IM(cc[ac+ido])), taui); | |
237 | |
238 RE(d2) = RE(c2) - IM(c3); | |
239 IM(d3) = IM(c2) - RE(c3); | |
240 RE(d3) = RE(c2) + IM(c3); | |
241 IM(d2) = IM(c2) + RE(c3); | |
10725 | 242 |
12527 | 243 #if 1 |
244 ComplexMult(&IM(ch[ah+l1*ido]), &RE(ch[ah+l1*ido]), | |
245 IM(d2), RE(d2), RE(wa1[i]), IM(wa1[i])); | |
246 ComplexMult(&IM(ch[ah+2*l1*ido]), &RE(ch[ah+2*l1*ido]), | |
247 IM(d3), RE(d3), RE(wa2[i]), IM(wa2[i])); | |
248 #else | |
249 ComplexMult(&RE(ch[ah+l1*ido]), &IM(ch[ah+l1*ido]), | |
250 RE(d2), IM(d2), RE(wa1[i]), IM(wa1[i])); | |
251 ComplexMult(&RE(ch[ah+2*l1*ido]), &IM(ch[ah+2*l1*ido]), | |
252 RE(d3), IM(d3), RE(wa2[i]), IM(wa2[i])); | |
253 #endif | |
254 } | |
255 } | |
256 } else { | |
257 for (k = 0; k < l1; k++) | |
258 { | |
259 for (i = 0; i < ido; i++) | |
260 { | |
261 ac = i + (3*k+1)*ido; | |
262 ah = i + k * ido; | |
10725 | 263 |
12527 | 264 RE(t2) = RE(cc[ac]) + RE(cc[ac+ido]); |
265 RE(c2) = RE(cc[ac-ido]) + MUL_F(RE(t2),taur); | |
266 IM(t2) = IM(cc[ac]) + IM(cc[ac+ido]); | |
267 IM(c2) = IM(cc[ac-ido]) + MUL_F(IM(t2),taur); | |
268 | |
269 RE(ch[ah]) = RE(cc[ac-ido]) + RE(t2); | |
270 IM(ch[ah]) = IM(cc[ac-ido]) + IM(t2); | |
271 | |
272 RE(c3) = MUL_F((RE(cc[ac]) - RE(cc[ac+ido])), taui); | |
273 IM(c3) = MUL_F((IM(cc[ac]) - IM(cc[ac+ido])), taui); | |
10725 | 274 |
12527 | 275 RE(d2) = RE(c2) + IM(c3); |
276 IM(d3) = IM(c2) + RE(c3); | |
277 RE(d3) = RE(c2) - IM(c3); | |
278 IM(d2) = IM(c2) - RE(c3); | |
279 | |
280 #if 1 | |
281 ComplexMult(&RE(ch[ah+l1*ido]), &IM(ch[ah+l1*ido]), | |
282 RE(d2), IM(d2), RE(wa1[i]), IM(wa1[i])); | |
283 ComplexMult(&RE(ch[ah+2*l1*ido]), &IM(ch[ah+2*l1*ido]), | |
284 RE(d3), IM(d3), RE(wa2[i]), IM(wa2[i])); | |
285 #else | |
286 ComplexMult(&IM(ch[ah+l1*ido]), &RE(ch[ah+l1*ido]), | |
287 IM(d2), RE(d2), RE(wa1[i]), IM(wa1[i])); | |
288 ComplexMult(&IM(ch[ah+2*l1*ido]), &RE(ch[ah+2*l1*ido]), | |
289 IM(d3), RE(d3), RE(wa2[i]), IM(wa2[i])); | |
290 #endif | |
291 } | |
10725 | 292 } |
293 } | |
294 } | |
295 } | |
296 | |
12527 | 297 |
298 static void passf4pos(const uint16_t ido, const uint16_t l1, const complex_t *cc, | |
299 complex_t *ch, const complex_t *wa1, const complex_t *wa2, | |
300 const complex_t *wa3) | |
10725 | 301 { |
302 uint16_t i, k, ac, ah; | |
303 | |
304 if (ido == 1) | |
305 { | |
306 for (k = 0; k < l1; k++) | |
307 { | |
10989 | 308 complex_t t1, t2, t3, t4; |
309 | |
10725 | 310 ac = 4*k; |
311 ah = k; | |
312 | |
12527 | 313 RE(t2) = RE(cc[ac]) + RE(cc[ac+2]); |
314 RE(t1) = RE(cc[ac]) - RE(cc[ac+2]); | |
10989 | 315 IM(t2) = IM(cc[ac]) + IM(cc[ac+2]); |
12527 | 316 IM(t1) = IM(cc[ac]) - IM(cc[ac+2]); |
10725 | 317 RE(t3) = RE(cc[ac+1]) + RE(cc[ac+3]); |
10989 | 318 IM(t4) = RE(cc[ac+1]) - RE(cc[ac+3]); |
319 IM(t3) = IM(cc[ac+3]) + IM(cc[ac+1]); | |
10725 | 320 RE(t4) = IM(cc[ac+3]) - IM(cc[ac+1]); |
321 | |
12527 | 322 RE(ch[ah]) = RE(t2) + RE(t3); |
10725 | 323 RE(ch[ah+2*l1]) = RE(t2) - RE(t3); |
10989 | 324 |
325 IM(ch[ah]) = IM(t2) + IM(t3); | |
10725 | 326 IM(ch[ah+2*l1]) = IM(t2) - IM(t3); |
10989 | 327 |
12527 | 328 RE(ch[ah+l1]) = RE(t1) + RE(t4); |
329 RE(ch[ah+3*l1]) = RE(t1) - RE(t4); | |
10989 | 330 |
12527 | 331 IM(ch[ah+l1]) = IM(t1) + IM(t4); |
332 IM(ch[ah+3*l1]) = IM(t1) - IM(t4); | |
10725 | 333 } |
334 } else { | |
335 for (k = 0; k < l1; k++) | |
336 { | |
10989 | 337 ac = 4*k*ido; |
338 ah = k*ido; | |
339 | |
10725 | 340 for (i = 0; i < ido; i++) |
341 { | |
10989 | 342 complex_t c2, c3, c4, t1, t2, t3, t4; |
10725 | 343 |
10989 | 344 RE(t2) = RE(cc[ac+i]) + RE(cc[ac+i+2*ido]); |
345 RE(t1) = RE(cc[ac+i]) - RE(cc[ac+i+2*ido]); | |
346 IM(t2) = IM(cc[ac+i]) + IM(cc[ac+i+2*ido]); | |
347 IM(t1) = IM(cc[ac+i]) - IM(cc[ac+i+2*ido]); | |
348 RE(t3) = RE(cc[ac+i+ido]) + RE(cc[ac+i+3*ido]); | |
349 IM(t4) = RE(cc[ac+i+ido]) - RE(cc[ac+i+3*ido]); | |
350 IM(t3) = IM(cc[ac+i+3*ido]) + IM(cc[ac+i+ido]); | |
351 RE(t4) = IM(cc[ac+i+3*ido]) - IM(cc[ac+i+ido]); | |
10725 | 352 |
12527 | 353 RE(c2) = RE(t1) + RE(t4); |
354 RE(c4) = RE(t1) - RE(t4); | |
10989 | 355 |
12527 | 356 IM(c2) = IM(t1) + IM(t4); |
357 IM(c4) = IM(t1) - IM(t4); | |
10725 | 358 |
10989 | 359 RE(ch[ah+i]) = RE(t2) + RE(t3); |
12527 | 360 RE(c3) = RE(t2) - RE(t3); |
10989 | 361 |
362 IM(ch[ah+i]) = IM(t2) + IM(t3); | |
12527 | 363 IM(c3) = IM(t2) - IM(t3); |
10989 | 364 |
12527 | 365 #if 1 |
366 ComplexMult(&IM(ch[ah+i+l1*ido]), &RE(ch[ah+i+l1*ido]), | |
367 IM(c2), RE(c2), RE(wa1[i]), IM(wa1[i])); | |
368 ComplexMult(&IM(ch[ah+i+2*l1*ido]), &RE(ch[ah+i+2*l1*ido]), | |
369 IM(c3), RE(c3), RE(wa2[i]), IM(wa2[i])); | |
370 ComplexMult(&IM(ch[ah+i+3*l1*ido]), &RE(ch[ah+i+3*l1*ido]), | |
371 IM(c4), RE(c4), RE(wa3[i]), IM(wa3[i])); | |
372 #else | |
373 ComplexMult(&RE(ch[ah+i+l1*ido]), &IM(ch[ah+i+l1*ido]), | |
374 RE(c2), IM(c2), RE(wa1[i]), IM(wa1[i])); | |
375 ComplexMult(&RE(ch[ah+i+2*l1*ido]), &IM(ch[ah+i+2*l1*ido]), | |
376 RE(c3), IM(c3), RE(wa2[i]), IM(wa2[i])); | |
377 ComplexMult(&RE(ch[ah+i+3*l1*ido]), &IM(ch[ah+i+3*l1*ido]), | |
378 RE(c4), IM(c4), RE(wa3[i]), IM(wa3[i])); | |
379 #endif | |
10725 | 380 } |
381 } | |
382 } | |
383 } | |
384 | |
12527 | 385 static void passf4neg(const uint16_t ido, const uint16_t l1, const complex_t *cc, |
386 complex_t *ch, const complex_t *wa1, const complex_t *wa2, | |
387 const complex_t *wa3) | |
10725 | 388 { |
389 uint16_t i, k, ac, ah; | |
390 | |
391 if (ido == 1) | |
392 { | |
393 for (k = 0; k < l1; k++) | |
394 { | |
12527 | 395 complex_t t1, t2, t3, t4; |
396 | |
397 ac = 4*k; | |
10725 | 398 ah = k; |
399 | |
12527 | 400 RE(t2) = RE(cc[ac]) + RE(cc[ac+2]); |
401 RE(t1) = RE(cc[ac]) - RE(cc[ac+2]); | |
402 IM(t2) = IM(cc[ac]) + IM(cc[ac+2]); | |
403 IM(t1) = IM(cc[ac]) - IM(cc[ac+2]); | |
404 RE(t3) = RE(cc[ac+1]) + RE(cc[ac+3]); | |
405 IM(t4) = RE(cc[ac+1]) - RE(cc[ac+3]); | |
406 IM(t3) = IM(cc[ac+3]) + IM(cc[ac+1]); | |
407 RE(t4) = IM(cc[ac+3]) - IM(cc[ac+1]); | |
10725 | 408 |
12527 | 409 RE(ch[ah]) = RE(t2) + RE(t3); |
410 RE(ch[ah+2*l1]) = RE(t2) - RE(t3); | |
411 | |
412 IM(ch[ah]) = IM(t2) + IM(t3); | |
413 IM(ch[ah+2*l1]) = IM(t2) - IM(t3); | |
10725 | 414 |
12527 | 415 RE(ch[ah+l1]) = RE(t1) - RE(t4); |
416 RE(ch[ah+3*l1]) = RE(t1) + RE(t4); | |
417 | |
418 IM(ch[ah+l1]) = IM(t1) - IM(t4); | |
419 IM(ch[ah+3*l1]) = IM(t1) + IM(t4); | |
10725 | 420 } |
421 } else { | |
422 for (k = 0; k < l1; k++) | |
423 { | |
12527 | 424 ac = 4*k*ido; |
425 ah = k*ido; | |
426 | |
10725 | 427 for (i = 0; i < ido; i++) |
428 { | |
12527 | 429 complex_t c2, c3, c4, t1, t2, t3, t4; |
430 | |
431 RE(t2) = RE(cc[ac+i]) + RE(cc[ac+i+2*ido]); | |
432 RE(t1) = RE(cc[ac+i]) - RE(cc[ac+i+2*ido]); | |
433 IM(t2) = IM(cc[ac+i]) + IM(cc[ac+i+2*ido]); | |
434 IM(t1) = IM(cc[ac+i]) - IM(cc[ac+i+2*ido]); | |
435 RE(t3) = RE(cc[ac+i+ido]) + RE(cc[ac+i+3*ido]); | |
436 IM(t4) = RE(cc[ac+i+ido]) - RE(cc[ac+i+3*ido]); | |
437 IM(t3) = IM(cc[ac+i+3*ido]) + IM(cc[ac+i+ido]); | |
438 RE(t4) = IM(cc[ac+i+3*ido]) - IM(cc[ac+i+ido]); | |
439 | |
440 RE(c2) = RE(t1) - RE(t4); | |
441 RE(c4) = RE(t1) + RE(t4); | |
442 | |
443 IM(c2) = IM(t1) - IM(t4); | |
444 IM(c4) = IM(t1) + IM(t4); | |
445 | |
446 RE(ch[ah+i]) = RE(t2) + RE(t3); | |
447 RE(c3) = RE(t2) - RE(t3); | |
448 | |
449 IM(ch[ah+i]) = IM(t2) + IM(t3); | |
450 IM(c3) = IM(t2) - IM(t3); | |
451 | |
452 #if 1 | |
453 ComplexMult(&RE(ch[ah+i+l1*ido]), &IM(ch[ah+i+l1*ido]), | |
454 RE(c2), IM(c2), RE(wa1[i]), IM(wa1[i])); | |
455 ComplexMult(&RE(ch[ah+i+2*l1*ido]), &IM(ch[ah+i+2*l1*ido]), | |
456 RE(c3), IM(c3), RE(wa2[i]), IM(wa2[i])); | |
457 ComplexMult(&RE(ch[ah+i+3*l1*ido]), &IM(ch[ah+i+3*l1*ido]), | |
458 RE(c4), IM(c4), RE(wa3[i]), IM(wa3[i])); | |
459 #else | |
460 ComplexMult(&IM(ch[ah+i+l1*ido]), &RE(ch[ah+i+l1*ido]), | |
461 IM(c2), RE(c2), RE(wa1[i]), IM(wa1[i])); | |
462 ComplexMult(&IM(ch[ah+i+2*l1*ido]), &RE(ch[ah+i+2*l1*ido]), | |
463 IM(c3), RE(c3), RE(wa2[i]), IM(wa2[i])); | |
464 ComplexMult(&IM(ch[ah+i+3*l1*ido]), &RE(ch[ah+i+3*l1*ido]), | |
465 IM(c4), RE(c4), RE(wa3[i]), IM(wa3[i])); | |
466 #endif | |
467 } | |
468 } | |
469 } | |
470 } | |
471 | |
472 static void passf5(const uint16_t ido, const uint16_t l1, const complex_t *cc, | |
473 complex_t *ch, const complex_t *wa1, const complex_t *wa2, const complex_t *wa3, | |
474 const complex_t *wa4, const int8_t isign) | |
475 { | |
476 static real_t tr11 = FRAC_CONST(0.309016994374947); | |
477 static real_t ti11 = FRAC_CONST(0.951056516295154); | |
478 static real_t tr12 = FRAC_CONST(-0.809016994374947); | |
479 static real_t ti12 = FRAC_CONST(0.587785252292473); | |
480 uint16_t i, k, ac, ah; | |
481 complex_t c2, c3, c4, c5, d3, d4, d5, d2, t2, t3, t4, t5; | |
10725 | 482 |
12527 | 483 if (ido == 1) |
484 { | |
485 if (isign == 1) | |
486 { | |
487 for (k = 0; k < l1; k++) | |
488 { | |
489 ac = 5*k + 1; | |
490 ah = k; | |
491 | |
492 RE(t2) = RE(cc[ac]) + RE(cc[ac+3]); | |
493 IM(t2) = IM(cc[ac]) + IM(cc[ac+3]); | |
494 RE(t3) = RE(cc[ac+1]) + RE(cc[ac+2]); | |
495 IM(t3) = IM(cc[ac+1]) + IM(cc[ac+2]); | |
496 RE(t4) = RE(cc[ac+1]) - RE(cc[ac+2]); | |
497 IM(t4) = IM(cc[ac+1]) - IM(cc[ac+2]); | |
498 RE(t5) = RE(cc[ac]) - RE(cc[ac+3]); | |
499 IM(t5) = IM(cc[ac]) - IM(cc[ac+3]); | |
500 | |
501 RE(ch[ah]) = RE(cc[ac-1]) + RE(t2) + RE(t3); | |
502 IM(ch[ah]) = IM(cc[ac-1]) + IM(t2) + IM(t3); | |
503 | |
504 RE(c2) = RE(cc[ac-1]) + MUL_F(RE(t2),tr11) + MUL_F(RE(t3),tr12); | |
505 IM(c2) = IM(cc[ac-1]) + MUL_F(IM(t2),tr11) + MUL_F(IM(t3),tr12); | |
506 RE(c3) = RE(cc[ac-1]) + MUL_F(RE(t2),tr12) + MUL_F(RE(t3),tr11); | |
507 IM(c3) = IM(cc[ac-1]) + MUL_F(IM(t2),tr12) + MUL_F(IM(t3),tr11); | |
508 | |
509 ComplexMult(&RE(c5), &RE(c4), | |
510 ti11, ti12, RE(t5), RE(t4)); | |
511 ComplexMult(&IM(c5), &IM(c4), | |
512 ti11, ti12, IM(t5), IM(t4)); | |
10725 | 513 |
12527 | 514 RE(ch[ah+l1]) = RE(c2) - IM(c5); |
515 IM(ch[ah+l1]) = IM(c2) + RE(c5); | |
516 RE(ch[ah+2*l1]) = RE(c3) - IM(c4); | |
517 IM(ch[ah+2*l1]) = IM(c3) + RE(c4); | |
518 RE(ch[ah+3*l1]) = RE(c3) + IM(c4); | |
519 IM(ch[ah+3*l1]) = IM(c3) - RE(c4); | |
520 RE(ch[ah+4*l1]) = RE(c2) + IM(c5); | |
521 IM(ch[ah+4*l1]) = IM(c2) - RE(c5); | |
522 } | |
523 } else { | |
524 for (k = 0; k < l1; k++) | |
525 { | |
526 ac = 5*k + 1; | |
527 ah = k; | |
528 | |
529 RE(t2) = RE(cc[ac]) + RE(cc[ac+3]); | |
530 IM(t2) = IM(cc[ac]) + IM(cc[ac+3]); | |
531 RE(t3) = RE(cc[ac+1]) + RE(cc[ac+2]); | |
532 IM(t3) = IM(cc[ac+1]) + IM(cc[ac+2]); | |
533 RE(t4) = RE(cc[ac+1]) - RE(cc[ac+2]); | |
534 IM(t4) = IM(cc[ac+1]) - IM(cc[ac+2]); | |
535 RE(t5) = RE(cc[ac]) - RE(cc[ac+3]); | |
536 IM(t5) = IM(cc[ac]) - IM(cc[ac+3]); | |
537 | |
538 RE(ch[ah]) = RE(cc[ac-1]) + RE(t2) + RE(t3); | |
539 IM(ch[ah]) = IM(cc[ac-1]) + IM(t2) + IM(t3); | |
540 | |
541 RE(c2) = RE(cc[ac-1]) + MUL_F(RE(t2),tr11) + MUL_F(RE(t3),tr12); | |
542 IM(c2) = IM(cc[ac-1]) + MUL_F(IM(t2),tr11) + MUL_F(IM(t3),tr12); | |
543 RE(c3) = RE(cc[ac-1]) + MUL_F(RE(t2),tr12) + MUL_F(RE(t3),tr11); | |
544 IM(c3) = IM(cc[ac-1]) + MUL_F(IM(t2),tr12) + MUL_F(IM(t3),tr11); | |
545 | |
546 ComplexMult(&RE(c4), &RE(c5), | |
547 ti12, ti11, RE(t5), RE(t4)); | |
548 ComplexMult(&IM(c4), &IM(c5), | |
549 ti12, ti12, IM(t5), IM(t4)); | |
10725 | 550 |
12527 | 551 RE(ch[ah+l1]) = RE(c2) + IM(c5); |
552 IM(ch[ah+l1]) = IM(c2) - RE(c5); | |
553 RE(ch[ah+2*l1]) = RE(c3) + IM(c4); | |
554 IM(ch[ah+2*l1]) = IM(c3) - RE(c4); | |
555 RE(ch[ah+3*l1]) = RE(c3) - IM(c4); | |
556 IM(ch[ah+3*l1]) = IM(c3) + RE(c4); | |
557 RE(ch[ah+4*l1]) = RE(c2) - IM(c5); | |
558 IM(ch[ah+4*l1]) = IM(c2) + RE(c5); | |
559 } | |
560 } | |
561 } else { | |
562 if (isign == 1) | |
563 { | |
564 for (k = 0; k < l1; k++) | |
565 { | |
566 for (i = 0; i < ido; i++) | |
567 { | |
568 ac = i + (k*5 + 1) * ido; | |
569 ah = i + k * ido; | |
570 | |
571 RE(t2) = RE(cc[ac]) + RE(cc[ac+3*ido]); | |
572 IM(t2) = IM(cc[ac]) + IM(cc[ac+3*ido]); | |
573 RE(t3) = RE(cc[ac+ido]) + RE(cc[ac+2*ido]); | |
574 IM(t3) = IM(cc[ac+ido]) + IM(cc[ac+2*ido]); | |
575 RE(t4) = RE(cc[ac+ido]) - RE(cc[ac+2*ido]); | |
576 IM(t4) = IM(cc[ac+ido]) - IM(cc[ac+2*ido]); | |
577 RE(t5) = RE(cc[ac]) - RE(cc[ac+3*ido]); | |
578 IM(t5) = IM(cc[ac]) - IM(cc[ac+3*ido]); | |
579 | |
580 RE(ch[ah]) = RE(cc[ac-ido]) + RE(t2) + RE(t3); | |
581 IM(ch[ah]) = IM(cc[ac-ido]) + IM(t2) + IM(t3); | |
582 | |
583 RE(c2) = RE(cc[ac-ido]) + MUL_F(RE(t2),tr11) + MUL_F(RE(t3),tr12); | |
584 IM(c2) = IM(cc[ac-ido]) + MUL_F(IM(t2),tr11) + MUL_F(IM(t3),tr12); | |
585 RE(c3) = RE(cc[ac-ido]) + MUL_F(RE(t2),tr12) + MUL_F(RE(t3),tr11); | |
586 IM(c3) = IM(cc[ac-ido]) + MUL_F(IM(t2),tr12) + MUL_F(IM(t3),tr11); | |
587 | |
588 ComplexMult(&RE(c5), &RE(c4), | |
589 ti11, ti12, RE(t5), RE(t4)); | |
590 ComplexMult(&IM(c5), &IM(c4), | |
591 ti11, ti12, IM(t5), IM(t4)); | |
592 | |
593 IM(d2) = IM(c2) + RE(c5); | |
594 IM(d3) = IM(c3) + RE(c4); | |
595 RE(d4) = RE(c3) + IM(c4); | |
596 RE(d5) = RE(c2) + IM(c5); | |
597 RE(d2) = RE(c2) - IM(c5); | |
598 IM(d5) = IM(c2) - RE(c5); | |
599 RE(d3) = RE(c3) - IM(c4); | |
600 IM(d4) = IM(c3) - RE(c4); | |
10725 | 601 |
12527 | 602 #if 1 |
603 ComplexMult(&IM(ch[ah+l1*ido]), &RE(ch[ah+l1*ido]), | |
604 IM(d2), RE(d2), RE(wa1[i]), IM(wa1[i])); | |
605 ComplexMult(&IM(ch[ah+2*l1*ido]), &RE(ch[ah+2*l1*ido]), | |
606 IM(d3), RE(d3), RE(wa2[i]), IM(wa2[i])); | |
607 ComplexMult(&IM(ch[ah+3*l1*ido]), &RE(ch[ah+3*l1*ido]), | |
608 IM(d4), RE(d4), RE(wa3[i]), IM(wa3[i])); | |
609 ComplexMult(&IM(ch[ah+4*l1*ido]), &RE(ch[ah+4*l1*ido]), | |
610 IM(d5), RE(d5), RE(wa4[i]), IM(wa4[i])); | |
611 #else | |
612 ComplexMult(&RE(ch[ah+l1*ido]), &IM(ch[ah+l1*ido]), | |
613 RE(d2), IM(d2), RE(wa1[i]), IM(wa1[i])); | |
614 ComplexMult(&RE(ch[ah+2*l1*ido]), &IM(ch[ah+2*l1*ido]), | |
615 RE(d3), IM(d3), RE(wa2[i]), IM(wa2[i])); | |
616 ComplexMult(&RE(ch[ah+3*l1*ido]), &IM(ch[ah+3*l1*ido]), | |
617 RE(d4), IM(d4), RE(wa3[i]), IM(wa3[i])); | |
618 ComplexMult(&RE(ch[ah+4*l1*ido]), &IM(ch[ah+4*l1*ido]), | |
619 RE(d5), IM(d5), RE(wa4[i]), IM(wa4[i])); | |
620 #endif | |
621 } | |
622 } | |
623 } else { | |
624 for (k = 0; k < l1; k++) | |
625 { | |
626 for (i = 0; i < ido; i++) | |
627 { | |
628 ac = i + (k*5 + 1) * ido; | |
629 ah = i + k * ido; | |
630 | |
631 RE(t2) = RE(cc[ac]) + RE(cc[ac+3*ido]); | |
632 IM(t2) = IM(cc[ac]) + IM(cc[ac+3*ido]); | |
633 RE(t3) = RE(cc[ac+ido]) + RE(cc[ac+2*ido]); | |
634 IM(t3) = IM(cc[ac+ido]) + IM(cc[ac+2*ido]); | |
635 RE(t4) = RE(cc[ac+ido]) - RE(cc[ac+2*ido]); | |
636 IM(t4) = IM(cc[ac+ido]) - IM(cc[ac+2*ido]); | |
637 RE(t5) = RE(cc[ac]) - RE(cc[ac+3*ido]); | |
638 IM(t5) = IM(cc[ac]) - IM(cc[ac+3*ido]); | |
10725 | 639 |
12527 | 640 RE(ch[ah]) = RE(cc[ac-ido]) + RE(t2) + RE(t3); |
641 IM(ch[ah]) = IM(cc[ac-ido]) + IM(t2) + IM(t3); | |
642 | |
643 RE(c2) = RE(cc[ac-ido]) + MUL_F(RE(t2),tr11) + MUL_F(RE(t3),tr12); | |
644 IM(c2) = IM(cc[ac-ido]) + MUL_F(IM(t2),tr11) + MUL_F(IM(t3),tr12); | |
645 RE(c3) = RE(cc[ac-ido]) + MUL_F(RE(t2),tr12) + MUL_F(RE(t3),tr11); | |
646 IM(c3) = IM(cc[ac-ido]) + MUL_F(IM(t2),tr12) + MUL_F(IM(t3),tr11); | |
647 | |
648 ComplexMult(&RE(c4), &RE(c5), | |
649 ti12, ti11, RE(t5), RE(t4)); | |
650 ComplexMult(&IM(c4), &IM(c5), | |
651 ti12, ti12, IM(t5), IM(t4)); | |
652 | |
653 IM(d2) = IM(c2) - RE(c5); | |
654 IM(d3) = IM(c3) - RE(c4); | |
655 RE(d4) = RE(c3) - IM(c4); | |
656 RE(d5) = RE(c2) - IM(c5); | |
657 RE(d2) = RE(c2) + IM(c5); | |
658 IM(d5) = IM(c2) + RE(c5); | |
659 RE(d3) = RE(c3) + IM(c4); | |
660 IM(d4) = IM(c3) + RE(c4); | |
661 | |
662 #if 1 | |
663 ComplexMult(&RE(ch[ah+l1*ido]), &IM(ch[ah+l1*ido]), | |
664 RE(d2), IM(d2), RE(wa1[i]), IM(wa1[i])); | |
665 ComplexMult(&RE(ch[ah+2*l1*ido]), &IM(ch[ah+2*l1*ido]), | |
666 RE(d3), IM(d3), RE(wa2[i]), IM(wa2[i])); | |
667 ComplexMult(&RE(ch[ah+3*l1*ido]), &IM(ch[ah+3*l1*ido]), | |
668 RE(d4), IM(d4), RE(wa3[i]), IM(wa3[i])); | |
669 ComplexMult(&RE(ch[ah+4*l1*ido]), &IM(ch[ah+4*l1*ido]), | |
670 RE(d5), IM(d5), RE(wa4[i]), IM(wa4[i])); | |
671 #else | |
672 ComplexMult(&IM(ch[ah+l1*ido]), &RE(ch[ah+l1*ido]), | |
673 IM(d2), RE(d2), RE(wa1[i]), IM(wa1[i])); | |
674 ComplexMult(&IM(ch[ah+2*l1*ido]), &RE(ch[ah+2*l1*ido]), | |
675 IM(d3), RE(d3), RE(wa2[i]), IM(wa2[i])); | |
676 ComplexMult(&IM(ch[ah+3*l1*ido]), &RE(ch[ah+3*l1*ido]), | |
677 IM(d4), RE(d4), RE(wa3[i]), IM(wa3[i])); | |
678 ComplexMult(&IM(ch[ah+4*l1*ido]), &RE(ch[ah+4*l1*ido]), | |
679 IM(d5), RE(d5), RE(wa4[i]), IM(wa4[i])); | |
680 #endif | |
681 } | |
10725 | 682 } |
683 } | |
684 } | |
685 } | |
686 | |
687 | |
688 /*---------------------------------------------------------------------- | |
689 cfftf1, cfftf, cfftb, cffti1, cffti. Complex FFTs. | |
690 ----------------------------------------------------------------------*/ | |
691 | |
12527 | 692 static INLINE void cfftf1pos(uint16_t n, complex_t *c, complex_t *ch, |
693 const uint16_t *ifac, const complex_t *wa, | |
694 const int8_t isign) | |
10725 | 695 { |
696 uint16_t i; | |
697 uint16_t k1, l1, l2; | |
698 uint16_t na, nf, ip, iw, ix2, ix3, ix4, ido, idl1; | |
699 | |
700 nf = ifac[1]; | |
701 na = 0; | |
702 l1 = 1; | |
703 iw = 0; | |
704 | |
705 for (k1 = 2; k1 <= nf+1; k1++) | |
706 { | |
707 ip = ifac[k1]; | |
708 l2 = ip*l1; | |
709 ido = n / l2; | |
710 idl1 = ido*l1; | |
711 | |
712 switch (ip) | |
713 { | |
10989 | 714 case 4: |
715 ix2 = iw + ido; | |
716 ix3 = ix2 + ido; | |
717 | |
718 if (na == 0) | |
12527 | 719 passf4pos((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw], &wa[ix2], &wa[ix3]); |
10989 | 720 else |
12527 | 721 passf4pos((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw], &wa[ix2], &wa[ix3]); |
10989 | 722 |
723 na = 1 - na; | |
724 break; | |
10725 | 725 case 2: |
726 if (na == 0) | |
12527 | 727 passf2pos((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw]); |
10725 | 728 else |
12527 | 729 passf2pos((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw]); |
10725 | 730 |
731 na = 1 - na; | |
732 break; | |
733 case 3: | |
734 ix2 = iw + ido; | |
735 | |
736 if (na == 0) | |
12527 | 737 passf3((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw], &wa[ix2], isign); |
10725 | 738 else |
12527 | 739 passf3((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw], &wa[ix2], isign); |
10725 | 740 |
741 na = 1 - na; | |
742 break; | |
743 case 5: | |
744 ix2 = iw + ido; | |
745 ix3 = ix2 + ido; | |
746 ix4 = ix3 + ido; | |
747 | |
748 if (na == 0) | |
12527 | 749 passf5((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4], isign); |
10725 | 750 else |
12527 | 751 passf5((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4], isign); |
752 | |
753 na = 1 - na; | |
754 break; | |
755 } | |
756 | |
757 l1 = l2; | |
758 iw += (ip-1) * ido; | |
759 } | |
760 | |
761 if (na == 0) | |
762 return; | |
763 | |
764 for (i = 0; i < n; i++) | |
765 { | |
766 RE(c[i]) = RE(ch[i]); | |
767 IM(c[i]) = IM(ch[i]); | |
768 } | |
769 } | |
770 | |
771 static INLINE void cfftf1neg(uint16_t n, complex_t *c, complex_t *ch, | |
772 const uint16_t *ifac, const complex_t *wa, | |
773 const int8_t isign) | |
774 { | |
775 uint16_t i; | |
776 uint16_t k1, l1, l2; | |
777 uint16_t na, nf, ip, iw, ix2, ix3, ix4, ido, idl1; | |
778 | |
779 nf = ifac[1]; | |
780 na = 0; | |
781 l1 = 1; | |
782 iw = 0; | |
783 | |
784 for (k1 = 2; k1 <= nf+1; k1++) | |
785 { | |
786 ip = ifac[k1]; | |
787 l2 = ip*l1; | |
788 ido = n / l2; | |
789 idl1 = ido*l1; | |
790 | |
791 switch (ip) | |
792 { | |
793 case 4: | |
794 ix2 = iw + ido; | |
795 ix3 = ix2 + ido; | |
796 | |
797 if (na == 0) | |
798 passf4neg((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw], &wa[ix2], &wa[ix3]); | |
799 else | |
800 passf4neg((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw], &wa[ix2], &wa[ix3]); | |
801 | |
802 na = 1 - na; | |
803 break; | |
804 case 2: | |
805 if (na == 0) | |
806 passf2neg((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw]); | |
807 else | |
808 passf2neg((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw]); | |
809 | |
810 na = 1 - na; | |
811 break; | |
812 case 3: | |
813 ix2 = iw + ido; | |
814 | |
815 if (na == 0) | |
816 passf3((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw], &wa[ix2], isign); | |
817 else | |
818 passf3((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw], &wa[ix2], isign); | |
819 | |
820 na = 1 - na; | |
821 break; | |
822 case 5: | |
823 ix2 = iw + ido; | |
824 ix3 = ix2 + ido; | |
825 ix4 = ix3 + ido; | |
826 | |
827 if (na == 0) | |
828 passf5((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4], isign); | |
829 else | |
830 passf5((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4], isign); | |
10725 | 831 |
832 na = 1 - na; | |
833 break; | |
834 } | |
835 | |
836 l1 = l2; | |
837 iw += (ip-1) * ido; | |
838 } | |
839 | |
840 if (na == 0) | |
841 return; | |
842 | |
843 for (i = 0; i < n; i++) | |
844 { | |
845 RE(c[i]) = RE(ch[i]); | |
846 IM(c[i]) = IM(ch[i]); | |
847 } | |
848 } | |
849 | |
850 void cfftf(cfft_info *cfft, complex_t *c) | |
851 { | |
12527 | 852 cfftf1neg(cfft->n, c, cfft->work, (const uint16_t*)cfft->ifac, (const complex_t*)cfft->tab, -1); |
10725 | 853 } |
854 | |
855 void cfftb(cfft_info *cfft, complex_t *c) | |
856 { | |
12527 | 857 cfftf1pos(cfft->n, c, cfft->work, (const uint16_t*)cfft->ifac, (const complex_t*)cfft->tab, +1); |
10725 | 858 } |
859 | |
860 static void cffti1(uint16_t n, complex_t *wa, uint16_t *ifac) | |
861 { | |
862 static uint16_t ntryh[4] = {3, 4, 2, 5}; | |
863 #ifndef FIXED_POINT | |
864 real_t arg, argh, argld, fi; | |
865 uint16_t ido, ipm; | |
866 uint16_t i1, k1, l1, l2; | |
867 uint16_t ld, ii, ip; | |
868 #endif | |
12527 | 869 uint16_t ntry = 0, i, j; |
10725 | 870 uint16_t ib; |
871 uint16_t nf, nl, nq, nr; | |
872 | |
873 nl = n; | |
874 nf = 0; | |
875 j = 0; | |
876 | |
877 startloop: | |
878 j++; | |
879 | |
880 if (j <= 4) | |
881 ntry = ntryh[j-1]; | |
882 else | |
883 ntry += 2; | |
884 | |
885 do | |
886 { | |
887 nq = nl / ntry; | |
888 nr = nl - ntry*nq; | |
889 | |
890 if (nr != 0) | |
891 goto startloop; | |
892 | |
893 nf++; | |
894 ifac[nf+1] = ntry; | |
895 nl = nq; | |
896 | |
897 if (ntry == 2 && nf != 1) | |
898 { | |
899 for (i = 2; i <= nf; i++) | |
900 { | |
901 ib = nf - i + 2; | |
902 ifac[ib+1] = ifac[ib]; | |
903 } | |
904 ifac[2] = 2; | |
905 } | |
906 } while (nl != 1); | |
907 | |
908 ifac[0] = n; | |
909 ifac[1] = nf; | |
910 | |
911 #ifndef FIXED_POINT | |
12527 | 912 argh = (real_t)2.0*(real_t)M_PI / (real_t)n; |
10725 | 913 i = 0; |
914 l1 = 1; | |
915 | |
916 for (k1 = 1; k1 <= nf; k1++) | |
917 { | |
918 ip = ifac[k1+1]; | |
919 ld = 0; | |
920 l2 = l1*ip; | |
921 ido = n / l2; | |
922 ipm = ip - 1; | |
923 | |
924 for (j = 0; j < ipm; j++) | |
925 { | |
926 i1 = i; | |
927 RE(wa[i]) = 1.0; | |
928 IM(wa[i]) = 0.0; | |
929 ld += l1; | |
930 fi = 0; | |
931 argld = ld*argh; | |
932 | |
933 for (ii = 0; ii < ido; ii++) | |
934 { | |
935 i++; | |
936 fi++; | |
937 arg = fi * argld; | |
10989 | 938 RE(wa[i]) = (real_t)cos(arg); |
12527 | 939 #if 1 |
10989 | 940 IM(wa[i]) = (real_t)sin(arg); |
12527 | 941 #else |
942 IM(wa[i]) = (real_t)-sin(arg); | |
943 #endif | |
10725 | 944 } |
945 | |
946 if (ip > 5) | |
947 { | |
948 RE(wa[i1]) = RE(wa[i]); | |
949 IM(wa[i1]) = IM(wa[i]); | |
950 } | |
951 } | |
952 l1 = l2; | |
953 } | |
954 #endif | |
955 } | |
956 | |
957 cfft_info *cffti(uint16_t n) | |
958 { | |
12527 | 959 cfft_info *cfft = (cfft_info*)faad_malloc(sizeof(cfft_info)); |
10725 | 960 |
961 cfft->n = n; | |
12527 | 962 cfft->work = (complex_t*)faad_malloc(n*sizeof(complex_t)); |
10725 | 963 |
964 #ifndef FIXED_POINT | |
12527 | 965 cfft->tab = (complex_t*)faad_malloc(n*sizeof(complex_t)); |
10725 | 966 |
967 cffti1(n, cfft->tab, cfft->ifac); | |
968 #else | |
969 cffti1(n, NULL, cfft->ifac); | |
970 | |
971 switch (n) | |
972 { | |
13453
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
973 case 64: cfft->tab = (complex_t*)cfft_tab_64; break; |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
974 case 512: cfft->tab = (complex_t*)cfft_tab_512; break; |
10725 | 975 #ifdef LD_DEC |
13453
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
976 case 256: cfft->tab = (complex_t*)cfft_tab_256; break; |
12527 | 977 #endif |
978 | |
979 #ifdef ALLOW_SMALL_FRAMELENGTH | |
13453
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
980 case 60: cfft->tab = (complex_t*)cfft_tab_60; break; |
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
981 case 480: cfft->tab = (complex_t*)cfft_tab_480; break; |
12527 | 982 #ifdef LD_DEC |
13453
6d50ef45a058
Update FAAD to a 2.1 beta CVS snapshot from 2004.07.12.
diego
parents:
12625
diff
changeset
|
983 case 240: cfft->tab = (complex_t*)cfft_tab_240; break; |
12527 | 984 #endif |
10725 | 985 #endif |
18141 | 986 case 128: cfft->tab = (complex_t*)cfft_tab_128; break; |
10725 | 987 } |
988 #endif | |
989 | |
990 return cfft; | |
991 } | |
992 | |
993 void cfftu(cfft_info *cfft) | |
994 { | |
12527 | 995 if (cfft->work) faad_free(cfft->work); |
10725 | 996 #ifndef FIXED_POINT |
12527 | 997 if (cfft->tab) faad_free(cfft->tab); |
10725 | 998 #endif |
999 | |
12527 | 1000 if (cfft) faad_free(cfft); |
10989 | 1001 } |
12527 | 1002 |