[转：http://xohit.blogchina.com/2228173.html]
实验用的头文件 MYFFT.H
作用：为帮助小虎子做实验，这个头文件提供了完整的一维与二维FFT算法，我想应改是够你折腾了吧！
#include <complex> // complex<float>
using namespace std;
typedef complex<float> Comp; // 复数类型定义
const float _2PI_ = 2.0f * 3.14159265f; // 常数2PI定义
const int MAX_N = 256; // 最大DFT点数
/*----*----*----*----*----*----*----*----*----*----*----*----*
FFT算法模块接口定义
*----*----*----*----*----*----*----*----*----*----*----*----*/
///////////////////////////////////////////
// Function name : BitReverse
// Description : 二进制倒序操作
// Return type : int
// Argument : int src 待倒读的数
// Argument : int size 二进制位数
int BitReverse(int src, int size)
{
int tmp = src;
int des = 0;
for (int i=size-1; i>=0; i--)
{
des = ((tmp & 0x1) << i) | des;
tmp = tmp >> 1;
}
return des;
}
//////////////////////////////////////////////////
// Function name : Reorder
// Description : 数据二进制整序
// Return type : void
// Argument : Comp x[MAX_N] 待整序数组
// Argument : int N FFT点数
// Argument : int M 点数的2的幂次
void Reorder(Comp x[MAX_N], int N, int M)
{
Comp new_x[MAX_N];
for (int i=0; i<N; i++)
new_x = x[BitReverse(i, M)];
// 重新存入原数据中(已经是二进制整序过了的数据)
for (i=0; i<N; i++)
x = new_x;
}
//////////////////////////////////////////////////
// Function name : CalcW
// Description : 计算旋转因子
// Return type : void
// Argument : Comp W[MAX_N] 存放因子的数组
// Argument : int N FFT的点数
// Argument : int flag 正逆变换标记
void CalcW(Comp W[MAX_N], int N, int flag)
{
for (int i=0; i<N/2; i++)
{
if (flag == 1)
W = exp(Comp(0, -_2PI_ * i / N)); // 正FFT变换
else
W = exp(Comp(0, _2PI_ * i / N)); // 逆FFT变换
}
}
/////////////////////////////////////////////////
// Function name : FFT_1D_Kernel
// Description : FFT算法核心
// Return type : void
// Argument : Comp* x 数据
// Argument : int M 幂次
// Argument : int flag 正逆变换标记
以下本应由自己完成。
void FFT_1D(Comp* x, int M, int flag)
{
int N = (1 << M);
// 二进制整序
Reorder(x, N, M);
// 旋转因子计算
Comp W[MAX_N];
CalcW(W, N, flag);
// 级内群数
int GroupNum = N/2; // 第一级的群数为N/2
// 群内蝶形单元数
int CellNum = 1; // 第一级的群内蝶形单元数为1
// 处理各级
for (int i=0; i<M; i++)
{
// 处理各群
for (int j=0; j<GroupNum; j++)
{
// 处理各蝶形单元
for (int k=0; k<CellNum; k++)
{
// (1) 计算出当前蝶形单元所含元素在数据数组中的位置
// 第一元素位置
int Pos1 = CellNum * j * 2 + k ; // 已经处理了前 j 群，每群有 CellNum 个单元，
每单元有 2 个元素
// 第二元素位置
int Pos2 = Pos1 + CellNum;
// (2) 计算旋转因子与单元的第二元素的复数乘积
Comp TMP = x[Pos2] * W[k * GroupNum] ;
// (3) 计算最终结果, 并存入到数组的原先位置
x[Pos2] = x[Pos1] - TMP ;
x[Pos1] = x[Pos1] + TMP ;
}
}
GroupNum >>= 1; // 级别增加, 则相应的群数减少一半
CellNum <<= 1; // 级别增加, 则相应的群内单元数增加一倍
}
// 如果是IFFT，各元素还要再除以N
if (flag != 1)
{
for (i=0; i<N; i++)
x /= N;
}
}
//////////////////////////////////////////////////////
// Function name : FFT_2D_Kernel
// Description : 2D FFT核心算法
// Return type : void
// Argument : Comp x[MAX_N][MAX_N] 二维数据
// Argument : int M 幂次
// Argument : int flag 正逆变换标记
以下本应由自己完成。
void FFT_2D(Comp x[MAX_N][MAX_N], int M, int flag)
{
int N = (1 << M);
// 先逐行进行 1D-FFT
for (int i=0; i<N; i++)
FFT_1D(x, M, flag); // <--- 计算结果再存入矩阵x中
// 再逐列进行 1D-FFT
Comp col[MAX_N];
for (int j=0; j<N; j++)
{
// 取得第j列的数据
for (int i=0; i<N; i++)
col = x[j];
// 对第j列数据进行 1D-FFT
FFT_1D(col, M, flag); // <--- 计算结果在数组col中
// 将结果放回矩阵第j列中
for (i=0; i<N; i++)
x[j] = col;
}
}
// <--- End of [FFT.H]