c-我如何才能获得与Matlab使用FFTW相同的希尔伯特变换结果

Question

我需要计算希尔伯特变换，并从我的信号中获得它的绝对数。我安装了FFTW，并遵循了这些YouTube教程，两者都很好。

然而，我像视频中那样实现了希尔伯特变换，我没有获得与Matlab计算的值相同的值。

我读过一些人通过使用FFT和其他计算来解决这个问题，但他们没有给出足够清楚的答案。有人能给我提供几行基于FFTW的代码吗？这样可以让结果与Matlab的计算结果相同吗？

这是希尔伯特变换代码，我跟随视频：

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

我要计算的主要程序：

pi16Buffer = (int16 *)pBuffer;
//int16* ch1Buffer; 
//int16* ch2Buffer;

double* ch1Buffer = NULL;
double* ch2Buffer = NULL;
fftw_complex result[N] ; // output array

// Allocate size to ch1 and ch2
ch1Buffer       = (double*)calloc(u32Size, sizeof(double));
ch2Buffer       = (double*)calloc(u32Size, sizeof(double));

//ch1ch2ch1ch2... fulfill the buffer
for (i = 0; i < u32Size/2; i++)
{
ch1Buffer[i] += (double)pi16Buffer[i*2];
ch2Buffer[i] += (double)pi16Buffer[i * 2 + 1];
}
// here hilbert on the whole ch2
hilbert(ch2Buffer, result); //hilbert(inputarray, outputarray)
for (i = 0; i < u32Size; i++)
{
if (ch1Buffer[i] > max1)  //Find max value in each segs of ch1 and ch2
max1 = ch1Buffer[i];
if (abs(result[i]) > max2)
max2 = abs(result[i]); // Calculate the absolute value of hilbert result;

}
Corrected = max2 / max1; //do the signal correction
free(ch1Buffer); //free buffer
free(ch2Buffer);


}
return Corrected;

Answer 1

我很难理解您的代码
这是我的测试代码，针对一个简单的N=4案例
它应该至少适用于所有偶数大小的输入。用奇数输入进行检查。

此代码计算分析信号。matlab的hilbert函数也做了同样的操作。

它的实数部分对应于输入的实数信号
希尔伯特变换对应于虚值。

其原理是确保其负频率的FFT值为零
这是通过在频域中进行简单的开窗来获得的。

#include <stdio.h>
#include <complex.h>
// Analytic signal calculation
// The real part of it corresponds to the input real signal
// The hilbert transform corresponds to the imaginary value
void print (complex *x, int n) {
for (int i = 0; i < n; ++i) {
printf ("(%lf, %lf) ", creal(x[i]), cimag(x[i]));
}
printf ("n");
}
void fft4 (complex *x, int n, int inversion) {
complex t0, t1, t2, t3;
t0 = x[0] + x[2];
t1 = x[1] + x[3];
t2 = x[0] - x[2];
if (!inversion) t3 = I * (x[1] - x[3]);
else t3 = -I * (x[1] - x[3]);
x[0] = t0 + t1;
x[2] = t0 - t1;
x[1] = t2 - t3;
x[3] = t2 + t3;
}
#define N 4
int main() {
const int n = N;
complex x[N] = {1, 2, 3, 4};
complex y[N] = {1, 2, 3, 4};

fft4 (y, n, 0); // direct FFT size 4

print (x, n);
print (y, n);

for (int i = 1; i < n/2; ++i) y[i] *= 2.0;
for (int i = n/2+1; i < n; ++i) y[i] = 0.0;

fft4 (y, n, 1);     // inverse FFT
for (int i = 0; i < n; ++i) y[i] /= n;  // normalisation
print (y, n);       // print analytical signal
return 0;
}