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.