我怎样才能让 FFTW 希尔伯特变换计算得更快？

问题描述

我正在使用来自 FFTW 源的希尔伯特变换函数。 因为我将在数据流模式下在我的 DAQ 中使用它。 该函数运行正常，但计算速度较慢，会导致 FIFO 溢出。 我听说将 fftw_plan() 从 hilbert() 移到外面以重复使用 plan 可能很有用，但是，一旦我这样做就会出错，在 Exception thrown at 0x0000000070769240 (libfftw3-3.dll) in CppFFTW.exe: 0xC0000005: Access violation reading location 0x0000000000000000. fftw_destroy_plan(plan);。有没有人有类似的经验或更好的解决方案来提高 hilbert() 计算？

这是我尝试过的（2020 12/26 编辑）：

#include <iostream>
#include <fftw3.h>
#include <time.h>


using namespace std;

//macros for real and imaginary parts
#define REAL 0
#define IMAG 1
//length of complex array
#define N 8

fftw_complex* out;
fftw_plan plan;


void hilbert(const double* in,fftw_complex* out)
{
    // copy the data to the complex array
    for (int i = 0; i < N; ++i) {
        out[i][REAL] = in[i];
        out[i][IMAG] = 0;
    }
    // creat a DFT plan and execute it
    //fftw_plan plan = fftw_plan_dft_1d(N,out,FFTW_FORWARD,FFTW_ESTIMATE);
    fftw_execute(plan);
    // destroy a plan to prevent memory leak
    //fftw_destroy_plan(plan);
    int hN = N >> 1; // half of the length (N/2)
    int numRem = hN; // the number of remaining elements
    // multiply the appropriate value by 2
    //(those should multiplied by 1 are left intact because they wouldn't change)
    for (int i = 1; i < hN; ++i) {
        out[i][REAL] *= 2;
        out[i][IMAG] *= 2;
    }
    // if the length is even,the number of the remaining elements decrease by 1
    if (N % 2 == 0)
        numRem--;
    else if (N > 1) {
        out[hN][REAL] *= 2;
        out[hN][IMAG] *= 2;
    }
    // set the remaining value to 0
    // (multiplying by 0 gives 0,so we don't care about the multiplicands)
    memset(&out[hN + 1][REAL],numRem * sizeof(fftw_complex));
    // creat a IDFT plan and execute it
    //plan = fftw_plan_dft_1d(N,FFTW_BACKWARD,FFTW_ESTIMATE);
    fftw_execute(plan);
    // do some cleaning
    //fftw_destroy_plan(plan);
    //fftw_cleanup();
    // scale the IDFT output
    for (int i = 0; i < N; ++i) {
        out[i][REAL] /= N;
        out[i][IMAG] /= N;
    }



   




}
/* displays complex numbers in the form a +/- bi. */
void displayComplex(fftw_complex* y)
{
    for (int i = 0; i < N; i++) {
        if (y[i][IMAG] < 0)
            cout << y[i][REAL] << "-" << abs(y[i][IMAG]) << "i" << endl;
        else
            cout << y[i][REAL] << "+" << y[i][IMAG] << "i" << endl;
    }
}



/* Test */
int main()
{
    fftw_plan plan = fftw_plan_dft_1d(N,FFTW_ESTIMATE);
    plan = fftw_plan_dft_1d(N,FFTW_ESTIMATE);
    

    
    // input array
    double x[N];
    // output array
    fftw_complex y[N];
    // fill the first of some numbers
    for (int i = 0; i < N; ++i) {
        x[i] = i + 1;  // i.e.{1 2 3 4 5 6 7 8}
    }
    //clock_t begin = clock();
    // compute the FFT of x and store the result in y.
    hilbert(x,y);
    // display the result
    cout << "Hilbert =" << endl;
    displayComplex(y);
    fftw_destroy_plan(plan);
    fftw_destroy_plan(plan);
}

这是原始代码：

#include <iostream>
#include <fftw3.h>


using namespace std;

//macros for real and imaginary parts
#define REAL 0
#define IMAG 1
//length of complex array
#define N 8

void hilbert(const double* in,fftw_complex* out)
{
    // copy the data to the complex array
    for (int i = 0; i < N; ++i) {
        out[i][REAL] = in[i];
        out[i][IMAG] = 0;
    }
    // creat a DFT plan and execute it
    fftw_plan plan = fftw_plan_dft_1d(N,FFTW_ESTIMATE);
    fftw_execute(plan);
    // destroy a plan to prevent memory leak
    fftw_destroy_plan(plan);
    int hN = N >> 1; // half of the length (N/2)
    int numRem = hN; // the number of remaining elements
    // multiply the appropriate value by 2
    //(those should multiplied by 1 are left intact because they wouldn't change)
    for (int i = 1; i < hN; ++i) {
        out[i][REAL] *= 2;
        out[i][IMAG] *= 2;
    }
    // if the length is even,numRem * sizeof(fftw_complex));
    // creat a IDFT plan and execute it
    plan = fftw_plan_dft_1d(N,FFTW_ESTIMATE);
    fftw_execute(plan);
    // do some cleaning
    fftw_destroy_plan(plan);
    fftw_cleanup();
    // scale the IDFT output
    for (int i = 0; i < N; ++i) {
        out[i][REAL] /= N;
        out[i][IMAG] /= N;
    }



   




}
/* displays complex numbers in the form a +/- bi. */
void displayComplex(fftw_complex* y)
{
    for (int i = 0; i < N; i++) {
        if (y[i][IMAG] < 0)
            cout << y[i][REAL] << "-" << abs(y[i][IMAG]) << "i" << endl;
        else
            cout << y[i][REAL] << "+" << y[i][IMAG] << "i" << endl;
    }
}



