Mercurial > mplayer.hg
annotate libfaad2/cfft.c @ 26843:7ba8584ff843
Remove incorrectly added formats
author | reimar |
---|---|
date | Sun, 25 May 2008 16:10:55 +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 |