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
```

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.

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

answered on Stack Overflow Jan 6, 2021 by Kevin

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
```

User contributions licensed under CC BY-SA 3.0