rtpse.c

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

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "rtdspc.h"

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

RTPSE.C - Real-Time Power spectral estimation using the FFT

This program does power spectral estimation on input samples.
The average power spectrum in each block is determined
and used to generate a series of outputs.

Length of each FFT snapshot: 64 points
Number of FFTs to average: 16 FFTs
Amount of overlap between each FFT: 60 points

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

/* FFT length must be a power of 2 */
#define FFT_LENGTH 64
#define M 6            /* must be log2(FFT_LENGTH) */

/* these variables global so they can be changed in real-time */
  int numav = 16;
  int ovlap = 60;

main()
{
  int            i,j,k;
  float          scale,tempflt;
  static float   mag[FFT_LENGTH], sig[FFT_LENGTH], hamw[FFT_LENGTH];
  static COMPLEX samp[FFT_LENGTH];

/* overall scale factor */
  scale = 1.0f/(float)FFT_LENGTH;
  scale *= scale/(float)numav;

/* calculate hamming window */
  tempflt = 8.0*atan(1.0)/(FFT_LENGTH-1);
  for(i = 0 ; i < FFT_LENGTH ; i++)
    hamw[i] = 0.54 - 0.46*cos(tempflt*i);
  
/* read in the first FFT_LENGTH samples, overlapped samples read in loop */
  for(i = 0 ; i < FFT_LENGTH ; i++) sig[i] = getinput();

  for(;;) {

    for (k=0; k<FFT_LENGTH; k++) mag[k] = 0;

    for (j=0; j<numav; j++){

      for (k=0; k<FFT_LENGTH; k++){
        samp[k].real = hamw[k]*sig[k];
        samp[k].imag = 0;
      }

      fft(samp,M);

      for (k=0; k<FFT_LENGTH; k++){
        tempflt  = samp[k].real * samp[k].real;
        tempflt += samp[k].imag * samp[k].imag;
        tempflt = scale*tempflt;
        mag[k] += tempflt;
      }

/* overlap the new samples with the old */
      for(k = 0 ; k < ovlap ; k++) sig[k] = sig[k+FFT_LENGTH-ovlap];
      for( ; k < FFT_LENGTH ; k++) sig[k] = getinput();
    }

/*  Take log after averaging the magnitudes.  */
    for (k=0; k<FFT_LENGTH/2; k++){
      tempflt = mag[k];
      if(tempflt < 1.e-10f) tempflt = 1.e-10f;
      sendout(10.0f*log10(tempflt));
    }
  }
}