基于C++任意点数的FFT／IFFT(时域和频域)实现

函数说明:更改主函数体中的N和length（=log2(N)）既可以实现任意点数(2的幂次)的FFT/ IFFT的实现，fft函数中flag标志位控制是正变换还是逆变换。

1.复数操作类

定义复数类，重载复数四则运算符号，重载输出运算符，重载赋值运算符。

/**********预编译文件头文件complex.h********/
#include"iostream"
using namespace std;
class complex
{
	double real,image;
public:
	complex(){real=0;image=0;}
	complex(float i,float j){real=i;image=j;}
	double getr(){return real;}
	double geti(){return image;}
	void show()
	{
		if(image>=0)
			cout<<real<<"+"<<image<<"j"<<"   ";
		else
			cout<<real<<image<<"j"<<"   ";
	}
	void setvalue(double i=0,double j=0)
	{
		real=abs(i)<0.0001?0:i;
		image=abs(j)<0.0001?0:j;
	}
	complex operator +(complex&);
	complex operator-(complex&);
	complex operator*(complex&);
	complex operator /(int n);
	void operator+=(complex&);
	void operator =(complex&);
	friend complex W(int m,int n, int flag);
};
complex complex::operator +(complex& c)
	{
		complex t;
		t.real=real+c.real;
		t.image=image+c.image;
		return t;
	}
complex complex::operator*(complex& c)
	{
		complex t;
		t.real=real*c.real-image*c.image ;
		t.image=real*c.image+image*c.real;
		return t;

	}
complex complex::operator/ (int n)
	{
		complex temp;
		temp.real=real/n;
		temp.image=image/n;
		return temp;
	 }
complex complex::operator -(complex& c)
	{
		complex t;
		t.real=real-c.real;
		t.image=image-c.image;
		return t;
	}
void complex::operator+=(complex&c)
	{
		real=real+c.real;
		image=image+c.image;
	}
void complex::operator =(complex&c)
	{
		real=abs(c.real)<0.00001?0:(c.real);
		image=abs(c.image)<0.00001?0:(c.image);
	}

2.主函数

/*******************主函数体***********************/
#include "stdafx.h"
#include"iostream"
#include"cmath"
#include"complex.h"
#define pi 3.141592657
#define N 16//需要运算的fft的点数
#define length 4//存储需要的二进制位,即length=log2(N)
using namespace std;
/*********************重新排列输入序列函数****************************/
void bit_reversal(complex a[],int n)
{
	int bitdesicion(unsigned);//用来重现排序原来的输入信号序列应该对应变换的信号的序号
	complex temp;
	int j;
	for(int i=0;i<n;i++)
	{
		j=(int)bitdesicion(i);
		if(i<j){temp=a[i];a[i]=a[j];a[j]=temp;}
	}
}
int bitdesicion(unsigned n)
{
    int j=0;
	double k;
	double result=0;
	int bit[length]={0};
	while(n){bit[j]=n%2;n=n/2;j++;}//将n分解成进制码存储于int指针类型bit中
	j=0;
	while(j<length)
	{k=length-1-j;result=result+bit[j]*pow(2.0,k);j++;}
     return (int)result;
}

/*********************显示信号数组的元素的函数*********************************/
void play(complex a[],int len)
{
	for(int i=0;i<len;i++)
	{
		a[i].show();
		if((i+1)%4==0) cout<<"  "<<endl;
	}
   cout<<'\n'<<endl;
}

