1 /*
2 * TwoLAME: an optimized MPEG Audio Layer Two encoder
3 *
4 * Copyright (C) 2001-2004 Michael Cheng
5 * Copyright (C) 2004-2018 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 */
22
23
24
25 /*
26 ** FFT and FHT routines
27 ** Copyright 1988, 1993; Ron Mayer
28 **
29 ** fht(fz,n);
30 ** Does a hartley transform of "n" points in the array "fz".
31 **
32 ** NOTE: This routine uses at least 2 patented algorithms, and may be
33 ** under the restrictions of a bunch of different organizations.
34 ** Although I wrote it completely myself; it is kind of a derivative
35 ** of a routine I once authored and released under the GPL, so it
36 ** may fall under the free software foundation's restrictions;
37 ** it was worked on as a Stanford Univ project, so they claim
38 ** some rights to it; it was further optimized at work here, so
39 ** I think this company claims parts of it. The patents are
40 ** held by R. Bracewell (the FHT algorithm) and O. Buneman (the
41 ** trig generator), both at Stanford Univ.
42 ** If it were up to me, I'd say go do whatever you want with it;
43 ** but it would be polite to give credit to the following people
44 ** if you use this anywhere:
45 ** Euler - probable inventor of the fourier transform.
46 ** Gauss - probable inventor of the FFT.
47 ** Hartley - probable inventor of the hartley transform.
48 ** Buneman - for a really cool trig generator
49 ** Mayer(me) - for authoring this particular version and
50 ** including all the optimizations in one package.
51 ** Thanks,
52 ** Ron Mayer; mayer@acuson.com
53 **
54 */
55
56 #include <stdio.h>
57 #include <math.h>
58
59 #include "twolame.h"
60 #include "common.h"
61 #include "fft.h"
62
63
64
65 #define SQRT2 1.4142135623730951454746218587388284504414
66
67
68 static const FLOAT costab[20] = {
69 .00000000000000000000000000000000000000000000000000,
70 .70710678118654752440084436210484903928483593768847,
71 .92387953251128675612818318939678828682241662586364,
72 .98078528040323044912618223613423903697393373089333,
73 .99518472667219688624483695310947992157547486872985,
74 .99879545620517239271477160475910069444320361470461,
75 .99969881869620422011576564966617219685006108125772,
76 .99992470183914454092164649119638322435060646880221,
77 .99998117528260114265699043772856771617391725094433,
78 .99999529380957617151158012570011989955298763362218,
79 .99999882345170190992902571017152601904826792288976,
80 .99999970586288221916022821773876567711626389934930,
81 .99999992646571785114473148070738785694820115568892,
82 .99999998161642929380834691540290971450507605124278,
83 .99999999540410731289097193313960614895889430318945,
84 .99999999885102682756267330779455410840053741619428
85 };
86 static const FLOAT sintab[20] = {
87 1.0000000000000000000000000000000000000000000000000,
88 .70710678118654752440084436210484903928483593768846,
89 .38268343236508977172845998403039886676134456248561,
90 .19509032201612826784828486847702224092769161775195,
91 .09801714032956060199419556388864184586113667316749,
92 .04906767432741801425495497694268265831474536302574,
93 .02454122852291228803173452945928292506546611923944,
94 .01227153828571992607940826195100321214037231959176,
95 .00613588464915447535964023459037258091705788631738,
96 .00306795676296597627014536549091984251894461021344,
97 .00153398018628476561230369715026407907995486457522,
98 .00076699031874270452693856835794857664314091945205,
99 .00038349518757139558907246168118138126339502603495,
100 .00019174759731070330743990956198900093346887403385,
101 .00009587379909597734587051721097647635118706561284,
102 .00004793689960306688454900399049465887274686668768
103 };
104
105 /* This is a simplified version for n an even power of 2 */
106 /* MFC: In the case of LayerII encoding, n==1024 always. */
107
fht(FLOAT * fz)108 static void fht(FLOAT * fz)
109 {
110 int i, k, k1, k2, k3, k4, kx;
111 FLOAT *fi, *fn, *gi;
112 FLOAT t_c, t_s;
113
114 FLOAT a;
115 static const struct {
116 unsigned short k1, k2;
117 } k1k2tab[8 * 62] = {
118 {
119 0x020, 0x010
120 }, {
121 0x040, 0x008
122 }, {
123 0x050, 0x028
124 }, {
125 0x060, 0x018
126 }, {
127 0x068, 0x058
128 }, {
129 0x070, 0x038
130 }, {
131 0x080, 0x004
132 }, {
133 0x088, 0x044
134 }, {
135 0x090, 0x024
136 }, {
137 0x098, 0x064
138 }, {
139 0x0a0, 0x014
140 }, {
141 0x0a4, 0x094
142 }, {
143 0x0a8, 0x054
144 }, {
145 0x0b0, 0x034
146 }, {
147 0x0b8, 0x074
148 }, {
149 0x0c0, 0x00c
150 }, {
151 0x0c4, 0x08c
152 }, {
153 0x0c8, 0x04c
154 }, {
155 0x0d0, 0x02c
156 }, {
157 0x0d4, 0x0ac
158 }, {
159 0x0d8, 0x06c
160 }, {
161 0x0e0, 0x01c
162 }, {
163 0x0e4, 0x09c
164 }, {
165 0x0e8, 0x05c
166 }, {
167 0x0ec, 0x0dc
168 }, {
169 0x0f0, 0x03c
170 }, {
171 0x0f4, 0x0bc
172 }, {
173 0x0f8, 0x07c
174 }, {
175 0x100, 0x002
176 }, {
177 0x104, 0x082
178 }, {
179 0x108, 0x042
180 }, {
181 0x10c, 0x0c2
182 }, {
183 0x110, 0x022
184 }, {
185 0x114, 0x0a2
186 }, {
187 0x118, 0x062
188 }, {
189 0x11c, 0x0e2
190 }, {
191 0x120, 0x012
192 }, {
193 0x122, 0x112
194 }, {
195 0x124, 0x092
196 }, {
197 0x128, 0x052
198 }, {
199 0x12c, 0x0d2
200 }, {
201 0x130, 0x032
202 }, {
203 0x134, 0x0b2
204 }, {
205 0x138, 0x072
206 }, {
207 0x13c, 0x0f2
208 }, {
209 0x140, 0x00a
210 }, {
211 0x142, 0x10a
212 }, {
213 0x144, 0x08a
214 }, {
215 0x148, 0x04a
216 }, {
217 0x14c, 0x0ca
218 }, {
219 0x150, 0x02a
220 }, {
221 0x152, 0x12a
222 }, {
223 0x154, 0x0aa
224 }, {
225 0x158, 0x06a
226 }, {
227 0x15c, 0x0ea
228 }, {
229 0x160, 0x01a
230 }, {
231 0x162, 0x11a
232 }, {
233 0x164, 0x09a
234 }, {
235 0x168, 0x05a
236 }, {
237 0x16a, 0x15a
238 }, {
239 0x16c, 0x0da
240 }, {
241 0x170, 0x03a
242 }, {
243 0x172, 0x13a
244 }, {
245 0x174, 0x0ba
246 }, {
247 0x178, 0x07a
248 }, {
249 0x17c, 0x0fa
250 }, {
251 0x180, 0x006
252 }, {
253 0x182, 0x106
254 }, {
255 0x184, 0x086
256 }, {
257 0x188, 0x046
258 }, {
259 0x18a, 0x146
260 }, {
261 0x18c, 0x0c6
262 }, {
263 0x190, 0x026
264 }, {
265 0x192, 0x126
266 }, {
267 0x194, 0x0a6
268 }, {
269 0x198, 0x066
270 }, {
271 0x19a, 0x166
272 }, {
273 0x19c, 0x0e6
274 }, {
275 0x1a0, 0x016
276 }, {
277 0x1a2, 0x116
278 }, {
279 0x1a4, 0x096
280 }, {
281 0x1a6, 0x196
282 }, {
283 0x1a8, 0x056
284 }, {
285 0x1aa, 0x156
286 }, {
287 0x1ac, 0x0d6
288 }, {
289 0x1b0, 0x036
290 }, {
291 0x1b2, 0x136
292 }, {
293 0x1b4, 0x0b6
294 }, {
295 0x1b8, 0x076
296 }, {
297 0x1ba, 0x176
298 }, {
299 0x1bc, 0x0f6
300 }, {
301 0x1c0, 0x00e
302 }, {
303 0x1c2, 0x10e
304 }, {
305 0x1c4, 0x08e
306 }, {
307 0x1c6, 0x18e
308 }, {
309 0x1c8, 0x04e
310 }, {
311 0x1ca, 0x14e
312 }, {
313 0x1cc, 0x0ce
314 }, {
315 0x1d0, 0x02e
316 }, {
317 0x1d2, 0x12e
318 }, {
319 0x1d4, 0x0ae
320 }, {
321 0x1d6, 0x1ae
322 }, {
323 0x1d8, 0x06e
324 }, {
325 0x1da, 0x16e
326 }, {
327 0x1dc, 0x0ee
328 }, {
329 0x1e0, 0x01e
330 }, {
331 0x1e2, 0x11e
332 }, {
333 0x1e4, 0x09e
334 }, {
335 0x1e6, 0x19e
336 }, {
337 0x1e8, 0x05e
338 }, {
339 0x1ea, 0x15e
340 }, {
341 0x1ec, 0x0de
342 }, {
343 0x1ee, 0x1de
344 }, {
345 0x1f0, 0x03e
346 }, {
347 0x1f2, 0x13e
348 }, {
349 0x1f4, 0x0be
350 }, {
351 0x1f6, 0x1be
352 }, {
353 0x1f8, 0x07e
354 }, {
355 0x1fa, 0x17e
356 }, {
357 0x1fc, 0x0fe
358 }, {
359 0x200, 0x001
360 }, {
361 0x202, 0x101
362 }, {
363 0x204, 0x081
364 }, {
365 0x206, 0x181
366 }, {
367 0x208, 0x041
368 }, {
369 0x20a, 0x141
370 }, {
371 0x20c, 0x0c1
372 }, {
373 0x20e, 0x1c1
374 }, {
375 0x210, 0x021
376 }, {
377 0x212, 0x121
378 }, {
379 0x214, 0x0a1
380 }, {
381 0x216, 0x1a1
382 }, {
383 0x218, 0x061
384 }, {
385 0x21a, 0x161
386 }, {
387 0x21c, 0x0e1
388 }, {
389 0x21e, 0x1e1
390 }, {
391 0x220, 0x011
392 }, {
393 0x221, 0x211
394 }, {
395 0x222, 0x111
396 }, {
397 0x224, 0x091
398 }, {
399 0x226, 0x191
400 }, {
401 0x228, 0x051
402 }, {
403 0x22a, 0x151
404 }, {
405 0x22c, 0x0d1
406 }, {
407 0x22e, 0x1d1
408 }, {
409 0x230, 0x031
410 }, {
411 0x232, 0x131
412 }, {
413 0x234, 0x0b1
414 }, {
415 0x236, 0x1b1
416 }, {
417 0x238, 0x071
418 }, {
419 0x23a, 0x171
420 }, {
421 0x23c, 0x0f1
422 }, {
423 0x23e, 0x1f1
424 }, {
425 0x240, 0x009
426 }, {
427 0x241, 0x209
428 }, {
429 0x242, 0x109
430 }, {
431 0x244, 0x089
432 }, {
433 0x246, 0x189
434 }, {
435 0x248, 0x049
436 }, {
437 0x24a, 0x149
438 }, {
439 0x24c, 0x0c9
440 }, {
441 0x24e, 0x1c9
442 }, {
443 0x250, 0x029
444 }, {
445 0x251, 0x229
446 }, {
447 0x252, 0x129
448 }, {
449 0x254, 0x0a9
450 }, {
451 0x256, 0x1a9
452 }, {
453 0x258, 0x069
454 }, {
455 0x25a, 0x169
456 }, {
457 0x25c, 0x0e9
458 }, {
459 0x25e, 0x1e9
460 }, {
461 0x260, 0x019
462 }, {
463 0x261, 0x219
464 }, {
465 0x262, 0x119
466 }, {
467 0x264, 0x099
468 }, {
469 0x266, 0x199
470 }, {
471 0x268, 0x059
472 }, {
473 0x269, 0x259
474 }, {
475 0x26a, 0x159
476 }, {
477 0x26c, 0x0d9
478 }, {
479 0x26e, 0x1d9
480 }, {
481 0x270, 0x039
482 }, {
483 0x271, 0x239
484 }, {
485 0x272, 0x139
486 }, {
487 0x274, 0x0b9
488 }, {
489 0x276, 0x1b9
490 }, {
491 0x278, 0x079
492 }, {
493 0x27a, 0x179
494 }, {
495 0x27c, 0x0f9
496 }, {
497 0x27e, 0x1f9
498 }, {
499 0x280, 0x005
500 }, {
501 0x281, 0x205
502 }, {
503 0x282, 0x105
504 }, {
505 0x284, 0x085
506 }, {
507 0x286, 0x185
508 }, {
509 0x288, 0x045
510 }, {
511 0x289, 0x245
512 }, {
513 0x28a, 0x145
514 }, {
515 0x28c, 0x0c5
516 }, {
517 0x28e, 0x1c5
518 }, {
519 0x290, 0x025
520 }, {
521 0x291, 0x225
522 }, {
523 0x292, 0x125
524 }, {
525 0x294, 0x0a5
526 }, {
527 0x296, 0x1a5
528 }, {
529 0x298, 0x065
530 }, {
531 0x299, 0x265
532 }, {
533 0x29a, 0x165
534 }, {
535 0x29c, 0x0e5
536 }, {
537 0x29e, 0x1e5
538 }, {
539 0x2a0, 0x015
540 }, {
541 0x2a1, 0x215
542 }, {
543 0x2a2, 0x115
544 }, {
545 0x2a4, 0x095
546 }, {
547 0x2a5, 0x295
548 }, {
549 0x2a6, 0x195
550 }, {
551 0x2a8, 0x055
552 }, {
553 0x2a9, 0x255
554 }, {
555 0x2aa, 0x155
556 }, {
557 0x2ac, 0x0d5
558 }, {
559 0x2ae, 0x1d5
560 }, {
561 0x2b0, 0x035
562 }, {
563 0x2b1, 0x235
564 }, {
565 0x2b2, 0x135
566 }, {
567 0x2b4, 0x0b5
568 }, {
569 0x2b6, 0x1b5
570 }, {
571 0x2b8, 0x075
572 }, {
573 0x2b9, 0x275
574 }, {
575 0x2ba, 0x175
576 }, {
577 0x2bc, 0x0f5
578 }, {
579 0x2be, 0x1f5
580 }, {
581 0x2c0, 0x00d
582 }, {
583 0x2c1, 0x20d
584 }, {
585 0x2c2, 0x10d
586 }, {
587 0x2c4, 0x08d
588 }, {
589 0x2c5, 0x28d
590 }, {
591 0x2c6, 0x18d
592 }, {
593 0x2c8, 0x04d
594 }, {
595 0x2c9, 0x24d
596 }, {
597 0x2ca, 0x14d
598 }, {
599 0x2cc, 0x0cd
600 }, {
601 0x2ce, 0x1cd
602 }, {
603 0x2d0, 0x02d
604 }, {
605 0x2d1, 0x22d
606 }, {
607 0x2d2, 0x12d
608 }, {
609 0x2d4, 0x0ad
610 }, {
611 0x2d5, 0x2ad
612 }, {
613 0x2d6, 0x1ad
614 }, {
615 0x2d8, 0x06d
616 }, {
617 0x2d9, 0x26d
618 }, {
619 0x2da, 0x16d
620 }, {
621 0x2dc, 0x0ed
622 }, {
623 0x2de, 0x1ed
624 }, {
625 0x2e0, 0x01d
626 }, {
627 0x2e1, 0x21d
628 }, {
629 0x2e2, 0x11d
630 }, {
631 0x2e4, 0x09d
632 }, {
633 0x2e5, 0x29d
634 }, {
635 0x2e6, 0x19d
636 }, {
637 0x2e8, 0x05d
638 }, {
639 0x2e9, 0x25d
640 }, {
641 0x2ea, 0x15d
642 }, {
643 0x2ec, 0x0dd
644 }, {
645 0x2ed, 0x2dd
646 }, {
647 0x2ee, 0x1dd
648 }, {
649 0x2f0, 0x03d
650 }, {
651 0x2f1, 0x23d
652 }, {
653 0x2f2, 0x13d
654 }, {
655 0x2f4, 0x0bd
656 }, {
657 0x2f5, 0x2bd
658 }, {
659 0x2f6, 0x1bd
660 }, {
661 0x2f8, 0x07d
662 }, {
663 0x2f9, 0x27d
664 }, {
665 0x2fa, 0x17d
666 }, {
667 0x2fc, 0x0fd
668 }, {
669 0x2fe, 0x1fd
670 }, {
671 0x300, 0x003
672 }, {
673 0x301, 0x203
674 }, {
675 0x302, 0x103
676 }, {
677 0x304, 0x083
678 }, {
679 0x305, 0x283
680 }, {
681 0x306, 0x183
682 }, {
683 0x308, 0x043
684 }, {
685 0x309, 0x243
686 }, {
687 0x30a, 0x143
688 }, {
689 0x30c, 0x0c3
690 }, {
691 0x30d, 0x2c3
692 }, {
693 0x30e, 0x1c3
694 }, {
695 0x310, 0x023
696 }, {
697 0x311, 0x223
698 }, {
699 0x312, 0x123
700 }, {
701 0x314, 0x0a3
702 }, {
703 0x315, 0x2a3
704 }, {
705 0x316, 0x1a3
706 }, {
707 0x318, 0x063
708 }, {
709 0x319, 0x263
710 }, {
711 0x31a, 0x163
712 }, {
713 0x31c, 0x0e3
714 }, {
715 0x31d, 0x2e3
716 }, {
717 0x31e, 0x1e3
718 }, {
719 0x320, 0x013
720 }, {
721 0x321, 0x213
722 }, {
723 0x322, 0x113
724 }, {
725 0x323, 0x313
726 }, {
727 0x324, 0x093
728 }, {
729 0x325, 0x293
730 }, {
731 0x326, 0x193
732 }, {
733 0x328, 0x053
734 }, {
735 0x329, 0x253
736 }, {
737 0x32a, 0x153
738 }, {
739 0x32c, 0x0d3
740 }, {
741 0x32d, 0x2d3
742 }, {
743 0x32e, 0x1d3
744 }, {
745 0x330, 0x033
746 }, {
747 0x331, 0x233
748 }, {
749 0x332, 0x133
750 }, {
751 0x334, 0x0b3
752 }, {
753 0x335, 0x2b3
754 }, {
755 0x336, 0x1b3
756 }, {
757 0x338, 0x073
758 }, {
759 0x339, 0x273
760 }, {
761 0x33a, 0x173
762 }, {
763 0x33c, 0x0f3
764 }, {
765 0x33d, 0x2f3
766 }, {
767 0x33e, 0x1f3
768 }, {
769 0x340, 0x00b
770 }, {
771 0x341, 0x20b
772 }, {
773 0x342, 0x10b
774 }, {
775 0x343, 0x30b
776 }, {
777 0x344, 0x08b
778 }, {
779 0x345, 0x28b
780 }, {
781 0x346, 0x18b
782 }, {
783 0x348, 0x04b
784 }, {
785 0x349, 0x24b
786 }, {
787 0x34a, 0x14b
788 }, {
789 0x34c, 0x0cb
790 }, {
791 0x34d, 0x2cb
792 }, {
793 0x34e, 0x1cb
794 }, {
795 0x350, 0x02b
796 }, {
797 0x351, 0x22b
798 }, {
799 0x352, 0x12b
800 }, {
801 0x353, 0x32b
802 }, {
803 0x354, 0x0ab
804 }, {
805 0x355, 0x2ab
806 }, {
807 0x356, 0x1ab
808 }, {
809 0x358, 0x06b
810 }, {
811 0x359, 0x26b
812 }, {
813 0x35a, 0x16b
814 }, {
815 0x35c, 0x0eb
816 }, {
817 0x35d, 0x2eb
818 }, {
819 0x35e, 0x1eb
820 }, {
821 0x360, 0x01b
822 }, {
823 0x361, 0x21b
824 }, {
825 0x362, 0x11b
826 }, {
827 0x363, 0x31b
828 }, {
829 0x364, 0x09b
830 }, {
831 0x365, 0x29b
832 }, {
833 0x366, 0x19b
834 }, {
835 0x368, 0x05b
836 }, {
837 0x369, 0x25b
838 }, {
839 0x36a, 0x15b
840 }, {
841 0x36b, 0x35b
842 }, {
843 0x36c, 0x0db
844 }, {
845 0x36d, 0x2db
846 }, {
847 0x36e, 0x1db
848 }, {
849 0x370, 0x03b
850 }, {
851 0x371, 0x23b
852 }, {
853 0x372, 0x13b
854 }, {
855 0x373, 0x33b
856 }, {
857 0x374, 0x0bb
858 }, {
859 0x375, 0x2bb
860 }, {
861 0x376, 0x1bb
862 }, {
863 0x378, 0x07b
864 }, {
865 0x379, 0x27b
866 }, {
867 0x37a, 0x17b
868 }, {
869 0x37c, 0x0fb
870 }, {
871 0x37d, 0x2fb
872 }, {
873 0x37e, 0x1fb
874 }, {
875 0x380, 0x007
876 }, {
877 0x381, 0x207
878 }, {
879 0x382, 0x107
880 }, {
881 0x383, 0x307
882 }, {
883 0x384, 0x087
884 }, {
885 0x385, 0x287
886 }, {
887 0x386, 0x187
888 }, {
889 0x388, 0x047
890 }, {
891 0x389, 0x247
892 }, {
893 0x38a, 0x147
894 }, {
895 0x38b, 0x347
896 }, {
897 0x38c, 0x0c7
898 }, {
899 0x38d, 0x2c7
900 }, {
901 0x38e, 0x1c7
902 }, {
903 0x390, 0x027
904 }, {
905 0x391, 0x227
906 }, {
907 0x392, 0x127
908 }, {
909 0x393, 0x327
910 }, {
911 0x394, 0x0a7
912 }, {
913 0x395, 0x2a7
914 }, {
915 0x396, 0x1a7
916 }, {
917 0x398, 0x067
918 }, {
919 0x399, 0x267
920 }, {
921 0x39a, 0x167
922 }, {
923 0x39b, 0x367
924 }, {
925 0x39c, 0x0e7
926 }, {
927 0x39d, 0x2e7
928 }, {
929 0x39e, 0x1e7
930 }, {
931 0x3a0, 0x017
932 }, {
933 0x3a1, 0x217
934 }, {
935 0x3a2, 0x117
936 }, {
937 0x3a3, 0x317
938 }, {
939 0x3a4, 0x097
940 }, {
941 0x3a5, 0x297
942 }, {
943 0x3a6, 0x197
944 }, {
945 0x3a7, 0x397
946 }, {
947 0x3a8, 0x057
948 }, {
949 0x3a9, 0x257
950 }, {
951 0x3aa, 0x157
952 }, {
953 0x3ab, 0x357
954 }, {
955 0x3ac, 0x0d7
956 }, {
957 0x3ad, 0x2d7
958 }, {
959 0x3ae, 0x1d7
960 }, {
961 0x3b0, 0x037
962 }, {
963 0x3b1, 0x237
964 }, {
965 0x3b2, 0x137
966 }, {
967 0x3b3, 0x337
968 }, {
969 0x3b4, 0x0b7
970 }, {
971 0x3b5, 0x2b7
972 }, {
973 0x3b6, 0x1b7
974 }, {
975 0x3b8, 0x077
976 }, {
977 0x3b9, 0x277
978 }, {
979 0x3ba, 0x177
980 }, {
981 0x3bb, 0x377
982 }, {
983 0x3bc, 0x0f7
984 }, {
985 0x3bd, 0x2f7
986 }, {
987 0x3be, 0x1f7
988 }, {
989 0x3c0, 0x00f
990 }, {
991 0x3c1, 0x20f
992 }, {
993 0x3c2, 0x10f
994 }, {
995 0x3c3, 0x30f
996 }, {
997 0x3c4, 0x08f
998 }, {
999 0x3c5, 0x28f
1000 }, {
1001 0x3c6, 0x18f
1002 }, {
1003 0x3c7, 0x38f
1004 }, {
1005 0x3c8, 0x04f
1006 }, {
1007 0x3c9, 0x24f
1008 }, {
1009 0x3ca, 0x14f
1010 }, {
1011 0x3cb, 0x34f
1012 }, {
1013 0x3cc, 0x0cf
1014 }, {
1015 0x3cd, 0x2cf
1016 }, {
1017 0x3ce, 0x1cf
1018 }, {
1019 0x3d0, 0x02f
1020 }, {
1021 0x3d1, 0x22f
1022 }, {
1023 0x3d2, 0x12f
1024 }, {
1025 0x3d3, 0x32f
1026 }, {
1027 0x3d4, 0x0af
1028 }, {
1029 0x3d5, 0x2af
1030 }, {
1031 0x3d6, 0x1af
1032 }, {
1033 0x3d7, 0x3af
1034 }, {
1035 0x3d8, 0x06f
1036 }, {
1037 0x3d9, 0x26f
1038 }, {
1039 0x3da, 0x16f
1040 }, {
1041 0x3db, 0x36f
1042 }, {
1043 0x3dc, 0x0ef
1044 }, {
1045 0x3dd, 0x2ef
1046 }, {
1047 0x3de, 0x1ef
1048 }, {
1049 0x3e0, 0x01f
1050 }, {
1051 0x3e1, 0x21f
1052 }, {
1053 0x3e2, 0x11f
1054 }, {
1055 0x3e3, 0x31f
1056 }, {
1057 0x3e4, 0x09f
1058 }, {
1059 0x3e5, 0x29f
1060 }, {
1061 0x3e6, 0x19f
1062 }, {
1063 0x3e7, 0x39f
1064 }, {
1065 0x3e8, 0x05f
1066 }, {
1067 0x3e9, 0x25f
1068 }, {
1069 0x3ea, 0x15f
1070 }, {
1071 0x3eb, 0x35f
1072 }, {
1073 0x3ec, 0x0df
1074 }, {
1075 0x3ed, 0x2df
1076 }, {
1077 0x3ee, 0x1df
1078 }, {
1079 0x3ef, 0x3df
1080 }, {
1081 0x3f0, 0x03f
1082 }, {
1083 0x3f1, 0x23f
1084 }, {
1085 0x3f2, 0x13f
1086 }, {
1087 0x3f3, 0x33f
1088 }, {
1089 0x3f4, 0x0bf
1090 }, {
1091 0x3f5, 0x2bf
1092 }, {
1093 0x3f6, 0x1bf
1094 }, {
1095 0x3f7, 0x3bf
1096 }, {
1097 0x3f8, 0x07f
1098 }, {
1099 0x3f9, 0x27f
1100 }, {
1101 0x3fa, 0x17f
1102 }, {
1103 0x3fb, 0x37f
1104 }, {
1105 0x3fc, 0x0ff
1106 }, {
1107 0x3fd, 0x2ff
1108 }, {
1109 0x3fe, 0x1ff
1110 }
1111 };
1112
1113
1114 {
1115 int i;
1116 for (i = 0; i < sizeof k1k2tab / sizeof k1k2tab[0]; ++i) {
1117 k1 = k1k2tab[i].k1;
1118 k2 = k1k2tab[i].k2;
1119 a = fz[k1];
1120 fz[k1] = fz[k2];
1121 fz[k2] = a;
1122 }
1123 }
1124
1125 for (fi = fz, fn = fz + 1024; fi < fn; fi += 4) {
1126 FLOAT f0, f1, f2, f3;
1127 f1 = fi[0] - fi[1];
1128 f0 = fi[0] + fi[1];
1129 f3 = fi[2] - fi[3];
1130 f2 = fi[2] + fi[3];
1131 fi[2] = (f0 - f2);
1132 fi[0] = (f0 + f2);
1133 fi[3] = (f1 - f3);
1134 fi[1] = (f1 + f3);
1135 }
1136
1137 k = 0;
1138 do {
1139 FLOAT s1, c1;
1140 k += 2;
1141 k1 = 1 << k;
1142 k2 = k1 << 1;
1143 k4 = k2 << 1;
1144 k3 = k2 + k1;
1145 kx = k1 >> 1;
1146 fi = fz;
1147 gi = fi + kx;
1148 fn = fz + 1024;
1149 do {
1150 FLOAT g0, f0, f1, g1, f2, g2, f3, g3;
1151 f1 = fi[0] - fi[k1];
1152 f0 = fi[0] + fi[k1];
1153 f3 = fi[k2] - fi[k3];
1154 f2 = fi[k2] + fi[k3];
1155 fi[k2] = f0 - f2;
1156 fi[0] = f0 + f2;
1157 fi[k3] = f1 - f3;
1158 fi[k1] = f1 + f3;
1159 g1 = gi[0] - gi[k1];
1160 g0 = gi[0] + gi[k1];
1161 g3 = SQRT2 * gi[k3];
1162 g2 = SQRT2 * gi[k2];
1163 gi[k2] = g0 - g2;
1164 gi[0] = g0 + g2;
1165 gi[k3] = g1 - g3;
1166 gi[k1] = g1 + g3;
1167 gi += k4;
1168 fi += k4;
1169 }
1170 while (fi < fn);
1171
1172 t_c = costab[k];
1173 t_s = sintab[k];
1174 c1 = 1;
1175 s1 = 0;
1176 for (i = 1; i < kx; i++) {
1177 FLOAT c2, s2;
1178 FLOAT t = c1;
1179 c1 = t * t_c - s1 * t_s;
1180 s1 = t * t_s + s1 * t_c;
1181 c2 = c1 * c1 - s1 * s1;
1182 s2 = 2 * (c1 * s1);
1183 fn = fz + 1024;
1184 fi = fz + i;
1185 gi = fz + k1 - i;
1186 do {
1187 FLOAT a, b, g0, f0, f1, g1, f2, g2, f3, g3;
1188 b = s2 * fi[k1] - c2 * gi[k1];
1189 a = c2 * fi[k1] + s2 * gi[k1];
1190 f1 = fi[0] - a;
1191 f0 = fi[0] + a;
1192 g1 = gi[0] - b;
1193 g0 = gi[0] + b;
1194 b = s2 * fi[k3] - c2 * gi[k3];
1195 a = c2 * fi[k3] + s2 * gi[k3];
1196 f3 = fi[k2] - a;
1197 f2 = fi[k2] + a;
1198 g3 = gi[k2] - b;
1199 g2 = gi[k2] + b;
1200 b = s1 * f2 - c1 * g3;
1201 a = c1 * f2 + s1 * g3;
1202 fi[k2] = f0 - a;
1203 fi[0] = f0 + a;
1204 gi[k3] = g1 - b;
1205 gi[k1] = g1 + b;
1206 b = c1 * g2 - s1 * f3;
1207 a = s1 * g2 + c1 * f3;
1208 gi[k2] = g0 - a;
1209 gi[0] = g0 + a;
1210 fi[k3] = f1 - b;
1211 fi[k1] = f1 + b;
1212 gi += k4;
1213 fi += k4;
1214 }
1215 while (fi < fn);
1216 }
1217 }
1218 while (k4 < 1024);
1219 }
1220
1221 #ifdef NEWATAN
1222 #define ATANSIZE 6000
1223 #define ATANSCALE 100.0
1224 /* Create a table of ATAN2 values.
1225 Valid for ratios of (y/x) from 0 to ATANSIZE/ATANSCALE (60)
1226 Largest absolute error in angle: 0.0167 radians i.e. ATANSCALE/ATANSIZE
1227 Depending on how you want to trade off speed/accuracy and mem usage, twiddle the defines
1228 MFC March 2003 */
1229 static FLOAT atan_t[ATANSIZE];
1230
atan_table(FLOAT y,FLOAT x)1231 static inline FLOAT atan_table(FLOAT y, FLOAT x)
1232 {
1233 int index;
1234
1235 index = (int) (ATANSCALE * fabs(y / x));
1236 if (index >= ATANSIZE)
1237 index = ATANSIZE - 1;
1238
1239 /* Have to work out the correct quadrant as well */
1240 if (y > 0 && x < 0)
1241 return (PI - atan_t[index]);
1242
1243 if (y < 0 && x > 0)
1244 return (-atan_t[index]);
1245
1246 if (y < 0 && x < 0)
1247 return (atan_t[index] - PI);
1248
1249 return (atan_t[index]);
1250 }
1251
atan_table_init(void)1252 static void atan_table_init(void)
1253 {
1254 int i;
1255 for (i = 0; i < ATANSIZE; i++)
1256 atan_t[i] = atan((FLOAT) i / ATANSCALE);
1257 }
1258
1259 #endif // NEWATAN
1260
1261 /* For variations on psycho model 2:
1262 N always equals 1024
1263 BUT in the returned values, no energy/phi is used at or above an index of 513 */
twolame_psycho_2_fft(FLOAT * x_real,FLOAT * energy,FLOAT * phi)1264 void twolame_psycho_2_fft(FLOAT * x_real, FLOAT * energy, FLOAT * phi)
1265 /* got rid of size "N" argument as it is always 1024 for layerII */
1266 {
1267 FLOAT imag, real;
1268 int i, j;
1269 #ifdef NEWATAN
1270 static int init = 0;
1271
1272 if (!init) {
1273 atan_table_init();
1274 init=1;
1275 }
1276 #endif
1277
1278
1279 fht(x_real);
1280
1281
1282 energy[0] = x_real[0] * x_real[0];
1283
1284 for (i = 1, j = 1023; i < 512; i++, j--) {
1285 imag = x_real[i];
1286 real = x_real[j];
1287 /* MFC FIXME Mar03 Why is this divided by 2.0? if a and b are the real and imaginary
1288 components then r = sqrt(a^2 + b^2), but, back in the psycho2 model, they calculate
1289 r=sqrt(energy), which, if you look at the original equation below is different */
1290 energy[i] = (imag * imag + real * real) / 2.0;
1291 if (energy[i] < 0.0005) {
1292 energy[i] = 0.0005;
1293 phi[i] = 0;
1294 } else
1295 #ifdef NEWATAN
1296 {
1297 phi[i] = atan_table(-imag, real) + PI / 4;
1298 }
1299 #else
1300 {
1301 phi[i] = atan2(-(FLOAT) imag, (FLOAT) real) + PI / 4;
1302 }
1303 #endif
1304 }
1305 energy[512] = x_real[512] * x_real[512];
1306 phi[512] = atan2(0.0, (FLOAT) x_real[512]);
1307 }
1308
1309
twolame_psycho_1_fft(FLOAT * x_real,FLOAT * energy,int N)1310 void twolame_psycho_1_fft(FLOAT * x_real, FLOAT * energy, int N)
1311 {
1312 FLOAT a, b;
1313 int i, j;
1314
1315 fht(x_real);
1316
1317 energy[0] = x_real[0] * x_real[0];
1318
1319 for (i = 1, j = N - 1; i < N / 2; i++, j--) {
1320 a = x_real[i];
1321 b = x_real[j];
1322 energy[i] = (a * a + b * b) / 2.0;
1323 }
1324 energy[N / 2] = x_real[N / 2] * x_real[N / 2];
1325 }
1326
1327
1328 // vim:ts=4:sw=4:nowrap:
1329