Hilbert Transform in Python?

Question

I am trying to write the Hilbert transform from scratch but not use any built in libraries other than for fft and ifft. I am not a mathematician by trade but I found these two algorithms online for the Hilbert transform, one in C and one in MATLAB. I tried to do implement them both but neither of them are giving me the same result as SciPy's Hilbert. I for sure have made something wrong in my implementation. Any insight would be greatly appreciated.

---

First implementation: (From MATLAB Website) Hilbert uses a four-step algorithm:

1. Calculate the FFT of the input sequence, storing the result in a vector x.

2. Create a vector h whose elements h(i) have the values:

   • 1 for i = 1, (n/2)+1

   • 2 for i = 2, 3, ... , (n/2)

   • 0 for i = (n/2)+2, ... , n

3. Calculate the element-wise product of x and h.

4. Calculate the inverse FFT of the sequence obtained in step 3 and returns the first n elements of the result.

My attempt:

def generate_array(n):
    a = np.hstack((np.full(n//2+1, 2), np.zeros(n//2-1)))
    a[[0, n//2]] = 1
    return a

def hilbert_from_scratch_2(u):
    fft_result = fft(u) #scipy fft

    n = len(u) 

    to_multiply = generate_array(n)

    result = np.multiply(n,to_multiply)

    return ifft(result) #scipy ifft

---

Second Implementation: (https://www.cfa.harvard.edu/~spaine/am/download/src/transform.c)

void hilbert(double *z, unsigned long n)
{
    double x;
    unsigned long i, n2;

    n2 = n << 1;
    /*
     * Compute the (bit-reversed) Fourier transform of z.
     */
    fft_dif(z, n);
    /*
     * Form the transform of the analytic sequence by zeroing
     * the transform for negative time, except for the (N/2)th.
     * element.  Since z is now in bit-reversed order, this means
     * zeroing every other complex element.  The array indices of
     * the elements to be zeroed are 6,7,10,11...etc. (The real
     * and imaginary parts of the (N/2)th element are in z[2] and
     * z[3], respectively.)
     */
    for (i = 6; i < n2; i += 4) {
        z[i] = 0.;
        z[i+1] = 0.;
    }
    /*
     * The 0th and (N/2)th elements get multiplied by 0.5.  Test
     * for the trivial 1-point transform, just in case.
     */
    z[0] *= 0.5;
    z[1] *= 0.5;
    if (n > 1) {
        z[2] *= 0.5;
        z[3] *= 0.5;
    }
    /*
     * Compute the inverse transform.
     */
    ifft_dit(z, n);
    /*
     * Normalize the array.  The factor of 2 is left over from
     * forming the transform in the time domain.
     */
    x = 2. / (double)n;
    for (i = 0; i < n2; ++i)
        z[i] *= x;
    return;
}   /* hilbert() */

My attempt:

def hilbert_from_scratch(signal):
    fast_ft = fft(signal) #scipy fft
    for i in range(6,len(signal),4):
        fast_ft[i] = 0
        fast_ft[i+1] = 0

    fast_ft[0] = fast_ft[0]*.5
    fast_ft[1] = fast_ft[1]*.5

    if(len(fast_ft) > 1):
        fast_ft[2] = fast_ft[2]*.5
        fast_ft[3] = fast_ft[3]*.5

    inverse_fft = ifft(fast_ft) #scipy ifft

    x = 2 / len(signal)

    for i in range(0,len(signal),1):
        inverse_fft[i] = inverse_fft[i]*x

    return inverse_fft

---

Any insight on why neither of them give the same result as SciPy's hilbert would be greatly appreciated.

Answer 1

I had a look at your code, made a couple edits, and compared it to the scipy and MATLAB Hilbert transforms. The function hilbert_from_scratch returns a complex sequence; the real components are the original signal and the complex components are the Hilbert transform. If you want just the Hilbert Transform, use np.imag on the returned array.

import math
from scipy.fftpack import *

def hilbert_from_scratch(u):
    # N : fft length
    # M : number of elements to zero out
    # U : DFT of u
    # v : IDFT of H(U)

    N = len(u)
    # take forward Fourier transform
    U = fft(u)
    M = N - N//2 - 1
    # zero out negative frequency components
    U[N//2+1:] = [0] * M
    # double fft energy except @ DC0
    U[1:N//2] = 2 * U[1:N//2]
    # take inverse Fourier transform
    v = ifft(U)
    return v


if __name__ == '__main__':
    N = 32
    f = 1
    dt = 1.0 / N
    y = []
    for n in range(N):
        x = 2*math.pi*f*dt*n
        y.append(math.sin(x))
    z1 = hilbert_from_scratch(y)
    z2 = hilbert(y)
    print(" n     y fromscratch scipy")
    for n in range(N):
        print('{:2d} {:+5.2f} {:+10.2f} {:+5.2f}'.format(n, y[n], z1[n], z2[n]))

Outputs:

 n     y fromscratch scipy
 0 +0.00 +0.00-1.00j +1.00
 1 +0.38 +0.38-0.92j +0.92
 2 +0.71 +0.71-0.71j +0.71
 3 +0.92 +0.92-0.38j +0.38
 4 +1.00 +1.00-0.00j +0.00
 5 +0.92 +0.92+0.38j -0.38
 6 +0.71 +0.71+0.71j -0.71
 7 +0.38 +0.38+0.92j -0.92
 8 +0.00 +0.00+1.00j -1.00
 9 -0.38 -0.38+0.92j -0.92
10 -0.71 -0.71+0.71j -0.71
11 -0.92 -0.92+0.38j -0.38
12 -1.00 -1.00+0.00j -0.00
13 -0.92 -0.92-0.38j +0.38
14 -0.71 -0.71-0.71j +0.71
15 -0.38 -0.38-0.92j +0.92

MATLAB:

>> y = sin(2*pi*linspace(0,1,17)); z = hilbert(y(1:end-1));
>> fprintf('%+2.2f %+2.2f\n',[real(z);imag(z)])
+0.00 -1.00
+0.38 -0.92
+0.71 -0.71
+0.92 -0.38
+1.00 -0.00
+0.92 +0.38
+0.71 +0.71
+0.38 +0.92
+0.00 +1.00
-0.38 +0.92
-0.71 +0.71
-0.92 +0.38
-1.00 +0.00
-0.92 -0.38
-0.71 -0.71
-0.38 -0.92