/* Test */
int main()
{
    
 
    // input array
    double x[N];
    // output array
    fftw_complex y[N];
    // fill the first of some numbers
    for (int i = 0; i < N; ++i) {
        x[i] = i + 1;  // i.e.{1 2 3 4 5 6 7 8}
    }
    // compute the FFT of x and store the result in y.
    hilbert(x,y);
    // display the result
    cout << "Hilbert =" << endl;
    displayComplex(y);
 
}

准确的输出

Hilbert =
1+3.82843i
2-1i
3-1i
4-1.82843i
5-1.82843i
6-1i
7-1i
8+3.82843i

解决方法

为了速度，您肯定希望尽可能重复使用 FFTW 计划。在多次调用 hilbert 时重用计划时，请删除 hilbert 中的 fftw_destroy_plan(plan) 行，否则该计划将无法在下一次调用中使用。相反，在程序结束时销毁计划。

编辑：我之前错过了您使用两个计划，一个用于向前转换，另一个用于向后转换。在 hilbert() 中，删除对 fftw_plan_dft_1d、fftw_destroy_plan 和 fftw_cleanup 的所有调用；唯一的 FFTW 调用应该是 fftw_execute。相反，在程序开始时预先创建两个计划，并在程序结束时销毁它们。

,

经过多次尝试，这是成功的 FFTW hilbert() 改进。 移出两个 fftw_plan 并让它们在 main() 中初始化，另外，fftw_destroy_plan 也被移走了。

#include <iostream>
#include <fftw3.h>
#include <time.h>


using namespace std;

//macros for real and imaginary parts
#define REAL 0
#define IMAG 1
//length of complex array
#define N 8
fftw_complex y[N];



void hilbert(const double* in,fftw_complex* out,fftw_plan plan_forward,fftw_plan plan_backward)
{
    // copy the data to the complex array
    for (int i = 0; i < N; ++i) {
        out[i][REAL] = in[i];
        out[i][IMAG] = 0;
    }
    // creat a DFT plan and execute it
    //fftw_plan plan = fftw_plan_dft_1d(N,out,FFTW_FORWARD,FFTW_ESTIMATE);
    fftw_execute(plan_forward);
    // destroy a plan to prevent memory leak
    //fftw_destroy_plan(plan_forward);
    int hN = N >> 1; // half of the length (N/2)
    int numRem = hN; // the number of remaining elements
    // multiply the appropriate value by 2
    //(those should multiplied by 1 are left intact because they wouldn't change)
    for (int i = 1; i < hN; ++i) {
        out[i][REAL] *= 2;
        out[i][IMAG] *= 2;
    }
    // if the length is even,the number of the remaining elements decrease by 1
    if (N % 2 == 0)
        numRem--;
    else if (N > 1) {
        out[hN][REAL] *= 2;
        out[hN][IMAG] *= 2;
    }
    // set the remaining value to 0
    // (multiplying by 0 gives 0,so we don't care about the multiplicands)
    memset(&out[hN + 1][REAL],numRem * sizeof(fftw_complex));
    // creat a IDFT plan and execute it
    //plan = fftw_plan_dft_1d(N,FFTW_BACKWARD,FFTW_ESTIMATE);
    fftw_execute(plan_backward);
    // do some cleaning
    //fftw_destroy_plan(plan_backward);
    //fftw_cleanup();
    // scale the IDFT output
    for (int i = 0; i < N; ++i) {
        out[i][REAL] /= N;
        out[i][IMAG] /= N;
    }








}
/* Displays complex numbers in the form a +/- bi. */
void displayComplex(fftw_complex* y)
{
    for (int i = 0; i < N; i++) {
        if (y[i][IMAG] < 0)
            cout << y[i][REAL] << "-" << abs(y[i][IMAG]) << "i" << endl;
        else
            cout << y[i][REAL] << "+" << y[i][IMAG] << "i" << endl;
    }
}



