C++中FFTW实现并行运算坑点记录

苦瓜大王_困版
452 阅读4分钟

问题定义

项目中有一个算法调用到了FFTW库，希望可以并行调用进行计算，大致长这样，fftw库在Line27-46使用：

void funcTest(const std::vector<double>& ad_data, double fs, std::vector<double>& VSD_1Hz, std::vector<double>& fre_1Hz)
{
    double U = 0.0;
    std::vector<double> data(ad_data.begin(), ad_data.end());
    std::vector<double> hann_window = hannWindow(data.size());
    for (size_t i = 0; i < data.size(); ++i) {
        data[i] *= hann_window[i];
        U += hann_window[i] * hann_window[i];
    }
    static fftw_complex* pFFTResult = nullptr;
    static unsigned int nFFTResultLength = 0;
    unsigned int nDataLen = data.size();
    unsigned int nDataHalfLen = nDataLen / 2 + 1;
    std::vector<double> vecPSDCache(nDataHalfLen, 0.0);
    int avg_points = 0;
    double f0 = 0.0;
    {
        int len = static_cast<int>(round(fs / 2) + 2);
        VSD_1Hz.resize(len, 0.0);
        fre_1Hz.resize(len, 0.0);

        double df = fs / 2 / nDataHalfLen;
        avg_points = static_cast<int>(round(1.0 / df));
        f0 = avg_points * df / 2;
    }
    {
        QMutexLock.lock();
        //执行FFT计算
        {
            fftw_plan plan;
            if (pFFTResult == nullptr || nDataHalfLen > nFFTResultLength) {
                if (pFFTResult) {
                    fftw_free(pFFTResult);
                }
                pFFTResult = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * nDataHalfLen);
                CHECK_RET_VOID(pFFTResult);
                nFFTResultLength = nDataHalfLen;
            }
            memset(pFFTResult, 0, sizeof(fftw_complex) * nFFTResultLength);
            plan = fftw_plan_dft_r2c_1d(nDataLen, data.data(), pFFTResult, FFTW_ESTIMATE);
            if (plan == nullptr) {
                fftw_free(pFFTResult);
                return;
            }
            fftw_execute(plan);
            fftw_destroy_plan(plan);
        }
        {
            CHECK_RET_VOID(pFFTResult);
            for (int i = 0; i < nDataHalfLen; ++i) {
                vecPSDCache[i] = std::pow(pFFTResult[i][0], 2) + std::pow(pFFTResult[i][1], 2);
            }
        }
        QMutexLock.unlock();
    }

会在多个并发线程中调用funcTest(m_calibrateData[chnId], fs, vsd, fre);，由于有线程锁的存在导致实际上是串行的。希望可以实现并发进行。

初步解决办法

初步解决方法是将funcTest放入一个类中实现，同时将需要竞争使用的资源：

    static fftw_complex* pFFTResult = nullptr;
    static unsigned int nFFTResultLength = 0;

改为类中的成员。然后每个线程都实例化一个该对象,每个线程都有自己的资源故而可以去掉线程锁QMutexLock。类的实现：

//.h
class STest {
private:
    std::vector<double> hannWindow(int size);

    fftw_complex *pFFTResult;//无需竞争
    unsigned int nFFTResultLength;//无需竞争
    fftw_plan fftPlan;
    int chnId;
    int count;

public:
    STest(int id): chnId(id), pFFTResult(nullptr), fftPlan(nullptr), nFFTResultLength(0) {}
    ~STest();
    void funcTest(const std::vector<double> &ad_data, double fs, std::vector<double> &VSD_1Hz, std::vector<double> &fre_1Hz);

};
//.cpp
double M_PI = 3.14159265358979323846;

STest::~STest() {
    if (pFFTResult) {
        fftw_free(pFFTResult);
    }
    if (fftPlan) {
        fftw_destroy_plan(fftPlan);
    }
}

// 计算汉宁窗
std::vector<double> STest::hannWindow(int size) {
    std::vector<double> window(size);
    for (int i = 0; i < size; ++i) {
        window[i] = 0.5 * (1 - std::cos(2 * M_PI *i / (size - 1)));
    }
    return window;
}

