Cinelerra: hvirtual/toolame-02l/fft.c Source File

00001 /*
00002 ** FFT and FHT routines
00003 **  Copyright 1988, 1993; Ron Mayer
00004 **  
00005 **  fht(fz,n);
00006 **      Does a hartley transform of "n" points in the array "fz".
00007 **      
00008 ** NOTE: This routine uses at least 2 patented algorithms, and may be
00009 **       under the restrictions of a bunch of different organizations.
00010 **       Although I wrote it completely myself; it is kind of a derivative
00011 **       of a routine I once authored and released under the GPL, so it
00012 **       may fall under the free software foundation's restrictions;
00013 **       it was worked on as a Stanford Univ project, so they claim
00014 **       some rights to it; it was further optimized at work here, so
00015 **       I think this company claims parts of it.  The patents are
00016 **       held by R. Bracewell (the FHT algorithm) and O. Buneman (the
00017 **       trig generator), both at Stanford Univ.
00018 **       If it were up to me, I'd say go do whatever you want with it;
00019 **       but it would be polite to give credit to the following people
00020 **       if you use this anywhere:
00021 **           Euler     - probable inventor of the fourier transform.
00022 **           Gauss     - probable inventor of the FFT.
00023 **           Hartley   - probable inventor of the hartley transform.
00024 **           Buneman   - for a really cool trig generator
00025 **           Mayer(me) - for authoring this particular version and
00026 **                       including all the optimizations in one package.
00027 **       Thanks,
00028 **       Ron Mayer; mayer@acuson.com
00029 **
00030 */
00031 #include <stdio.h>
00032 #include <math.h>
00033 #include "common.h"
00034 #include "fft.h"
00035 #define         SQRT2                   1.4142135623730951454746218587388284504414
00036 
00037 
00038 static FLOAT costab[20] = {
00039   .00000000000000000000000000000000000000000000000000,
00040   .70710678118654752440084436210484903928483593768847,
00041   .92387953251128675612818318939678828682241662586364,
00042   .98078528040323044912618223613423903697393373089333,
00043   .99518472667219688624483695310947992157547486872985,
00044   .99879545620517239271477160475910069444320361470461,
00045   .99969881869620422011576564966617219685006108125772,
00046   .99992470183914454092164649119638322435060646880221,
00047   .99998117528260114265699043772856771617391725094433,
00048   .99999529380957617151158012570011989955298763362218,
00049   .99999882345170190992902571017152601904826792288976,
00050   .99999970586288221916022821773876567711626389934930,
00051   .99999992646571785114473148070738785694820115568892,
00052   .99999998161642929380834691540290971450507605124278,
00053   .99999999540410731289097193313960614895889430318945,
00054   .99999999885102682756267330779455410840053741619428
00055 };
00056 static FLOAT sintab[20] = {
00057   1.0000000000000000000000000000000000000000000000000,
00058   .70710678118654752440084436210484903928483593768846,
00059   .38268343236508977172845998403039886676134456248561,
00060   .19509032201612826784828486847702224092769161775195,
00061   .09801714032956060199419556388864184586113667316749,
00062   .04906767432741801425495497694268265831474536302574,
00063   .02454122852291228803173452945928292506546611923944,
00064   .01227153828571992607940826195100321214037231959176,
00065   .00613588464915447535964023459037258091705788631738,
00066   .00306795676296597627014536549091984251894461021344,
00067   .00153398018628476561230369715026407907995486457522,
00068   .00076699031874270452693856835794857664314091945205,
00069   .00038349518757139558907246168118138126339502603495,
00070   .00019174759731070330743990956198900093346887403385,
00071   .00009587379909597734587051721097647635118706561284,
00072   .00004793689960306688454900399049465887274686668768
00073 };
00074 
00075 /* This is a simplified version for n an even power of 2 */
00076 /* MFC: In the case of LayerII encoding, n==1024 always. */
00077 
00078 static void fht (FLOAT * fz)
00079 {
00080   int i, k, k1, k2, k3, k4, kx;
00081   FLOAT *fi, *fn, *gi;
00082   FLOAT t_c, t_s;
00083 
00084   FLOAT a;
00085   static const struct {
00086     unsigned short k1, k2;
00087   } k1k2tab[8 * 62] = {
00088     {
00089     0x020, 0x010}
00090     , {
00091     0x040, 0x008}
00092     , {
00093     0x050, 0x028}
00094     , {
00095     0x060, 0x018}
00096     , {
00097     0x068, 0x058}
00098     , {
00099     0x070, 0x038}
00100     , {
00101     0x080, 0x004}
00102     , {
00103     0x088, 0x044}
00104     , {
00105     0x090, 0x024}
00106     , {
00107     0x098, 0x064}
00108     , {
00109     0x0a0, 0x014}
00110     , {
00111     0x0a4, 0x094}
00112     , {
00113     0x0a8, 0x054}
00114     , {
00115     0x0b0, 0x034}
00116     , {
00117     0x0b8, 0x074}
00118     , {
00119     0x0c0, 0x00c}
00120     , {
00121     0x0c4, 0x08c}
00122     , {
00123     0x0c8, 0x04c}
00124     , {
00125     0x0d0, 0x02c}
00126     , {
00127     0x0d4, 0x0ac}
00128     , {
00129     0x0d8, 0x06c}
00130     , {
00131     0x0e0, 0x01c}
00132     , {
00133     0x0e4, 0x09c}
00134     , {
00135     0x0e8, 0x05c}
00136     , {
00137     0x0ec, 0x0dc}
00138     , {
00139     0x0f0, 0x03c}
00140     , {
00141     0x0f4, 0x0bc}
00142     , {
00143     0x0f8, 0x07c}
00144     , {
00145     0x100, 0x002}
00146     , {
00147     0x104, 0x082}
00148     , {
00149     0x108, 0x042}
00150     , {
00151     0x10c, 0x0c2}
00152     , {
00153     0x110, 0x022}
00154     , {
00155     0x114, 0x0a2}
00156     , {
00157     0x118, 0x062}
00158     , {
00159     0x11c, 0x0e2}
00160     , {
00161     0x120, 0x012}
00162     , {
00163     0x122, 0x112}
00164     , {
00165     0x124, 0x092}
00166     , {
00167     0x128, 0x052}
00168     , {
00169     0x12c, 0x0d2}
00170     , {
00171     0x130, 0x032}
00172     , {
00173     0x134, 0x0b2}
00174     , {
00175     0x138, 0x072}
00176     , {
00177     0x13c, 0x0f2}
00178     , {
00179     0x140, 0x00a}
00180     , {
00181     0x142, 0x10a}
00182     , {
00183     0x144, 0x08a}
00184     , {
00185     0x148, 0x04a}
00186     , {
00187     0x14c, 0x0ca}
00188     , {
00189     0x150, 0x02a}
00190     , {
00191     0x152, 0x12a}
00192     , {
00193     0x154, 0x0aa}
00194     , {
00195     0x158, 0x06a}
00196     , {
00197     0x15c, 0x0ea}
00198     , {
00199     0x160, 0x01a}
00200     , {
00201     0x162, 0x11a}
00202     , {
00203     0x164, 0x09a}
00204     , {
00205     0x168, 0x05a}
00206     , {
00207     0x16a, 0x15a}
00208     , {
00209     0x16c, 0x0da}
00210     , {
00211     0x170, 0x03a}
00212     , {
00213     0x172, 0x13a}
00214     , {
00215     0x174, 0x0ba}
00216     , {
00217     0x178, 0x07a}
00218     , {
00219     0x17c, 0x0fa}
00220     , {
00221     0x180, 0x006}
00222     , {
00223     0x182, 0x106}
00224     , {
00225     0x184, 0x086}
00226     , {
00227     0x188, 0x046}
00228     , {
00229     0x18a, 0x146}
00230     , {
00231     0x18c, 0x0c6}
00232     , {
00233     0x190, 0x026}
00234     , {
00235     0x192, 0x126}
00236     , {
00237     0x194, 0x0a6}
00238     , {
00239     0x198, 0x066}
00240     , {
00241     0x19a, 0x166}
00242     , {
00243     0x19c, 0x0e6}
00244     , {
00245     0x1a0, 0x016}
00246     , {
00247     0x1a2, 0x116}
00248     , {
00249     0x1a4, 0x096}
00250     , {
00251     0x1a6, 0x196}
00252     , {
00253     0x1a8, 0x056}
00254     , {
00255     0x1aa, 0x156}
00256     , {
00257     0x1ac, 0x0d6}
00258     , {
00259     0x1b0, 0x036}
00260     , {
00261     0x1b2, 0x136}
00262     , {
00263     0x1b4, 0x0b6}
00264     , {
00265     0x1b8, 0x076}
00266     , {
00267     0x1ba, 0x176}
00268     , {
00269     0x1bc, 0x0f6}
00270     , {
00271     0x1c0, 0x00e}
00272     , {
00273     0x1c2, 0x10e}
00274     , {
00275     0x1c4, 0x08e}
00276     , {
00277     0x1c6, 0x18e}
00278     , {
00279     0x1c8, 0x04e}
00280     , {
00281     0x1ca, 0x14e}
00282     , {
00283     0x1cc, 0x0ce}
00284     , {
00285     0x1d0, 0x02e}
00286     , {
00287     0x1d2, 0x12e}
00288     , {
00289     0x1d4, 0x0ae}
00290     , {
00291     0x1d6, 0x1ae}
00292     , {
00293     0x1d8, 0x06e}
00294     , {
00295     0x1da, 0x16e}
00296     , {
00297     0x1dc, 0x0ee}
00298     , {
00299     0x1e0, 0x01e}
00300     , {
00301     0x1e2, 0x11e}
00302     , {
00303     0x1e4, 0x09e}
00304     , {
00305     0x1e6, 0x19e}
00306     , {
00307     0x1e8, 0x05e}
00308     , {
00309     0x1ea, 0x15e}
00310     , {
00311     0x1ec, 0x0de}
00312     , {
00313     0x1ee, 0x1de}
00314     , {
00315     0x1f0, 0x03e}
00316     , {
00317     0x1f2, 0x13e}
00318     , {
00319     0x1f4, 0x0be}
00320     , {
00321     0x1f6, 0x1be}
00322     , {
00323     0x1f8, 0x07e}
00324     , {
00325     0x1fa, 0x17e}
00326     , {
00327     0x1fc, 0x0fe}
00328     , {
00329     0x200, 0x001}
00330     , {
00331     0x202, 0x101}
00332     , {
00333     0x204, 0x081}
00334     , {
00335     0x206, 0x181}
00336     , {
00337     0x208, 0x041}
00338     , {
00339     0x20a, 0x141}
00340     , {
00341     0x20c, 0x0c1}
00342     , {
00343     0x20e, 0x1c1}
00344     , {
00345     0x210, 0x021}
00346     , {
00347     0x212, 0x121}
00348     , {
00349     0x214, 0x0a1}
00350     , {
00351     0x216, 0x1a1}
00352     , {
00353     0x218, 0x061}
00354     , {
00355     0x21a, 0x161}
00356     , {
00357     0x21c, 0x0e1}
00358     , {
00359     0x21e, 0x1e1}
00360     , {
00361     0x220, 0x011}
00362     , {
00363     0x221, 0x211}
00364     , {
00365     0x222, 0x111}
00366     , {
00367     0x224, 0x091}
00368     , {
00369     0x226, 0x191}
00370     , {
00371     0x228, 0x051}
00372     , {
00373     0x22a, 0x151}
00374     , {
00375     0x22c, 0x0d1}
00376     , {
00377     0x22e, 0x1d1}
00378     , {
00379     0x230, 0x031}
00380     , {
00381     0x232, 0x131}
00382     , {
00383     0x234, 0x0b1}
00384     , {
00385     0x236, 0x1b1}
00386     , {
00387     0x238, 0x071}
00388     , {
00389     0x23a, 0x171}
00390     , {
00391     0x23c, 0x0f1}
00392     , {
00393     0x23e, 0x1f1}
00394     , {
00395     0x240, 0x009}
00396     , {
00397     0x241, 0x209}
00398     , {
00399     0x242, 0x109}
00400     , {
00401     0x244, 0x089}
00402     , {
00403     0x246, 0x189}
00404     , {
00405     0x248, 0x049}
00406     , {
00407     0x24a, 0x149}
00408     , {
00409     0x24c, 0x0c9}
00410     , {
00411     0x24e, 0x1c9}
00412     , {
00413     0x250, 0x029}
00414     , {
00415     0x251, 0x229}
00416     , {
00417     0x252, 0x129}
00418     , {
00419     0x254, 0x0a9}
00420     , {
00421     0x256, 0x1a9}
00422     , {
00423     0x258, 0x069}
00424     , {
00425     0x25a, 0x169}
00426     , {
00427     0x25c, 0x0e9}
00428     , {
00429     0x25e, 0x1e9}
00430     , {
00431     0x260, 0x019}
00432     , {
00433     0x261, 0x219}
00434     , {
00435     0x262, 0x119}
00436     , {
00437     0x264, 0x099}
00438     , {
00439     0x266, 0x199}
00440     , {
00441     0x268, 0x059}
00442     , {
00443     0x269, 0x259}
00444     , {
00445     0x26a, 0x159}
00446     , {
00447     0x26c, 0x0d9}
00448     , {
00449     0x26e, 0x1d9}
00450     , {
00451     0x270, 0x039}
00452     , {
00453     0x271, 0x239}
00454     , {
00455     0x272, 0x139}
00456     , {
00457     0x274, 0x0b9}
00458     , {
00459     0x276, 0x1b9}
00460     , {
00461     0x278, 0x079}
00462     , {
00463     0x27a, 0x179}
00464     , {
00465     0x27c, 0x0f9}
00466     , {
00467     0x27e, 0x1f9}
00468     , {
00469     0x280, 0x005}
00470     , {
00471     0x281, 0x205}
00472     , {
00473     0x282, 0x105}
00474     , {
00475     0x284, 0x085}
00476     , {
00477     0x286, 0x185}
00478     , {
00479     0x288, 0x045}
00480     , {
00481     0x289, 0x245}
00482     , {
00483     0x28a, 0x145}
00484     , {
00485     0x28c, 0x0c5}
00486     , {
00487     0x28e, 0x1c5}
00488     , {
00489     0x290, 0x025}
00490     , {
00491     0x291, 0x225}
00492     , {
00493     0x292, 0x125}
00494     , {
00495     0x294, 0x0a5}
00496     , {
00497     0x296, 0x1a5}
00498     , {
00499     0x298, 0x065}
00500     , {
00501     0x299, 0x265}
00502     , {
00503     0x29a, 0x165}
00504     , {
00505     0x29c, 0x0e5}
00506     , {
00507     0x29e, 0x1e5}
00508     , {
00509     0x2a0, 0x015}
00510     , {
00511     0x2a1, 0x215}
00512     , {
00513     0x2a2, 0x115}
00514     , {
00515     0x2a4, 0x095}
00516     , {
00517     0x2a5, 0x295}
00518     , {
00519     0x2a6, 0x195}
00520     , {
00521     0x2a8, 0x055}
00522     , {
00523     0x2a9, 0x255}
00524     , {
00525     0x2aa, 0x155}
00526     , {
00527     0x2ac, 0x0d5}
00528     , {
00529     0x2ae, 0x1d5}
00530     , {
00531     0x2b0, 0x035}
00532     , {
00533     0x2b1, 0x235}
00534     , {
00535     0x2b2, 0x135}
00536     , {
00537     0x2b4, 0x0b5}
00538     , {
00539     0x2b6, 0x1b5}
00540     , {
00541     0x2b8, 0x075}
00542     , {
00543     0x2b9, 0x275}
00544     , {
00545     0x2ba, 0x175}
00546     , {
00547     0x2bc, 0x0f5}
00548     , {
00549     0x2be, 0x1f5}
00550     , {
00551     0x2c0, 0x00d}
00552     , {
00553     0x2c1, 0x20d}
00554     , {
00555     0x2c2, 0x10d}
00556     , {
00557     0x2c4, 0x08d}
00558     , {
00559     0x2c5, 0x28d}
00560     , {
00561     0x2c6, 0x18d}
00562     , {
00563     0x2c8, 0x04d}
00564     , {
00565     0x2c9, 0x24d}
00566     , {
00567     0x2ca, 0x14d}
00568     , {
00569     0x2cc, 0x0cd}
00570     , {
00571     0x2ce, 0x1cd}
00572     , {
00573     0x2d0, 0x02d}
00574     , {
00575     0x2d1, 0x22d}
00576     , {
00577     0x2d2, 0x12d}
00578     , {
00579     0x2d4, 0x0ad}
00580     , {
00581     0x2d5, 0x2ad}
00582     , {
00583     0x2d6, 0x1ad}
00584     , {
00585     0x2d8, 0x06d}
00586     , {
00587     0x2d9, 0x26d}
00588     , {
00589     0x2da, 0x16d}
00590     , {
00591     0x2dc, 0x0ed}
00592     , {
00593     0x2de, 0x1ed}
00594     , {
00595     0x2e0, 0x01d}
00596     , {
00597     0x2e1, 0x21d}
00598     , {
00599     0x2e2, 0x11d}
00600     , {
00601     0x2e4, 0x09d}
00602     , {
00603     0x2e5, 0x29d}
00604     , {
00605     0x2e6, 0x19d}
00606     , {
00607     0x2e8, 0x05d}
00608     , {
00609     0x2e9, 0x25d}
00610     , {
00611     0x2ea, 0x15d}
00612     , {
00613     0x2ec, 0x0dd}
00614     , {
00615     0x2ed, 0x2dd}
00616     , {
00617     0x2ee, 0x1dd}
00618     , {
00619     0x2f0, 0x03d}
00620     , {
00621     0x2f1, 0x23d}
00622     , {
00623     0x2f2, 0x13d}
00624     , {
00625     0x2f4, 0x0bd}
00626     , {
00627     0x2f5, 0x2bd}
00628     , {
00629     0x2f6, 0x1bd}
00630     , {
00631     0x2f8, 0x07d}
00632     , {
00633     0x2f9, 0x27d}
00634     , {
00635     0x2fa, 0x17d}
00636     , {
00637     0x2fc, 0x0fd}
00638     , {
00639     0x2fe, 0x1fd}
00640     , {
00641     0x300, 0x003}
00642     , {
00643     0x301, 0x203}
00644     , {
00645     0x302, 0x103}
00646     , {
00647     0x304, 0x083}
00648     , {
00649     0x305, 0x283}
00650     , {
00651     0x306, 0x183}
00652     , {
00653     0x308, 0x043}
00654     , {
00655     0x309, 0x243}
00656     , {
00657     0x30a, 0x143}
00658     , {
00659     0x30c, 0x0c3}
00660     , {
00661     0x30d, 0x2c3}
00662     , {
00663     0x30e, 0x1c3}
00664     , {
00665     0x310, 0x023}
00666     , {
00667     0x311, 0x223}
00668     , {
00669     0x312, 0x123}
00670     , {
00671     0x314, 0x0a3}
00672     , {
00673     0x315, 0x2a3}
00674     , {
00675     0x316, 0x1a3}
00676     , {
00677     0x318, 0x063}
00678     , {
00679     0x319, 0x263}
00680     , {
00681     0x31a, 0x163}
00682     , {
00683     0x31c, 0x0e3}
00684     , {
00685     0x31d, 0x2e3}
00686     , {
00687     0x31e, 0x1e3}
00688     , {
00689     0x320, 0x013}
00690     , {
00691     0x321, 0x213}
00692     , {
00693     0x322, 0x113}
00694     , {
00695     0x323, 0x313}
00696     , {
00697     0x324, 0x093}
00698     , {
00699     0x325, 0x293}
00700     , {
00701     0x326, 0x193}
00702     , {
00703     0x328, 0x053}
00704     , {
00705     0x329, 0x253}
00706     , {
00707     0x32a, 0x153}
00708     , {
00709     0x32c, 0x0d3}
00710     , {
00711     0x32d, 0x2d3}
00712     , {
00713     0x32e, 0x1d3}
00714     , {
00715     0x330, 0x033}
00716     , {
00717     0x331, 0x233}
00718     , {
00719     0x332, 0x133}
00720     , {
00721     0x334, 0x0b3}
00722     , {
00723     0x335, 0x2b3}
00724     , {
00725     0x336, 0x1b3}
00726     , {
00727     0x338, 0x073}
00728     , {
00729     0x339, 0x273}
00730     , {
00731     0x33a, 0x173}
00732     , {
00733     0x33c, 0x0f3}
00734     , {
00735     0x33d, 0x2f3}
00736     , {
00737     0x33e, 0x1f3}
00738     , {
00739     0x340, 0x00b}
00740     , {
00741     0x341, 0x20b}
00742     , {
00743     0x342, 0x10b}
00744     , {
00745     0x343, 0x30b}
00746     , {
00747     0x344, 0x08b}
00748     , {
00749     0x345, 0x28b}
00750     , {
00751     0x346, 0x18b}
00752     , {
00753     0x348, 0x04b}
00754     , {
00755     0x349, 0x24b}
00756     , {
00757     0x34a, 0x14b}
00758     , {
00759     0x34c, 0x0cb}
00760     , {
00761     0x34d, 0x2cb}
00762     , {
00763     0x34e, 0x1cb}
00764     , {
00765     0x350, 0x02b}
00766     , {
00767     0x351, 0x22b}
00768     , {
00769     0x352, 0x12b}
00770     , {
00771     0x353, 0x32b}
00772     , {
00773     0x354, 0x0ab}
00774     , {
00775     0x355, 0x2ab}
00776     , {
00777     0x356, 0x1ab}
00778     , {
00779     0x358, 0x06b}
00780     , {
00781     0x359, 0x26b}
00782     , {
00783     0x35a, 0x16b}
00784     , {
00785     0x35c, 0x0eb}
00786     , {
00787     0x35d, 0x2eb}
00788     , {
00789     0x35e, 0x1eb}
00790     , {
00791     0x360, 0x01b}
00792     , {
00793     0x361, 0x21b}
00794     , {
00795     0x362, 0x11b}
00796     , {
00797     0x363, 0x31b}
00798     , {
00799     0x364, 0x09b}
00800     , {
00801     0x365, 0x29b}
00802     , {
00803     0x366, 0x19b}
00804     , {
00805     0x368, 0x05b}
00806     , {
00807     0x369, 0x25b}
00808     , {
00809     0x36a, 0x15b}
00810     , {
00811     0x36b, 0x35b}
00812     , {
00813     0x36c, 0x0db}
00814     , {
00815     0x36d, 0x2db}
00816     , {
00817     0x36e, 0x1db}
00818     , {
00819     0x370, 0x03b}
00820     , {
00821     0x371, 0x23b}
00822     , {
00823     0x372, 0x13b}
00824     , {
00825     0x373, 0x33b}
00826     , {
00827     0x374, 0x0bb}
00828     , {
00829     0x375, 0x2bb}
00830     , {
00831     0x376, 0x1bb}
00832     , {
00833     0x378, 0x07b}
00834     , {
00835     0x379, 0x27b}
00836     , {
00837     0x37a, 0x17b}
00838     , {
00839     0x37c, 0x0fb}
00840     , {
00841     0x37d, 0x2fb}
00842     , {
00843     0x37e, 0x1fb}
00844     , {
00845     0x380, 0x007}
00846     , {
00847     0x381, 0x207}
00848     , {
00849     0x382, 0x107}
00850     , {
00851     0x383, 0x307}
00852     , {
00853     0x384, 0x087}
00854     , {
00855     0x385, 0x287}
00856     , {
00857     0x386, 0x187}
00858     , {
00859     0x388, 0x047}
00860     , {
00861     0x389, 0x247}
00862     , {
00863     0x38a, 0x147}
00864     , {
00865     0x38b, 0x347}
00866     , {
00867     0x38c, 0x0c7}
00868     , {
00869     0x38d, 0x2c7}
00870     , {
00871     0x38e, 0x1c7}
00872     , {
00873     0x390, 0x027}
00874     , {
00875     0x391, 0x227}
00876     , {
00877     0x392, 0x127}
00878     , {
00879     0x393, 0x327}
00880     , {
00881     0x394, 0x0a7}
00882     , {
00883     0x395, 0x2a7}
00884     , {
00885     0x396, 0x1a7}
00886     , {
00887     0x398, 0x067}
00888     , {
00889     0x399, 0x267}
00890     , {
00891     0x39a, 0x167}
00892     , {
00893     0x39b, 0x367}
00894     , {
00895     0x39c, 0x0e7}
00896     , {
00897     0x39d, 0x2e7}
00898     , {
00899     0x39e, 0x1e7}
00900     , {
00901     0x3a0, 0x017}
00902     , {
00903     0x3a1, 0x217}
00904     , {
00905     0x3a2, 0x117}
00906     , {
00907     0x3a3, 0x317}
00908     , {
00909     0x3a4, 0x097}
00910     , {
00911     0x3a5, 0x297}
00912     , {
00913     0x3a6, 0x197}
00914     , {
00915     0x3a7, 0x397}
00916     , {
00917     0x3a8, 0x057}
00918     , {
00919     0x3a9, 0x257}
00920     , {
00921     0x3aa, 0x157}
00922     , {
00923     0x3ab, 0x357}
00924     , {
00925     0x3ac, 0x0d7}
00926     , {
00927     0x3ad, 0x2d7}
00928     , {
00929     0x3ae, 0x1d7}
00930     , {
00931     0x3b0, 0x037}
00932     , {
00933     0x3b1, 0x237}
00934     , {
00935     0x3b2, 0x137}
00936     , {
00937     0x3b3, 0x337}
00938     , {
00939     0x3b4, 0x0b7}
00940     , {
00941     0x3b5, 0x2b7}
00942     , {
00943     0x3b6, 0x1b7}
00944     , {
00945     0x3b8, 0x077}
00946     , {
00947     0x3b9, 0x277}
00948     , {
00949     0x3ba, 0x177}
00950     , {
00951     0x3bb, 0x377}
00952     , {
00953     0x3bc, 0x0f7}
00954     , {
00955     0x3bd, 0x2f7}
00956     , {
00957     0x3be, 0x1f7}
00958     , {
00959     0x3c0, 0x00f}
00960     , {
00961     0x3c1, 0x20f}
00962     , {
00963     0x3c2, 0x10f}
00964     , {
00965     0x3c3, 0x30f}
00966     , {
00967     0x3c4, 0x08f}
00968     , {
00969     0x3c5, 0x28f}
00970     , {
00971     0x3c6, 0x18f}
00972     , {
00973     0x3c7, 0x38f}
00974     , {
00975     0x3c8, 0x04f}
00976     , {
00977     0x3c9, 0x24f}
00978     , {
00979     0x3ca, 0x14f}
00980     , {
00981     0x3cb, 0x34f}
00982     , {
00983     0x3cc, 0x0cf}
00984     , {
00985     0x3cd, 0x2cf}
00986     , {
00987     0x3ce, 0x1cf}
00988     , {
00989     0x3d0, 0x02f}
00990     , {
00991     0x3d1, 0x22f}
00992     , {
00993     0x3d2, 0x12f}
00994     , {
00995     0x3d3, 0x32f}
00996     , {
00997     0x3d4, 0x0af}
00998     , {
00999     0x3d5, 0x2af}
01000     , {
01001     0x3d6, 0x1af}
01002     , {
01003     0x3d7, 0x3af}
01004     , {
01005     0x3d8, 0x06f}
01006     , {
01007     0x3d9, 0x26f}
01008     , {
01009     0x3da, 0x16f}
01010     , {
01011     0x3db, 0x36f}
01012     , {
01013     0x3dc, 0x0ef}
01014     , {
01015     0x3dd, 0x2ef}
01016     , {
01017     0x3de, 0x1ef}
01018     , {
01019     0x3e0, 0x01f}
01020     , {
01021     0x3e1, 0x21f}
01022     , {
01023     0x3e2, 0x11f}
01024     , {
01025     0x3e3, 0x31f}
01026     , {
01027     0x3e4, 0x09f}
01028     , {
01029     0x3e5, 0x29f}
01030     , {
01031     0x3e6, 0x19f}
01032     , {
01033     0x3e7, 0x39f}
01034     , {
01035     0x3e8, 0x05f}
01036     , {
01037     0x3e9, 0x25f}
01038     , {
01039     0x3ea, 0x15f}
01040     , {
01041     0x3eb, 0x35f}
01042     , {
01043     0x3ec, 0x0df}
01044     , {
01045     0x3ed, 0x2df}
01046     , {
01047     0x3ee, 0x1df}
01048     , {
01049     0x3ef, 0x3df}
01050     , {
01051     0x3f0, 0x03f}
01052     , {
01053     0x3f1, 0x23f}
01054     , {
01055     0x3f2, 0x13f}
01056     , {
01057     0x3f3, 0x33f}
01058     , {
01059     0x3f4, 0x0bf}
01060     , {
01061     0x3f5, 0x2bf}
01062     , {
01063     0x3f6, 0x1bf}
01064     , {
01065     0x3f7, 0x3bf}
01066     , {
01067     0x3f8, 0x07f}
01068     , {
01069     0x3f9, 0x27f}
01070     , {
01071     0x3fa, 0x17f}
01072     , {
01073     0x3fb, 0x37f}
01074     , {
01075     0x3fc, 0x0ff}
01076     , {
01077     0x3fd, 0x2ff}
01078     , {
01079     0x3fe, 0x1ff}
01080   };
01081   {
01082     int i;
01083     for (i = 0; i < sizeof k1k2tab / sizeof k1k2tab[0]; ++i) {
01084       k1 = k1k2tab[i].k1;
01085       k2 = k1k2tab[i].k2;
01086       a = fz[k1];
01087       fz[k1] = fz[k2];
01088       fz[k2] = a;
01089     }
01090   }
01091 
01092   for (fi = fz, fn = fz + 1024; fi < fn; fi += 4) {
01093     FLOAT f0, f1, f2, f3;
01094     f1 = fi[0] - fi[1];
01095     f0 = fi[0] + fi[1];
01096     f3 = fi[2] - fi[3];
01097     f2 = fi[2] + fi[3];
01098     fi[2] = (f0 - f2);
01099     fi[0] = (f0 + f2);
01100     fi[3] = (f1 - f3);
01101     fi[1] = (f1 + f3);
01102   }
01103 
01104   k = 0;
01105   do {
01106     FLOAT s1, c1;
01107     k += 2;
01108     k1 = 1 << k;
01109     k2 = k1 << 1;
01110     k4 = k2 << 1;
01111     k3 = k2 + k1;
01112     kx = k1 >> 1;
01113     fi = fz;
01114     gi = fi + kx;
01115     fn = fz + 1024;
01116     do {
01117       FLOAT g0, f0, f1, g1, f2, g2, f3, g3;
01118       f1 = fi[0] - fi[k1];
01119       f0 = fi[0] + fi[k1];
01120       f3 = fi[k2] - fi[k3];
01121       f2 = fi[k2] + fi[k3];
01122       fi[k2] = f0 - f2;
01123       fi[0] = f0 + f2;
01124       fi[k3] = f1 - f3;
01125       fi[k1] = f1 + f3;
01126       g1 = gi[0] - gi[k1];
01127       g0 = gi[0] + gi[k1];
01128       g3 = SQRT2 * gi[k3];
01129       g2 = SQRT2 * gi[k2];
01130       gi[k2] = g0 - g2;
01131       gi[0] = g0 + g2;
01132       gi[k3] = g1 - g3;
01133       gi[k1] = g1 + g3;
01134       gi += k4;
01135       fi += k4;
01136     }
01137     while (fi < fn);
01138     t_c = costab[k];
01139     t_s = sintab[k];
01140     c1 = 1;
01141     s1 = 0;
01142     for (i = 1; i < kx; i++) {
01143       FLOAT c2, s2;
01144       FLOAT t = c1;
01145       c1 = t * t_c - s1 * t_s;
01146       s1 = t * t_s + s1 * t_c;
01147       c2 = c1 * c1 - s1 * s1;
01148       s2 = 2 * (c1 * s1);
01149       fn = fz + 1024;
01150       fi = fz + i;
01151       gi = fz + k1 - i;
01152       do {
01153         FLOAT a, b, g0, f0, f1, g1, f2, g2, f3, g3;
01154         b = s2 * fi[k1] - c2 * gi[k1];
01155         a = c2 * fi[k1] + s2 * gi[k1];
01156         f1 = fi[0] - a;
01157         f0 = fi[0] + a;
01158         g1 = gi[0] - b;
01159         g0 = gi[0] + b;
01160         b = s2 * fi[k3] - c2 * gi[k3];
01161         a = c2 * fi[k3] + s2 * gi[k3];
01162         f3 = fi[k2] - a;
01163         f2 = fi[k2] + a;
01164         g3 = gi[k2] - b;
01165         g2 = gi[k2] + b;
01166         b = s1 * f2 - c1 * g3;
01167         a = c1 * f2 + s1 * g3;
01168         fi[k2] = f0 - a;
01169         fi[0] = f0 + a;
01170         gi[k3] = g1 - b;
01171         gi[k1] = g1 + b;
01172         b = c1 * g2 - s1 * f3;
01173         a = s1 * g2 + c1 * f3;
01174         gi[k2] = g0 - a;
01175         gi[0] = g0 + a;
01176         fi[k3] = f1 - b;
01177         fi[k1] = f1 + b;
01178         gi += k4;
01179         fi += k4;
01180       }
01181       while (fi < fn);
01182     }
01183   }
01184   while (k4 < 1024);
01185 }
01186 
01187 #ifdef NEWATAN
01188 #define ATANSIZE 2000
01189 #define ATANSCALE 50.0
01190   static FLOAT atan_t[ATANSIZE];
01191 
01192 INLINE FLOAT atan_table(FLOAT y, FLOAT x) {
01193   int index;
01194 
01195   index = (int)(ATANSCALE * fabs(y/x));
01196   if (index>=ATANSIZE)
01197     index = ATANSIZE-1;
01198 
01199   if (y>0 && x<0)
01200     return( PI - atan_t[index] );
01201 
01202   if (y<0 && x>0)
01203     return( -atan_t[index] );
01204 
01205   if (y<0 && x<0)
01206     return( atan_t[index] - PI );
01207 
01208   return(atan_t[index]);
01209 }
01210 
01211 void atan_table_init(void) {
01212   int i;
01213   for (i=0;i<ATANSIZE;i++)
01214     atan_t[i] = atan((double)i/ATANSCALE);
01215 }
01216 
01217 #endif //NEWATAN
01218 
01219 /* For variations on psycho model 2:
01220    N always equals 1024
01221    BUT in the returned values, no energy/phi is used at or above an index of 513 */
01222 void psycho_2_fft (FLOAT * x_real, FLOAT * energy, FLOAT * phi) 
01223      /* got rid of size "N" argument as it is always 1024 for layerII */
01224 {
01225   FLOAT a, b;
01226   int i, j;
01227 #ifdef NEWATAN
01228   static int init=0;
01229 
01230   if (!init) {
01231     atan_table_init();
01232     init++;
01233   }
01234 #endif
01235 
01236 
01237   fht (x_real);
01238 
01239 
01240   energy[0] = x_real[0] * x_real[0];
01241 
01242   for (i = 1, j = 1023; i < 512; i++, j--) {
01243     a = x_real[i];
01244     b = x_real[j];
01245     /* MFC FIXME Mar03 Why is this divided by 2.0?
01246        if a and b are the real and imaginary components then
01247        r = sqrt(a^2 + b^2),
01248        but, back in the psycho2 model, they calculate r=sqrt(energy), 
01249        which, if you look at the original equation below is different */
01250     energy[i] = (a * a + b * b) / 2.0;
01251     if (energy[i] < 0.0005) {
01252       energy[i] = 0.0005;
01253       phi[i] = 0;
01254     } else      
01255 #ifdef NEWATAN
01256       {         
01257         phi[i] = atan_table(-a, b) + PI/4;
01258       }
01259 #else
01260       {
01261       phi[i] = atan2(-(double)a, (double)b) + PI/4;
01262       }
01263 #endif
01264   }
01265   energy[512] = x_real[512] * x_real[512];
01266   phi[512] = atan2 (0.0, (double) x_real[512]);
01267 }
01268 
01269 
01270 void psycho_1_fft (FLOAT * x_real, FLOAT * energy, int N)
01271 {
01272   FLOAT a, b;
01273   int i, j;
01274 
01275   fht (x_real);
01276 
01277   energy[0] = x_real[0] * x_real[0];
01278 
01279   for (i = 1, j = N - 1; i < N / 2; i++, j--) {
01280     a = x_real[i];
01281     b = x_real[j];
01282     energy[i] = (a * a + b * b) / 2.0;
01283   }
01284   energy[N / 2] = x_real[N / 2] * x_real[N / 2];
01285 }
01286 
01287 
01288