C语言实现FFT or IFFT

matlab实现

close all;
clear all;
clc;

count=2;
y=[1+2i,2+3i];
y_FFT=fft(y,count);  % 傅里叶变换
P=abs(y_FFT)/count;  %  幅度谱
u=max(P);            % 得到波束1接收到的信号的幅度值

仿真结果

C语言实现

/************FFT***********/
#include   <stdio.h>
#include   <math.h>
#include   <stdlib.h>

#define PI 3.1415926535897931
#define   N   1000
  
typedef   struct//结构体
 {
	double   real;
	double   img;
 }complex;
  
complex   x[N],  *W;//输出序列的值
int   size_x=0;     //输入序列的长度，只限2的N次方


/*****************************
            子函数
*****************************/
  void   fft(); /*快速傅里叶变换*/
  void   ifft(); /*快速傅里叶逆变换*/
  void   initW();
  void   change();
  void   add(complex   ,complex   ,complex   *);   /*复数加法*/   
  void   mul(complex   ,complex   ,complex   *);   /*复数乘法*/   
  void   sub(complex   ,complex   ,complex   *);   /*复数减法*/   
  void   divi(complex   ,complex   ,complex   *);/*复数除法*/   
  void   output(); /*输出结果*/




/*****************************
             主函数
*****************************/
int   main()
 {
	int i,method;
	printf("Please   input   the   size   of   x:\n");/*输入序列的长度*/
	scanf("%d",&size_x);
	printf("Please   input   the   data   in   x[N]:(such as:5 6)\n"); /*输入序列对应的值*/
	for(i=0;i<size_x;i++)
	{
		scanf("%lf %lf",&x[i].real,&x[i].img);
		initW();   
	}
	printf("Use   FFT(0)   or   IFFT(1)?\n");   /*选择FFT或逆FFT运算*/
	scanf("%d",&method);   
	if(method==0)   
		fft();   
	else   
		ifft();   
	output();   
	return 0;   
 }   
 
  
/*****************************
       基-2 FFT运算
*****************************/
void   fft()   
{   
  int   i=0,j=0,k=0,l=0;   
  complex   up,down,product;   
  change();   
  for(i=0;i<log(size_x)/log(2) ;i++)
  {       
	l=1<<i;   
	for(j=0;j<size_x;j+=2*l)
	{                           
		for(k=0;k<l;k++)
		{                 
			mul(x[j+k+l],W[size_x*k/2/l],&product);   
			add(x[j+k],product,&up);   
			sub(x[j+k],product,&down);   
			x[j+k]=up;   
			x[j+k+l]=down;   
		}   
	}   
  }   
}   
    
/*****************************
         IFFT运算
*****************************/
void   ifft()   
{   
  int   i=0,j=0,k=0,l=size_x;   
  complex   up,down;   
  for(i=0;i<   (int)(   log(size_x)/log(2)   );i++) /*蝶形运算*/  
  {     
	l/=2;   
	for(j=0;j<size_x;j+=   2*l   )  
	{                         
		for(k=0;k<l;k++)
		{               
			add(x[j+k],x[j+k+l],&up);   
			up.real/=2;up.img/=2;   
			sub(x[j+k],x[j+k+l],&down);   
			down.real/=2;down.img/=2;   
			divi(down,W[size_x*k/2/l],&down);   
			x[j+k]=up;   
			x[j+k+l]=down;   
		}   
	}   
  }   
  change();   
}   
    
/*****************************
         
*****************************/
void   initW()   
 {   
	int   i;   
	W=(complex   *)malloc(sizeof(complex)   *   size_x);   
	for(i=0;i<size_x;i++)   
	{   
		W[i].real=cos(2*PI/size_x*i);   
		W[i].img=-1*sin(2*PI/size_x*i);   
	}   
 }   
 
   
/*****************************
         
*****************************/
void   change()   
{   
	complex   temp;   
	unsigned   short   i=0,j=0,k=0;   
	double   t;   
	for(i=0;i<size_x;i++)   
	{   
		k=i;j=0;   
		t=(log(size_x)/log(2));   
		while(   (t--)>0   )   
		{   
			j=j<<1;   
			j|=(k   &   1);   
			k=k>>1;   
		}   
		if(j>i)   
		{   
			temp=x[i];   
			x[i]=x[j];   
			x[j]=temp;   
		}   
	}   
}   
    
/*****************************
         输出结果
*****************************/
void   output()   
{   
	int   i;   
	printf("The   result   are   as   follows\n");   
	for(i=0;i<size_x;i++)   
	{   
		printf("%.4f",x[i].real);   
		if(x[i].img>=0.0001)   
			printf("+%.4fj\n",x[i].img);   
		else   if(fabs(x[i].img)<0.0001)   
			printf("\n");   
		else     
			printf("%.4fj\n",x[i].img);   
	}   
} 


/*****************************
         复数加法
*****************************/  
void   add(complex   a,complex   b,complex   *c)   
{   
	c->real=a.real+b.real;   
	c->img=a.img+b.img;   
 }   
  
/*****************************
         复数乘法
*****************************/   
void   mul(complex   a,complex   b,complex   *c)   
{   
	c->real=a.real*b.real   -   a.img*b.img;   
	c->img=a.real*b.img   +   a.img*b.real;   
}
  
/*****************************
         复数减法
*****************************/ 
void   sub(complex   a,complex   b,complex   *c)   
{   
	c->real=a.real-b.real;   
	c->img=a.img-b.img;   
}
  
/*****************************
         复数除法
*****************************/ 
void   divi(complex   a,complex   b,complex   *c)   
{   
	c->real=(   a.real*b.real+a.img*b.img   )/( b.real*b.real+b.img*b.img);    
	c->img=(   a.img*b.real-a.real*b.img)/(b.real*b.real+b.img*b.img);   
}

仿真结果