当前位置：首页 > java >正文

FFT实现（Cooley-Tukey算法）

java 2025/5/5 6:45:17

傅里叶变换理论基础

傅里叶变换是一种将信号从时域转换到频域的数学工具。对于一个离散的序列 (x[n])（(n = 0,1,\cdots,N - 1)），其离散傅里叶变换（DFT）定义为：

[X[k]=\sum_{n = 0}^{N - 1}x[n]e^{-j\frac{2\pi}{N}kn}, \quad k = 0,1,\cdots,N - 1]

其中，(X[k]) 是频域中的第 (k) 个分量，(x[n]) 是时域中的第 (n) 个样本，(N) 是序列的长度，(j=\sqrt{-1})。

直接计算 DFT 的时间复杂度是 (O(N^2))，当 (N) 很大时，计算量非常大。而 Cooley - Tukey 算法是一种快速傅里叶变换（FFT）算法，它通过利用旋转因子 (W_N^{kn}=e{-j\frac{2\pi}{N}kn}) 的周期性和对称性，将计算复杂度降低到 (O(N\log N))。

代码与傅里叶变换公式的结合分析

1. 函数定义与变量初始化

vector<complex<float>> FFT(const vector<float>& frame) {const int N = frame.size();vector<complex<float>> x(N);

FFT 函数接收一个长度为 (N) 的时域信号 frame，并创建一个长度为 (N) 的复数向量 x 用于存储中间结果和最终的频域信号。这里的 (N) 对应 DFT 公式中的序列长度。

2. 位反转排序

// 位反转排序
for(int i=0; i<N; ++i){int rev = bit_reverse(i, log2(N));x[i] = frame[rev];
}

Cooley - Tukey 算法采用分治法，将原始序列按照奇偶分成子序列。位反转排序是为了将序列重新排列，使得后续的蝶形运算能够正确进行。在频域中，这样的排序可以将不同频率分量正确地组合起来。从数学角度看，它是将序列的索引进行了重新映射，以满足分治算法的要求。

3. 蝶形运算

// 蝶形运算
for(int s=1; s<=log2(N); ++s){int m = 1 << s;complex<float> wm = exp(-2if * M_PI / m);for(int k=0; k<N; k+=m){complex<float> w = 1;for(int j=0; j<m/2; ++j){complex<float> t = w * x[k+j+m/2];x[k+j+m/2] = x[k+j] - t;x[k+j] += t;w *= wm;}}
}

外层循环：for(int s = 1; s <= log2(N); ++s) 控制蝶形运算的级数。对于长度为 (N = 2^n) 的序列，总共需要进行 (n=\log_2 N) 级蝶形运算。每一级的分组大小 (m = 2^s)。
旋转因子 (W_m)：complex<float> wm = exp(-2if * M_PI / m) 计算当前级的旋转因子 (W_m = e^{-j\frac{2\pi}{m}})。旋转因子是傅里叶变换中的关键元素，它体现了不同频率分量之间的相位关系。
中层循环：for(int k = 0; k < N; k += m) 遍历所有的分组。每个分组的长度为 (m)。
内层循环：for(int j = 0; j < m/2; ++j) 对每个分组内的元素进行蝶形运算。蝶形运算的公式如下：

设 (X_1[k]) 和 (X_2[k]) 是两个子序列的频域表示，(W_m^j) 是旋转因子，则蝶形运算可以表示为：

[
\begin{cases}
X[k]=X_1[k]+W_m^jX_2[k]\
X[k + \frac{m}{2}]=X_1[k]-W_m^jX_2[k]
\end{cases}
]

在代码中，t = w * x[k + j + m / 2] 对应 (W_m^jX_2[k])，x[k + j + m / 2] = x[k + j] - t 和 x[k + j] += t 分别对应上述公式中的两个式子。

4. 返回结果

return x;

最终返回的向量 x 就是经过 FFT 变换后的频域信号 (X[k])，满足离散傅里叶变换的定义。

代码存在的问题及修正

代码中 2if 是不合法的 C++ 语法，应该替换为 std::complex<float>(0, 2)。同时，bit_reverse 函数需要自行实现，以下是修正后的完整代码示例：

#include <cmath>
#include <complex>
#include <vector>// 位反转函数
int bit_reverse(int num, int bits) {int reversed = 0;for (int i = 0; i < bits; ++i) {reversed = (reversed << 1) | ((num >> i) & 1);}return reversed;
}std::vector<std::complex<float>> FFT(const std::vector<float>& frame) {const int N = frame.size();std::vector<std::complex<float>> x(N);// 位反转排序for(int i = 0; i < N; ++i) {int rev = bit_reverse(i, static_cast<int>(std::log2(N)));x[i] = frame[rev];}// 蝶形运算for(int s = 1; s <= std::log2(N); ++s) {int m = 1 << s;std::complex<float> wm = std::exp(-std::complex<float>(0, 2) * static_cast<float>(M_PI) / m);for(int k = 0; k < N; k += m) {std::complex<float> w = 1;for(int j = 0; j < m / 2; ++j) {std::complex<float> t = w * x[k + j + m / 2];x[k + j + m / 2] = x[k + j] - t;x[k + j] += t;w *= wm;}}}return x;
}