基于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);
}
/*******************主函数体***********************/
#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[N];complex output[N]; 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; }
