如何获得与 Matlab 使用 FFTW 相同的希尔伯特变换结果？

Question

我需要计算希尔伯特变换并从我的信号中得到它的绝对数。 我安装了 FFTW 并关注了these YouTube tutorials，都运行良好。

但是，我按照视频的方式实现了希尔伯特变换，我没有得到与 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;
}