Mercurial > audlegacy

comparison Plugins/Input/aac/libfaad2/cfft.c @ 61:fa848bd484d8 trunk

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[svn] Move plugins to Plugins/
author nenolod
date Fri, 28 Oct 2005 22:58:11 -0700
parents
children 0a2ad94e8607
comparison
equal deleted inserted replaced
60:1771f253e1b2 61:fa848bd484d8
1 /*
2 ** FAAD2 - Freeware Advanced Audio (AAC) Decoder including SBR decoding
3 ** Copyright (C) 2003 M. Bakker, Ahead Software AG, http://www.nero.com
4 **
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.
9 **
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.
14 **
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
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 **
25 ** $Id: cfft.c,v 1.19 2003/11/12 20:47:57 menno Exp $
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
46 /*----------------------------------------------------------------------
47 passf2, passf3, passf4, passf5. Complex FFT passes fwd and bwd.
48 ----------------------------------------------------------------------*/
49
50 static void passf2(const uint16_t ido, const uint16_t l1, const complex_t *cc,
51 complex_t *ch, const complex_t *wa, const int8_t isign)
52 {
53 uint16_t i, k, ah, ac;
54
55 if (ido == 1)
56 {
57 for (k = 0; k < l1; k++)
58 {
59 ah = 2*k;
60 ac = 4*k;
61
62 RE(ch[ah]) = RE(cc[ac]) + RE(cc[ac+1]);
63 RE(ch[ah+l1]) = RE(cc[ac]) - RE(cc[ac+1]);
64 IM(ch[ah]) = IM(cc[ac]) + IM(cc[ac+1]);
65 IM(ch[ah+l1]) = IM(cc[ac]) - IM(cc[ac+1]);
66 }
67 } else {
68 if (isign == 1)
69 {
70 for (k = 0; k < l1; k++)
71 {
72 ah = k*ido;
73 ac = 2*k*ido;
74
75 for (i = 0; i < ido; i++)
76 {
77 complex_t t2;
78
79 RE(ch[ah+i]) = RE(cc[ac+i]) + RE(cc[ac+i+ido]);
80 RE(t2) = RE(cc[ac+i]) - RE(cc[ac+i+ido]);
81
82 IM(ch[ah+i]) = IM(cc[ac+i]) + IM(cc[ac+i+ido]);
83 IM(t2) = IM(cc[ac+i]) - IM(cc[ac+i+ido]);
84
85 ComplexMult(&IM(ch[ah+i+l1*ido]), &RE(ch[ah+i+l1*ido]),
86 IM(t2), RE(t2), RE(wa[i]), IM(wa[i]));
87 }
88 }
89 } else {
90 for (k = 0; k < l1; k++)
91 {
92 ah = k*ido;
93 ac = 2*k*ido;
94
95 for (i = 0; i < ido; i++)
96 {
97 complex_t t2;
98
99 RE(ch[ah+i]) = RE(cc[ac+i]) + RE(cc[ac+i+ido]);
100 RE(t2) = RE(cc[ac+i]) - RE(cc[ac+i+ido]);
101
102 IM(ch[ah+i]) = IM(cc[ac+i]) + IM(cc[ac+i+ido]);
103 IM(t2) = IM(cc[ac+i]) - IM(cc[ac+i+ido]);
104
105 ComplexMult(&RE(ch[ah+i+l1*ido]), &IM(ch[ah+i+l1*ido]),
106 RE(t2), IM(t2), RE(wa[i]), IM(wa[i]));
107 }
108 }
109 }
110 }
111 }
112
113
114 static void passf3(const uint16_t ido, const uint16_t l1, const complex_t *cc,
115 complex_t *ch, const complex_t *wa1, const complex_t *wa2,
116 const int8_t isign)
117 {
118 static real_t taur = FRAC_CONST(-0.5);
119 static real_t taui = FRAC_CONST(0.866025403784439);
120 uint16_t i, k, ac, ah;
121 complex_t c2, c3, d2, d3, t2;
122
123 if (ido == 1)
124 {
125 if (isign == 1)
126 {
127 for (k = 0; k < l1; k++)
128 {
129 ac = 3*k+1;
130 ah = k;
131
132 RE(t2) = RE(cc[ac]) + RE(cc[ac+1]);
133 IM(t2) = IM(cc[ac]) + IM(cc[ac+1]);
134 RE(c2) = RE(cc[ac-1]) + MUL_F(RE(t2),taur);
135 IM(c2) = IM(cc[ac-1]) + MUL_F(IM(t2),taur);
136
137 RE(ch[ah]) = RE(cc[ac-1]) + RE(t2);
138 IM(ch[ah]) = IM(cc[ac-1]) + IM(t2);
139
140 RE(c3) = MUL_F((RE(cc[ac]) - RE(cc[ac+1])), taui);
141 IM(c3) = MUL_F((IM(cc[ac]) - IM(cc[ac+1])), taui);
142
143 RE(ch[ah+l1]) = RE(c2) - IM(c3);
144 IM(ch[ah+l1]) = IM(c2) + RE(c3);
145 RE(ch[ah+2*l1]) = RE(c2) + IM(c3);
146 IM(ch[ah+2*l1]) = IM(c2) - RE(c3);
147 }
148 } else {
149 for (k = 0; k < l1; k++)
150 {
151 ac = 3*k+1;
152 ah = k;
153
154 RE(t2) = RE(cc[ac]) + RE(cc[ac+1]);
155 IM(t2) = IM(cc[ac]) + IM(cc[ac+1]);
156 RE(c2) = RE(cc[ac-1]) + MUL_F(RE(t2),taur);
157 IM(c2) = IM(cc[ac-1]) + MUL_F(IM(t2),taur);
158
159 RE(ch[ah]) = RE(cc[ac-1]) + RE(t2);
160 IM(ch[ah]) = IM(cc[ac-1]) + IM(t2);
161
162 RE(c3) = MUL_F((RE(cc[ac]) - RE(cc[ac+1])), taui);
163 IM(c3) = MUL_F((IM(cc[ac]) - IM(cc[ac+1])), taui);
164
165 RE(ch[ah+l1]) = RE(c2) + IM(c3);
166 IM(ch[ah+l1]) = IM(c2) - RE(c3);
167 RE(ch[ah+2*l1]) = RE(c2) - IM(c3);
168 IM(ch[ah+2*l1]) = IM(c2) + RE(c3);
169 }
170 }
171 } else {
172 if (isign == 1)
173 {
174 for (k = 0; k < l1; k++)
175 {
176 for (i = 0; i < ido; i++)
177 {
178 ac = i + (3*k+1)*ido;
179 ah = i + k * ido;
180
181 RE(t2) = RE(cc[ac]) + RE(cc[ac+ido]);
182 RE(c2) = RE(cc[ac-ido]) + MUL_F(RE(t2),taur);
183 IM(t2) = IM(cc[ac]) + IM(cc[ac+ido]);
184 IM(c2) = IM(cc[ac-ido]) + MUL_F(IM(t2),taur);
185
186 RE(ch[ah]) = RE(cc[ac-ido]) + RE(t2);
187 IM(ch[ah]) = IM(cc[ac-ido]) + IM(t2);
188
189 RE(c3) = MUL_F((RE(cc[ac]) - RE(cc[ac+ido])), taui);
190 IM(c3) = MUL_F((IM(cc[ac]) - IM(cc[ac+ido])), taui);
191
192 RE(d2) = RE(c2) - IM(c3);
193 IM(d3) = IM(c2) - RE(c3);
194 RE(d3) = RE(c2) + IM(c3);
195 IM(d2) = IM(c2) + RE(c3);
196
197 ComplexMult(&IM(ch[ah+l1*ido]), &RE(ch[ah+l1*ido]),
198 IM(d2), RE(d2), RE(wa1[i]), IM(wa1[i]));
199 ComplexMult(&IM(ch[ah+2*l1*ido]), &RE(ch[ah+2*l1*ido]),
200 IM(d3), RE(d3), RE(wa2[i]), IM(wa2[i]));
201 }
202 }
203 } else {
204 for (k = 0; k < l1; k++)
205 {
206 for (i = 0; i < ido; i++)
207 {
208 ac = i + (3*k+1)*ido;
209 ah = i + k * ido;
210
211 RE(t2) = RE(cc[ac]) + RE(cc[ac+ido]);
212 RE(c2) = RE(cc[ac-ido]) + MUL_F(RE(t2),taur);
213 IM(t2) = IM(cc[ac]) + IM(cc[ac+ido]);
214 IM(c2) = IM(cc[ac-ido]) + MUL_F(IM(t2),taur);
215
216 RE(ch[ah]) = RE(cc[ac-ido]) + RE(t2);
217 IM(ch[ah]) = IM(cc[ac-ido]) + IM(t2);
218
219 RE(c3) = MUL_F((RE(cc[ac]) - RE(cc[ac+ido])), taui);
220 IM(c3) = MUL_F((IM(cc[ac]) - IM(cc[ac+ido])), taui);
221
222 RE(d2) = RE(c2) + IM(c3);
223 IM(d3) = IM(c2) + RE(c3);
224 RE(d3) = RE(c2) - IM(c3);
225 IM(d2) = IM(c2) - RE(c3);
226
227 ComplexMult(&RE(ch[ah+l1*ido]), &IM(ch[ah+l1*ido]),
228 RE(d2), IM(d2), RE(wa1[i]), IM(wa1[i]));
229 ComplexMult(&RE(ch[ah+2*l1*ido]), &IM(ch[ah+2*l1*ido]),
230 RE(d3), IM(d3), RE(wa2[i]), IM(wa2[i]));
231 }
232 }
233 }
234 }
235 }
236
237 static void passf4(const uint16_t ido, const uint16_t l1, const complex_t *cc,
238 complex_t *ch, const complex_t *wa1, const complex_t *wa2,
239 const complex_t *wa3, const int8_t isign)
240 {
241 uint16_t i, k, ac, ah;
242
243 if (ido == 1)
244 {
245 if (isign == 1)
246 {
247 for (k = 0; k < l1; k++)
248 {
249 complex_t t1, t2, t3, t4;
250
251 ac = 4*k;
252 ah = k;
253
254 RE(t2) = RE(cc[ac]) + RE(cc[ac+2]);
255 RE(t1) = RE(cc[ac]) - RE(cc[ac+2]);
256 IM(t2) = IM(cc[ac]) + IM(cc[ac+2]);
257 IM(t1) = IM(cc[ac]) - IM(cc[ac+2]);
258 RE(t3) = RE(cc[ac+1]) + RE(cc[ac+3]);
259 IM(t4) = RE(cc[ac+1]) - RE(cc[ac+3]);
260 IM(t3) = IM(cc[ac+3]) + IM(cc[ac+1]);
261 RE(t4) = IM(cc[ac+3]) - IM(cc[ac+1]);
262
263 RE(ch[ah]) = RE(t2) + RE(t3);
264 RE(ch[ah+2*l1]) = RE(t2) - RE(t3);
265
266 IM(ch[ah]) = IM(t2) + IM(t3);
267 IM(ch[ah+2*l1]) = IM(t2) - IM(t3);
268
269 RE(ch[ah+l1]) = RE(t1) + RE(t4);
270 RE(ch[ah+3*l1]) = RE(t1) - RE(t4);
271
272 IM(ch[ah+l1]) = IM(t1) + IM(t4);
273 IM(ch[ah+3*l1]) = IM(t1) - IM(t4);
274 }
275 } else {
276 for (k = 0; k < l1; k++)
277 {
278 complex_t t1, t2, t3, t4;
279
280 ac = 4*k;
281 ah = k;
282
283 RE(t2) = RE(cc[ac]) + RE(cc[ac+2]);
284 RE(t1) = RE(cc[ac]) - RE(cc[ac+2]);
285 IM(t2) = IM(cc[ac]) + IM(cc[ac+2]);
286 IM(t1) = IM(cc[ac]) - IM(cc[ac+2]);
287 RE(t3) = RE(cc[ac+1]) + RE(cc[ac+3]);
288 IM(t4) = RE(cc[ac+1]) - RE(cc[ac+3]);
289 IM(t3) = IM(cc[ac+3]) + IM(cc[ac+1]);
290 RE(t4) = IM(cc[ac+3]) - IM(cc[ac+1]);
291
292 RE(ch[ah]) = RE(t2) + RE(t3);
293 RE(ch[ah+2*l1]) = RE(t2) - RE(t3);
294
295 IM(ch[ah]) = IM(t2) + IM(t3);
296 IM(ch[ah+2*l1]) = IM(t2) - IM(t3);
297
298 RE(ch[ah+l1]) = RE(t1) - RE(t4);
299 RE(ch[ah+3*l1]) = RE(t1) + RE(t4);
300
301 IM(ch[ah+l1]) = IM(t1) - IM(t4);
302 IM(ch[ah+3*l1]) = IM(t1) + IM(t4);
303 }
304 }
305 } else {
306 if (isign == 1)
307 {
308 for (k = 0; k < l1; k++)
309 {
310 ac = 4*k*ido;
311 ah = k*ido;
312
313 for (i = 0; i < ido; i++)
314 {
315 complex_t c2, c3, c4, t1, t2, t3, t4;
316
317 RE(t2) = RE(cc[ac+i]) + RE(cc[ac+i+2*ido]);
318 RE(t1) = RE(cc[ac+i]) - RE(cc[ac+i+2*ido]);
319 IM(t2) = IM(cc[ac+i]) + IM(cc[ac+i+2*ido]);
320 IM(t1) = IM(cc[ac+i]) - IM(cc[ac+i+2*ido]);
321 RE(t3) = RE(cc[ac+i+ido]) + RE(cc[ac+i+3*ido]);
322 IM(t4) = RE(cc[ac+i+ido]) - RE(cc[ac+i+3*ido]);
323 IM(t3) = IM(cc[ac+i+3*ido]) + IM(cc[ac+i+ido]);
324 RE(t4) = IM(cc[ac+i+3*ido]) - IM(cc[ac+i+ido]);
325
326 RE(c2) = RE(t1) + RE(t4);
327 RE(c4) = RE(t1) - RE(t4);
328
329 IM(c2) = IM(t1) + IM(t4);
330 IM(c4) = IM(t1) - IM(t4);
331
332 RE(ch[ah+i]) = RE(t2) + RE(t3);
333 RE(c3) = RE(t2) - RE(t3);
334
335 IM(ch[ah+i]) = IM(t2) + IM(t3);
336 IM(c3) = IM(t2) - IM(t3);
337
338 ComplexMult(&IM(ch[ah+i+l1*ido]), &RE(ch[ah+i+l1*ido]),
339 IM(c2), RE(c2), RE(wa1[i]), IM(wa1[i]));
340 ComplexMult(&IM(ch[ah+i+2*l1*ido]), &RE(ch[ah+i+2*l1*ido]),
341 IM(c3), RE(c3), RE(wa2[i]), IM(wa2[i]));
342 ComplexMult(&IM(ch[ah+i+3*l1*ido]), &RE(ch[ah+i+3*l1*ido]),
343 IM(c4), RE(c4), RE(wa3[i]), IM(wa3[i]));
344 }
345 }
346 } else {
347 for (k = 0; k < l1; k++)
348 {
349 ac = 4*k*ido;
350 ah = k*ido;
351
352 for (i = 0; i < ido; i++)
353 {
354 complex_t c2, c3, c4, t1, t2, t3, t4;
355
356 RE(t2) = RE(cc[ac+i]) + RE(cc[ac+i+2*ido]);
357 RE(t1) = RE(cc[ac+i]) - RE(cc[ac+i+2*ido]);
358 IM(t2) = IM(cc[ac+i]) + IM(cc[ac+i+2*ido]);
359 IM(t1) = IM(cc[ac+i]) - IM(cc[ac+i+2*ido]);
360 RE(t3) = RE(cc[ac+i+ido]) + RE(cc[ac+i+3*ido]);
361 IM(t4) = RE(cc[ac+i+ido]) - RE(cc[ac+i+3*ido]);
362 IM(t3) = IM(cc[ac+i+3*ido]) + IM(cc[ac+i+ido]);
363 RE(t4) = IM(cc[ac+i+3*ido]) - IM(cc[ac+i+ido]);
364
365 RE(c2) = RE(t1) - RE(t4);
366 RE(c4) = RE(t1) + RE(t4);
367
368 IM(c2) = IM(t1) - IM(t4);
369 IM(c4) = IM(t1) + IM(t4);
370
371 RE(ch[ah+i]) = RE(t2) + RE(t3);
372 RE(c3) = RE(t2) - RE(t3);
373
374 IM(ch[ah+i]) = IM(t2) + IM(t3);
375 IM(c3) = IM(t2) - IM(t3);
376
377 ComplexMult(&RE(ch[ah+i+l1*ido]), &IM(ch[ah+i+l1*ido]),
378 RE(c2), IM(c2), RE(wa1[i]), IM(wa1[i]));
379 ComplexMult(&RE(ch[ah+i+2*l1*ido]), &IM(ch[ah+i+2*l1*ido]),
380 RE(c3), IM(c3), RE(wa2[i]), IM(wa2[i]));
381 ComplexMult(&RE(ch[ah+i+3*l1*ido]), &IM(ch[ah+i+3*l1*ido]),
382 RE(c4), IM(c4), RE(wa3[i]), IM(wa3[i]));
383 }
384 }
385 }
386 }
387 }
388
389 static void passf5(const uint16_t ido, const uint16_t l1, const complex_t *cc,
390 complex_t *ch, const complex_t *wa1, const complex_t *wa2, const complex_t *wa3,
391 const complex_t *wa4, const int8_t isign)
392 {
393 static real_t tr11 = FRAC_CONST(0.309016994374947);
394 static real_t ti11 = FRAC_CONST(0.951056516295154);
395 static real_t tr12 = FRAC_CONST(-0.809016994374947);
396 static real_t ti12 = FRAC_CONST(0.587785252292473);
397 uint16_t i, k, ac, ah;
398 complex_t c2, c3, c4, c5, d3, d4, d5, d2, t2, t3, t4, t5;
399
400 if (ido == 1)
401 {
402 if (isign == 1)
403 {
404 for (k = 0; k < l1; k++)
405 {
406 ac = 5*k + 1;
407 ah = k;
408
409 RE(t2) = RE(cc[ac]) + RE(cc[ac+3]);
410 IM(t2) = IM(cc[ac]) + IM(cc[ac+3]);
411 RE(t3) = RE(cc[ac+1]) + RE(cc[ac+2]);
412 IM(t3) = IM(cc[ac+1]) + IM(cc[ac+2]);
413 RE(t4) = RE(cc[ac+1]) - RE(cc[ac+2]);
414 IM(t4) = IM(cc[ac+1]) - IM(cc[ac+2]);
415 RE(t5) = RE(cc[ac]) - RE(cc[ac+3]);
416 IM(t5) = IM(cc[ac]) - IM(cc[ac+3]);
417
418 RE(ch[ah]) = RE(cc[ac-1]) + RE(t2) + RE(t3);
419 IM(ch[ah]) = IM(cc[ac-1]) + IM(t2) + IM(t3);
420
421 RE(c2) = RE(cc[ac-1]) + MUL_F(RE(t2),tr11) + MUL_F(RE(t3),tr12);
422 IM(c2) = IM(cc[ac-1]) + MUL_F(IM(t2),tr11) + MUL_F(IM(t3),tr12);
423 RE(c3) = RE(cc[ac-1]) + MUL_F(RE(t2),tr12) + MUL_F(RE(t3),tr11);
424 IM(c3) = IM(cc[ac-1]) + MUL_F(IM(t2),tr12) + MUL_F(IM(t3),tr11);
425
426 ComplexMult(&RE(c5), &RE(c4),
427 ti11, ti12, RE(t5), RE(t4));
428 ComplexMult(&IM(c5), &IM(c4),
429 ti11, ti12, IM(t5), IM(t4));
430
431 RE(ch[ah+l1]) = RE(c2) - IM(c5);
432 IM(ch[ah+l1]) = IM(c2) + RE(c5);
433 RE(ch[ah+2*l1]) = RE(c3) - IM(c4);
434 IM(ch[ah+2*l1]) = IM(c3) + RE(c4);
435 RE(ch[ah+3*l1]) = RE(c3) + IM(c4);
436 IM(ch[ah+3*l1]) = IM(c3) - RE(c4);
437 RE(ch[ah+4*l1]) = RE(c2) + IM(c5);
438 IM(ch[ah+4*l1]) = IM(c2) - RE(c5);
439 }
440 } else {
441 for (k = 0; k < l1; k++)
442 {
443 ac = 5*k + 1;
444 ah = k;
445
446 RE(t2) = RE(cc[ac]) + RE(cc[ac+3]);
447 IM(t2) = IM(cc[ac]) + IM(cc[ac+3]);
448 RE(t3) = RE(cc[ac+1]) + RE(cc[ac+2]);
449 IM(t3) = IM(cc[ac+1]) + IM(cc[ac+2]);
450 RE(t4) = RE(cc[ac+1]) - RE(cc[ac+2]);
451 IM(t4) = IM(cc[ac+1]) - IM(cc[ac+2]);
452 RE(t5) = RE(cc[ac]) - RE(cc[ac+3]);
453 IM(t5) = IM(cc[ac]) - IM(cc[ac+3]);
454
455 RE(ch[ah]) = RE(cc[ac-1]) + RE(t2) + RE(t3);
456 IM(ch[ah]) = IM(cc[ac-1]) + IM(t2) + IM(t3);
457
458 RE(c2) = RE(cc[ac-1]) + MUL_F(RE(t2),tr11) + MUL_F(RE(t3),tr12);
459 IM(c2) = IM(cc[ac-1]) + MUL_F(IM(t2),tr11) + MUL_F(IM(t3),tr12);
460 RE(c3) = RE(cc[ac-1]) + MUL_F(RE(t2),tr12) + MUL_F(RE(t3),tr11);
461 IM(c3) = IM(cc[ac-1]) + MUL_F(IM(t2),tr12) + MUL_F(IM(t3),tr11);
462
463 ComplexMult(&RE(c4), &RE(c5),
464 ti12, ti11, RE(t5), RE(t4));
465 ComplexMult(&IM(c4), &IM(c5),
466 ti12, ti12, IM(t5), IM(t4));
467
468 RE(ch[ah+l1]) = RE(c2) + IM(c5);
469 IM(ch[ah+l1]) = IM(c2) - RE(c5);
470 RE(ch[ah+2*l1]) = RE(c3) + IM(c4);
471 IM(ch[ah+2*l1]) = IM(c3) - RE(c4);
472 RE(ch[ah+3*l1]) = RE(c3) - IM(c4);
473 IM(ch[ah+3*l1]) = IM(c3) + RE(c4);
474 RE(ch[ah+4*l1]) = RE(c2) - IM(c5);
475 IM(ch[ah+4*l1]) = IM(c2) + RE(c5);
476 }
477 }
478 } else {
479 if (isign == 1)
480 {
481 for (k = 0; k < l1; k++)
482 {
483 for (i = 0; i < ido; i++)
484 {
485 ac = i + (k*5 + 1) * ido;
486 ah = i + k * ido;
487
488 RE(t2) = RE(cc[ac]) + RE(cc[ac+3*ido]);
489 IM(t2) = IM(cc[ac]) + IM(cc[ac+3*ido]);
490 RE(t3) = RE(cc[ac+ido]) + RE(cc[ac+2*ido]);
491 IM(t3) = IM(cc[ac+ido]) + IM(cc[ac+2*ido]);
492 RE(t4) = RE(cc[ac+ido]) - RE(cc[ac+2*ido]);
493 IM(t4) = IM(cc[ac+ido]) - IM(cc[ac+2*ido]);
494 RE(t5) = RE(cc[ac]) - RE(cc[ac+3*ido]);
495 IM(t5) = IM(cc[ac]) - IM(cc[ac+3*ido]);
496
497 RE(ch[ah]) = RE(cc[ac-ido]) + RE(t2) + RE(t3);
498 IM(ch[ah]) = IM(cc[ac-ido]) + IM(t2) + IM(t3);
499
500 RE(c2) = RE(cc[ac-ido]) + MUL_F(RE(t2),tr11) + MUL_F(RE(t3),tr12);
501 IM(c2) = IM(cc[ac-ido]) + MUL_F(IM(t2),tr11) + MUL_F(IM(t3),tr12);
502 RE(c3) = RE(cc[ac-ido]) + MUL_F(RE(t2),tr12) + MUL_F(RE(t3),tr11);
503 IM(c3) = IM(cc[ac-ido]) + MUL_F(IM(t2),tr12) + MUL_F(IM(t3),tr11);
504
505 ComplexMult(&RE(c5), &RE(c4),
506 ti11, ti12, RE(t5), RE(t4));
507 ComplexMult(&IM(c5), &IM(c4),
508 ti11, ti12, IM(t5), IM(t4));
509
510 IM(d2) = IM(c2) + RE(c5);
511 IM(d3) = IM(c3) + RE(c4);
512 RE(d4) = RE(c3) + IM(c4);
513 RE(d5) = RE(c2) + IM(c5);
514 RE(d2) = RE(c2) - IM(c5);
515 IM(d5) = IM(c2) - RE(c5);
516 RE(d3) = RE(c3) - IM(c4);
517 IM(d4) = IM(c3) - RE(c4);
518
519 ComplexMult(&IM(ch[ah+l1*ido]), &RE(ch[ah+l1*ido]),
520 IM(d2), RE(d2), RE(wa1[i]), IM(wa1[i]));
521 ComplexMult(&IM(ch[ah+2*l1*ido]), &RE(ch[ah+2*l1*ido]),
522 IM(d3), RE(d3), RE(wa2[i]), IM(wa2[i]));
523 ComplexMult(&IM(ch[ah+3*l1*ido]), &RE(ch[ah+3*l1*ido]),
524 IM(d4), RE(d4), RE(wa3[i]), IM(wa3[i]));
525 ComplexMult(&IM(ch[ah+4*l1*ido]), &RE(ch[ah+4*l1*ido]),
526 IM(d5), RE(d5), RE(wa4[i]), IM(wa4[i]));
527 }
528 }
529 } else {
530 for (k = 0; k < l1; k++)
531 {
532 for (i = 0; i < ido; i++)
533 {
534 ac = i + (k*5 + 1) * ido;
535 ah = i + k * ido;
536
537 RE(t2) = RE(cc[ac]) + RE(cc[ac+3*ido]);
538 IM(t2) = IM(cc[ac]) + IM(cc[ac+3*ido]);
539 RE(t3) = RE(cc[ac+ido]) + RE(cc[ac+2*ido]);
540 IM(t3) = IM(cc[ac+ido]) + IM(cc[ac+2*ido]);
541 RE(t4) = RE(cc[ac+ido]) - RE(cc[ac+2*ido]);
542 IM(t4) = IM(cc[ac+ido]) - IM(cc[ac+2*ido]);
543 RE(t5) = RE(cc[ac]) - RE(cc[ac+3*ido]);
544 IM(t5) = IM(cc[ac]) - IM(cc[ac+3*ido]);
545
546 RE(ch[ah]) = RE(cc[ac-ido]) + RE(t2) + RE(t3);
547 IM(ch[ah]) = IM(cc[ac-ido]) + IM(t2) + IM(t3);
548
549 RE(c2) = RE(cc[ac-ido]) + MUL_F(RE(t2),tr11) + MUL_F(RE(t3),tr12);
550 IM(c2) = IM(cc[ac-ido]) + MUL_F(IM(t2),tr11) + MUL_F(IM(t3),tr12);
551 RE(c3) = RE(cc[ac-ido]) + MUL_F(RE(t2),tr12) + MUL_F(RE(t3),tr11);
552 IM(c3) = IM(cc[ac-ido]) + MUL_F(IM(t2),tr12) + MUL_F(IM(t3),tr11);
553
554 ComplexMult(&RE(c4), &RE(c5),
555 ti12, ti11, RE(t5), RE(t4));
556 ComplexMult(&IM(c4), &IM(c5),
557 ti12, ti12, IM(t5), IM(t4));
558
559 IM(d2) = IM(c2) - RE(c5);
560 IM(d3) = IM(c3) - RE(c4);
561 RE(d4) = RE(c3) - IM(c4);
562 RE(d5) = RE(c2) - IM(c5);
563 RE(d2) = RE(c2) + IM(c5);
564 IM(d5) = IM(c2) + RE(c5);
565 RE(d3) = RE(c3) + IM(c4);
566 IM(d4) = IM(c3) + RE(c4);
567
568 ComplexMult(&RE(ch[ah+l1*ido]), &IM(ch[ah+l1*ido]),
569 RE(d2), IM(d2), RE(wa1[i]), IM(wa1[i]));
570 ComplexMult(&RE(ch[ah+2*l1*ido]), &IM(ch[ah+2*l1*ido]),
571 RE(d3), IM(d3), RE(wa2[i]), IM(wa2[i]));
572 ComplexMult(&RE(ch[ah+3*l1*ido]), &IM(ch[ah+3*l1*ido]),
573 RE(d4), IM(d4), RE(wa3[i]), IM(wa3[i]));
574 ComplexMult(&RE(ch[ah+4*l1*ido]), &IM(ch[ah+4*l1*ido]),
575 RE(d5), IM(d5), RE(wa4[i]), IM(wa4[i]));
576 }
577 }
578 }
579 }
580 }
581
582
583 /*----------------------------------------------------------------------
584 cfftf1, cfftf, cfftb, cffti1, cffti. Complex FFTs.
585 ----------------------------------------------------------------------*/
586
587 INLINE void cfftf1(uint16_t n, complex_t *c, complex_t *ch,
588 uint16_t *ifac, complex_t *wa, int8_t isign)
589 {
590 uint16_t i;
591 uint16_t k1, l1, l2;
592 uint16_t na, nf, ip, iw, ix2, ix3, ix4, ido, idl1;
593
594 nf = ifac[1];
595 na = 0;
596 l1 = 1;
597 iw = 0;
598
599 for (k1 = 2; k1 <= nf+1; k1++)
600 {
601 ip = ifac[k1];
602 l2 = ip*l1;
603 ido = n / l2;
604 idl1 = ido*l1;
605
606 switch (ip)
607 {
608 case 4:
609 ix2 = iw + ido;
610 ix3 = ix2 + ido;
611
612 if (na == 0)
613 passf4(ido, l1, c, ch, &wa[iw], &wa[ix2], &wa[ix3], isign);
614 else
615 passf4(ido, l1, ch, c, &wa[iw], &wa[ix2], &wa[ix3], isign);
616
617 na = 1 - na;
618 break;
619 case 2:
620 if (na == 0)
621 passf2(ido, l1, c, ch, &wa[iw], isign);
622 else
623 passf2(ido, l1, ch, c, &wa[iw], isign);
624
625 na = 1 - na;
626 break;
627 case 3:
628 ix2 = iw + ido;
629
630 if (na == 0)
631 passf3(ido, l1, c, ch, &wa[iw], &wa[ix2], isign);
632 else
633 passf3(ido, l1, ch, c, &wa[iw], &wa[ix2], isign);
634
635 na = 1 - na;
636 break;
637 case 5:
638 ix2 = iw + ido;
639 ix3 = ix2 + ido;
640 ix4 = ix3 + ido;
641
642 if (na == 0)
643 passf5(ido, l1, c, ch, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4], isign);
644 else
645 passf5(ido, l1, ch, c, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4], isign);
646
647 na = 1 - na;
648 break;
649 }
650
651 l1 = l2;
652 iw += (ip-1) * ido;
653 }
654
655 if (na == 0)
656 return;
657
658 for (i = 0; i < n; i++)
659 {
660 RE(c[i]) = RE(ch[i]);
661 IM(c[i]) = IM(ch[i]);
662 }
663 }
664
665 void cfftf(cfft_info *cfft, complex_t *c)
666 {
667 cfftf1(cfft->n, c, cfft->work, cfft->ifac, cfft->tab, -1);
668 }
669
670 void cfftb(cfft_info *cfft, complex_t *c)
671 {
672 cfftf1(cfft->n, c, cfft->work, cfft->ifac, cfft->tab, +1);
673 }
674
675 static void cffti1(uint16_t n, complex_t *wa, uint16_t *ifac)
676 {
677 static uint16_t ntryh[4] = {3, 4, 2, 5};
678 #ifndef FIXED_POINT
679 real_t arg, argh, argld, fi;
680 uint16_t ido, ipm;
681 uint16_t i1, k1, l1, l2;
682 uint16_t ld, ii, ip;
683 #endif
684 uint16_t ntry, i, j;
685 uint16_t ib;
686 uint16_t nf, nl, nq, nr;
687
688 nl = n;
689 nf = 0;
690 j = 0;
691
692 startloop:
693 j++;
694
695 if (j <= 4)
696 ntry = ntryh[j-1];
697 else
698 ntry += 2;
699
700 do
701 {
702 nq = nl / ntry;
703 nr = nl - ntry*nq;
704
705 if (nr != 0)
706 goto startloop;
707
708 nf++;
709 ifac[nf+1] = ntry;
710 nl = nq;
711
712 if (ntry == 2 && nf != 1)
713 {
714 for (i = 2; i <= nf; i++)
715 {
716 ib = nf - i + 2;
717 ifac[ib+1] = ifac[ib];
718 }
719 ifac[2] = 2;
720 }
721 } while (nl != 1);
722
723 ifac[0] = n;
724 ifac[1] = nf;
725
726 #ifndef FIXED_POINT
727 argh = (real_t)2.0*M_PI / (real_t)n;
728 i = 0;
729 l1 = 1;
730
731 for (k1 = 1; k1 <= nf; k1++)
732 {
733 ip = ifac[k1+1];
734 ld = 0;
735 l2 = l1*ip;
736 ido = n / l2;
737 ipm = ip - 1;
738
739 for (j = 0; j < ipm; j++)
740 {
741 i1 = i;
742 RE(wa[i]) = 1.0;
743 IM(wa[i]) = 0.0;
744 ld += l1;
745 fi = 0;
746 argld = ld*argh;
747
748 for (ii = 0; ii < ido; ii++)
749 {
750 i++;
751 fi++;
752 arg = fi * argld;
753 RE(wa[i]) = (real_t)cos(arg);
754 IM(wa[i]) = (real_t)sin(arg);
755 }
756
757 if (ip > 5)
758 {
759 RE(wa[i1]) = RE(wa[i]);
760 IM(wa[i1]) = IM(wa[i]);
761 }
762 }
763 l1 = l2;
764 }
765 #endif
766 }
767
768 cfft_info *cffti(uint16_t n)
769 {
770 cfft_info *cfft = (cfft_info*)malloc(sizeof(cfft_info));
771
772 cfft->n = n;
773 cfft->work = (complex_t*)malloc(n*sizeof(complex_t));
774
775 #ifndef FIXED_POINT
776 cfft->tab = (complex_t*)malloc(n*sizeof(complex_t));
777
778 cffti1(n, cfft->tab, cfft->ifac);
779 #else
780 cffti1(n, NULL, cfft->ifac);
781
782 switch (n)
783 {
784 case 64: cfft->tab = cfft_tab_64; break;
785 case 512: cfft->tab = cfft_tab_512; break;
786 #ifdef LD_DEC
787 case 256: cfft->tab = cfft_tab_256; break;
788 #endif
789
790 #ifdef ALLOW_SMALL_FRAMELENGTH
791 case 60: cfft->tab = cfft_tab_60; break;
792 case 480: cfft->tab = cfft_tab_480; break;
793 #ifdef LD_DEC
794 case 240: cfft->tab = cfft_tab_240; break;
795 #endif
796 #endif
797 }
798 #endif
799
800 return cfft;
801 }
802
803 void cfftu(cfft_info *cfft)
804 {
805 if (cfft->work) free(cfft->work);
806 #ifndef FIXED_POINT
807 if (cfft->tab) free(cfft->tab);
808 #endif
809
810 if (cfft) free(cfft);
811 }