什么是快速傅里叶变换?数学公式表达和C++代码实现
FFT 简介
快速傅里叶变换(Fast Fourier Transform,FFT)是一种高效的计算离散傅里叶变换(Discrete Fourier Transform,DFT)的算法。DFT 是一种将时间域信号转换为频域信号的方法,常用于信号处理、图像处理、模式识别等领域。
FFT 算法的基本思想是将 DFT 的计算复杂度从 O ( N 2 ) O(N^2) O(N2) 降低到 O ( N log N ) O(N\log N) O(NlogN),其中 N N N 表示信号的长度。FFT 算法的具体实现基于分治法和蝴蝶算法,将 DFT 的计算分解为多个子问题,通过递归和迭代的方式计算。
FFT 数学公式表达
数学公式表达如下:
设 x n x_n xn 表示长度为 N N N 的信号, X k X_k Xk 表示其在频域中的第 k k k 个分量,则 DFT 的计算公式为:
X k = ∑ n = 0 N − 1 x n e − 2 π i n k / N X_k = \sum_{n=0}^{N-1} x_n e^{-2\pi i nk/N} Xk=n=0∑N−1xne−2πink/N
其中, i i i 表示虚数单位。FFT 算法通过将 N N N 个点的 DFT 分解为两个 N / 2 N/2 N/2 个点的 DFT,从而实现了计算复杂度的降低。
FFT 算法的基本思路是将 N N N 个点的 DFT 表示为两个 N / 2 N/2 N/2 个点的 DFT 的和差形式:
X k = ∑ n = 0 N − 1 x n e − 2 π i n k / N = ∑ m = 0 N / 2 − 1 x 2 m e − 2 π i k ( 2 m ) / N + ∑ m = 0 N / 2 − 1 x 2 m + 1 e − 2 π i k ( 2 m + 1 ) / N = ∑ m = 0 N / 2 − 1 x 2 m e − 2 π i m k / ( N / 2 ) + e − 2 π i k / N ∑ m = 0 N / 2 − 1 x 2 m + 1 e − 2 π i m k / ( N / 2 ) = X k ( e v e n ) + e − 2 π i k / N X k ( o d d ) \begin{aligned} X_k &= \sum_{n=0}^{N-1} x_n e^{-2\pi i nk/N} \ &= \sum_{m=0}^{N/2-1} x_{2m} e^{-2\pi i k(2m)/N} + \sum_{m=0}^{N/2-1} x_{2m+1} e^{-2\pi i k(2m+1)/N} \ &= \sum_{m=0}^{N/2-1} x_{2m} e^{-2\pi i mk/(N/2)} + e^{-2\pi ik/N} \sum_{m=0}^{N/2-1} x_{2m+1} e^{-2\pi imk/(N/2)} \ &= X_k^{(even)} + e^{-2\pi ik/N} X_k^{(odd)} \end{aligned} Xk=n=0∑N−1xne−2πink/N =m=0∑N/2−1x2me−2πik(2m)/N+m=0∑N/2−1x2m+1e−2πik(2m+1)/N =m=0∑N/2−1x2me−2πimk/(N/2)+e−2πik/Nm=0∑N/2−1x2m+1e−2πimk/(N/2) =Xk(even)+e−2πik/NXk(odd)
其中, X k ( e v e n ) X_k^{(even)} Xk(even) 和 X k ( o d d ) X_k^{(odd)} Xk(odd) 分别表示 x 2 m x_{2m} x2m 和 x 2 m + 1 x_{2m+1} x2m+1 的 DFT。
上式表明, N N N 个点的 DFT 可以通过两个 N / 2 N/2 N/2 个点的 DFT 计算得到,从而实现了分治的思想。
C++ 代码实现
C++ 代码实现如下:
#include <iostream>
#include <cmath>
#include <complex>
using namespace std;
const double PI = acos(-1);
// FFT算法
void fft(complex<double>* x, int n, int inv) {
if (n == 1) {
return;
}
complex<double> *even = new complex<double>[n/2];
complex<double> *odd = new complex<double>[n/2];
for (int i = 0; i < n/2; ++i) {
even[i] = x[2*i];
odd[i] = x[2*i+1];
}
fft(even, n/2, inv);
fft(odd, n/2, inv);
for (int k = 0; k < n/2; ++k) {
complex<double> t = polar(1.0, -inv*2*PI*k/n) * odd[k];
x[k] = even[k] + t;
x[k+n/2] = even[k] - t;
}
delete[] even;
delete[] odd;
}
// DFT算法
void dft(complex<double>* x, int n) {
complex<double> *y = new complex<double>[n];
for (int k = 0; k < n; ++k) {
y[k] = 0;
for (int j = 0; j < n; ++j) {
y[k] += x[j] * polar(1.0, -2PIk*j/n);
}
}
for (int k = 0; k < n; ++k) {
x[k] = y[k];
}
delete[] y;
}
// 主函数
int main() {
int n = 8;
complex<double> *x = new complex<double>[n];
for (int i = 0; i < n; ++i) {
x[i] = complex<double>(i+1, 0);
}
// DFT算法
dft(x, n);
for (int i = 0; i < n; ++i) {
cout << x[i] << " ";
}
cout << endl;
// FFT算法
fft(x, n, 1); // 计算正变换
for (int i = 0; i < n; ++i) {
cout << x[i] << " ";
}
cout << endl;
fft(x, n, -1); // 计算逆变换
for (int i = 0; i < n; ++i) {
cout << x[i] / n << " "; // 除以n是为了归一化
}
cout << endl;
delete[] x;
return 0;
}
上述代码实现了 FFT 算法和 DFT 算法(在主函数中分别调用)。其中,FFT 算法通过递归和迭代实现了 DFT 的计算,并且可以计算正变换和逆变换。需要注意的是,FFT 算法计算的结果需要进行归一化,即除以信号长度 N N N。