0
$\begingroup$

I am playing around with the Hilbert transform using Python's scipy.signal.hilbert package to better understand what the various instantaneous properties are given different input signals. One input signal I analyzed was a 10Hz sine wave. Here is the code which generates this signal:

    self.duration = 10
    self.fs = 1600
    samples = int(self.fs * self.duration)
    self.time = np.arange(samples) / self.fs
    self.signal = np.sin(10 * 2.0 * np.pi * self.time)

I go ahead and generate the Hilbert properties (instantaneous phase, frequency, amplitude, regenerated carrier wave):

demean = self.signal - np.mean(self.signal)
hil = hilbert(demean)
inst_amplitude = np.abs(hil)
inst_phase = np.unwrap(np.angle(hil))
inst_frequency = np.diff(inst_phase) / (2 * np.pi) * sampling_frequency
inst_frequency = np.insert(inst_frequency, 0, 0)
regen_carrier = np.cos(inst_phase)

Then I plot these properties: enter image description here

Why does instantaneous amplitude increase with respect to time? I would think it should remain constant?

If I change the duration of my signal to, say, 40sec, the instantaneous amplitude also increases with respect to time. Between the two plots, it appears that the instantaneous amplitude increases in time increments of a factor of 2. So it increases at 0.25sec, 0.5sec, 1sec, 2sec, 4sec, 8sec, 16sec, 32sec, etc. Why is this happening? The code is so simple to produce instantaneous amplitude, too, so I don't understand where I am going wrong.

enter image description here

Here's the scipy.signal.hilbert code:

def hilbert(x, N=None, axis=-1):
"""
Compute the analytic signal, using the Hilbert transform.

The transformation is done along the last axis by default.

Parameters
----------
x : array_like
    Signal data.  Must be real.
N : int, optional
    Number of Fourier components.  Default: ``x.shape[axis]``
axis : int, optional
    Axis along which to do the transformation.  Default: -1.

Returns
-------
xa : ndarray
    Analytic signal of `x`, of each 1-D array along `axis`

Notes
-----
The analytic signal ``x_a(t)`` of signal ``x(t)`` is:

.. math:: x_a = F^{-1}(F(x) 2U) = x + i y

where `F` is the Fourier transform, `U` the unit step function,
and `y` the Hilbert transform of `x`. [1]_

In other words, the negative half of the frequency spectrum is zeroed
out, turning the real-valued signal into a complex signal.  The Hilbert
transformed signal can be obtained from ``np.imag(hilbert(x))``, and the
original signal from ``np.real(hilbert(x))``.

Examples
---------
In this example we use the Hilbert transform to determine the amplitude
envelope and instantaneous frequency of an amplitude-modulated signal.

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from scipy.signal import hilbert, chirp

>>> duration = 1.0
>>> fs = 400.0
>>> samples = int(fs*duration)
>>> t = np.arange(samples) / fs

We create a chirp of which the frequency increases from 20 Hz to 100 Hz and
apply an amplitude modulation.

>>> signal = chirp(t, 20.0, t[-1], 100.0)
>>> signal *= (1.0 + 0.5 * np.sin(2.0*np.pi*3.0*t) )

The amplitude envelope is given by magnitude of the analytic signal. The
instantaneous frequency can be obtained by differentiating the
instantaneous phase in respect to time. The instantaneous phase corresponds
to the phase angle of the analytic signal.

>>> analytic_signal = hilbert(signal)
>>> amplitude_envelope = np.abs(analytic_signal)
>>> instantaneous_phase = np.unwrap(np.angle(analytic_signal))
>>> instantaneous_frequency = (np.diff(instantaneous_phase) /
...                            (2.0*np.pi) * fs)

>>> fig = plt.figure()
>>> ax0 = fig.add_subplot(211)
>>> ax0.plot(t, signal, label='signal')
>>> ax0.plot(t, amplitude_envelope, label='envelope')
>>> ax0.set_xlabel("time in seconds")
>>> ax0.legend()
>>> ax1 = fig.add_subplot(212)
>>> ax1.plot(t[1:], instantaneous_frequency)
>>> ax1.set_xlabel("time in seconds")
>>> ax1.set_ylim(0.0, 120.0)

References
----------
.. [1] Wikipedia, "Analytic signal".
       https://en.wikipedia.org/wiki/Analytic_signal
.. [2] Leon Cohen, "Time-Frequency Analysis", 1995. Chapter 2.
.. [3] Alan V. Oppenheim, Ronald W. Schafer. Discrete-Time Signal
       Processing, Third Edition, 2009. Chapter 12.
       ISBN 13: 978-1292-02572-8

"""
x = np.asarray(x)
if np.iscomplexobj(x):
    raise ValueError("x must be real.")
if N is None:
    N = x.shape[axis]
if N <= 0:
    raise ValueError("N must be positive.")

Xf = sp_fft.fft(x, N, axis=axis)
h = np.zeros(N)
if N % 2 == 0:
    h[0] = h[N // 2] = 1
    h[1:N // 2] = 2
else:
    h[0] = 1
    h[1:(N + 1) // 2] = 2

if x.ndim > 1:
    ind = [np.newaxis] * x.ndim
    ind[axis] = slice(None)
    h = h[tuple(ind)]
x = sp_fft.ifft(Xf * h, axis=axis)
return x
$\endgroup$
2
  • 1
    $\begingroup$ Did you notice the scaling is 1E-31? It appears to be floating point error but is essentially amplitude 1 throughout. $\endgroup$ – Dan Boschen Mar 2 at 23:39
  • $\begingroup$ Yes, I did notice this. What would cause the floating point error to change like this in multiples of 2? Is this a scipy.signal.hilbert thing? or more of a computer thing? $\endgroup$ – curious individual Mar 3 at 14:58

Your Answer

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

Browse other questions tagged or ask your own question.