How could I make FFTW Hilbert transform calculate faster?

Question

I'm using the Hilbert transform function from the FFTW source. Because I am going to use it in my DAQ with data streaming mode. The function is working fine but the calculation speed is slow which will cause the FIFO overflow. I've heard that move the fftw_plan() outside from the hilbert() for reuse the plan might be useful, however, it's an error once I did that, saying Exception thrown at 0x0000000070769240 (libfftw3-3.dll) in CppFFTW.exe: 0xC0000005: Access violation reading location 0x0000000000000000. at the fftw_destroy_plan(plan);. Does anyone has similar experiences or even better solution to boost up the hilbert() calculation?

Here is what I've tried (2020 12/26 edited):

#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, 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], 0, numRem * sizeof(fftw_complex));
    // creat a IDFT plan and execute it
    //plan = fftw_plan_dft_1d(N, out, out, 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, out, out, FFTW_FORWARD, FFTW_ESTIMATE);
    plan = fftw_plan_dft_1d(N, out, out, FFTW_BACKWARD, 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);
}

Here is the original code:

#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, out, 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], 0, numRem * sizeof(fftw_complex));
    // creat a IDFT plan and execute it
    plan = fftw_plan_dft_1d(N, out, out, 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()
{
    
 
    // 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);
 
}

The exact output

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

Answer 1

After several tried, here is the successful FFTW hilbert() improvement. Moved out the two fftw_plan and let them initialized in the main(), plus, the fftw_destroy_plan were moved too.

#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, 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], 0, numRem * sizeof(fftw_complex));
    // creat a IDFT plan and execute it
    //plan = fftw_plan_dft_1d(N, out, out, 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, y, FFTW_FORWARD, FFTW_ESTIMATE);
    fftw_plan plan_backward = fftw_plan_dft_1d(N, y, y, FFTW_BACKWARD, 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, y, plan_forward, plan_backward);
    
    // display the result
    cout << "Hilbert =" << endl;
    
    displayComplex(y);
    fftw_destroy_plan(plan_forward);
    fftw_destroy_plan(plan_backward);
   
}

Answer 2

For speed, you definitely want to reuse FFTW plans if possible. When reusing a plan across multiple calls to hilbert, remove the fftw_destroy_plan(plan) line within hilbert, otherwise the plan won't be valid to use in the next call. Instead, destroy the plan at the end of the program.

Edit: I had missed previously that you use two plans, one for forward transformation, another for backward. In hilbert(), remove all calls to fftw_plan_dft_1d, fftw_destroy_plan, and fftw_cleanup; the only FFTW call should be fftw_execute. Instead, create both plans up front once at the beginning of the program, and destroy them at the end of the program.

Answer 3

Prefaced by my top comments ...

Okay ... After a few false starts ...

Your second version had some issues.

Notably, you were not allocating out.

Also, to be safe, I believe that you need two fftw_plan structs, one for forward and one for backward.

These should be allocated/initialized once in main.

And, remove any allocation/destruction calls in hilbert. That way, it can be called multiple times by just changing the data placed in out.

The cleanup code can go to the bottom of main.

I've created a cleaned up and working version of your second version. I've shown the difference of old vs new code with cpp conditionals:

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

So, here it is:

#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, out, out, FFTW_FORWARD, 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, so we don't care about the multiplicands)
    memset(&out[hN + 1][REAL], 0, numRem * sizeof(fftw_complex));
    // creat a IDFT plan and execute it
#if 0
    plan_bwd = fftw_plan_dft_1d(N, out, out, FFTW_BACKWARD, 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, out, out, FFTW_FORWARD, FFTW_ESTIMATE);
    plan_bwd = fftw_plan_dft_1d(N, out, out, FFTW_BACKWARD, 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();
}

---

Here is the program output:

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