上一篇随笔，简要写了一下FFT中数组重新排序的算法。现在把完整的FFT代码分享给大家(有比较详细的注释)。
/*2015年11月10日于河北工业大学*/
#include <complex>
#include <iostream.h>
#include <math.h>
#include <stdlib.h>
const int N=8; //数组的长度
const double PI=3.141592653589793; //圆周率
const double ZERO= 1.0E-14; //当绝对值小于这个数的时候认为是0
typedef std::complex<double> complex;
//函数原型声明
void Reverse(complex src[],int n);
void FFT(complex src[],int n);
void Display(complex f[],int n);
/* 主函数 */
int main()
{
complex src[8]; //原始数组
for(int i=0;i<8;i++)
{
src[i]=complex( double(i), 0.0 );
}
cout<<"原始数组为："<<endl;
Display(src,N);
//一维快速傅里叶变换
FFT(src,N);
cout<<"变换后的数组为："<<endl;
Display(src,N);
return 0;
}
/***********************************************
***** 一维快速离散傅里叶变换 *************
***** 输入：数组长度n *************
***** 虚部 imag *************
***** 实部real *************
************************************************/
void FFT(complex src[],int n)
{
//原始数组按位码倒序排列
Reverse(src,n);
cout<<"倒置后的数组为："<<endl;
Display(src,n);
//计算M，M=log(2,n)表示共有M级蝶形运算
for(int i=n,M=1;(i=i/2)!=1;M++);
for(int L=1;L<=M;L++)
{
// int arrayNum=n/pow(2,L); //被分成的小的数组的个数
int elementNum=pow(2,L); //每个小数组的元素的个数
int uniteachArr=elementNum/2; //每个小数组中所包含的蝶形运算单元的个数
int offset=pow(2,L-1); //每个蝶形运算节点的距离
//cout<<"偏移量为："<<offset<<endl;
for(int i=0;i<n;i=i+elementNum)
{
for(int j=i;j<=i+uniteachArr-1;j++)
{
double r=(j-i)*pow(2,M-L); //W(r,N)中的r
// complex W=complex(cos(2*PI*r/N),sin(2*PI*r/N));
complex temp=src[j+offset]*complex(cos(2*PI*r/double(n)),-sin(2*PI*r/double(n)));
src[j+offset]=src[j]-temp;
src[j]=src[j]+temp;
}
}
}
}
/***********************************************
***** 函数功能：实现原始数组重新排序 *******
***** 主要实现按位码倒序的排列 *****
***** 函数输入：原数组 ****
************************************************/
void Reverse(complex src[],int n)
{
/*
p1和p2表示要进行交换的两个元素的索引
其中p1从第二个元素到倒数第二个元素
扫描，p2则是通过算法确定的与p1进行
交换的另一个元素的索引
*/
int p1;
int p2;
int middle; //表示中心点
int offset; //表示偏移量
/* 中心点和偏移量是计算p2的重要参数 */
p2=n/2; //与p1=1时交换的元素的索引为N/2
for(p1=1;p1<=n-2;/*第一个元素和最后一个元素不用交换*/p1++)
{
/*当p1<p2时交换，避免重复交换*/
if(p1<p2)
{
/*交换*/
complex temp=src[p1];
src[p1]=src[p2];
src[p2]=temp;
}
/*下面开始计算下一个p2*/
middle=n/2; //初始的中心点
offset=p2; //初试的偏移量
//当中心点小于等于偏移量时要更新
while(middle<=offset)
{
offset=offset-middle;
middle=middle/2;
}
p2=middle+offset;//p2就等于新的中心点和偏移量的和
}
}
/*显示数据*/
void Display(complex f[],int n)
{
for(int i=0;i<n;i++)
{
cout.width(9);
cout.setf(ios::fixed);
cout.precision(4);
cout<<f[i].real();
cout.flags(ios::showpos);
cout.width(9);
cout.setf(ios::fixed);
cout.precision(4);
cout<<f[i].imag()<<'i'<<' ';
cout.unsetf(ios::showpos);
if((i+1)%3==0) cout<<endl;
}
cout<<endl;
}