toolame-02l/fft.c

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