mfccfft.cpp

这是我本人编写的一个32位定点小数运算的函数库。对于没有浮点运算器的场合

// MfccFFT.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include <FFT.h>
#include <math.h>
#include <conio.h>

//#define PI   3.14159265358979
#define TPI  6.28318530717959

void fft(double	s[], int invert, int len)
{
   int ii,jj,n,nn,limit,m,j,inc,i;
   double wx,wr,wpr,wpi,wi,theta;
   double xre,xri,x;
   
   n = len;
   nn=n / 2; j = 1;
   for (ii=1;ii<=nn;ii++) {
      i = 2 * ii - 1;
      if (j>i) {
         xre = s[j]; xri = s[j + 1];
         s[j] = s[i];  s[j + 1] = s[i + 1];
         s[i] = xre; s[i + 1] = xri;
      }
      m = n / 2;
      while (m >= 2  && j > m) {
         j -= m; m /= 2;
      }
      j += m;
   };
   limit = 2;
   while (limit < n) {
      inc = 2 * limit; theta = TPI / limit;
      if (invert) theta = -theta;
      x = sin(0.5 * theta);
      wpr = -2.0 * x * x; wpi = sin(theta); 
      wr = 1.0; wi = 0.0;
      for (ii=1; ii<=limit/2; ii++) {
         m = 2 * ii - 1;
         for (jj = 0; jj<=(n - m) / inc;jj++) {
            i = m + jj * inc;
            j = i + limit;
            xre = wr * s[j] - wi * s[j + 1];
            xri = wr * s[j + 1] + wi * s[j];
            s[j] = s[i] - xre; s[j + 1] = s[i + 1] - xri;
            s[i] = s[i] + xre; s[i + 1] = s[i + 1] + xri;
         }
         wx = wr;
         wr = wr * wpr - wi * wpi + wr;
         wi = wi * wpr + wx * wpi + wi;
      }
      limit = inc;
   }
   if (invert)
      for (i = 1;i<=n;i++) 
         s[i] = s[i] / nn;
   
   return;
}

#define N 1024

int main(int argc, char* argv[])
{
	double tdata[2*N+1];

	CComplex<double> temp[N];
	tdata[0]=0;
	for (int i=0; i<400; i++){
		double d1=(double)rand()/RAND_MAX,d2=(double)rand()/RAND_MAX;
		tdata[i*2+1]=d1; tdata[i*2+2]=d2;
		temp[i].real()=d2; temp[i].imag()=d1;
	}
	for (i=400;i<N;i++) {
		tdata[i*2+1]=0; tdata[i*2+2]=0;
		temp[i].real()=0; temp[i].imag()=0;
	}

	FFT(temp,N); fft(tdata,0,N*2);
	double sum=0;
	for (i=0;i<N;i++) {
		double d=fabs(temp[i].real()-tdata[i*2+2])+fabs(temp[i].imag()-tdata[i*2+1]);
		sum+=d;
	}
	printf("Delta=%lf\n",sum/N); getch();
	return 0;
}