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