递归FFT处理大数据集时程序中途终止问题求助

问题详情 #

研究FFT算法时，用超过65536条数据测试效率，输入200000条数据时程序无报错直接终止，递归流程无法完成。程序终止时输出如下：

PS C:\Users\user\proj> ./a.exe
Loading data
Calculating fft
200000
100000
50000
25000
12500
6250
3125
1562
781
PS C:\Users\user\proj> 

用到的代码 #

FFT实现程序 #

#include <complex>
#include <iostream>
#include <vector>
#include <cmath>
#include <chrono>
#include <fstream>
#include <sstream>

using namespace std;
using namespace std::chrono;

const double PI = 3.141592653589793238460;

typedef complex<double> cd;
typedef vector<cd> vcd;

void fft(vcd& a) {
    int n = a.size();
    cout << a.size() << endl;
    if (n <= 1) return;

    vcd a0(n / 2), a1(n / 2);
    for (int i = 0; 2 * i < n; i++) {
        a0[i] = a[2*i];
        a1[i] = a[2*i + 1];
    }

    fft(a0);
    fft(a1);

    double ang = 2 * PI / n;
    cd w(1), wn(cos(ang), sin(ang));
    for (int i = 0; 2 * i < n; i++) {
        a[i] = a0[i] + w * a1[i];
        a[i + n/2] = a0[i] - w * a1[i];
        w *= wn;
    }
}

int main() {
    ifstream dataFile("data.txt");
    string line;
    vcd data;

    if (!dataFile.is_open()) {
        cout << "Couldn't open the file" << endl;
        return 1;
    }
    cout << "Loading data" << endl;
    while (getline(dataFile, line)) {
        stringstream ss(line);
        double val;
        ss >> val;
        data.push_back(val);
    }

    dataFile.close();

    auto start = high_resolution_clock::now();
    cout << "Calculating fft" << endl;
    fft(data);

    auto stop = high_resolution_clock::now();
    auto duration = duration_cast<microseconds>(stop - start);

    cout << "FFT result: " << endl;
    // for (int i = 0; i < data.size(); i++)
        // cout << data[i] << endl;

    cout << "Time taken FFT: " << duration.count() << endl;

    return 0;
}

输入数据生成程序 #

#include <iostream>
#include <fstream>
#include <cmath>

using namespace std;

const double PI = 3.141592653589793238460;

int main() {
    // int N = 200000; // stops
    int N = 65536; // works
    ofstream dataFile("data.txt");

    if (!dataFile.is_open()) {
        return 1;
    }

    for (int i = 0; i < N; i++) {
        double dataPoint = 0.5 * sin(2 * PI * 5 * i / N) + 0.5 * sin(2 * PI * 10 * i / N) + 0.3 * sin(2 * PI * 20 * i / N);
        dataFile << dataPoint << endl;
    }

    dataFile.close();
    cout << "Data saved" << endl;

    return 0;
}

问题原因 #

你的递归FFT实现仅支持长度为2的整数次幂的输入数据，但200000并不是2的幂：

• 递归过程中，数据长度不断除以2，到781时（200000÷256=781.25，取整后为781），此时n=781，n/2=390。
• 循环for (int i = 0; 2 * i < n; i++)会执行到i=390，此时2*i+1=781，而原数组a的长度是781，索引范围是0~780，数组越界访问直接导致程序崩溃终止。
• 65536是2^16，递归过程中每次长度都是偶数，不会触发越界问题，因此可以正常运行。

解决方案 #

方案1：补零到最近的2的幂次 #

在FFT计算前，将输入数据补零，使其长度变为大于等于原长度的最小2的幂：

// 在main函数中调用fft(data)前添加以下代码
int n = 1;
while (n < data.size()) n <<= 1;
data.resize(n, cd(0));

处理后输入长度会被调整为262144（2^18），可正常完成递归FFT。

方案2：改用迭代版FFT #

递归版FFT处理大尺寸数据时还可能遇到栈溢出问题（递归深度过大），迭代版FFT（蝶形算法）可避免该问题，同时性能通常更优。

方案3：实现支持任意长度的FFT算法 #

如果需要处理非2幂长度的数据，可实现混合基FFT或Bluestein算法，但实现复杂度更高，多数场景下补零到2的幂是更简单高效的选择。

内容的提问来源于stack exchange，提问作者zarky