How could I make FFTW Hilbert transform calculate faster? [closed]

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/30 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



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);
   
}

The exact output value

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

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);
 
}
```

Answer 1

Updated Code Review for 2020-12-30 Code

Some of the original code review items have been addressed. Good work! Here are some remaining ideas:

• I assume your purpose here is to make the hilbert function faster, because you will call it multiple times. You've moved the FFT plan creation outside of hilbert, which is probably a huge speed increase. However, you still call fftw_destroy_plan within hilbert. This is a problem if you call hilbert more than once. You should move fftw_destroy_plan outside of hilbert, so it would be called from main instead.
• You have three different for loops inside hilbert. One copies data from input vector to output vector. One scales half the output vector by a factor of 2. One scales the output vector by 1/N. These loops can all be combined into a single loop at the beginning of hilbert. This works because an FFT is linear, so scaling the input by a constant factor is the same as scaling the output by a constant factor. This new loop would look something like:

      auto scale = 2.0 / N;
      for (int i = 0; i < N; ++i)
      {
          out[i][REAL] = in[i] * scale;
          out[i][IMAG] = 0.0;
      }
  

  Then you would get rid of the other for loops. You would also get rid of the code which multiplies out[hN] by 2 (but keep the call to memset which zeros the second half of the output). Finally, you need to fix out[0], which is now too large by a factor of two:

      out[0][REAL] /= 2.0;
      out[0][IMAG] /= 2.0;
  

  However, the speed increase from combining these loops is probably small compared to the FFT time.
• It may take a little more work, but you could change the forward FFT so that it takes its input directly from the real-valued input vector. This is probably close to twice as fast as the complex FFT you are currently doing (but this is just for the forward FFT, so your overall speed increase might be up to 25%). Still, I think it is likely that your complex FFT is already fast enough now that you've removed the creation/destruction of the plan outside of hilbert.

Original Code Review for Original 2020-12-26 Code

• In general, using global variables is a bad idea (see below for my reasoning for this statement). You declare a global variable out. It is never used, because your hilbert function takes a parameter named out. So just delete the global out variable.

• You also declare the FFT plan as a global variable. This plan is what gets used by the hilbert function. But when you initialize the plan in main, you declare a local plan variable. The local plan variable gets initialized by fft_plan_dft_1d, but that doesn't initialize the global variable. The local plan variable is never used. Better would be to get rid of the global plan variable, and use only the one declared in main. Pass it as a parameter into the hilbert function. (This may be one cause of your exception, since you're currently using an uninitialized plan variable inside hilbert.)

• You only have one plan variable, but it looks like you're trying to make it hold two plans. That won't work. It will remember only the last call to fft_plan_dft_1d, so both times you call fft_execute, you are getting an inverse FFT. And when you call fft_destroy_plan twice, the second one may cause an exception.

  To fix this, you'd want two different plan variables, maybe call them plan_forward and plan_inverse or something like that. Pass them both as parameters into the hilbert function.

• Having using namespace std; is generally not a good idea. Instead, when you want to use something from that namespace, put std:: in front of the name. For example, use std::cout instead of cout.

• You have int hN = N >> 1;. I personally would prefer int hN = N / 2, as it seems clearer.

• You copy the input array into the output array and then do an in-place FFT. I think it should be possible to have fftw do a non-in-place FFT where the input vector is real-valued and the output is complex. This should be slightly more efficient. But I'm too lazy to look up the details of how to do it.

Thoughts on Global Variables

My original review just said "In general, using global variables is a bad idea", without any explanation or context. Some of the comments on this answer either disagreed or wanted more explanation, so here's my reasoning.

As with most things related to coding style, reasonable people may have differing opinions. Furthermore, nothing is absolute, so I certainly acknowledge that there are cases where global variables are useful or even necessary. Perhaps I should have said "use global variables only when really necessary".

• What do commonly-used coding guidelines say?

  There are lots of coding guidelines on the internet. Many of them in some way discourage the use of global variables. Some examples:

  • The C++ Core Guidelines, maintained by Bjarne Stroustrup and Herb Sutter, say to avoid non-const global variables.
  • The Google C++ Style Guide forbids objects with static storage duration unless they are trivially destructable, and also discourages static storage duration if the object uses dynamic initialization. (Note: global variables have static storage duration.)
  • In Joint Strike Fighter C++ Coding Standards, based on the MISRA C++ coding standards, AV Rule 98 says "Every nonlocal name, except main(), should be placed in some namespace." (This doesn't prohibit global variables, but does help avoid some of the problems caused by global variables.)
  • In University of Michigan C++ Coding Standards, global variables are acceptable only in certain limited circumstances.
  • Some random code guideline found on the internet discourages use of global variables, suggesting that you instead use variables declared in functions.

• What can go wrong when you use global variables

  • Initialization issues. If global variables are declared in two different files, there is no specification of the order in which they are initialized. If the initialization of one of them depends on the value of the other, you have undefined behavior. It may appear to work sometimes, but it can't be relied upon.
  • Multi-thread issues. If multiple threads simultaneously try to use or modify a global variable, you can get unexpected or undefined behavior. Making this work correctly can be tricky, sometimes requiring the use of mutexes or similar constructs. Often, these constructs slow down the code substantially.
  • Name conflicts within a file. You can have a function-local variable and a global variable with the same name. They are different variables, and this can lead to confusion. Maybe you meant the local variable to access the global variable, but accidentally declared it locally. Or maybe you meant for the local variable to be different than the global variable, but a subsequent maintainer may not realize that and think they are the same variable. These kind of issues can lead to bugs.
  • Name conflicts with external libraries. If you have an external global variable, it can conflict with other external global variables with the same name. These kind of conflicts can make it difficult to integrate separate software libraries.

In many cases, global variables offer no benefit in exchange for the above potential problems (the original post is one example of that). Therefore, I think it is best to avoid global variables unless they are really necessary.