fft.cpp.txt

FFT 算法的C语言班版描述。本人亲自编写

#include "stdio.h"
#include <math.h>       //调用源（头）文件
#define  N   64

//------/*定义一个复数结构*/-----------------------------------------
struct    compx          
  {float  real;
   float  imag;
  };
//-------------------------------------------------------------------
//*********************************************************************
//********程序接口！由ad采集的数据放入 结构体S的数组中！***************
//********注意字节的大小和长度！                        ***************
//*********************************************************************
struct compx   s[N];         //FFT输入输出均从是s[1]开始存入  
struct  compx temp; 
struct compx   EE (struct compx ,struct compx); //定义复数相乘结构

 void  FFT (struct compx  *xin ,int n)        ;  /*定义FFT函数*/


//----两复数相乘的程序------------------------------------------------

struct compx EE(struct compx a1,struct compx b2) 

{struct compx b3;                          //b3保存两复数间的结果
    b3.real=a1.real*b2.real-a1.imag*b2.imag; //两复数间的运算
    b3.imag=a1.real*b2.imag+a1.imag*b2.real;
return(b3);    /*返回结果*/
}
//----------------------------------------------------
//----------------------------------------------------
//------- /*FFT函数体*/

void  FFT(struct compx *xin ,int n)  
 {
       /*定义变量*/
     int      f ;
	 int	  m ;
	 int     nv2;
	 int     nm1;
	 int       i; 
     int       k; 
	 int     j=0;                  /*j 存是i的倒位序二进制数*/
	 int    leng;                        
  struct  compx   t;          //定义结构变量*/v, w,
     nv2=N/2;                     /*最高位值的权值*/
     f=n;                        /*f为中间变量*/
    for(m=1;(f=f/2)!=1;m++)
	{;}                         /*求级数m*/
   nm1=n-1;       /*nm1为数组长度*/
 //--对数组进行预处理的的程序-----
//************************************************************

  for(i=0;i<n;i++)                   /*倒序*/
  {
	 if(i<j)                        /*j应该是i的倒位序二进制数。。*/                  
	 {       t=xin[j];
        xin[j]=xin[i];
        xin[i]=t;
     }                             /*i<j则换位*/
     k=nv2;                       /*k为倒序中相应位置的权值*/
     while(k<=j) 
     {j=j-k;
      k=k/2;  
     if(k==0 ) break;//     || j==0
	 }                   /*k<j时最高为变为0*/
    j=j+k;              /* j为数组中的位数，是一个十进制数*/ 
   }
 
//***************************************************************
//***********蝶形运算**************自己编的****************
     int    le ,lei,u,er;
     int    xu  ;                         // 每一级的角度旋转倍乘因子
     double  pi =2*3.14159265	;	 
     for  (leng=1;leng<=m;leng++)                /*leng控制级数     */
	 {    le  =pow(2,leng);                       /*le等于2的leng次方*/
          lei =le/2;                             /*蝶形两节点间的距离*/
	      xu  =n/le;
		  for( i=0;i<n;i=i+le)            //  计算每一个蝴蝶
		  {
			  for(j=i+lei,u=0;j<i+le;j++,u++ )           //先对蝴蝶的下面的尾巴做变换；
			  { 
			     temp.real= cos(pi*xu*u/n);
				 temp.imag=(-sin(pi*xu*u/n));
				 xin[j]=EE(temp,xin[j]); 		  
			  }
			  for(j=i;j<i+lei;j++)                
			  {   er=j+lei;
				  temp              =xin[j];	
				 xin[j].real	    =xin[j].real+xin[er].real;
			     xin[j].imag	    =xin[j].imag+xin[er].imag;
			     xin[er].real       =temp.real  -xin[er].real;
				 xin[er].imag       =temp.imag  -xin[er].imag;
			  }  
		  
		  }	 
	 }

}

//----------------------------------------------------
//---*定义主函数*/------------------------------------
void    main()                 
 {  
    float  qq ;	
    int i;     //变量初始化，N为总点数，i为每点数
   for(i=0;i<10;i++)
   {s[i].real=0.0;s[i].imag=0.0;}
   for(i=10;i<N;i++)
   {s[i].real=exp(-(i-10)/15.0)*sin(6.283153*(i-10)/16.0);s[i].imag=0.0;}
//----------------------------------------------------  
    FFT(s,N);               /*调用FFt*/
//----------------------------------------------------
 for(i=1 ; i<=N ;i++)                 /*输出*/
 {    qq = (s[i].real * s[i].real) + (s[i].imag+s[i].imag);
      qq=  sqrt(qq);
	   printf("%10.7f",s[ i ].real);
       printf(" + ");
       printf("%10.7f",s[ i ].imag);
       printf("j  q= %10.7f   ",qq);
      printf("\n");
    }
//---------------------------------------------------- 
 }