fftdll.cpp

实现正弦信号的快速傅里叶变换

#include "fftdll.h"
#include <math.h>
#include <afxwin.h>
#define M_PI 3.1416



void fft(double *xr,double *xi,int n)
 {
   int i,j,a,b,k;
   int ep,arg,mt,s0,s1;
   double pr,pi,ph,t1,t2;
   double *s,*c;
   c=(double *)calloc(n,sizeof(double));
   if(c==NULL) exit(1);
   s=(double *)calloc(n,sizeof(double));
   if(s==NULL) exit(1);
   j=0;
   
   for(i=0;i<n-1;i++)
   {
     if(i<j)
     {

		 t1=xr[i];xr[i]=xr[j];xr[j]=t1;
         t2=xi[i];xi[i]=xi[j];xi[j]=t2;

     }
     k=n/2;
     while(k<=j)
     {
	j=j-k;
	k=n/2;
     }
     j=j+k;
   }
   ep=0;
   i=n;
   while(i!=1)
   {
      ep=ep+1;
      i=i/2;
   }
   ph=2*M_PI/n;
   for(i=0;i<n;i++)
   {
     s[i]=sin(ph*i);
     c[i]=cos(ph*i);
   }
   a=2;
   b=1;
   for(mt=1;mt<=ep;mt++)
   {
      s0=n/a;
      s1=0;
      for(k=0;k<b;k++)
      {
	i=k;
	while(i<n)
	{
	  arg=i+b;
	  if(k==0)
	  {
	    pr=xr[arg];
	    pi=xi[arg];
	  }
	  else
	  {
	    pr=xr[arg]*c[s1]-xi[arg]*s[s1];
	    pi=xr[arg]*s[s1]+xi[arg]*c[s1];
	  }
	  xr[arg]=xr[i]-pr;
	  xi[arg]=xi[i]-pi;
	  xr[i]=xr[i]+pr;
	  xi[i]=xi[i]+pi;
	  i=i+a;
	}
	s1=s1+s0;
      }
      a=2*a;
      b=b*2;
   }
   free(c);
   free(s);
 }