/************获得加权系数W,flag标记为正变换或者逆变换系数*******************/
void getW(complex w[],int len,int flag)
{
	if(flag)
	for(int i=0 ; i<len ; i++ )
			w[i].setvalue(cos(pi*2*i/len),-sin(pi*2*i/len));
	else
	for(int j=0 ; j<len ; j++ )
			w[j].setvalue(cos(pi*2*j/len),sin(pi*2*j/len) );
}
/****************************基2时域FFT算法***********************************/
void fft_temporal(complex input[],int len,int flag)
{
	int i,j;int L,k;
	void bit_reversal(complex a[],int);
	bit_reversal(input,len);//len 代表输入数据的长度
	cout<<"重新排序后输入信号为:"<<endl;play(input,len);
	complex* w=new complex[len];getW(w,len,flag);//获得快速傅里叶的系数,flag为标志位，代表获得正变换系数，代表获得逆变换系数
	if(flag)
	    cout<<"FFT系数如下"<<endl;
	else 
		cout<<"IFFT系数如下"<<endl;
    play(w,len);
	complex* tempdata=new complex[len];
	complex* tempw=new complex[len/2];

/*********************最核心--蝶形算法**************************************/
for(i=1;i<=length;i++)//i代表第i级蝶形,length表示log2(N)
	{		
		L=(int)pow(2.0,i);
		for(j=0;j<len;j+=L)//L为相应单个蝶形运算包含的节点数也是每组蝶形运算的间隔数
		{ 
			for(k=0;k<L/2;k++)  
				tempw[k]=w[k*len/L]*input[j+L/2+k]; //傅里叶加权系数和对应输入信号的乘积事先准备好
			for(k=0;k<L/2;k++)//L/2将每一个蝶形算法依据是加还是减够成两个部分，每个部分是蝶形点间隔数的一半
			{
				tempdata[j+k]=input[j+k]+tempw[k];
			    tempdata[j+L/2+k]=input[j+k]-tempw[k];
			} 
		}
	       for(j=0;j<len;j++) input[j]=tempdata[j];
	 }
 	if(!flag) for(j=0;j<len;j++) input[j]=tempdata[j]/len;//逆变换必须乘加权系数/N才能得到正确的结果
        delete[]w;
        delete[]tempw;
	delete[]tempdata;
}

/**************************************基2频域FFT算法***********************************/
void fft_frequency(complex input[],int len,int flag)
{
	int i,j;int L,k;
	void bit_reversal(complex a[],int);
	complex* w=new complex[len];getW(w,len,flag);//获得快速傅里叶的系数,flag为标志位，代表获得正变换系数，代表逆变换
	complex* tempdata=new complex[len];
	if(flag)
	    cout<<"FFT系数如下"<<endl;
	else 
		cout<<"IFFT系数如下"<<endl;
    play(w,len);
/********************************最核心--蝶形算法*****************************************/
for(i=1;i<=length;i++)//i代表第i级蝶形,length表示log2(N)
	{		
		L=(int)pow(2.0,length+1-i);
		for(j=0;j<len;j+=L)//L为相应单个蝶形运算包含的节点数也是每组蝶形运算的间隔数
		{ 
			for(k=0;k<L/2;k++)//L/2将每一个蝶形算法依据是加还是减够成两个部分，每个部分是蝶形点间隔数的一半
			{
				tempdata[j+k]=input[j+k]+input[j+L/2+k];
			    tempdata[j+L/2+k]=w[k*len/L]*(input[j+k]-input[j+L/2+k]);
			}
		    for(k=0;k<len;k++) input[k]=tempdata[k];
		}
	 }
 	if(!flag) for(j=0;j<len;j++) input[j]=input[j]/len;//逆变换必须乘加权系数/N才能得到正确的结果
	bit_reversal(input,len);//len 代表输入数据的长度，频率的FFT在最后要重新排序
    delete[]w;delete[]tempdata;
}

/************************主函数********************************************/
int _tmain(int argc, _TCHAR* argv[])
{
/***********************函数声明*****************************/
   void play(complex a[],int);//显示整个数组<pre name="code" class="cpp">   void getW(complex w[],int len,int flag);//获得正变换(flag=1)或者逆变换系数(flag=0)
   void fft_temporal(complex input[],int len,int flag);//基时间的fft,flag控制为正变换或者逆变换
   void fft_frequency(complex input[],int len,int flag);//基频率的fft
/************************************************************/
   complex input
;

complex output
;
   for(int i=0;i<N;i++)  input[i].setvalue(i,0);
   cout<<"输入信号为:"<<endl;play(input,N);
   fft_temporal(input,N,1);//"1控制为正变换"
   cout<<"基2时域FFT输出信号为:"<<endl;play(input,N);
   for(int j=0;j<N;j++) 

output[j]=input[j];
   fft_temporal(output,N,0);//"0"控制为逆变换
   cout<<"基2时域IFFT输出信号为:"<<endl;play(output,N);
   fft_frequency(output,N,1);
   cout<<"基2频率FFT输出信号为:"<<endl;play(output,N);
   fft_frequency(output,N,0);
   cout<<"基2频率IFFT输出信号为:"<<endl;play(output,N);
   return 0;
}