我怎样才能使 FFTW Hilbert 变换计算得更快？
How could I make FFTW Hilbert transform calculate faster?
我正在使用 FFTW 源中的希尔伯特变换函数。
因为我将在我的 DAQ 中以数据流模式使用它。
函数运行正常，但计算速度慢，会导致FIFO溢出。
我听说将 fftw_plan() 从 hilbert() 移到外面以重用 plan 可能会有用，但是，一旦我这样做，它就会出错，在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, 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);
}
原代码如下：
#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);
}
准确输出
Hilbert =
1+3.82843i
2-1i
3-1i
4-1.82843i
5-1.82843i
6-1i
7-1i
8+3.82843i
以我的热门评论开头...
好吧...经过几次错误的开始...
你的第二个版本有一些问题。
值得注意的是，您没有分配out。
此外，为了安全起见，我认为您需要 两个 fftw_plan 结构，一个用于前向，一个用于后向。
这些应该是 allocated/initialized 一次 在 main.
并且，删除 hilbert 中的任何 allocation/destruction 个调用。这样，只需更改 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, 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();
}
程序输出如下：
Hilbert =
0.125+0i
0.5+0i
0.75+0i
1+0i
0.625+0i
0+0i
0+0i
0+0i
0.000171
为了速度，如果可能的话，您绝对希望重用 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, 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);
}
我正在使用 FFTW 源中的希尔伯特变换函数。
因为我将在我的 DAQ 中以数据流模式使用它。
函数运行正常，但计算速度慢，会导致FIFO溢出。
我听说将 fftw_plan() 从 hilbert() 移到外面以重用 plan 可能会有用，但是，一旦我这样做，它就会出错，在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, 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);
}
原代码如下：
#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);
}
准确输出
Hilbert =
1+3.82843i
2-1i
3-1i
4-1.82843i
5-1.82843i
6-1i
7-1i
8+3.82843i
以我的热门评论开头...
好吧...经过几次错误的开始...
你的第二个版本有一些问题。
值得注意的是，您没有分配out。
此外，为了安全起见，我认为您需要 两个 fftw_plan 结构，一个用于前向，一个用于后向。
这些应该是 allocated/initialized 一次 在 main.
并且，删除 hilbert 中的任何 allocation/destruction 个调用。这样，只需更改 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, 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();
}
程序输出如下：
Hilbert =
0.125+0i
0.5+0i
0.75+0i
1+0i
0.625+0i
0+0i
0+0i
0+0i
0.000171
为了速度，如果可能的话，您绝对希望重用 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, 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);
}