fft.c

2N点实序列FFT_C程序 上《数字信号处理》课程时编的

#include  <math.h>
#include  <stdlib.h>
# define pi 3.14159265
struct compx { double  real, imag;};/****定义结构体类型****/
struct compx *fft(struct compx *,  int);/**FFT函数声明**/

/******定义复数乘法运算函数******/
struct compx  EE(struct compx b1,struct compx b2)
{struct compx  b3;
b3.real=b1.real*b2.real-b1.imag*b2.imag;
b3.imag=b1.real*b2.imag+b1.imag*b2.real;
return  (b3);
}

/**FFT函数定义，本函数参数如下：
输入：
  xin  为结构体compx  数组的首地址，相当于需要处理的复数信号的首地址；
  N    为结构体compx  数组的元素的个数，即信号的点数；
输出：
  xin  为结构体compx  数组的首地址,原址输出。  **/

struct compx *fft(struct compx * xin , int N)
{
  int m,LH,i,k,j,l,f;
  struct compx * y;
  struct compx  v,w,t ;

      /**开始计算倒序值**/
  f=N;
  for (m=1; (f=f/2)!=1; m++) {;}
  y=(struct compx *)calloc(N,sizeof(struct compx));
  LH=N/2;   j=0;
  for (i=1;i<=N-2;i++)  
		{    if(i<j)  {t=xin[j];  xin[j]=xin[i];  xin[i]=t;	}
             k=LH; 
	         while  (k<=j)  {j=j-k;  k=k/2;}
	         j=j+k; 	
	         xin[i]=y[j];
		}

      /**开始计算FFT算法**/

	  for(l=1;l<=m;l++)
	  { int le,B,ip; /*ip 相当于算法中的k+B  */
		le=(int) pow(2,l);  /*le 相当于算法中 2的L次方 */
	    B=le/2;
		v.real=1.0; v.imag=0.0;
        w.real=cos(pi/B);  w.imag=-sin(pi/B); /*W 相当于算法中的相当于旋转因子 */

		for(j=0;j<=B-1;j++)
		{ 
			if(j=0){ v.real=1.0; v.imag=0.0;}  /*第一级旋转因子为数值为1*/
			else {v=EE(v,w); }    /* 求旋转因子的P次方 P 等于J 乘以2的M—L次方*/

			for (k=j;i<=N-1;i=i+le)
			{ip=i+B;
			t=EE(xin[ip],v);   /* t 即算法中的  X(k+B)*W_N_P */ 

			xin[ip].real=xin[i].real-t.real;  
			xin[ip].imag=xin[i].imag-t.imag;   /*X(k+B)=X(k)-X(k+B)*W_N_P*/ 

            xin[i].real=xin[i].real+t.real;
			xin[i].imag=xin[i].imag-t.imag;  /*X(k)=X(k)+X(k+B)*W_N_P*/ 
			}	
		}
	  }  

return  xin;

}