rfast32.c

dsp AD公司ADSP21的代码,里面有FFT FIR IIR EQULIZER G722_21F 等可以在项目中直接应用的代码.此代码的来源是ADI公司自己出版的书籍,此书在美国购得

#include <math.h>
#include <libap.h>
#include "rtdspc.h"
#include "filter.h"

/***********************************************************************

RFAST.C - Realtime fast convolution using the FFT

This program performs fast convolution using the FFT.  It performs
the convolution required to implement a 35 point FIR filter
(stored in variable fir_lpf35) on an
arbitrary length realtime input.  The filter is
a LPF with 40 dB out of band rejection.  The 3 dB point is at a
relative frequency of approximately .25*fs.

The filter will switch in and out every 32 blocks (32*FFT_LENGTH).

************************************************************************/

/* FFT length must be a power of 2 */
#define FFT_LENGTH 1024
#define M 10         /* must be log2(FFT_LENGTH) */
#define FILTER_LENGTH 35

    float input_save[FILTER_LENGTH];
    COMPLEX samp[FFT_LENGTH], filt[FFT_LENGTH];
    int bypass;
    int count;

void main()
{
  register int          i, j;
  register float        tempflt,inv_len;

  bypass = 0;

  int_disable();

/* copy the filter into complex array and scale by 1/N for inverse FFT */
  inv_len = 1.0/FFT_LENGTH;  
  for(i = 0 ; i < FILTER_LENGTH ; i++) {
      filt[i].real = inv_len*fir_lpf35[i];
      filt[i].imag = 0.0;
  }
  for( ; i < FFT_LENGTH ; i++) {
      filt[i].real = 0.0;
      filt[i].imag = 0.0;
  }

/* FFT the zero filled filter impulse response */
  ffta(FFT_LENGTH,M,filt);

/* read in one FFT worth of samples to start */
  for(i = 0 ; i < FFT_LENGTH-FILTER_LENGTH ; i++) {
    samp[i].real = getinput();
    samp[i].imag = 0.0;
  }

/* save the last FILTER_LENGTH points for next time */
  for(j = 0 ; j < FILTER_LENGTH ; j++, i++) {
    input_save[j] = samp[i].real = getinput();
    samp[i].imag = 0.0;
  }

  count = 0;
  while(1) {

/* do FFT of samples */
    ffta(FFT_LENGTH,M,samp);

/* Multiply the two transformed sequences */
/* swap the real and imag outputs to allow a forward FFT instead of
inverse FFT */
    if(!bypass) {
     for(i = 0 ; i < FFT_LENGTH ; i++) {
      tempflt = samp[i].real * filt[i].real
                             - samp[i].imag * filt[i].imag;
      samp[i].real = samp[i].real * filt[i].imag
                             + samp[i].imag * filt[i].real;
      samp[i].imag = tempflt;
     }
    }
    else {
     for(i = 0 ; i < FFT_LENGTH ; i++) {
      tempflt = inv_len*samp[i].real;
      samp[i].real = inv_len*samp[i].imag;
      samp[i].imag = tempflt;
     }
    }

/* Inverse fft the multiplied sequences */
    ffta(FFT_LENGTH,M,samp);

/* Write the result out */
/* because a forward FFT was used for the inverse FFT,
the output is in the imag part */
    for(i = FILTER_LENGTH ; i < FFT_LENGTH ; i++) sendout(samp[i].imag);

/* overlap the last FILTER_LENGTH-1 input data points in the next FFT */
    for(i = 0; i < FILTER_LENGTH ; i++) {
      samp[i].real = input_save[i];
      samp[i].imag = 0.0;
    }

    for( ; i < FFT_LENGTH-FILTER_LENGTH ; i++) {
      samp[i].real = getinput();
      samp[i].imag = 0.0;
    }

/* save the last FILTER_LENGTH points for next time */
    for(j = 0 ; j < FILTER_LENGTH ; j++, i++) {
      input_save[j] = samp[i].real = getinput();
      samp[i].imag = 0.0;
    }
/* toggle bypass switch every 32 blocks */
    if(((count++) & 31) == 0) bypass ^= 1;
  }
}