fft.cpp

序执行时要求输入序列的长度N_pre

#include "stdio.h"
#include "stdlib.h"
#include "math.h"
void main()
{   
#define  pi  3.14159
	int N_pre;
	int N;
	int M;  //级数
	int G;  //组数
	int i,j,k,mid;
	float *X_re,*X_im;
	float exchange_re,exchange_im;
	float W_re,W_im;
	float temp_p_re,temp_p_im;
    float temp_q_re,temp_q_im;
	float Cmidmult_re,Cmidmult_im;

	printf("input the size of input array:\n");
	scanf("%d",&N_pre);
	printf("N_pre=%d\n",N_pre); 
	
    for(N=1,i=0;i<N_pre;i++)
	{
		N=N<<1;
		if(N_pre>N)  continue;
		else  break;
	}
        
	M=i+1;  //M是FFT级数 2^M=N;  N是补零调整满足2的整数幂的最小的数

  
 printf("N=%d,M=%d\n",N,M);

 X_re=(float *)calloc(sizeof(float),N);
 X_im=(float *)calloc(sizeof(float),N);

 /*****************输入x的值，并判断补零**************/
 for(i=0;i<N_pre;i++)
 {
   printf("input the value of X_re[%d]:\n",i);
   scanf("%f",&X_re[i]);
   printf("input the value of X_im[%d]:\n",i);
   scanf("%f",&X_im[i]);
 } 

  for(;i<N;i++)
  {
	  X_re[i]=0;
      X_im[i]=0;
  }


/**************************换位操作************************/
for(j=i=0;i<N-1;i++)
{
	if(i<j)
	{
	  exchange_re=X_re[i];
	  exchange_im=X_im[i];
      X_re[i]=X_re[j];
      X_im[i]=X_im[j];
      X_re[j]=exchange_re;
	  X_im[j]=exchange_im;
	}
	 mid=N>>1;
	 while(mid<=j)
	 {
		 j=j-mid;
	     mid=mid>>1;
	 }
	  j=j+mid;
}



/*******************变换开始****************************/

/***************级（K）循环开始*************************/ 


for(k=0;k<M;k++)
{
	G=N>>(k+1);  //组数确定

	/*****************组（G ）循环开始**************************/
	for(j=0;j<G;j++)  
	{ 
		for(i=0;i<(1<<k);i++)//每组中有2^k个蝶形运算
		{
          W_re=(float)cos((2*pi*i)/(2<<k)); //除以2^k+1;
		  W_im=(float)(0.0-sin((2*pi*i)/(2<<k))); //除以2^k+1;
             

          temp_p_re=X_re[j*N/G+i];
		  temp_p_im=X_im[j*N/G+i];
		 
		  temp_q_re=X_re[j*N/G+i+(1<<k)];
		  temp_q_im=X_im[j*N/G+i+(1<<k)];
          
		  Cmidmult_re=W_re*temp_q_re-W_im*temp_q_im;
		  Cmidmult_im=W_re*temp_q_im+W_im*temp_q_re;

          X_re[j*N/G+i]=temp_p_re+Cmidmult_re;
          X_im[j*N/G+i]=temp_p_im+Cmidmult_im;

          X_re[j*N/G+i+(1<<k)]=temp_p_re-Cmidmult_re;
		  X_im[j*N/G+i+(1<<k)]=temp_p_im-Cmidmult_im;
         
		}
		
	}
   
} 
/*********************输出FFT结果***********************/	  


for(i=0;i<N;i++)
{   if(X_im[i]>=0.0)
	printf(" the FFT value of Xk[%d] is %f +%f i:\n",i,X_re[i],X_im[i]);	
    else
    printf(" the FFT value of Xk[%d] is %f %f i:\n",i,X_re[i],X_im[i]);
//  printf(" the FFT value of Xk[%d] is %f:\n",i,);
}
 	  
	  
	  free(X_re);
      free(X_im);
}