/* Test */
int main()
{
    // input array
    double x[N];
    // output array
    //fftw_complex y[N];
    fftw_plan plan_forward = fftw_plan_dft_1d(N,y,FFTW_ESTIMATE);
    fftw_plan plan_backward = fftw_plan_dft_1d(N,FFTW_ESTIMATE);
    



   
    // fill the first of some numbers
    for (int i = 0; i < N; ++i) {
        x[i] = i + 1;  // i.e.{1 2 3 4 5 6 7 8}
    }
    // compute the FFT of x and store the result in y.
    hilbert(x,plan_forward,plan_backward);
    
    // display the result
    cout << "Hilbert =" << endl;
    
    displayComplex(y);
    fftw_destroy_plan(plan_forward);
    fftw_destroy_plan(plan_backward);
   
}

,

以我最喜欢的评论开头...

好的...经过几次错误的启动...

您的第二个版本有一些问题。

值得注意的是，您没有分配out。

另外，为了安全起见，我相信您需要两个 fftw_plan 结构，一个用于向前，一个用于向后。

这些应该在main中一次分配/初始化。

并且，删除 hilbert 中的所有分配/销毁调用。这样，只需更改放置在 out 中的数据即可多次调用它。

清理代码可以到main的底部。

我已经创建了您的第二个版本的清理和工作版本。我已经用 cpp 条件展示了新旧代码的区别：

#if 0
// old code ...
#else
// new code ...
#endif

所以，这里是：

#include <iostream>
#include <cstring>
#include <fftw3.h>
#include <time.h>

using namespace std;

//macros for real and imaginary parts
#define REAL 0
#define IMAG 1

//length of complex array
#define N 8

fftw_plan plan_fwd;
fftw_plan plan_bwd;
fftw_complex *out;

void
hilbert(const double *in,fftw_complex *out)
{
    // copy the data to the complex array
    for (int i = 0; i < N; ++i) {
        out[i][REAL] = in[i];
        out[i][IMAG] = 0;
    }

    // creat a DFT plan and execute it
    // fftw_plan plan = fftw_plan_dft_1d(N,FFTW_ESTIMATE); // This line has been moved to the main function

    fftw_execute(plan_fwd);

    // destroy a plan to prevent memory leak
#if 0
    fftw_destroy_plan(plan_fwd);
#endif
    int hN = N >> 1;                    // half of the length (N/2)
    int numRem = hN;                    // the number of remaining elements

    // multiply the appropriate value by 2
    // (those should multiplied by 1 are left intact because they wouldn't change)
    for (int i = 1; i < hN; ++i) {
        out[i][REAL] *= 2;
        out[i][IMAG] *= 2;
    }

    // if the length is even,the number of the remaining elements decrease by 1
    if (N % 2 == 0)
        numRem--;
    else if (N > 1) {
        out[hN][REAL] *= 2;
        out[hN][IMAG] *= 2;
    }

    // set the remaining value to 0
    // (multiplying by 0 gives 0,numRem * sizeof(fftw_complex));
    // creat a IDFT plan and execute it
#if 0
    plan_bwd = fftw_plan_dft_1d(N,FFTW_ESTIMATE);
#endif
    fftw_execute(plan_bwd);
    // do some cleaning
#if 0
    fftw_destroy_plan(plan_bwd);
    fftw_cleanup();
#endif

    // scale the IDFT output
    for (int i = 0; i < N; ++i) {
        out[i][REAL] /= N;
        out[i][IMAG] /= N;
    }
}

/* Displays complex numbers in the form a +/- bi. */
void
displayComplex(fftw_complex * y)
{
    for (int i = 0; i < N; i++) {
        if (y[i][IMAG] < 0)
            cout << y[i][REAL] << "-" << abs(y[i][IMAG]) << "i" << endl;
        else
            cout << y[i][REAL] << "+" << y[i][IMAG] << "i" << endl;
    }
}

/* Test */
int
main()
{

#if 1
    out = (fftw_complex *) malloc(sizeof(fftw_complex) * N);
#endif

    plan_fwd = fftw_plan_dft_1d(N,FFTW_ESTIMATE);
    plan_bwd = fftw_plan_dft_1d(N,FFTW_ESTIMATE);

    // input array
    double x[N];

    // output array
    fftw_complex y[N];

    // fill the first of some numbers
    for (int i = 0; i < N; ++i) {
        x[i] = i + 1;                   // i.e.{1 2 3 4 5 6 7 8}
    }
    clock_t begin = clock();

    // compute the FFT of x and store the result in y.
    hilbert(x,y);
    // display the result
    cout << "Hilbert =" << endl;
    displayComplex(y);
    clock_t end = clock();
    double time_spent = (double) (end - begin) / CLOCKS_PER_SEC;

    printf("%f",time_spent);

    fftw_destroy_plan(plan_fwd);
    fftw_destroy_plan(plan_bwd);
    fftw_cleanup();
}

---

这是程序输出：

Hilbert =
0.125+0i
0.5+0i
0.75+0i
1+0i
0.625+0i
0+0i
0+0i
0+0i
0.000171