Cooley-Tukey FFT算法的非递归实现
一维情况
#include <iostream> #include <complex> #include <cmath> const double PI = 3.14159265358979323846; // 交换位置 void swap(std::complex<double>& a, std::complex<double>& b) { std::complex<double> temp = a; a = b; b = temp; } // 频率倒位 void bitReverse(std::complex<double>* point, int N) { int j = 0; for (int i = 1; i < N - 1; ++i) { int k = N / 2; while (j >= k) { j -= k; k /= 2; } j += k; if (i < j) { swap(point[i], point[j]); } } } // Cooley-Tukey算法的非递归实现 void FFT(std::complex<double>* point, int N) { bitReverse(point, N); for (int s = 2; s <= N; s *= 2) { // 蝶形运算的Stage循环 double angle = -2 * PI / s; std::complex<double> w(1, 0); std::complex<double> wn(cos(angle), sin(angle)); for (int k = 0; k < N; k += s) { // 一组蝶形运算group，注：每个group的蝶形因子乘数不同 w = std::complex<double>(1, 0); for (int j = 0; j < s / 2; ++j) { // 一个蝶形运算 注：group内的蝶形运算 std::complex<double> u = point[k + j]; std::complex<double> v = point[k + j + s / 2] * w; point[k + j] = u + v; point[k + j + s / 2] = u - v; w *= wn; } } } } // FFT逆变换 void IFFT(std::complex<double>* point, int N) { // 共轭变换(即：实部不变，虚部取相反数) for (int i = 0; i < N; ++i) { point[i] = std::conj(point[i]); } // 应用FFT FFT(point, N); // 共轭变换(即：实部不变，虚部取相反数)并归一化 for (int i = 0; i < N; ++i) { point[i] = std::conj(point[i]); point[i] /= N; } } // FFT变换测试 void FFTTest() { const int N = 8; // 输入序列的长度 std::complex<double> input[N] = { 0, 1, 2, 3, 4, 5, 6, 7 }; // 输入序列 // 将输入序列扩展为2的幂次方长度 int paddedLength = (int)pow(2, ceil(log2(N))); std::complex<double>* paddedInput = new std::complex<double>[paddedLength]; for (int i = 0; i < N; ++i) { paddedInput[i] = input[i]; } // 补零到2的n次幂 for (int i = N; i < paddedLength; ++i) { paddedInput[i] = 0; } // 执行FFT变换 FFT(paddedInput, paddedLength); std::cout << "========== FFT变换 ========== " << std::endl; // 输出原始数据 std::cout << "原始数据:----------" << std::endl; for (int k = 0; k < N; ++k) { std::cout << "f[" << k << "] = " << input[k] << std::endl; } // 输出变换结果 std::cout << "变换结果:----------" << std::endl; for (int k = 0; k < paddedLength; ++k) { std::cout << "F[" << k << "] = " << paddedInput[k] << std::endl; } // 输出幅值结果 std::cout << "幅值结果:----------" << std::endl; for (int k = 0; k < paddedLength; ++k) { std::cout << "Length[" << k << "] = " << std::abs(paddedInput[k]) << std::endl; // 计算复数的模 } delete[] paddedInput; } // FFT逆变换测试 void IFFTTest() { const int N = 8; // 输入序列的长度 std::complex<double> input[N] = { {28,0}, {-4,9.65685}, {-4,4}, {-4,1.65685}, {-4,0}, {-4,-1.65685}, {-4,-4}, {-4,-9.65685} }; // 输入序列 // 将输入序列扩展为2的幂次方长度 int paddedLength = (int)pow(2, ceil(log2(N))); std::complex<double>* paddedInput = new std::complex<double>[paddedLength]; for (int i = 0; i < N; ++i) { paddedInput[i] = input[i]; } // 补零到2的n次幂 for (int i = N; i < paddedLength; ++i) { paddedInput[i] = 0; } // 执行FFT逆变换 IFFT(paddedInput, paddedLength); std::cout << "========== FFT逆变换 ========== " << std::endl; // 输出原始数据 std::cout << "原始数据:----------" << std::endl; for (int k = 0; k < N; ++k) { std::cout << "F[" << k << "] = " << input[k] << std::endl; } // 输出变换结果 std::cout << "变换结果:----------" << std::endl; for (int k = 0; k < paddedLength; ++k) { std::cout << "f[" << k << "] = " << paddedInput[k] << std::endl; } delete[] paddedInput; } int main() { // 执行FFT变换测试 FFTTest(); std::cout << std::endl; // 执行FFT逆变换测试 IFFTTest(); return 0; }
输出结果如下：
========== FFT变换 ========== 原始数据:---------- f[0] = (0,0) f[1] = (1,0) f[2] = (2,0) f[3] = (3,0) f[4] = (4,0) f[5] = (5,0) f[6] = (6,0) f[7] = (7,0) 变换结果:---------- F[0] = (28,0) F[1] = (-4,9.65685) F[2] = (-4,4) F[3] = (-4,1.65685) F[4] = (-4,0) F[5] = (-4,-1.65685) F[6] = (-4,-4) F[7] = (-4,-9.65685) 幅值结果:---------- Length[0] = 28 Length[1] = 10.4525 Length[2] = 5.65685 Length[3] = 4.32957 Length[4] = 4 Length[5] = 4.32957 Length[6] = 5.65685 Length[7] = 10.4525 ========== FFT逆变换 ========== 原始数据:---------- F[0] = (28,0) F[1] = (-4,9.65685) F[2] = (-4,4) F[3] = (-4,1.65685) F[4] = (-4,0) F[5] = (-4,-1.65685) F[6] = (-4,-4) F[7] = (-4,-9.65685) 变换结果:---------- f[0] = (0,-0) f[1] = (1,1.72255e-16) f[2] = (2,4.44089e-16) f[3] = (3,3.82857e-16) f[4] = (4,-0) f[5] = (5,-4.979e-17) f[6] = (6,-4.44089e-16) f[7] = (7,-5.05322e-16)
二维情况
二维FFT是在一维FFT基础上实现，实现过程为：
1.对二维输入数据的每一行进行FFT，变换结果仍然按行存入二维数组中。
2.在1的结果基础上，对每一列进行FFT，再存入原来数组，即得到二维FFT结果。
#include <iostream> #include <complex> #include <cmath> const double PI = 3.14159265358979323846; // 交换位置 void swap(std::complex<double>& a, std::complex<double>& b) { std::complex<double> temp = a; a = b; b = temp; } // 频率倒位 void bitReverse(std::complex<double>* point, int N) { int j = 0; for (int i = 1; i < N - 1; ++i) { int k = N / 2; while (j >= k) { j -= k; k /= 2; } j += k; if (i < j) { swap(point[i], point[j]); } } } // 一维Cooley-Tukey算法的非递归实现 void FFT(std::complex<double>* point, int N) { bitReverse(point, N); for (int s = 2; s <= N; s *= 2) { double angle = -2 * PI / s; std::complex<double> w(1, 0); std::complex<double> wn(cos(angle), sin(angle)); for (int k = 0; k < N; k += s) { w = std::complex<double>(1, 0); for (int j = 0; j < s / 2; ++j) { std::complex<double> u = point[k + j]; std::complex<double> v = point[k + j + s / 2] * w; point[k + j] = u + v; point[k + j + s / 2] = u - v; w *= wn; } } } } // 一维FFT逆变换 void IFFT(std::complex<double>* point, int N) { // 共轭变换(即：实部不变，虚部取相反数) for (int i = 0; i < N; ++i) { point[i] = std::conj(point[i]); } // 应用FFT FFT(point, N); // 共轭变换(即：实部不变，虚部取相反数)并归一化 for (int i = 0; i < N; ++i) { point[i] = std::conj(point[i]); point[i] /= N; } } // 二维Cooley-Tukey算法的非递归实现 void FFT2D(std::complex<double>* point, int rows, int cols) { // 对每一行进行FFT变换 for (int i = 0; i < rows; ++i) { FFT(&point[i * cols], cols); } // 对每一列进行FFT变换 std::complex<double>* column = new std::complex<double>[rows]; for (int j = 0; j < cols; ++j) { for (int i = 0; i < rows; ++i) { column[i] = point[i * cols + j]; } FFT(column, rows); for (int i = 0; i < rows; ++i) { point[i * cols + j] = column[i]; } } delete[] column; } // 二维FFT逆变换 void IFFT2D(std::complex<double>* point, int rows, int cols) { // 对每一行进行逆变换 for (int i = 0; i < rows; ++i) { IFFT(&point[i * cols], cols); } // 对每一列进行逆变换 std::complex<double>* column = new std::complex<double>[rows]; for (int j = 0; j < cols; ++j) { for (int i = 0; i < rows; ++i) { column[i] = point[i * cols + j]; } IFFT(column, rows); for (int i = 0; i < rows; ++i) { point[i * cols + j] = column[i]; } } delete[] column; } // FFT变换测试 void FFTTest() { const int rows = 4; // 输入矩阵的行数 const int cols = 4; // 输入矩阵的列数 std::complex<double> input[rows * cols] = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16 }; // 输入矩阵 // 将输入矩阵的行数和列数扩展为2的幂次方长度 int paddedRows = (int)pow(2, ceil(log2(rows))); int paddedCols = (int)pow(2, ceil(log2(cols))); // 将长宽补零到2的n次幂 std::complex<double>* paddedInput = new std::complex<double>[paddedRows * paddedCols]; for (int i = 0; i < rows; ++i) { for (int j = 0; j < cols; ++j) { paddedInput[i * paddedCols + j] = input[i * cols + j]; } for (int j = cols; j < paddedCols; ++j) { paddedInput[i * paddedCols + j] = 0; } } for (int i = rows; i < paddedRows; ++i) { for (int j = 0; j < paddedCols; ++j) { paddedInput[i * paddedCols + j] = 0; } } // 执行2D FFT变换 FFT2D(paddedInput, paddedRows, paddedCols); std::cout << "========== FFT变换 ========== " << std::endl; // 输出原始数据 std::cout << "原始数据:----------" << std::endl; for (int i = 0; i < rows; ++i) { for (int j = 0; j < cols; ++j) { std::cout << "f[" << i << "][" << j << "] = " << input[i * cols + j] << std::endl; } } // 输出变换结果 std::cout << "变换结果:----------" << std::endl; for (int i = 0; i < paddedRows; ++i) { for (int j = 0; j < paddedCols; ++j) { std::cout << "F[" << i << "][" << j << "] = " << paddedInput[i * paddedCols + j] << std::endl; } } // 输出幅值结果 std::cout << "幅值结果:----------" << std::endl; for (int i = 0; i < paddedRows; ++i) { for (int j = 0; j < paddedCols; ++j) { std::cout << "Length[" << i << "][" << j << "] = " << std::abs(paddedInput[i * paddedCols + j]) << std::endl; // 计算复数的模 } } delete[] paddedInput; } // FFT逆变换测试 void IFFTTest() { const int rows = 4; // 输入矩阵的行数 const int cols = 4; // 输入矩阵的列数 std::complex<double> input[rows * cols] = { {136,0}, {-8,8}, {-8,0}, {-8,-8}, {-32,32}, {0,0}, {0,0}, {0,0}, {-32,0}, {0,0}, {0,0}, {0,0}, {-32,-32}, {0,0}, {0,0}, {0,0} }; // 输入矩阵 // 将输入矩阵的行数和列数扩展为2的幂次方长度 int paddedRows = (int)pow(2, ceil(log2(rows))); int paddedCols = (int)pow(2, ceil(log2(cols))); // 将长宽补零到2的n次幂 std::complex<double>* paddedInput = new std::complex<double>[paddedRows * paddedCols]; for (int i = 0; i < rows; ++i) { for (int j = 0; j < cols; ++j) { paddedInput[i * paddedCols + j] = input[i * cols + j]; } for (int j = cols; j < paddedCols; ++j) { paddedInput[i * paddedCols + j] = 0; } } for (int i = rows; i < paddedRows; ++i) { for (int j = 0; j < paddedCols; ++j) { paddedInput[i * paddedCols + j] = 0; } } // 执行2D FFT变换 IFFT2D(paddedInput, paddedRows, paddedCols); std::cout << "========== 执行FFT逆变换 ========== " << std::endl; // 输出原始数据 std::cout << "原始数据:----------" << std::endl; for (int i = 0; i < rows; ++i) { for (int j = 0; j < cols; ++j) { std::cout << "F[" << i << "][" << j << "] = " << input[i * cols + j] << std::endl; } } // 输出变换结果 std::cout << "变换结果:----------" << std::endl; for (int i = 0; i < paddedRows; ++i) { for (int j = 0; j < paddedCols; ++j) { std::cout << "f[" << i << "][" << j << "] = " << paddedInput[i * paddedCols + j] << std::endl; } } delete[] paddedInput; } int main() { // 执行FFT变换测试 FFTTest(); std::cout << std::endl; // 执行FFT逆变换测试 IFFTTest(); return 0; }
输出结果如下：
========== FFT变换 ========== 原始数据:---------- f[0][0] = (1,0) f[0][1] = (2,0) f[0][2] = (3,0) f[0][3] = (4,0) f[1][0] = (5,0) f[1][1] = (6,0) f[1][2] = (7,0) f[1][3] = (8,0) f[2][0] = (9,0) f[2][1] = (10,0) f[2][2] = (11,0) f[2][3] = (12,0) f[3][0] = (13,0) f[3][1] = (14,0) f[3][2] = (15,0) f[3][3] = (16,0) 变换结果:---------- F[0][0] = (136,0) F[0][1] = (-8,8) F[0][2] = (-8,0) F[0][3] = (-8,-8) F[1][0] = (-32,32) F[1][1] = (0,0) F[1][2] = (0,0) F[1][3] = (0,0) F[2][0] = (-32,0) F[2][1] = (0,0) F[2][2] = (0,0) F[2][3] = (0,0) F[3][0] = (-32,-32) F[3][1] = (0,0) F[3][2] = (0,0) F[3][3] = (0,0) 幅值结果:---------- Length[0][0] = 136 Length[0][1] = 11.3137 Length[0][2] = 8 Length[0][3] = 11.3137 Length[1][0] = 45.2548 Length[1][1] = 0 Length[1][2] = 0 Length[1][3] = 0 Length[2][0] = 32 Length[2][1] = 0 Length[2][2] = 0 Length[2][3] = 0 Length[3][0] = 45.2548 Length[3][1] = 0 Length[3][2] = 0 Length[3][3] = 0 ========== 执行FFT逆变换 ========== 原始数据:---------- F[0][0] = (136,0) F[0][1] = (-8,8) F[0][2] = (-8,0) F[0][3] = (-8,-8) F[1][0] = (-32,32) F[1][1] = (0,0) F[1][2] = (0,0) F[1][3] = (0,0) F[2][0] = (-32,0) F[2][1] = (0,0) F[2][2] = (0,0) F[2][3] = (0,0) F[3][0] = (-32,-32) F[3][1] = (0,0) F[3][2] = (0,0) F[3][3] = (0,0) 变换结果:---------- f[0][0] = (1,-0) f[0][1] = (2,6.12323e-17) f[0][2] = (3,-0) f[0][3] = (4,-6.12323e-17) f[1][0] = (5,2.44929e-16) f[1][1] = (6,3.06162e-16) f[1][2] = (7,2.44929e-16) f[1][3] = (8,1.83697e-16) f[2][0] = (9,-0) f[2][1] = (10,6.12323e-17) f[2][2] = (11,-0) f[2][3] = (12,-6.12323e-17) f[3][0] = (13,-2.44929e-16) f[3][1] = (14,-1.83697e-16) f[3][2] = (15,-2.44929e-16) f[3][3] = (16,-3.06162e-16)
扩展阅读
