5
$\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
3
$\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 at 11:03
  • $\begingroup$ @DrSokoban The link to the signal is in the original question $\endgroup$ – endolith Jul 25 at 14:31
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

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.