【多项式】快速傅里叶变换（FFT）

1. 傅里叶变换基础

1.1 傅里叶级数与傅里叶变换

傅里叶变换的核心思想是：任何周期函数都可以表示为一系列正弦和余弦函数的和。对于非周期函数，可以将其视为周期无限大的周期函数。

连续傅里叶变换（Continuous Fourier Transform）：
$$F(\omega )={\int }_{-\infty }^{\infty }f(t)e^{-i\omega t}dt$$

离散傅里叶变换（Discrete Fourier Transform, DFT）：
$$X[k]=\sum _{n=0}^{N-1}x[n]e^{-i\frac{2\pi }{N}kn},k=0,1,\dots ,N-1$$

1.2 为什么需要FFT？

直接计算DFT的时间复杂度为O(N²)，当N很大时计算量巨大。FFT通过分治策略将复杂度降低到O(N log N)，是信号处理领域的革命性算法。

2. FFT核心思想：分治法

2.1 蝶形运算（Butterfly Operation）

FFT的核心是将长度为N的DFT分解为两个长度为N/2的DFT，然后合并结果。

设序列x[n]长度为N（N为2的幂），将其分为偶数项和奇数项：

• 偶数项：$x_e[m]=x[2m],m=0,1,\dots ,N/2-1$
• 奇数项：$x_o[m]=x[2m+1],m=0,1,\dots ,N/2-1$

则DFT可表示为：
$$X[k]=\sum _{m=0}^{N/2-1}{x}_e[m]W_N^{2km}+W_N^k\sum _{m=0}^{N/2-1}{x}_o[m]W_N^{2km}$$

其中$W_N=e^{-i\frac{2\pi }{N}}$是旋转因子（Twiddle Factor）。

注意到$W_N^2=W_{N/2}$，所以：
$$X[k]=X_e[k]+W_N^k{X}_o[k]$$
$$X[k+N/2]=X_e[k]-W_N^k{X}_o[k]$$

2.2 递归实现

#include <iostream>
#include <vector>
#include <complex>
#include <cmath>

using namespace std;

const double PI = acos(-1.0);

// 复数类型
typedef complex<double> Complex;

// 递归FFT实现
vector<Complex> recursive_fft(const vector<Complex>& a) {
    int n = a.size();
    
    // 基本情况：长度为1时直接返回
    if (n == 1) {
        return vector<Complex>(1, a[0]);
    }
    
    // 分离偶数和奇数索引
    vector<Complex> a_even(n/2), a_odd(n/2);
    for (int i = 0; i < n/2; i++) {
        a_even[i] = a[2*i];
        a_odd[i] = a[2*i + 1];
    }
    
    // 递归计算
    vector<Complex> y_even = recursive_fft(a_even);
    vector<Complex> y_odd = recursive_fft(a_odd);
    
    // 合并结果
    vector<Complex> y(n);
    for (int k = 0; k < n/2; k++) {
        Complex w = Complex(cos(2*PI*k/n), -sin(2*PI*k/n)); // W_n^k
        Complex t = w * y_odd[k];
        y[k] = y_even[k] + t;
        y[k + n/2] = y_even[k] - t;
    }
    
    return y;
}

3. 迭代实现

递归实现虽然直观，但存在函数调用开销。迭代版本通过位逆序重排（Bit-reversal permutation）避免递归。

3.1 位逆序重排

观察递归过程，输入序列的最终顺序是按位逆序排列的。例如N=8时：

• 0: 000 → 000: 0
• 1: 001 → 100: 4
• 2: 010 → 010: 2
• 3: 011 → 110: 6
• 4: 100 → 001: 1
• 5: 101 → 101: 5
• 6: 110 → 011: 3
• 7: 111 → 111: 7

// 位逆序重排
void bit_reverse(vector<Complex>& a) {
    int n = a.size();
    int log_n = 0;
    while ((1 << log_n) < n) log_n++;
    
    for (int i = 0; i < n; i++) {
        int rev = 0;
        for (int j = 0; j < log_n; j++) {
            if (i & (1 << j)) {
                rev |= (1 << (log_n - 1 - j));
            }
        }
        if (i < rev) {
            swap(a[i], a[rev]);
        }
    }
}

3.2 迭代FFT实现

// 迭代FFT实现
vector<Complex> iterative_fft(vector<Complex> a, bool invert = false) {
    int n = a.size();
    
    // 位逆序重排
    bit_reverse(a);
    
    // 迭代计算
    for (int len = 2; len <= n; len <<= 1) { // 当前子问题长度
        double angle = 2 * PI / len * (invert ? 1 : -1);
        Complex wlen(cos(angle), sin(angle));
        
        for (int i = 0; i < n; i += len) { // 处理每个子问题
            Complex w(1);
            for (int j = 0; j < len/2; j++) {
                Complex u = a[i + j];
                Complex v = a[i + j + len/2] * w;
                a[i + j] = u + v;
                a[i + j + len/2] = u - v;
                w *= wlen;
            }
        }
    }
    
    // 逆变换需要除以n
    if (invert) {
        for (int i = 0; i < n; i++) {
            a[i] /= n;
        }
    }
    
    return a;
}

4. 完整的FFT类实现

class FFT {
private:
    static const double PI;
    
