home *** CD-ROM | disk | FTP | other *** search
/ Audio 4.94 - Over 11,000 Files / audio-11000.iso / msdos / sndbords / proaudio / freq3 / realfft.c < prev    next >
C/C++ Source or Header  |  1993-09-07  |  4KB  |  161 lines
  1. /*
  2.  *     Program: REALFFT.C
  3.  *      Author: Philip VanBaren
  4.  *        Date: 2 September 1993
  5.  *
  6.  * Description: These routines perform an FFT on real data.
  7.  *              On a 486/33 compiled using Borland C++ 3.1 with full
  8.  *              speed optimization and a small memory model, a 1024 point 
  9.  *              FFT takes about 16ms.
  10.  *              This code is for integer data, but could be converted
  11.  *              to float or double simply by changing the data types
  12.  *              and getting rid of the bit-shifting necessary to prevent
  13.  *              overflow/underflow in fixed-point calculations.
  14.  *
  15.  *  Note: Output is BIT-REVERSED! so you must use the BitReversed to
  16.  *        get legible output, (i.e. Real_i = buffer[ BitReversed[i] ]
  17.  *                                  Imag_i = buffer[ BitReversed[i]+1 ] )
  18.  *        Input is in normal order.
  19.  */
  20.  
  21. #include <math.h>
  22. #include <stdlib.h>
  23. #include <stdio.h>
  24.  
  25. int *BitReversed;
  26. int *SinTable;
  27. int Points = 0;
  28.  
  29. /*
  30.  *  Initialize the Sine table and Twiddle pointers (bit-reversed pointers)
  31.  *  for the FFT routine.
  32.  */
  33. void InitializeFFT(int fftlen)
  34. {
  35.    int i;
  36.    int temp;
  37.    int mask;
  38.  
  39.    /*
  40.     *  FFT size is only half the number of data points
  41.     *  The full FFT output can be reconstructed from this FFT's output.
  42.     *  (This optimization can be made since the data is real.)
  43.     */
  44.    Points = fftlen/2;
  45.  
  46.    if((SinTable=(int *)malloc(2*Points*sizeof(int)))==NULL)
  47.    {
  48.       puts("Error allocating memory for Sine table.");
  49.       exit(1);
  50.    }
  51.    if((BitReversed=(int *)malloc(Points*sizeof(int)))==NULL)
  52.    {
  53.       puts("Error allocating memory for BitReversed.");
  54.       exit(1);
  55.    }
  56.  
  57.    for(i=0;i<Points;i++)
  58.    {
  59.       temp=0;
  60.       for(mask=Points/2;mask>0;mask >>= 1)
  61.          temp=(temp >> 1) + (i&mask ? Points : 0);
  62.  
  63.       BitReversed[i]=temp;
  64.    }
  65.  
  66.    for(i=0;i<Points;i++)
  67.    {
  68.       SinTable[BitReversed[i]  ]=floor(-32768*sin(2*M_PI*i/(2*Points))+0.5);
  69.       SinTable[BitReversed[i]+1]=floor(-32768*cos(2*M_PI*i/(2*Points))+0.5);
  70.    }
  71. }
  72.  
  73. /*
  74.  *  Free up the memory allotted for Sin table and Twiddle Pointers
  75.  */
  76. void EndFFT()
  77. {
  78.    free(BitReversed);
  79.    free(SinTable);
  80.    Points=0;
  81. }
  82.  
  83. int *A,*B;
  84. int *sptr;
  85. int *endptr1,*endptr2;
  86. int *br1,*br2;
  87. long HRplus,HRminus,HIplus,HIminus;
  88.  
  89. /*
  90.  *  Actual FFT routine.  Must call InitializeFFT(fftlen) first!
  91.  */
  92. void realfft(int *buffer)
  93. {
  94.    int ButterfliesPerGroup=Points/2;
  95.  
  96.    endptr1=buffer+Points*2;
  97.  
  98.    /*
  99.     *  Butterfly:
  100.     *     Ain-----Aout
  101.     *         \ /
  102.     *         / \
  103.     *     Bin-----Bout
  104.     */
  105.  
  106.    while(ButterfliesPerGroup>0)
  107.    {
  108.       A=buffer;
  109.       B=buffer+ButterfliesPerGroup*2;
  110.       sptr=SinTable;
  111.  
  112.       while(A<endptr1)
  113.       {
  114.          register int sin=*sptr;
  115.          register int cos=*(sptr+1);
  116.          endptr2=B;
  117.          while(A<endptr2)
  118.          {
  119.             long v1=((long)*B*cos + (long)*(B+1)*sin) >> 15;
  120.             long v2=((long)*B*sin - (long)*(B+1)*cos) >> 15;
  121.             *(A++)=(*(B++)=(*A+v1)>>1)-v1;
  122.             *(A++)=(*(B++)=(*A-v2)>>1)+v2;
  123.          }
  124.          A=B;
  125.          B+=ButterfliesPerGroup*2;
  126.          sptr+=2;
  127.       }
  128.       ButterfliesPerGroup >>= 1;
  129.    }
  130.    /*
  131.     *  Massage output to get the output for a real input sequence.
  132.     */
  133.    br1=BitReversed+1;
  134.    br2=BitReversed+Points-1;
  135.  
  136.    while(br1<br2)
  137.    {
  138.       register long temp1;
  139.       register long temp2;
  140.       int sin=SinTable[*br1];
  141.       int cos=SinTable[*br1+1];
  142.       A=buffer+*br1;
  143.       B=buffer+*br2;
  144.       HRplus = (HRminus = *A     - *B    ) + (*B << 1);
  145.       HIplus = (HIminus = *(A+1) - *(B+1)) + (*(B+1) << 1);
  146.       temp1  = ((long)sin*HRminus - (long)cos*HIplus) >> 15;
  147.       temp2  = ((long)cos*HRminus + (long)sin*HIplus) >> 15;
  148.       *B     = (*A     = (HRplus  + temp1) >> 1) - temp1;
  149.       *(B+1) = (*(A+1) = (HIminus + temp2) >> 1) - HIminus;
  150.  
  151.       br1++;
  152.       br2--;
  153.    }
  154.    /*
  155.     *  Handle DC bin separately
  156.     */
  157.    buffer[0]+=buffer[1];
  158.    buffer[1]=0;
  159. }
  160.  
  161.