home *** CD-ROM | disk | FTP | other *** search
/ AmigActive 4 / AACD04.ISO / AACD / Sound / LAME / Source / subs.c < prev    next >
C/C++ Source or Header  |  1999-06-02  |  7KB  |  249 lines
  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.  
  32. #include <math.h>
  33. #include "common.h"
  34.  
  35. static FLOAT costab[20]=
  36.     {
  37.      .00000000000000000000000000000000000000000000000000,
  38.      .70710678118654752440084436210484903928483593768847,
  39.      .92387953251128675612818318939678828682241662586364,
  40.      .98078528040323044912618223613423903697393373089333,
  41.      .99518472667219688624483695310947992157547486872985,
  42.      .99879545620517239271477160475910069444320361470461,
  43.      .99969881869620422011576564966617219685006108125772,
  44.      .99992470183914454092164649119638322435060646880221,
  45.      .99998117528260114265699043772856771617391725094433,
  46.      .99999529380957617151158012570011989955298763362218,
  47.      .99999882345170190992902571017152601904826792288976,
  48.      .99999970586288221916022821773876567711626389934930,
  49.      .99999992646571785114473148070738785694820115568892,
  50.      .99999998161642929380834691540290971450507605124278,
  51.      .99999999540410731289097193313960614895889430318945,
  52.      .99999999885102682756267330779455410840053741619428
  53.     };
  54. static FLOAT sintab[20]=
  55.     {
  56.      1.0000000000000000000000000000000000000000000000000,
  57.      .70710678118654752440084436210484903928483593768846,
  58.      .38268343236508977172845998403039886676134456248561,
  59.      .19509032201612826784828486847702224092769161775195,
  60.      .09801714032956060199419556388864184586113667316749,
  61.      .04906767432741801425495497694268265831474536302574,
  62.      .02454122852291228803173452945928292506546611923944,
  63.      .01227153828571992607940826195100321214037231959176,
  64.      .00613588464915447535964023459037258091705788631738,
  65.      .00306795676296597627014536549091984251894461021344,
  66.      .00153398018628476561230369715026407907995486457522,
  67.      .00076699031874270452693856835794857664314091945205,
  68.      .00038349518757139558907246168118138126339502603495,
  69.      .00019174759731070330743990956198900093346887403385,
  70.      .00009587379909597734587051721097647635118706561284,
  71.      .00004793689960306688454900399049465887274686668768
  72.     };
  73.  
  74. #define SQRT2   2*0.70710678118654752440084436210484
  75.  
  76. /* This is a simplified version for n an even power of 2 */
  77.  
  78. static void fht(FLOAT *fz, int n)
  79. {
  80.  int i,k,k1,k2,k3,k4,kx;
  81.  FLOAT *fi,*fn,*gi;
  82.  FLOAT t_c,t_s;
  83.  
  84.  for (k1=1,k2=0;k1<n;k1++)
  85.     {
  86.      FLOAT a;
  87.      for (k=n>>1; (!((k2^=k)&k)); k>>=1);
  88.      if (k1>k2)
  89.     {
  90.          a=fz[k1];fz[k1]=fz[k2];fz[k2]=a;
  91.     }
  92.     }
  93.   for (fi=fz,fn=fz+n;fi<fn;fi+=4)
  94.      {
  95.       FLOAT f0,f1,f2,f3;
  96.       f1     = fi[0 ]-fi[1 ];
  97.       f0     = fi[0 ]+fi[1 ];
  98.       f3     = fi[2 ]-fi[3 ];
  99.       f2     = fi[2 ]+fi[3 ];
  100.       fi[2 ] = (f0-f2);    
  101.       fi[0 ] = (f0+f2);
  102.       fi[3 ] = (f1-f3);    
  103.       fi[1 ] = (f1+f3);
  104.      }
  105.  
  106.  k=0;
  107.  do
  108.     {
  109.      FLOAT s1,c1;
  110.      k  += 2;
  111.      k1  = 1  << k;
  112.      k2  = k1 << 1;
  113.      k4  = k2 << 1;
  114.      k3  = k2 + k1;
  115.      kx  = k1 >> 1;
  116.      fi  = fz;
  117.      gi  = fi + kx;
  118.      fn  = fz + n;
  119.      do
  120.         {
  121.          FLOAT g0,f0,f1,g1,f2,g2,f3,g3;
  122.          f1      = fi[0 ] - fi[k1];
  123.          f0      = fi[0 ] + fi[k1];
  124.          f3      = fi[k2] - fi[k3];
  125.          f2      = fi[k2] + fi[k3];
  126.          fi[k2]  = f0      - f2;
  127.          fi[0 ]  = f0      + f2;
  128.          fi[k3]  = f1      - f3;
  129.          fi[k1]  = f1      + f3;
  130.          g1      = gi[0 ] - gi[k1];
  131.          g0      = gi[0 ] + gi[k1];
  132.          g3      = SQRT2  * gi[k3];
  133.          g2      = SQRT2  * gi[k2];
  134.          gi[k2]  = g0      - g2;
  135.          gi[0 ]  = g0      + g2;
  136.          gi[k3]  = g1      - g3;
  137.          gi[k1]  = g1      + g3;
  138.          gi     += k4;
  139.          fi     += k4;
  140.         } while (fi<fn);
  141.      t_c = costab[k];
  142.      t_s = sintab[k];
  143.      c1 = 1;
  144.      s1 = 0;
  145.      for (i=1;i<kx;i++)
  146.         {
  147.      FLOAT c2,s2;
  148.          FLOAT t = c1;
  149.          c1 = t*t_c - s1*t_s;
  150.          s1 = t*t_s + s1*t_c;
  151.          c2 = c1*c1 - s1*s1;
  152.          s2 = 2*(c1*s1);
  153.          fn = fz + n;
  154.          fi = fz +i;
  155.          gi = fz +k1-i;
  156.          do
  157.         {
  158.          FLOAT a,b,g0,f0,f1,g1,f2,g2,f3,g3;
  159.          b       = s2*fi[k1] - c2*gi[k1];
  160.          a       = c2*fi[k1] + s2*gi[k1];
  161.          f1      = fi[0 ]    - a;
  162.          f0      = fi[0 ]    + a;
  163.          g1      = gi[0 ]    - b;
  164.          g0      = gi[0 ]    + b;
  165.          b       = s2*fi[k3] - c2*gi[k3];
  166.          a       = c2*fi[k3] + s2*gi[k3];
  167.          f3      = fi[k2]    - a;
  168.          f2      = fi[k2]    + a;
  169.          g3      = gi[k2]    - b;
  170.          g2      = gi[k2]    + b;
  171.          b       = s1*f2     - c1*g3;
  172.          a       = c1*f2     + s1*g3;
  173.          fi[k2]  = f0        - a;
  174.          fi[0 ]  = f0        + a;
  175.          gi[k3]  = g1        - b;
  176.          gi[k1]  = g1        + b;
  177.          b       = c1*g2     - s1*f3;
  178.          a       = s1*g2     + c1*f3;
  179.          gi[k2]  = g0        - a;
  180.          gi[0 ]  = g0        + a;
  181.          fi[k3]  = f1        - b;
  182.          fi[k1]  = f1        + b;
  183.          gi     += k4;
  184.          fi     += k4;
  185.         } while (fi<fn);
  186.         }
  187.     } while (k4<n);
  188. }
  189.  
  190. void fft(FLOAT *x_real, FLOAT *energy, FLOAT *ax, FLOAT *bx, int N)
  191. {
  192.  FLOAT a,b;
  193.  int i,j;
  194.  
  195.  fht(x_real,N);
  196.  
  197.  
  198.  energy[0] = x_real[0] * x_real[0];
  199.  ax[0] = bx[0] = x_real[0];
  200.  
  201.  for (i=1,j=N-1;i<N/2;i++,j--) {
  202.    a = ax[i] = x_real[i];
  203.    b = bx[i] = x_real[j];
  204.    energy[i]=(a*a + b*b)/2;
  205.    
  206.    if (energy[i] < 0.0005) {
  207.      energy[i] = 0.0005;
  208.      ax[i] = bx[i] = sqrt(0.0005);
  209.    }
  210.  }
  211.  energy[N/2] = x_real[N/2] * x_real[N/2];
  212.  ax[N/2] = bx[N/2] = x_real[N/2];
  213. }
  214.  
  215.  
  216.  
  217.  
  218. double fft_side( FLOAT in[2][1024], int s)
  219. {
  220.  double energy,t;
  221.  FLOAT a,b;
  222.  int i,j;
  223.  
  224.  if(!s)
  225.  {
  226.   energy = (in[0][0] - in[1][0]) * (in[0][0] - in[1][0])/2;
  227.   s++;
  228.  }
  229.  else
  230.   energy=0.0;
  231.  
  232.  for (i=s,j=1024-s;i<512;i++,j--) {
  233.   a = in[0][i] - in[1][i];
  234.   b = in[0][j] - in[1][j];
  235.  
  236.   t=(a*a + b*b)/4;
  237.  
  238.      if (t < 0.0005) {
  239.        t = 0.0005;
  240.      }
  241.    energy += t;
  242.  }
  243.  return energy + (in[0][512] - in[1][512]) * (in[0][512] - in[1][512])/2;
  244.  
  245.  
  246. }
  247.  
  248.  
  249.