fft.c

通过读取已有的由matlab产生的数据文件进行1024点浮点fft运算的c程序.其中,radix2()实现基2算法,ChangeOrder(）实现变址运算.

#define PI 3.14159265358979
#define M 10
#define N 1024
#include <stdio.h>
#include <math.h>
/*声明函数*/
void radix2(float *xr, float *xi, float *wr, float *wi);
void ChangeOrder(float *xr, float *xi);
/*定义全局变量*/
float xr[N],xi[N];       //定义输入数列的实部和虚部数组
float x[2*N];            //定义输入数列数组    
float xm[N];             //定义输入数列的模的数组
float x1r[N],x1i[N];     //定义输出数列的实部和虚部数组
float w1[N],w2[N];       //定义旋转因子的实部和虚部数组
float y[2*N];            //定义输出数列的数组
float ym[N];             //定义输出数列的模的数组
int nx1;                 //定义FFT计算的点数
int i;     
 
/*主函数开始*/
main()
{
  float delta;
  FILE *fp;

   nx1=N;
  
  fp=fopen("sine.dat","r");
  if(!fp) 
   { 
     printf("Can not open file!\n"); 
     exit(0); 
     } 
  for (i=0; i<N; i++)
   	 fscanf(fp,"%f\n",&xr[i]);
  fclose(fp); 
   	
   	
  /*init the fft coefficients*/
  delta=(float)(2*PI/nx1);
  for (i=0; i<nx1/2; i++)
      {
        w1[i]=(float)(cos(i*delta));
        w2[i]=(float)(sin(i*delta));
      }
/*构造要做fft的序列*/
  for (i=0; i<nx1; i++)
      {
        x1r[i]=xr[i];
        xi[i]=x1i[i]=0.0;
        x[2*i]=x1r[i];     //构造要做FFT的序列，按实部0、虚部0、实部1、虚部1……顺序构造
        x[2*i+1]=x1i[i];
      }
  for (i=0; i<nx1; i++)
      {
        xm[i]=(float)(sqrt((xr[i]*xr[i])+(xi[i]*xi[i])));    //计算输入数组的模
      }
  ChangeOrder(x1r,x1i);                             //调用倒序子程序进行排序
/*fft*/ 
  radix2(x1r,x1i,w1,w2);                            //调用FFT算法子程序
  for (i=0; i<nx1; i++)
      {
        float sum=0.0;
        sum=(float)(sqrt((x1r[i]*x1r[i])+(x1i[i]*x1i[i])));  //计算输出序列的模，放在ym数组中
        ym[i]=sum;
        y[2*i]=(float)(sqrt((x1r[i])*(x1r[i])));         
        y[2*i+1]=(float)(sqrt((x1i[i])*(x1i[i])));
      }    
}

/*基2fft算法*/
  void radix2(float xr[], float xi[], float wr[], float wi[])
  {
    int L,B,J,P,k;
    float rPartKB,iPartKB;
    for (L=1; L<=M; L++)
        {
          B=(int)(pow(2,(L-1)));
          for (J=0; J<=B-1; J++)
              {  
                P=J*(int)(pow(2,(M-L)));
                for (k=J; k<=N-1; k+=(int)(pow(2,L)))
                    {
                      rPartKB=xr[k+B]*wr[P]+xi[k+B]*wi[P];
                      iPartKB=xi[k+B]*wr[P]-xr[k+B]*wi[P];                  
                      xr[k+B]=xr[k]-rPartKB;
                      xi[k+B]=xi[k]-iPartKB;
                      xr[k]=xr[k]+rPartKB;
                      xi[k]=xi[k]+iPartKB;
                    }
              }
        }
  }
  
  /*倒序子程序*/
  void ChangeOrder(float xr[], float xi[])
  {  
    int LH,N1,I,J,K;
    float T;
    LH=N/2; J=LH; N1=N-2;
    for (I=1; I<=N1; I++)
       {
         if(I<J)
         {
           T=xr[I]; xr[I]=xr[J]; xr[J]=T;
           T=xi[I]; xi[I]=xi[J]; xi[J]=T;
         }
         K=LH;
         while(J>=K)
         {
           J=J-K;
           K=K/2;
         }
         J=J+K;
       }
  }