home *** CD-ROM | disk | FTP | other *** search
/ Da Capo / da_capo_vol1.bin / programs / amiga / misc / mpegaudio / subs.c < prev    next >
C/C++ Source or Header  |  1994-03-21  |  6KB  |  163 lines
  1. /**********************************************************************
  2. Copyright (c) 1991 MPEG/audio software simulation group, All Rights Reserved
  3. subs.c
  4. **********************************************************************/
  5. /**********************************************************************
  6.  * MPEG/audio coding/decoding software, work in progress              *
  7.  *   NOT for public distribution until verified and approved by the   *
  8.  *   MPEG/audio committee.  For further information, please contact   *
  9.  *   Davis Pan, 508-493-2241, e-mail: pan@3d.enet.dec.com             *
  10.  *                                                                    *
  11.  * VERSION 3.9                                                        *
  12.  *   changes made since last update:                                  *
  13.  *   date   programmers         comment                               *
  14.  * 2/25/91  Davis Pan           start of version 1.0 records          *
  15.  * 5/10/91  W. Joseph Carter    Ported to Macintosh and Unix.         *
  16.  * 7/10/91  Earle Jennings      Ported to MsDos from Macintosh        *
  17.  *                              Replacement of one float with FLOAT   *
  18.  * 2/11/92  W. Joseph Carter    Added type casting to memset() args.  *
  19.  * 4/27/92  Masahiro Iwadare    Added 256 point version for Layer III *
  20.  **********************************************************************/
  21.  
  22. #include "common.h"
  23. #include "encoder.h"
  24.  
  25. /*****************************************************************************
  26.  ************************** Start of Subroutines *****************************
  27.  *****************************************************************************/
  28.  
  29. /*****************************************************************************
  30.  * FFT computes fast fourier transform of BLKSIZE samples of data            *
  31.  *   uses decimation-in-frequency algorithm described in "Digital            *
  32.  *   Signal Processing" by Oppenheim and Schafer, refer to pages 304         *
  33.  *   (flow graph) and 330-332 (Fortran program in problem 5)                 *
  34.  *   to get the inverse fft, change line 20 from                             *
  35.  *                 w_imag[L] = -sin(PI/le1);                                 *
  36.  *                          to                                               *
  37.  *                 w_imag[L] = sin(PI/le1);                                  *
  38.  *                                                                           *
  39.  *   required constants:                                                     *
  40.  *         #define      PI          3.14159265358979                         *
  41.  *         #define      BLKSIZE     1024                                     *
  42.  *         #define      LOGBLKSIZE  10                                       *
  43.  *         #define      BLKSIZE_S   256                                      *
  44.  *         #define      LOGBLKSIZE_S 8                                       *
  45.  *                                                                           *
  46.  *****************************************************************************/
  47. #define      BLKSIZE_S   256
  48. #define      LOGBLKSIZE_S 8
  49.  
  50. void fft(x_real,x_imag, energy, phi, N)
  51. FLOAT x_real[BLKSIZE], x_imag[BLKSIZE], energy[BLKSIZE], phi[BLKSIZE];
  52. int    N;
  53. {
  54.  int     M,MM1;
  55.  static int     init=0;
  56.  int     NV2, NM1, MP;
  57.  static double  w_real[2][LOGBLKSIZE], w_imag[2][LOGBLKSIZE];
  58.  int            i,j,k,L;
  59.  int            ip, le,le1;
  60.  double         t_real, t_imag, u_real, u_imag;
  61.  
  62.  if(init==0) {
  63.     memset((char *) w_real, 0, sizeof(w_real));  /* preset statics to 0 */
  64.     memset((char *) w_imag, 0, sizeof(w_imag));  /* preset statics to 0 */
  65.     M = LOGBLKSIZE;
  66.     for(L=0; L<M; L++){
  67.        le = 1 << (M-L);
  68.        le1 = le >> 1;
  69.        w_real[0][L] = cos(PI/le1);
  70.        w_imag[0][L] = -sin(PI/le1);
  71.     }          
  72.     M = LOGBLKSIZE_S;
  73.     for(L=0; L<M; L++){
  74.        le = 1 << (M-L);
  75.        le1 = le >> 1;
  76.        w_real[1][L] = cos(PI/le1);
  77.        w_imag[1][L] = -sin(PI/le1);
  78.     }          
  79.     init++;
  80.  }
  81.  switch(N) {
  82.     case BLKSIZE:
  83.             M = LOGBLKSIZE;
  84.             MP = 0;
  85.             break;
  86.     case BLKSIZE_S:
  87.             M = LOGBLKSIZE_S;
  88.             MP = 1;
  89.             break;
  90.     default:    printf("Error: Bad FFT Size in subs.c\n");
  91.             exit(-1);
  92.  }
  93.  MM1 = M-1;
  94.  NV2 = N >> 1;
  95.  NM1 = N - 1;
  96.  for(L=0; L<MM1; L++){
  97.     le = 1 << (M-L);
  98.     le1 = le >> 1;
  99.     u_real = 1;
  100.     u_imag = 0;
  101.     for(j=0; j<le1; j++){
  102.        for(i=j; i<N; i+=le){
  103.           ip = i + le1;
  104.           t_real = x_real[i] + x_real[ip];
  105.           t_imag = x_imag[i] + x_imag[ip];
  106.           x_real[ip] = x_real[i] - x_real[ip];
  107.           x_imag[ip] = x_imag[i] - x_imag[ip];
  108.           x_real[i] = t_real;
  109.           x_imag[i] = t_imag;
  110.           t_real = x_real[ip];
  111.           x_real[ip] = x_real[ip]*u_real - x_imag[ip]*u_imag;
  112.           x_imag[ip] = x_imag[ip]*u_real + t_real*u_imag;
  113.        }
  114.        t_real = u_real;
  115.        u_real = u_real*w_real[MP][L] - u_imag*w_imag[MP][L];
  116.        u_imag = u_imag*w_real[MP][L] + t_real*w_imag[MP][L];
  117.     }
  118.  }
  119.  /* special case: L = M-1; all Wn = 1 */
  120.  for(i=0; i<N; i+=2){
  121.     ip = i + 1;
  122.     t_real = x_real[i] + x_real[ip];
  123.     t_imag = x_imag[i] + x_imag[ip];
  124.     x_real[ip] = x_real[i] - x_real[ip];
  125.     x_imag[ip] = x_imag[i] - x_imag[ip];
  126.     x_real[i] = t_real;
  127.     x_imag[i] = t_imag;
  128.     energy[i] = x_real[i]*x_real[i] + x_imag[i]*x_imag[i];
  129.     if(energy[i] <= 0.0005){phi[i] = 0;energy[i] = 0.0005;}
  130.     else phi[i] = atan2((double) x_imag[i],(double) x_real[i]);
  131.     energy[ip] = x_real[ip]*x_real[ip] + x_imag[ip]*x_imag[ip];
  132.     if(energy[ip] == 0)phi[ip] = 0;
  133.     else phi[ip] = atan2((double) x_imag[ip],(double) x_real[ip]);
  134.  }
  135.  /* this section reorders the data to the correct ordering */
  136.  j = 0;
  137.  for(i=0; i<NM1; i++){
  138.     if(i<j){
  139. /* use this section only if you need the FFT in complex number form *
  140.  * (and in the correct ordering)                                    */
  141.        t_real = x_real[j];
  142.        t_imag = x_imag[j];
  143.        x_real[j] = x_real[i];
  144.        x_imag[j] = x_imag[i];
  145.        x_real[i] = t_real;
  146.        x_imag[i] = t_imag;
  147. /* reorder the energy and phase, phi                                        */
  148.        t_real = energy[j];
  149.        energy[j] = energy[i];
  150.        energy[i] = t_real;
  151.        t_real = phi[j];
  152.        phi[j] = phi[i];
  153.        phi[i] = t_real;
  154.     }
  155.     k=NV2;
  156.     while(k<=j){
  157.        j = j-k;
  158.        k = k >> 1;
  159.     }
  160.     j = j+k;
  161.  }
  162. }
  163.