6
$\begingroup$

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: Raw signal.

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

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 implementation

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?

$\endgroup$
  • $\begingroup$ It seems like there some kind of low-pass filtering in the Matlab code? Care to show us your Matlab code? $\endgroup$ – Ben Jan 10 '18 at 20:58
  • $\begingroup$ @Ben I am just using the default matlab code hilbert(), I dont think I am using any low pass filter. The edits contain the matlab code $\endgroup$ – Spandyie Jan 10 '18 at 21:54
  • 2
    $\begingroup$ For starters, doesn't the hilbert function in Matlab return complex values? If so, I would expect that you plot the module of the hilbert transform and not the complex values. $\endgroup$ – Ben Jan 10 '18 at 22:16
  • $\begingroup$ @Ben Its absolute of hilbert. There was a typo in matlab code $\endgroup$ – Spandyie Jan 10 '18 at 22:25
  • $\begingroup$ It still seems like you have some kind of smoothing in Matlab. Why not do it with a simple sine wave in both Python and Matlab ? You should get a constant enveloppe $\endgroup$ – Ben Jan 11 '18 at 14:27
5
$\begingroup$

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)

enter image description here

enter image description here

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

$\endgroup$
  • $\begingroup$ could you share your signal.txt, because I am having the same problem as OP. $\endgroup$ – Dr Sokoban Jul 25 '19 at 11:03
  • $\begingroup$ @DrSokoban The link to the signal is in the original question $\endgroup$ – endolith Jul 25 '19 at 14:31
2
$\begingroup$

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.

$\endgroup$
1
$\begingroup$

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))."

$\endgroup$
0
$\begingroup$

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()
$\endgroup$
  • 1
    $\begingroup$ There is no lowpass filter in the matlab code. It's just the magnitude of the analytical signal, which should produce an envelope in both cases. $\endgroup$ – endolith Apr 9 '18 at 19:39
0
$\begingroup$

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)

enter image description here

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

enter image description here

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

enter image description here

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...

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.