How could I obtain the same Hilbert transform result as Matlab does with FFTW?

Question

I need to calculate the Hilbert transform and get its absolute number from my signal. I installed the FFTW and followed these YouTube tutorials, both worked good.

However, I implemented the Hilbert Transform as the video did, I didn't obtain the value as same as Matlab calculated.

I've read some of the people who solved this problem by just using the FFT and doing other calculations, but they didn't give the answer clear enough. Could anyone provide me several lines of code which based on the FFTW so that can let result become same to Matlab calculated result?

Here is the Hilbert transform code which I followed the video:

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

My main program to calculate:

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

I have some difficulty to understand your code.
Here is my test code, on a simple N=4 case.
It should work at least for all even sizes of the input. To be checked with odd input.

This code calculates the analytic signal. The hilbert function of matlab does the same.

The real part of it corresponds to the input real signal.
The Hilbert transform corresponds to the imaginary value.

The principle is to insure that the value of its FFT for negative frequencies is zero.
This is obtained by a simple windowing in the frequency domain.

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