【发布时间】：2020-12-22 21:11:40
【问题描述】：
我正在使用来自 FFTW 源的希尔伯特变换函数。
因为我打算在数据流模式下在我的 DAQ 中使用它。
该函数工作正常，但计算速度慢，会导致 FIFO 溢出。
我听说将fftw_plan() 从hilbert() 移到外部以重用plan 可能很有用，但是，一旦我这样做，就会出错，在fftw_destroy_plan(plan); 上说Exception thrown at 0x0000000070769240 (libfftw3-3.dll) in CppFFTW.exe: 0xC0000005: Access violation reading location 0x0000000000000000.。有没有人有类似的经验或更好的解决方案来提高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
【问题讨论】：
-
也许您可以通过在您的陈述中将问题重新集中在这一部分上来回答这个问题：“但是，一旦我这样做了，它就会出现问题。”。 FFTW 的文档确切地说要重新使用 fftw_plan。那么，您尝试过什么？您看到的“错误”效果是什么？
-
已编辑，之所以这么说是因为我不确定我的方法是否正确。
-
为什么要覆盖一个干净的结构？代替
out[i][REAL]，执行：out[i].re。我什至不确定它们是否等效 [yours might be UB] 并且.re是结构的定义方式和 lib 的期望。同样适用于out[i].im -
改几个字母有用吗...？
-
out是一个指针，因此通过out[i][anything]访问它是错误的。我很惊讶它编译了，因为没有办法告诉编译器行数。它 is UB 并且您正在访问 超出 数组的末尾。当你这样做时：out[i][1]你不设置out[i]的.im部分。你在破坏out[i + 1]。最终，这会导致您在 [stack] 数组的末尾乱涂乱画——这就是 UB [未定义的行为]。
标签： c++ signal-processing fftw