1 /*
2 * TwoLAME: an optimized MPEG Audio Layer Two encoder
3 *
4 * Copyright (C) 2001-2004 Michael Cheng
5 * Copyright (C) 2004-2006 The TwoLAME Project
6 *
7 * This library is free software; you can redistribute it and/or
8 * modify it under the terms of the GNU Lesser General Public
9 * License as published by the Free Software Foundation; either
10 * version 2.1 of the License, or (at your option) any later version.
11 *
12 * This library is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 * Lesser General Public License for more details.
16 *
17 * You should have received a copy of the GNU Lesser General Public
18 * License along with this library; if not, write to the Free Software
19 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 *
21 * $Id$
22 *
23 */
24
25
26
27 /*
28 ** FFT and FHT routines
29 ** Copyright 1988, 1993; Ron Mayer
30 **
31 ** fht(fz,n);
32 ** Does a hartley transform of "n" points in the array "fz".
33 **
34 ** NOTE: This routine uses at least 2 patented algorithms, and may be
35 ** under the restrictions of a bunch of different organizations.
36 ** Although I wrote it completely myself; it is kind of a derivative
37 ** of a routine I once authored and released under the GPL, so it
38 ** may fall under the free software foundation's restrictions;
39 ** it was worked on as a Stanford Univ project, so they claim
40 ** some rights to it; it was further optimized at work here, so
41 ** I think this company claims parts of it. The patents are
42 ** held by R. Bracewell (the FHT algorithm) and O. Buneman (the
43 ** trig generator), both at Stanford Univ.
44 ** If it were up to me, I'd say go do whatever you want with it;
45 ** but it would be polite to give credit to the following people
46 ** if you use this anywhere:
47 ** Euler - probable inventor of the fourier transform.
48 ** Gauss - probable inventor of the FFT.
49 ** Hartley - probable inventor of the hartley transform.
50 ** Buneman - for a really cool trig generator
51 ** Mayer(me) - for authoring this particular version and
52 ** including all the optimizations in one package.
53 ** Thanks,
54 ** Ron Mayer; mayer@acuson.com
55 **
56 */
57
58 #include <stdio.h>
59 #include <math.h>
60
61 #include "twolame.h"
62 #include "common.h"
63 #include "fft.h"
64
65
66
67 #define SQRT2 1.4142135623730951454746218587388284504414
68
69
70 static const FLOAT costab[20] = {
71 .00000000000000000000000000000000000000000000000000,
72 .70710678118654752440084436210484903928483593768847,
73 .92387953251128675612818318939678828682241662586364,
74 .98078528040323044912618223613423903697393373089333,
75 .99518472667219688624483695310947992157547486872985,
76 .99879545620517239271477160475910069444320361470461,
77 .99969881869620422011576564966617219685006108125772,
78 .99992470183914454092164649119638322435060646880221,
79 .99998117528260114265699043772856771617391725094433,
80 .99999529380957617151158012570011989955298763362218,
81 .99999882345170190992902571017152601904826792288976,
82 .99999970586288221916022821773876567711626389934930,
83 .99999992646571785114473148070738785694820115568892,
84 .99999998161642929380834691540290971450507605124278,
85 .99999999540410731289097193313960614895889430318945,
86 .99999999885102682756267330779455410840053741619428
87 };
88 static const FLOAT sintab[20] = {
89 1.0000000000000000000000000000000000000000000000000,
90 .70710678118654752440084436210484903928483593768846,
91 .38268343236508977172845998403039886676134456248561,
92 .19509032201612826784828486847702224092769161775195,
93 .09801714032956060199419556388864184586113667316749,
94 .04906767432741801425495497694268265831474536302574,
95 .02454122852291228803173452945928292506546611923944,
96 .01227153828571992607940826195100321214037231959176,
97 .00613588464915447535964023459037258091705788631738,
98 .00306795676296597627014536549091984251894461021344,
99 .00153398018628476561230369715026407907995486457522,
100 .00076699031874270452693856835794857664314091945205,
101 .00038349518757139558907246168118138126339502603495,
102 .00019174759731070330743990956198900093346887403385,
103 .00009587379909597734587051721097647635118706561284,
104 .00004793689960306688454900399049465887274686668768
105 };
106
107 /* This is a simplified version for n an even power of 2 */
108 /* MFC: In the case of LayerII encoding, n==1024 always. */
109
fht(FLOAT * fz)110 static void fht(FLOAT * fz)
111 {
112 int i, k, k1, k2, k3, k4, kx;
113 FLOAT *fi, *fn, *gi;
114 FLOAT t_c, t_s;
115
116 FLOAT a;
117 static const struct {
118 unsigned short k1, k2;
119 } k1k2tab[8 * 62] = {
120 {
121 0x020, 0x010}, {
122 0x040, 0x008}, {
123 0x050, 0x028}, {
124 0x060, 0x018}, {
125 0x068, 0x058}, {
126 0x070, 0x038}, {
127 0x080, 0x004}, {
128 0x088, 0x044}, {
129 0x090, 0x024}, {
130 0x098, 0x064}, {
131 0x0a0, 0x014}, {
132 0x0a4, 0x094}, {
133 0x0a8, 0x054}, {
134 0x0b0, 0x034}, {
135 0x0b8, 0x074}, {
136 0x0c0, 0x00c}, {
137 0x0c4, 0x08c}, {
138 0x0c8, 0x04c}, {
139 0x0d0, 0x02c}, {
140 0x0d4, 0x0ac}, {
141 0x0d8, 0x06c}, {
142 0x0e0, 0x01c}, {
143 0x0e4, 0x09c}, {
144 0x0e8, 0x05c}, {
145 0x0ec, 0x0dc}, {
146 0x0f0, 0x03c}, {
147 0x0f4, 0x0bc}, {
148 0x0f8, 0x07c}, {
149 0x100, 0x002}, {
150 0x104, 0x082}, {
151 0x108, 0x042}, {
152 0x10c, 0x0c2}, {
153 0x110, 0x022}, {
154 0x114, 0x0a2}, {
155 0x118, 0x062}, {
156 0x11c, 0x0e2}, {
157 0x120, 0x012}, {
158 0x122, 0x112}, {
159 0x124, 0x092}, {
160 0x128, 0x052}, {
161 0x12c, 0x0d2}, {
162 0x130, 0x032}, {
163 0x134, 0x0b2}, {
164 0x138, 0x072}, {
165 0x13c, 0x0f2}, {
166 0x140, 0x00a}, {
167 0x142, 0x10a}, {
168 0x144, 0x08a}, {
169 0x148, 0x04a}, {
170 0x14c, 0x0ca}, {
171 0x150, 0x02a}, {
172 0x152, 0x12a}, {
173 0x154, 0x0aa}, {
174 0x158, 0x06a}, {
175 0x15c, 0x0ea}, {
176 0x160, 0x01a}, {
177 0x162, 0x11a}, {
178 0x164, 0x09a}, {
179 0x168, 0x05a}, {
180 0x16a, 0x15a}, {
181 0x16c, 0x0da}, {
182 0x170, 0x03a}, {
183 0x172, 0x13a}, {
184 0x174, 0x0ba}, {
185 0x178, 0x07a}, {
186 0x17c, 0x0fa}, {
187 0x180, 0x006}, {
188 0x182, 0x106}, {
189 0x184, 0x086}, {
190 0x188, 0x046}, {
191 0x18a, 0x146}, {
192 0x18c, 0x0c6}, {
193 0x190, 0x026}, {
194 0x192, 0x126}, {
195 0x194, 0x0a6}, {
196 0x198, 0x066}, {
197 0x19a, 0x166}, {
198 0x19c, 0x0e6}, {
199 0x1a0, 0x016}, {
200 0x1a2, 0x116}, {
201 0x1a4, 0x096}, {
202 0x1a6, 0x196}, {
203 0x1a8, 0x056}, {
204 0x1aa, 0x156}, {
205 0x1ac, 0x0d6}, {
206 0x1b0, 0x036}, {
207 0x1b2, 0x136}, {
208 0x1b4, 0x0b6}, {
209 0x1b8, 0x076}, {
210 0x1ba, 0x176}, {
211 0x1bc, 0x0f6}, {
212 0x1c0, 0x00e}, {
213 0x1c2, 0x10e}, {
214 0x1c4, 0x08e}, {
215 0x1c6, 0x18e}, {
216 0x1c8, 0x04e}, {
217 0x1ca, 0x14e}, {
218 0x1cc, 0x0ce}, {
219 0x1d0, 0x02e}, {
220 0x1d2, 0x12e}, {
221 0x1d4, 0x0ae}, {
222 0x1d6, 0x1ae}, {
223 0x1d8, 0x06e}, {
224 0x1da, 0x16e}, {
225 0x1dc, 0x0ee}, {
226 0x1e0, 0x01e}, {
227 0x1e2, 0x11e}, {
228 0x1e4, 0x09e}, {
229 0x1e6, 0x19e}, {
230 0x1e8, 0x05e}, {
231 0x1ea, 0x15e}, {
232 0x1ec, 0x0de}, {
233 0x1ee, 0x1de}, {
234 0x1f0, 0x03e}, {
235 0x1f2, 0x13e}, {
236 0x1f4, 0x0be}, {
237 0x1f6, 0x1be}, {
238 0x1f8, 0x07e}, {
239 0x1fa, 0x17e}, {
240 0x1fc, 0x0fe}, {
241 0x200, 0x001}, {
242 0x202, 0x101}, {
243 0x204, 0x081}, {
244 0x206, 0x181}, {
245 0x208, 0x041}, {
246 0x20a, 0x141}, {
247 0x20c, 0x0c1}, {
248 0x20e, 0x1c1}, {
249 0x210, 0x021}, {
250 0x212, 0x121}, {
251 0x214, 0x0a1}, {
252 0x216, 0x1a1}, {
253 0x218, 0x061}, {
254 0x21a, 0x161}, {
255 0x21c, 0x0e1}, {
256 0x21e, 0x1e1}, {
257 0x220, 0x011}, {
258 0x221, 0x211}, {
259 0x222, 0x111}, {
260 0x224, 0x091}, {
261 0x226, 0x191}, {
262 0x228, 0x051}, {
263 0x22a, 0x151}, {
264 0x22c, 0x0d1}, {
265 0x22e, 0x1d1}, {
266 0x230, 0x031}, {
267 0x232, 0x131}, {
268 0x234, 0x0b1}, {
269 0x236, 0x1b1}, {
270 0x238, 0x071}, {
271 0x23a, 0x171}, {
272 0x23c, 0x0f1}, {
273 0x23e, 0x1f1}, {
274 0x240, 0x009}, {
275 0x241, 0x209}, {
276 0x242, 0x109}, {
277 0x244, 0x089}, {
278 0x246, 0x189}, {
279 0x248, 0x049}, {
280 0x24a, 0x149}, {
281 0x24c, 0x0c9}, {
282 0x24e, 0x1c9}, {
283 0x250, 0x029}, {
284 0x251, 0x229}, {
285 0x252, 0x129}, {
286 0x254, 0x0a9}, {
287 0x256, 0x1a9}, {
288 0x258, 0x069}, {
289 0x25a, 0x169}, {
290 0x25c, 0x0e9}, {
291 0x25e, 0x1e9}, {
292 0x260, 0x019}, {
293 0x261, 0x219}, {
294 0x262, 0x119}, {
295 0x264, 0x099}, {
296 0x266, 0x199}, {
297 0x268, 0x059}, {
298 0x269, 0x259}, {
299 0x26a, 0x159}, {
300 0x26c, 0x0d9}, {
301 0x26e, 0x1d9}, {
302 0x270, 0x039}, {
303 0x271, 0x239}, {
304 0x272, 0x139}, {
305 0x274, 0x0b9}, {
306 0x276, 0x1b9}, {
307 0x278, 0x079}, {
308 0x27a, 0x179}, {
309 0x27c, 0x0f9}, {
310 0x27e, 0x1f9}, {
311 0x280, 0x005}, {
312 0x281, 0x205}, {
313 0x282, 0x105}, {
314 0x284, 0x085}, {
315 0x286, 0x185}, {
316 0x288, 0x045}, {
317 0x289, 0x245}, {
318 0x28a, 0x145}, {
319 0x28c, 0x0c5}, {
320 0x28e, 0x1c5}, {
321 0x290, 0x025}, {
322 0x291, 0x225}, {
323 0x292, 0x125}, {
324 0x294, 0x0a5}, {
325 0x296, 0x1a5}, {
326 0x298, 0x065}, {
327 0x299, 0x265}, {
328 0x29a, 0x165}, {
329 0x29c, 0x0e5}, {
330 0x29e, 0x1e5}, {
331 0x2a0, 0x015}, {
332 0x2a1, 0x215}, {
333 0x2a2, 0x115}, {
334 0x2a4, 0x095}, {
335 0x2a5, 0x295}, {
336 0x2a6, 0x195}, {
337 0x2a8, 0x055}, {
338 0x2a9, 0x255}, {
339 0x2aa, 0x155}, {
340 0x2ac, 0x0d5}, {
341 0x2ae, 0x1d5}, {
342 0x2b0, 0x035}, {
343 0x2b1, 0x235}, {
344 0x2b2, 0x135}, {
345 0x2b4, 0x0b5}, {
346 0x2b6, 0x1b5}, {
347 0x2b8, 0x075}, {
348 0x2b9, 0x275}, {
349 0x2ba, 0x175}, {
350 0x2bc, 0x0f5}, {
351 0x2be, 0x1f5}, {
352 0x2c0, 0x00d}, {
353 0x2c1, 0x20d}, {
354 0x2c2, 0x10d}, {
355 0x2c4, 0x08d}, {
356 0x2c5, 0x28d}, {
357 0x2c6, 0x18d}, {
358 0x2c8, 0x04d}, {
359 0x2c9, 0x24d}, {
360 0x2ca, 0x14d}, {
361 0x2cc, 0x0cd}, {
362 0x2ce, 0x1cd}, {
363 0x2d0, 0x02d}, {
364 0x2d1, 0x22d}, {
365 0x2d2, 0x12d}, {
366 0x2d4, 0x0ad}, {
367 0x2d5, 0x2ad}, {
368 0x2d6, 0x1ad}, {
369 0x2d8, 0x06d}, {
370 0x2d9, 0x26d}, {
371 0x2da, 0x16d}, {
372 0x2dc, 0x0ed}, {
373 0x2de, 0x1ed}, {
374 0x2e0, 0x01d}, {
375 0x2e1, 0x21d}, {
376 0x2e2, 0x11d}, {
377 0x2e4, 0x09d}, {
378 0x2e5, 0x29d}, {
379 0x2e6, 0x19d}, {
380 0x2e8, 0x05d}, {
381 0x2e9, 0x25d}, {
382 0x2ea, 0x15d}, {
383 0x2ec, 0x0dd}, {
384 0x2ed, 0x2dd}, {
385 0x2ee, 0x1dd}, {
386 0x2f0, 0x03d}, {
387 0x2f1, 0x23d}, {
388 0x2f2, 0x13d}, {
389 0x2f4, 0x0bd}, {
390 0x2f5, 0x2bd}, {
391 0x2f6, 0x1bd}, {
392 0x2f8, 0x07d}, {
393 0x2f9, 0x27d}, {
394 0x2fa, 0x17d}, {
395 0x2fc, 0x0fd}, {
396 0x2fe, 0x1fd}, {
397 0x300, 0x003}, {
398 0x301, 0x203}, {
399 0x302, 0x103}, {
400 0x304, 0x083}, {
401 0x305, 0x283}, {
402 0x306, 0x183}, {
403 0x308, 0x043}, {
404 0x309, 0x243}, {
405 0x30a, 0x143}, {
406 0x30c, 0x0c3}, {
407 0x30d, 0x2c3}, {
408 0x30e, 0x1c3}, {
409 0x310, 0x023}, {
410 0x311, 0x223}, {
411 0x312, 0x123}, {
412 0x314, 0x0a3}, {
413 0x315, 0x2a3}, {
414 0x316, 0x1a3}, {
415 0x318, 0x063}, {
416 0x319, 0x263}, {
417 0x31a, 0x163}, {
418 0x31c, 0x0e3}, {
419 0x31d, 0x2e3}, {
420 0x31e, 0x1e3}, {
421 0x320, 0x013}, {
422 0x321, 0x213}, {
423 0x322, 0x113}, {
424 0x323, 0x313}, {
425 0x324, 0x093}, {
426 0x325, 0x293}, {
427 0x326, 0x193}, {
428 0x328, 0x053}, {
429 0x329, 0x253}, {
430 0x32a, 0x153}, {
431 0x32c, 0x0d3}, {
432 0x32d, 0x2d3}, {
433 0x32e, 0x1d3}, {
434 0x330, 0x033}, {
435 0x331, 0x233}, {
436 0x332, 0x133}, {
437 0x334, 0x0b3}, {
438 0x335, 0x2b3}, {
439 0x336, 0x1b3}, {
440 0x338, 0x073}, {
441 0x339, 0x273}, {
442 0x33a, 0x173}, {
443 0x33c, 0x0f3}, {
444 0x33d, 0x2f3}, {
445 0x33e, 0x1f3}, {
446 0x340, 0x00b}, {
447 0x341, 0x20b}, {
448 0x342, 0x10b}, {
449 0x343, 0x30b}, {
450 0x344, 0x08b}, {
451 0x345, 0x28b}, {
452 0x346, 0x18b}, {
453 0x348, 0x04b}, {
454 0x349, 0x24b}, {
455 0x34a, 0x14b}, {
456 0x34c, 0x0cb}, {
457 0x34d, 0x2cb}, {
458 0x34e, 0x1cb}, {
459 0x350, 0x02b}, {
460 0x351, 0x22b}, {
461 0x352, 0x12b}, {
462 0x353, 0x32b}, {
463 0x354, 0x0ab}, {
464 0x355, 0x2ab}, {
465 0x356, 0x1ab}, {
466 0x358, 0x06b}, {
467 0x359, 0x26b}, {
468 0x35a, 0x16b}, {
469 0x35c, 0x0eb}, {
470 0x35d, 0x2eb}, {
471 0x35e, 0x1eb}, {
472 0x360, 0x01b}, {
473 0x361, 0x21b}, {
474 0x362, 0x11b}, {
475 0x363, 0x31b}, {
476 0x364, 0x09b}, {
477 0x365, 0x29b}, {
478 0x366, 0x19b}, {
479 0x368, 0x05b}, {
480 0x369, 0x25b}, {
481 0x36a, 0x15b}, {
482 0x36b, 0x35b}, {
483 0x36c, 0x0db}, {
484 0x36d, 0x2db}, {
485 0x36e, 0x1db}, {
486 0x370, 0x03b}, {
487 0x371, 0x23b}, {
488 0x372, 0x13b}, {
489 0x373, 0x33b}, {
490 0x374, 0x0bb}, {
491 0x375, 0x2bb}, {
492 0x376, 0x1bb}, {
493 0x378, 0x07b}, {
494 0x379, 0x27b}, {
495 0x37a, 0x17b}, {
496 0x37c, 0x0fb}, {
497 0x37d, 0x2fb}, {
498 0x37e, 0x1fb}, {
499 0x380, 0x007}, {
500 0x381, 0x207}, {
501 0x382, 0x107}, {
502 0x383, 0x307}, {
503 0x384, 0x087}, {
504 0x385, 0x287}, {
505 0x386, 0x187}, {
506 0x388, 0x047}, {
507 0x389, 0x247}, {
508 0x38a, 0x147}, {
509 0x38b, 0x347}, {
510 0x38c, 0x0c7}, {
511 0x38d, 0x2c7}, {
512 0x38e, 0x1c7}, {
513 0x390, 0x027}, {
514 0x391, 0x227}, {
515 0x392, 0x127}, {
516 0x393, 0x327}, {
517 0x394, 0x0a7}, {
518 0x395, 0x2a7}, {
519 0x396, 0x1a7}, {
520 0x398, 0x067}, {
521 0x399, 0x267}, {
522 0x39a, 0x167}, {
523 0x39b, 0x367}, {
524 0x39c, 0x0e7}, {
525 0x39d, 0x2e7}, {
526 0x39e, 0x1e7}, {
527 0x3a0, 0x017}, {
528 0x3a1, 0x217}, {
529 0x3a2, 0x117}, {
530 0x3a3, 0x317}, {
531 0x3a4, 0x097}, {
532 0x3a5, 0x297}, {
533 0x3a6, 0x197}, {
534 0x3a7, 0x397}, {
535 0x3a8, 0x057}, {
536 0x3a9, 0x257}, {
537 0x3aa, 0x157}, {
538 0x3ab, 0x357}, {
539 0x3ac, 0x0d7}, {
540 0x3ad, 0x2d7}, {
541 0x3ae, 0x1d7}, {
542 0x3b0, 0x037}, {
543 0x3b1, 0x237}, {
544 0x3b2, 0x137}, {
545 0x3b3, 0x337}, {
546 0x3b4, 0x0b7}, {
547 0x3b5, 0x2b7}, {
548 0x3b6, 0x1b7}, {
549 0x3b8, 0x077}, {
550 0x3b9, 0x277}, {
551 0x3ba, 0x177}, {
552 0x3bb, 0x377}, {
553 0x3bc, 0x0f7}, {
554 0x3bd, 0x2f7}, {
555 0x3be, 0x1f7}, {
556 0x3c0, 0x00f}, {
557 0x3c1, 0x20f}, {
558 0x3c2, 0x10f}, {
559 0x3c3, 0x30f}, {
560 0x3c4, 0x08f}, {
561 0x3c5, 0x28f}, {
562 0x3c6, 0x18f}, {
563 0x3c7, 0x38f}, {
564 0x3c8, 0x04f}, {
565 0x3c9, 0x24f}, {
566 0x3ca, 0x14f}, {
567 0x3cb, 0x34f}, {
568 0x3cc, 0x0cf}, {
569 0x3cd, 0x2cf}, {
570 0x3ce, 0x1cf}, {
571 0x3d0, 0x02f}, {
572 0x3d1, 0x22f}, {
573 0x3d2, 0x12f}, {
574 0x3d3, 0x32f}, {
575 0x3d4, 0x0af}, {
576 0x3d5, 0x2af}, {
577 0x3d6, 0x1af}, {
578 0x3d7, 0x3af}, {
579 0x3d8, 0x06f}, {
580 0x3d9, 0x26f}, {
581 0x3da, 0x16f}, {
582 0x3db, 0x36f}, {
583 0x3dc, 0x0ef}, {
584 0x3dd, 0x2ef}, {
585 0x3de, 0x1ef}, {
586 0x3e0, 0x01f}, {
587 0x3e1, 0x21f}, {
588 0x3e2, 0x11f}, {
589 0x3e3, 0x31f}, {
590 0x3e4, 0x09f}, {
591 0x3e5, 0x29f}, {
592 0x3e6, 0x19f}, {
593 0x3e7, 0x39f}, {
594 0x3e8, 0x05f}, {
595 0x3e9, 0x25f}, {
596 0x3ea, 0x15f}, {
597 0x3eb, 0x35f}, {
598 0x3ec, 0x0df}, {
599 0x3ed, 0x2df}, {
600 0x3ee, 0x1df}, {
601 0x3ef, 0x3df}, {
602 0x3f0, 0x03f}, {
603 0x3f1, 0x23f}, {
604 0x3f2, 0x13f}, {
605 0x3f3, 0x33f}, {
606 0x3f4, 0x0bf}, {
607 0x3f5, 0x2bf}, {
608 0x3f6, 0x1bf}, {
609 0x3f7, 0x3bf}, {
610 0x3f8, 0x07f}, {
611 0x3f9, 0x27f}, {
612 0x3fa, 0x17f}, {
613 0x3fb, 0x37f}, {
614 0x3fc, 0x0ff}, {
615 0x3fd, 0x2ff}, {
616 0x3fe, 0x1ff}
617 };
618
619
620 {
621 int i;
622 for (i = 0; i < sizeof k1k2tab / sizeof k1k2tab[0]; ++i) {
623 k1 = k1k2tab[i].k1;
624 k2 = k1k2tab[i].k2;
625 a = fz[k1];
626 fz[k1] = fz[k2];
627 fz[k2] = a;
628 }
629 }
630
631 for (fi = fz, fn = fz + 1024; fi < fn; fi += 4) {
632 FLOAT f0, f1, f2, f3;
633 f1 = fi[0] - fi[1];
634 f0 = fi[0] + fi[1];
635 f3 = fi[2] - fi[3];
636 f2 = fi[2] + fi[3];
637 fi[2] = (f0 - f2);
638 fi[0] = (f0 + f2);
639 fi[3] = (f1 - f3);
640 fi[1] = (f1 + f3);
641 }
642
643 k = 0;
644 do {
645 FLOAT s1, c1;
646 k += 2;
647 k1 = 1 << k;
648 k2 = k1 << 1;
649 k4 = k2 << 1;
650 k3 = k2 + k1;
651 kx = k1 >> 1;
652 fi = fz;
653 gi = fi + kx;
654 fn = fz + 1024;
655 do {
656 FLOAT g0, f0, f1, g1, f2, g2, f3, g3;
657 f1 = fi[0] - fi[k1];
658 f0 = fi[0] + fi[k1];
659 f3 = fi[k2] - fi[k3];
660 f2 = fi[k2] + fi[k3];
661 fi[k2] = f0 - f2;
662 fi[0] = f0 + f2;
663 fi[k3] = f1 - f3;
664 fi[k1] = f1 + f3;
665 g1 = gi[0] - gi[k1];
666 g0 = gi[0] + gi[k1];
667 g3 = SQRT2 * gi[k3];
668 g2 = SQRT2 * gi[k2];
669 gi[k2] = g0 - g2;
670 gi[0] = g0 + g2;
671 gi[k3] = g1 - g3;
672 gi[k1] = g1 + g3;
673 gi += k4;
674 fi += k4;
675 }
676 while (fi < fn);
677
678 t_c = costab[k];
679 t_s = sintab[k];
680 c1 = 1;
681 s1 = 0;
682 for (i = 1; i < kx; i++) {
683 FLOAT c2, s2;
684 FLOAT t = c1;
685 c1 = t * t_c - s1 * t_s;
686 s1 = t * t_s + s1 * t_c;
687 c2 = c1 * c1 - s1 * s1;
688 s2 = 2 * (c1 * s1);
689 fn = fz + 1024;
690 fi = fz + i;
691 gi = fz + k1 - i;
692 do {
693 FLOAT a, b, g0, f0, f1, g1, f2, g2, f3, g3;
694 b = s2 * fi[k1] - c2 * gi[k1];
695 a = c2 * fi[k1] + s2 * gi[k1];
696 f1 = fi[0] - a;
697 f0 = fi[0] + a;
698 g1 = gi[0] - b;
699 g0 = gi[0] + b;
700 b = s2 * fi[k3] - c2 * gi[k3];
701 a = c2 * fi[k3] + s2 * gi[k3];
702 f3 = fi[k2] - a;
703 f2 = fi[k2] + a;
704 g3 = gi[k2] - b;
705 g2 = gi[k2] + b;
706 b = s1 * f2 - c1 * g3;
707 a = c1 * f2 + s1 * g3;
708 fi[k2] = f0 - a;
709 fi[0] = f0 + a;
710 gi[k3] = g1 - b;
711 gi[k1] = g1 + b;
712 b = c1 * g2 - s1 * f3;
713 a = s1 * g2 + c1 * f3;
714 gi[k2] = g0 - a;
715 gi[0] = g0 + a;
716 fi[k3] = f1 - b;
717 fi[k1] = f1 + b;
718 gi += k4;
719 fi += k4;
720 }
721 while (fi < fn);
722 }
723 }
724 while (k4 < 1024);
725 }
726
727 #ifdef NEWATAN
728 #define ATANSIZE 6000
729 #define ATANSCALE 100.0
730 /* Create a table of ATAN2 values.
731 Valid for ratios of (y/x) from 0 to ATANSIZE/ATANSCALE (60)
732 Largest absolute error in angle: 0.0167 radians i.e. ATANSCALE/ATANSIZE
733 Depending on how you want to trade off speed/accuracy and mem usage, twiddle the defines
734 MFC March 2003 */
735 static FLOAT atan_t[ATANSIZE];
736
atan_table(FLOAT y,FLOAT x)737 static inline FLOAT atan_table(FLOAT y, FLOAT x)
738 {
739 int index;
740
741 index = (int) (ATANSCALE * fabs(y / x));
742 if (index >= ATANSIZE)
743 index = ATANSIZE - 1;
744
745 /* Have to work out the correct quadrant as well */
746 if (y > 0 && x < 0)
747 return (PI - atan_t[index]);
748
749 if (y < 0 && x > 0)
750 return (-atan_t[index]);
751
752 if (y < 0 && x < 0)
753 return (atan_t[index] - PI);
754
755 return (atan_t[index]);
756 }
757
atan_table_init(void)758 static void atan_table_init(void)
759 {
760 int i;
761 for (i = 0; i < ATANSIZE; i++)
762 atan_t[i] = atan((FLOAT) i / ATANSCALE);
763 }
764
765 #endif // NEWATAN
766
767 /* For variations on psycho model 2:
768 N always equals 1024
769 BUT in the returned values, no energy/phi is used at or above an index of 513 */
psycho_2_fft(FLOAT * x_real,FLOAT * energy,FLOAT * phi)770 void psycho_2_fft(FLOAT * x_real, FLOAT * energy, FLOAT * phi)
771 /* got rid of size "N" argument as it is always 1024 for layerII */
772 {
773 FLOAT imag, real;
774 int i, j;
775 #ifdef NEWATAN
776 static int init = 0;
777
778 if (!init) {
779 atan_table_init();
780 init++;
781 }
782 #endif
783
784
785 fht(x_real);
786
787
788 energy[0] = x_real[0] * x_real[0];
789
790 for (i = 1, j = 1023; i < 512; i++, j--) {
791 imag = x_real[i];
792 real = x_real[j];
793 /* MFC FIXME Mar03 Why is this divided by 2.0? if a and b are the real and imaginary
794 components then r = sqrt(a^2 + b^2), but, back in the psycho2 model, they calculate
795 r=sqrt(energy), which, if you look at the original equation below is different */
796 energy[i] = (imag * imag + real * real) / 2.0;
797 if (energy[i] < 0.0005) {
798 energy[i] = 0.0005;
799 phi[i] = 0;
800 } else
801 #ifdef NEWATAN
802 {
803 phi[i] = atan_table(-imag, real) + PI / 4;
804 }
805 #else
806 {
807 phi[i] = atan2(-(FLOAT) imag, (FLOAT) real) + PI / 4;
808 }
809 #endif
810 }
811 energy[512] = x_real[512] * x_real[512];
812 phi[512] = atan2(0.0, (FLOAT) x_real[512]);
813 }
814
815
psycho_1_fft(FLOAT * x_real,FLOAT * energy,int N)816 void psycho_1_fft(FLOAT * x_real, FLOAT * energy, int N)
817 {
818 FLOAT a, b;
819 int i, j;
820
821 fht(x_real);
822
823 energy[0] = x_real[0] * x_real[0];
824
825 for (i = 1, j = N - 1; i < N / 2; i++, j--) {
826 a = x_real[i];
827 b = x_real[j];
828 energy[i] = (a * a + b * b) / 2.0;
829 }
830 energy[N / 2] = x_real[N / 2] * x_real[N / 2];
831 }
832
833
834 // vim:ts=4:sw=4:nowrap:
835