Addressing Hilbert transform edge effects

Question

I'm experimenting with the Hilbert transform in Python and just now understanding how potentially severe the edge effects are. Here is the code I'm using:

import numpy as np
import scipy.signal
import matplotlib.pyplot as plt
option = 1
if option == 1:
    x = np.sin(np.arange(100))
elif option == 2:
    x = np.sin(np.arange(100)*np.pi/10)
elif option == 3:
    x = np.sin(np.arange(100)/4)
hilbert_data = scipy.signal.hilbert(x, axis=0)
real = hilbert_data.real
imag = hilbert_data.imag
amp_data = np.absolute(hilbert_data)
phase_data = np.unwrap(np.angle(hilbert_data), axis=0)
freq_data = np.diff(phase_data, axis=0)
fig, axs = plt.subplots(2, 2, figsize=(12,5))
axs[0, 0].plot(x)
axs[0, 0].set_title('Signal')
axs[0, 1].plot(amp_data, 'tab:orange')
axs[0, 1].set_title('Amplitude')
axs[0, 1].set_ylim([0.8,1.2])
axs[1, 0].plot(phase_data, 'tab:green')
axs[1, 0].set_title('Phase')
axs[1, 1].plot(freq_data, 'tab:red')
axs[1, 1].set_title('Frequency')

These are the outputs when I run it with options 1, 2, and 3, respectively:

---

---

What's noticeable is how different the edge effects are in each example. With option 2, there are essentially no edge effects. With option 1, the edge effects on amplitude are oscillatory and frequency spikes up at the ends. With option 3, it's the opposite - amplitude spikes and frequency is oscillatory.

Assuming that there is reason to believe that a signal is relatively consistent (e.g., the example above is based on signal with uniform amplitude/frequency over time), is there a way to use the outputs in the middle to correct the outputs on the ends? In the example above, it might make sense to just take the amplitude/frequency from the the centermost point, however, for real signals, that would likely be inaccurate. The difficulty is that the edge effects appear to be highly sensitive to small changes in the input as with the above.

Are there ways to determine the instantaneous amplitude, frequency, and phase other than the Hilbert transform? Any thoughts are appreciated.

Answer 1

The Hilbert Transform as implemented is a filter with a long time domain response (and on its own is non-causal), so requires time history and a matching delay in order to produce a result. Given one set of samples, the Hilbert Transform and the Analytic Signal (the signal as the real component and the Hilbert Transform as the imaginary component), there is no way that it can produce an instantaneous result of the signal's envelope based on just that one sample. Further the more history we have, the closer we can approximate a true Hilbert impulse response and then improve the accuracy of the result, but using this to predict the instantaneous phase/freq/amplitude would require a consistency on those parameters over the memory of the filter (which then puts constraints on the total bandwidth of the signal).

Further it is to be noted that the hilbert function in MATLAB/Octave and Python's scipy.signal are due to the simplistic (and fast!) approach of simply nulling the negative frequency axis bins (the upper half of the FFT bins) that those algorithms use. That approach of determining the time domain impulse response for filters from the IFFT is called the "Frequency Sampling" method of filter design. It is the most intuitive when a specific frequency domain response is desired (just IFFT the response), but it also has the worst performance compared to other methods such as windowing or using the optimized algorithms known as "least-squares" and "equiripple" (Parks-McLellan) that are used in the 'firls' and 'firpm' or 'remez' commands; this is because the Frequency Sampling method will result in a response that is exact at the bin frequencies used, but have more error compared to the other methods at all other frequencies that are in between the FFT bins, and more error in the time domain responses as well. However this will not completely resolve the time requirement issue for the OP which will occur when using the Hilbert; an averaging interval is required before the first "valid" sample can be used.

To minimize errors, consider creating the Hilbert Transform based on windowing the ideal time domain Discrete Hilbert impulse response. Alternatively the Hilbert can also be derived from the coefficients of a half-band filter which are determined using the optimized algorithms, and the python `'remez' function includes an option for implementing the Hilbert directly.

Below are the results of recent simulations I did directly comparing the frequency response achievable with both odd and even samples (using 91 and 92 length samples) for the hilbert function compared to hilbert filters created using Kaiser windowed ideal Discrete Hilbert impulse response.

The following shows the resulting amplitude ripple for the frequency response of a Hilbert impulse response created with 91 and 92 samples using the hilbert function. Notice the significant difference between using even and odd samples. This is a direct effect of the time domain aliasing with the zeroing of FFT bins approach the hilbert function uses.

For an apples to apples comparison, below is the frequency response for the better case of 92 samples for the hilbert derived response compared to using the Kaiser windowed Discrete Hilbert Impulse response with the same number of samples and overall frequency coverage:

For Hilbert filters, odd length impulse responses may be desired since as linear phase solutions will have an integer sample delay.

Below shows a comparison of what can be achieved with different filter design approaches where each has been matched to the same frequency coverage (which is narrower than the hilbert comparison given above, but includes another Kaiser windowed solution):

To determine the best estimate of instantaneous amplitude, frequency and phase; consider instead sampling the waveform in quadrature with two ADCs rather than creating the quadrature channel from one waveform, and correcting for the time offset between the samples as I detail in this post.