    // 位逆序重排
    static void bit_reverse(vector<Complex>& a) {
        int n = a.size();
        int log_n = 0;
        while ((1 << log_n) < n) log_n++;
        
        for (int i = 0; i < n; i++) {
            int rev = 0;
            for (int j = 0; j < log_n; j++) {
                if (i & (1 << j)) {
                    rev |= (1 << (log_n - 1 - j));
                }
            }
            if (i < rev) {
                swap(a[i], a[rev]);
            }
        }
    }
    
public:
    FFT() = default;
    
    // 快速傅里叶变换
    static vector<Complex> fft(const vector<Complex>& input) {
        vector<Complex> a = input;
        int n = a.size();
        
        // 确保长度是2的幂
        if ((n & (n-1)) != 0) {
            int new_n = 1;
            while (new_n < n) new_n <<= 1;
            a.resize(new_n, Complex(0, 0));
        }
        
        return iterative_fft(a, false);
    }
    
    // 快速傅里叶逆变换
    static vector<Complex> ifft(const vector<Complex>& input) {
        vector<Complex> a = input;
        int n = a.size();
        
        // 确保长度是2的幂
        if ((n & (n-1)) != 0) {
            int new_n = 1;
            while (new_n < n) new_n <<= 1;
            a.resize(new_n, Complex(0, 0));
        }
        
        return iterative_fft(a, true);
    }
    
    // 卷积计算
    static vector<Complex> convolution(const vector<Complex>& a, const vector<Complex>& b) {
        vector<Complex> fa = fft(a);
        vector<Complex> fb = fft(b);
        
        int n = fa.size();
        vector<Complex> fc(n);
        for (int i = 0; i < n; i++) {
            fc[i] = fa[i] * fb[i];
        }
        
        return ifft(fc);
    }
    
    // 实数序列的FFT
    static vector<Complex> fft_real(const vector<double>& input) {
        vector<Complex> complex_input(input.size());
        for (int i = 0; i < input.size(); i++) {
            complex_input[i] = Complex(input[i], 0);
        }
        return fft(complex_input);
    }
};

const double FFT::PI = acos(-1.0);

5. 优化技巧

5.1 预计算旋转因子

// 预计算常用旋转因子
class FFTOptimizer {
private:
    static unordered_map<int, vector<Complex>> precomputed_w;
    
public:
    static const vector<Complex>& get_w(int n) {
        if (precomputed_w.find(n) == precomputed_w.end()) {
            vector<Complex> w(n/2);
            double angle = 2 * PI / n;
            for (int i = 0; i < n/2; i++) {
                w[i] = Complex(cos(angle * i), -sin(angle * i));
            }
            precomputed_w[n] = w;
        }
        return precomputed_w[n];
    }
};

5.2 原地计算优化

// 原地FFT实现
void in_place_fft(vector<Complex>& a, bool invert = false) {
    int n = a.size();
    bit_reverse(a);
    
    for (int len = 2; len <= n; len <<= 1) {
        double angle = 2 * PI / len * (invert ? 1 : -1);
        Complex wlen(cos(angle), sin(angle));
        
        for (int i = 0; i < n; i += len) {
            Complex w(1);
            for (int j = 0; j < len/2; j++) {
                Complex u = a[i + j];
                Complex v = a[i + j + len/2] * w;
                a[i + j] = u + v;
                a[i + j + len/2] = u - v;
                w *= wlen;
            }
        }
    }
    
    if (invert) {
        for (auto& x : a) x /= n;
    }
}

6. 应用示例

6.1 信号频谱分析

void signal_analysis() {
    // 生成测试信号：10Hz + 20Hz正弦波
    const int N = 1024;
    vector<double> signal(N);
    for (int i = 0; i < N; i++) {
        double t = i / 100.0; // 采样率100Hz
        signal[i] = sin(2*PI*10*t) + 0.5*sin(2*PI*20*t);
    }
    
    // 计算FFT
    vector<Complex> spectrum = FFT::fft_real(signal);
    
    // 输出频谱
    cout << "Frequency\tMagnitude" << endl;
    for (int i = 0; i < N/2; i++) {
        double freq = i * 100.0 / N; // 频率
        double magnitude = abs(spectrum[i]);
        if (magnitude > 0.1) { // 只显示显著成分
            cout << freq << "\t\t" << magnitude << endl;
        }
    }
}

6.2 多项式乘法

vector<double> polynomial_multiply(const vector<double>& a, const vector<double>& b) {
    vector<Complex> ca(a.begin(), a.end());
    vector<Complex> cb(b.begin(), b.end());
    
    vector<Complex> result = FFT::convolution(ca, cb);
    
    vector<double> real_result(result.size());
    for (int i = 0; i < result.size(); i++) {
        real_result[i] = result[i].real();
    }
    
    return real_result;
}

7. 复杂度分析

• 时间复杂度：O(N log N)
• 空间复杂度：O(N)
• 精度：受浮点数精度限制，大N时可能有累积误差

8. 注意事项

1. 长度要求：标准FFT要求序列长度为2的幂，可通过补零处理
2. 精度问题：浮点数计算存在舍入误差
3. 内存访问：位逆序重排可能影响缓存性能
4. 边界情况：处理长度为1的特殊情况

9. 总结

FFT是数字信号处理的核心算法，通过分治策略将DFT的O(N²)复杂度降低到O(N log N)。迭代实现比递归实现更高效，实际应用中还需考虑精度、内存和特定硬件的优化。