fft.rar_ c++ FFT_FFT 单片机_fft_二维 傅立叶_二维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]