Why is Scipy implementation of Hilbert() function different from Matlab implementation of the function?

Question

I am trying to fit Hilbert envelop to a high frequency ultrasonic signal of frequency 250 KHZ and sampling rate 12000000. Raw signal looks like below: .

I used hilbert() function from scipy.signal package in python this is what it looks like .

The python code looks like below

from scipy.signal import hilbert
import numpy as np
 def Hilbert(self,i=0):
    analytical_signal = hilbert(self.sensor["s"+str(i)])
    amplitude_envelope = np.abs(analytical_signal)
    return amplitude_envelope

The Matlab implementation looks like this

Matlab code is as follows:

figure;
plot(abs(hilbert(signal)),'r');
hold on; 
plot(signal,'b');

The line data is as follows Signal. I am wondering which one is correct?

Answer 1

It works fine for me:

from scipy.signal import hilbert
import numpy as np
from matplotlib.pyplot import plot

sensor = np.loadtxt('signal.txt')
plot(sensor)

analytical_signal = hilbert(sensor)
plot(analytical_signal.real)
plot(analytical_signal.imag)

amplitude_envelope = np.abs(analytical_signal)
plot(amplitude_envelope)

What are you doing differently? Maybe you're throwing away the imaginary part of analytical_signal somehow?

Answer 2

You have to detrend the signal. If you have a DC component the behauviour changes completely.

dt=0.0001
t = np.arange(0,0.1,dt)
x = (1+np.cos(2*np.pi*50*t))*np.cos(2*np.pi*1000*t)
plt.plot(x)

h=abs(hilbert(x))    
plt.plot(h)

x=x+1
h=abs(hilbert(x))
plt.plot(h)

This is because hilbert(x) returns de analytical function xr(t)+jxh(t), where xh is the Hilbert's Transform and xr is x(t), the original signal. xh is the same for x(t) than for x'=x(t)+k (Hilbert's transform of a constant is zero). So, when you are calculating the envelope as abs(hilbert(xr)) what you get sqrt(xr^2+xh^2). If you use x'=x(t)+k you get sqrt((xr+k)^2+xh^2).

In Matlab, the envelope() function substracts the mean before applying hibert(). If you use abs(hilbert()) you should obtain the same result than with Python.

PD: After writting this I download your signal and I test the abs(hilbert()) on it. In Python I obtain the envelope...

Answer 3

I had the same issue. As https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.hilbert.html says @Notes, you need to need to take the imaginary part of the hilbet transform result... "The Hilbert transformed signal can be obtained from np.imag(hilbert(x)), and the original signal from np.real(hilbert(x))."

Answer 4

I don't know python. But theoretically,

Hilbert transformation is done by:

1. Real part of the signal
2. Rotating the phase of the signal by 90°
3. Analytical signal = real + i*(rotated signal).
4. Envelope is a distance function. It's the distance between the center of the analytic signal to the amplitude of the sample.
5. Instantaneous frequency is the angle.

So, I guess from the plot that you forget to consider the imaginary part of the signal. That's why you don't see the envelope but the absolute value of the real signal.

Answer 5

You should apply a LPF after applying the hilbert method. Here is an example.

import numpy
from scipy.signal import butter, filtfilt, hilbert
import matplotlib.pyplot as plt

def FilteredSignal(signal, fs, cutoff):
    B, A = butter(1, cutoff / (fs / 2), btype='low')
    filtered_signal = filtfilt(B, A, signal, axis=0)
    return filtered_signal

fs = 10000.
T = .1
time = numpy.arange(0., T, 1 / fs)
frequency = 1000
noise = numpy.random.normal(0, 1, int(fs/10))
signal = numpy.sin(2 * numpy.pi * frequency * time) + 0.2 * noise
analyticSignal = hilbert(signal)
amplitudeEvelope = numpy.abs(analyticSignal)
cutoff = 1000
filteredSignal = FilteredSignal(amplitudeEvelope, fs, cutoff)

fig2, ax2 = plt.subplots(1, 1)
ax2.plot(time, signal)
ax2.plot(time, filteredSignal)
ax2.set_xlabel('Tiempo')
ax2.set_ylabel('Amplitud')
ax2.grid(True)
plt.show()