rfast.c

来自「dsp AD公司ADSP21的代码,里面有FFT FIR IIR EQULIZE」· C语言 代码 · 共 103 行

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.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.

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

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

void main()
{
  int          i, j;
  float        tempflt;
  COMPLEX      *samp, *filt;
  static float input_save[FILTER_LENGTH];

/* power of 2 length of FFT and complex allocation */
  samp = (COMPLEX *) calloc(FFT_LENGTH, sizeof(COMPLEX));
  if(!samp){
    exit(1);
  }

/* Zero fill the filter to the sequence length */
  filt = (COMPLEX *) calloc(FFT_LENGTH, sizeof(COMPLEX));
  if(!filt){
    exit(1);
  }

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

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

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

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

  while(1) {

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

/* Multiply the two transformed sequences */
/* swap the real and imag outputs to allow a forward FFT instead of
inverse FFT */
    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;
    }

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

/* 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;
    }
  }
}