void STest::funcTest(const std::vector<double> &ad_data, double fs, std::vector<double> &VSD_1Hz, std::vector<double> &fre_1Hz) {
   double U = 0.0; // 复制数据到FFT输入数组 
   std::vector<double> data(ad_data.begin(), ad_data.end()); // 应用汉宁窗 
   std::vector<double> hann_window = hannWindow(data.size()); 
   for (size_t i = 0; i < data.size(); ++i) { 
       data[i] *= hann_window[i]; 
       U += hann_window[i] * hann_window[i]; 
   }    
    fftw_complex *pFFTResult = nullptr; 
    unsigned int nFFTResultLength = 0;
    unsigned int nDataLen = data.size(); 
    unsigned int nDataHalfLen = nDataLen / 2 + 1; // 计算功率谱密度和1Hz平均-准备工作
    std::vector<double> vecPSDCache(nDataHalfLen, 0.0); 
    int avg_points = 0; double f0 = 0.0; 
    { 
        int len = static_cast<int>(round(fs / 2) + 2); 
        VSD_1Hz.resize(len, 0.0); 
        fre_1Hz.resize(len, 0.0); 
        double df = fs / 2 / nDataHalfLen; 
        avg_points = static_cast<int>(round(1.0 / df)); 
        f0 = avg_points *df / 2; 
    } 
    { 
        //执行FFT计算 
        { 
            if (pFFTResult == nullptr || nDataHalfLen > nFFTResultLength){ 
                if (pFFTResult) { 
                    fftw_free(pFFTResult); 
                } 
                pFFTResult = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * nDataHalfLen);
                CHECK_RET_VOID(pFFTResult); 
                nFFTResultLength = nDataHalfLen; 
            } 
            memset(pFFTResult, 0, sizeof(fftw_complex) * nFFTResultLength + 1); 
            qDebug() << chnId << "-- plan前，第" << count << "次"; 
            fftPlan = fftw_plan_dft_r2c_1d(nDataLen, data.data(), pFFTResult, FFTW_ESTIMATE); 
            if (fftPlan == nullptr) { 
                fftw_free(pFFTResult); 
                return; 
            } 
            fftw_execute(fftPlan); 
            if(fftPlan != nullptr) { 
                fftw_destroy_plan(fftPlan); 
            } 
            qDebug() << chnId << "-- plan结束，第" << count << "次"; 
        }
        { 
            CHECK_RET_VOID(pFFTResult); 
            for (int i = 0; i < nDataHalfLen; ++i) { 
                vecPSDCache[i] = std::pow(pFFTResult[i][0], 2) + std::pow(pFFTResult[i][1], 2); 
            } 
        } 
    }
    count++;
}

//多线程中调用
m_STest[chnId]->funcTest(m_calibrateData[chnId], fs, vsd, fre);

再次遇到bug

使用上述代码，在debug模式运行下报错：打印结果如下：

0 -- plan前，第 0 次 
4 -- plan前，第 0 次 
8 -- plan前，第 0 次 
12 -- plan前，第 0 次 
Critical error detected c0000374

程序是在Line73fftPlan = fftw_plan_dft_r2c_1d(nDataLen, data.data(), pFFTResult, FFTW_ESTIMATE); 处报错。

这块卡了我2-3天，反复在找是不是还有哪个成员是竞争关系，哪个指针是没有实例化就调用，哪个指针是重复删除了（百度搜索'Critical error detected c0000374'看到的解决方案）。

直到我看到了这个问答：如何在多线程中使用FFTW？-腾讯云开发者社区-腾讯云 (tencent.com)，和我的使用场景类似。下面的评论给了我启示：

FFTW人员为线程安全主题这里提供了一个很好的总结。包装:除了fftw_execute之外，没有什么是线程安全的，所以您必须注意，例如，只有一个线程创建计划。然而，并行执行它们应该是没有问题的。

也就是说Line73fftPlan = fftw_plan_dft_r2c_1d(nDataLen, data.data(), pFFTResult, FFTW_ESTIMATE);在fftw库中调用fftw_plan_dft_r2c_1d方法来获得fftPlan是线程不安全的，因此这行代码是不能多线程并发使用的。需要给这行代码加上线程锁。

最终代码

//.h
class STest {
private:
    std::vector<double> hannWindow(int size);

    fftw_complex *pFFTResult;//无需竞争
    unsigned int nFFTResultLength;//无需竞争
    fftw_plan fftPlan;
    int chnId;
    static std::mutex psdMutex;
    
public:
    STest(int id): chnId(id), pFFTResult(nullptr), fftPlan(nullptr), nFFTResultLength(0) {}
    ~STest();
    void funcTest(const std::vector<double> &ad_data, double fs, std::vector<double> &VSD_1Hz, std::vector<double> &fre_1Hz);

};
//.cpp
std::mutex STest::psdMutex;

void STest::funcTest(const std::vector<double> &ad_data, double fs, std::vector<double> &VSD_1Hz, std::vector<double> &fre_1Hz) {
            //······
            psdMutex.lock(); // fftw_plan_dft_r2c_1d时加上锁
            fftPlan = fftw_plan_dft_r2c_1d(nDataLen, data.data(), pFFTResult, FFTW_ESTIMATE); 
            psdMutex.unlock();
            //······
}

总结一下：就是那句：除了fftw_execute之外，没有什么是线程安全的，所以想要多线程执行fftw，只有一个线程创建所有的计划。然后，并行执行它们。或者创建计划时加锁

avatar
文章
阅读
粉丝