C++循环调用FFTW3执行FFT计算时结果偶发异常问题求助

问题描述 #

• 开发场景：基于C++结合Eigen库、fftw3（3.3.10版本）实现FFT功能，读取CSV存储的测量数据完成频谱分析
• 异常现象：循环重复执行分析逻辑时，FFT结果偶发异常，部分循环返回正确结果、部分返回错误结果，一致性无法保证
• 初步怀疑方向：内存泄漏、变量类型转换错误

复现代码 #

#define EIGEN_FFTW_DEFAULT
#include <iostream>
#include <string>
#include <vector>
#include <cmath>
#include <fstream>
#include <sstream>
#include "Eigen/Dense"
#include <fftw3.h>
using namespace Eigen;

// define functions
template <typename T>
T readCSV(const std::string &path);
int nextpow2(int n);
VectorXd offsetData(VectorXd v);
void fftw_test(VectorXd x);

// calculate exponent of next higher power of 2
int nextpow2(int n)
{
    if (n < 0) // n must be int
        return 0;
    if (n == 1) // n
        return 1;

    return (int)floor(log2(n - 1)) + 1.0;
};

// Read Measurement Data
template <typename T>
T readCSV(const std::string &path)
{
    std::ifstream file;
    std::string line;
    std::string cell;
    std::vector<double> row;
    uint rows = 0;

    file.open(path);
    std::cout << "Opend file: " << path << std::endl;

    std::getline(file, line); // skip the first header line

    while (std::getline(file, line))
    {
        std::stringstream lineStream(line);

        while (std::getline(lineStream, cell, ','))
        {
            row.push_back(std::stod(cell)); // insert value as double
        }
        ++rows;
    }
    return Map<const Matrix<typename T::Scalar, T ::RowsAtCompileTime, T::ColsAtCompileTime, RowMajor> >(row.data(), rows, row.size() / rows);
};

void fftw_test(VectorXd x)
{
    // Convert data unit
    x = x * 980.665 * 10;                 // Unit conversion：[G] to [cm/sec^2] to [mm/sec^2]
    int ns = x.size();                    // number of samples
    int nfft = std::pow(2, nextpow2(ns)); // number of fft

    // Zero padding to array
    VectorXd xpad;
    int npad = nfft - ns;

    if (npad > 0)
    {
        xpad = VectorXd(nfft);

        for (int i = 0; i < ns; ++i)
        {
            xpad(i) = x(i);
        }
    }
    else
    {
        xpad = x;
    }

    int N = nfft;
    fftw_complex *in, *out, *in2;
    in = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * N);
    out = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * N);
    in2 = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * N);

    fftw_plan p, q;
    for (int i = 0; i < N; i++)
    {
        in[i][0] = (double)xpad(i);
        in[i][1] = 0;
    }

    p = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE | FFTW_PRESERVE_INPUT);
    fftw_execute(p);

    for (int i = 0; i < 10; i++)
    {
        printf("in: %3d %+9.5f %+9.5f I\n", i, in[i][0], in[i][1]);
    }

    for (int i = 0; i < 10; i++)
    {
        printf("freq: %3d %+9.5f %+9.5f I\n", i, out[i][0], out[i][1]);
    }

    fftw_destroy_plan(p);
    fftw_free(in);
    fftw_free(out);
    fftw_free(in2);
    fftw_cleanup();
};

VectorXd offsetData(VectorXd v)
{
    // Offset by mean values
    int ns = v.size(); // number of samples
    VectorXd ones = MatrixXd::Ones(ns, 1);
    v = v - v.mean() * ones;
    return v;
};

int main()
{
    // Read measured data from csv file
    MatrixXd measuredData = readCSV<MatrixXd>("./sampleCsv/20220208-134655_A351AU.csv");

    // Extract a vertical acceleration column
    VectorXd Acc = measuredData.col(4);
    VectorXd Acc_offset = offsetData(Acc / 1000);

    for (int i = 0; i < 100; ++i)
    {
        // fftw bug test
        printf("loop: %ith \n", i);
        fftw_test(Acc_offset);
    };

    return 0;
}

根因分析 #

核心问题是CSV读取函数存在悬空内存引用，触发未定义行为，这是结果随机异常的直接原因：

1. readCSV函数中std::vector<double> row是栈上局部变量，函数返回时通过Eigen::Map直接映射row.data()的内存构造返回值，但Map不会主动拷贝源数据
2. 函数执行结束后局部变量row被自动析构，其占用的内存被释放归还给系统，返回的Eigen矩阵实际指向已经失效的内存地址
3. 后续循环访问该矩阵数据时，失效内存的内容是不确定的：如果内存未被其他变量覆盖就会得到正确结果，一旦被其他操作写入新值就会得到错误结果，完全符合偶发异常的特征

其余次要问题：

• 零填充逻辑存在缺陷：新建xpad向量时Eigen不会自动初始化元素为0，仅填充了前ns个有效采样点，末尾npad个填充位是随机垃圾值，会导致FFT结果固定偏移
• 冗余资源操作：每次调用fftw_test都执行fftw_cleanup()，会清空FFTW全局缓存的优化计划，反复创建销毁计划反而会降低性能，无需每次调用都执行清理
• 冗余变量：定义了未使用的fftw_plan q和fftw_complex *in2，属于无效代码

修复方案 #

1. 修复CSV读取的悬空引用问题：将Map映射的内存拷贝到独立的矩阵对象中再返回，保证返回值持有独立的内存空间，修改后的readCSV返回逻辑如下：

// 原返回逻辑直接返回Map，改为先构造临时矩阵再返回，触发数据拷贝
return Matrix<typename T::Scalar, T::RowsAtCompileTime, T::ColsAtCompileTime, RowMajor>(
    Map<const Matrix<typename T::Scalar, T ::RowsAtCompileTime, T::ColsAtCompileTime, RowMajor> >(
        row.data(), rows, row.size() / rows
    )
);

2. 修复零填充逻辑：创建xpad后先将所有元素初始化为0，再填充有效采样点：

if (npad > 0)
{
    xpad = VectorXd::Zero(nfft); // 初始化为全0
    for (int i = 0; i < ns; ++i)
    {
        xpad(i) = x(i);
    }
}

3. 移除不必要的fftw_cleanup()调用：该函数仅需在程序完全结束FFT操作前调用一次即可，无需每次执行FFT都调用
4. 删除未使用的冗余变量q和in2，避免无意义的内存申请释放

正确结果参考 #

freq:   0  -0.00000  +0.00000 I
freq:   1 +320.64441 -83.56961 I
freq:   2 -113.66004 -195.80680 I
freq:   3 -28.57778 -13.57046 I
freq:   4 -47.71908 +185.43538 I
freq:   5 +381.01770 +92.18739 I
freq:   6 +430.73267 -348.16464 I
freq:   7 -111.55714 -796.10333 I
freq:   8 -810.79331 -273.42916 I
freq:   9 -624.83461 +607.38775 I

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