ex4_3.c

c语言和c++是目前应用的比较流行的语言,它简单扼要,易理解,很有前途

/* Program for FFT */

#include <math.h>
#include <stdio.h>


const float pi=3.1415926;
int N;							/* FFT点数 */
float x_re[300], x_im[300];	 	/* 输入信号序列 */
float y_re[300], y_im[300];		/* 输出频谱序列 */
float w_re, w_im;				/* 蝶形因子 */
int m;							/* 蝶形运算的级数,即Log2(N) */
float t_re, t_im, v_re, v_im;  	/* 临时变量 */
int j,i,k,f,n;   				
int a, b, c;


main()
{ 	
	N=256;

	/* 初始化数据空间 */
	for(i=0; i<300; i++)
	{
		x_re[i]=0;
		x_im[i]=0;
	}
	
	/* 设定输入信号序列为单边指数函数 */
	/* 考虑到实际衰减很快，因此可以忽略后面大半部分数值 */
	for(i=0;i<=N;i++)
	{
		 x_re[i]=exp(-i);
		 x_im[i]=0;		
	}
	
	/* 复制到输出数组 */
	for(i=0; i<300; i++)
	{
		y_re[i]=x_re[i];
		y_im[i]=x_im[i];
	}

	/* 用雷德算法对输入信号序列进行倒序重排 */
  	j=0;
  	for(i=0;i<N;i++)
  	{	
  		if(i<j)
  		{ 
  		  t_re=y_re[j]; 
  		  t_im=y_im[j]; 
  		  y_re[j]=y_re[i]; 
  		  y_im[j]=y_im[i];
  		  y_re[i]=t_re;
		  y_im[i]=t_im;
  		}
  		k=N/2;
  		while((k<=j)&(k>0))
  		{ 
  		  j=j-k; 
  		  k=k/2;
  		}
  		j=j+k;
    }
    
    /* 计算蝶形运算的级数log2(N) */
	f=N;
	for(m=1; (f=f/2)!=1; m++);
	
    
    /*** FFT ***/
    for(n=1; n<=m; n++)
    {   
    	a=1;				/* a=2的n次方 */
    	for(i=0;i<n;i++)
    		a=a*2;
    	
    	b=a/2;
    	
    	v_re=1.0; 			/* 蝶形因子 */
    	v_im=0.0;
    	w_re=cos(pi/b);		
    	w_im=-sin(pi/b);
    	
    	for(j=0;j<b;j++)	/* 蝶形运算 */
    	{
    		for(i=j;i<N;i=i+a)
    		{
    			c=i+b;
   
       			t_re=y_re[c]*v_re-y_im[c]*v_im;
    			t_im=y_re[c]*v_im+y_im[c]*v_re;
    			
    			y_re[c]=y_re[i]-t_re;
    			y_im[c]=y_im[i]-t_im;
    			
    			y_re[i]=y_re[i]+t_re;
    			y_im[i]=y_im[i]+t_im;
    		}
    		
    		t_re=v_re*w_re-v_im*w_im;
    		t_im=v_re*w_im+v_im*w_re;
    		
    		v_re=t_re;
    		v_im=t_im;
    	}
    }
    
}