Discrepancies between FFT-based Hilbert Transform and FIR filter results

Question

I compared the results of Hilbert Transform as provided in scipy.signal.hilbert() (which is uses I/FFT as you can see in the source code) with its approximation implemented as FIR filter using coeffs that I found in this paper. I constructed a quick example where the input signal x0 is a sine with a linearly decreasing frequency w0. The Hilbert transform h is calculated using the scipy function[1], the quadrature Q is calculated using the FIR filter with coeffs from the code in paper [3] and the in-phase component I is a copy of the input signal - it should be the input signal delayed by 3, I know, but for the sake of the example I wanted to remove the lag, because I wanted to see how well the reconstructed signal x2 is matching the original x0. Further the amplitude and phase series A1, ph1 where computed from hand A2, ph2 where computed from I and Q. Reconstructed signals x1 and x2 were computed from the corresponding amplitude and phase series, then finally the instantaneous frequencies w1 and w2 were calculated from ph1 and ph2 respectively. Here is the Python code with all the calculations and plotting:

import numpy as np
from scipy.signal import lfilter
import matplotlib.pyplot as plt

# Construct the input signal
idx = np.arange(512)
w0 = [2*np.pi/(80 + x/8) for x in idx]
x0 = np.sin(w0*idx)
# Calculate the transforms
h = hilbert(x0)
I = x0 # normally it should be np.roll(x0, 3)
Q = lfilter([0.25, 0, 0.75, 0, -0.25, 0, -0.75], [1], x0)
# Calculate amplitude arrays
A1 = np.abs(h)
A2 = np.sqrt(I**2 + Q**2)
# Calcuclate phase arrays
ph1 = np.angle(h)
ph2 = -np.arctan2(Q, I)
# Calculate instant. frequency
w1 = np.diff(np.unwrap(ph1))
w2 = np.diff(np.unwrap(ph2))

# Plot the results
fig = plt.figure()
ax = plt.subplot(311)
plt.plot(x0, lw=3, color='b', label='signal')
plt.plot(A1*np.cos(ph1), lw=2, color='g', label='x1')
plt.plot(A1, lw=2, ls='dotted', color='g', label='A1')
plt.plot(A2*np.cos(ph2), lw=1, color='k', label='x2')
plt.plot(A2, lw=1, ls='dotted', color='k', label='A2')
plt.legend(bbox_to_anchor=(1, 1), loc=2, borderaxespad=0.)
plt.subplot(312, sharex=ax)
plt.plot(ph1, lw=2, color='g', label='ph1')
plt.plot(ph2, lw=1, color='k', label='ph2')
plt.legend(bbox_to_anchor=(1, 1), loc=2, borderaxespad=0.)
plt.subplot(313, sharex=ax)
plt.plot(w0, lw=3, color='b', label='w0')
plt.plot(w1, lw=2, color='g', label='w1')
plt.plot(w2, lw=1, color='k', label='w2')
plt.legend(bbox_to_anchor=(1, 1), loc=2, borderaxespad=0.)
fig.show()

And this is how the figure looks like:

The reconstructed signals x1 and x2 seem to be tightly matching the input signal despite of the discrepancies between the amplitudes A1 and A2, and phases ph1 and ph2. The latter further affect the instantaneous frequencies w1 and w2. The one calculated from scipy.signal.hilbert() is at least a fairly straight line (displaced by an offset which should probably be included in the calculations at some point), but w2 has sudden peaks in places where A2 has dips and it does not follow w0 that well.

Hilbert Transform is widely used for instantaneous frequency estimation, so considering the above example I have been wondering whether there are some steps or details missing from my calculations or is that how accurate the approximation should be? I have been using scipy.signal.hilbert for some time and am aware of the edge effects, but the FIR method seems to be out of line throughout. Could that be due to the truncated coeffs? I did not try longer filters, because I do not know exactly how to calculate the coeffs, but the paper [3] promised that the proposed approximation was good enough to be used to successfully predict turns in gold prices during high volatility periods, so I took the authors' word for that. I actually saw similar phase shifting filters in other books or papers and one thing that really surprised me was that every time authors used exactly the same arrays of coeffs for detrending filters as they used to calculate the quadrature component ([3] is no exception: detrender and Q1 are results of exactly the same filters). I am still new to DSP, so if I am missing something important here please let me know.

Answer 1

There are several reasons why the two results don't match:

1. the coefficients of the FIR Hilbert transformer are wrong
2. the FIR Hilbert transformer is too short to even come close to the performance of the FFT-based implementation
3. the frequency of the input signal is too low for the FIR Hilbert transformer to perform properly. A FIR Hilbert transformer always performs poorly at very low frequencies; you can shift the frequency range towards low frequencies by increasing the filter length, but there will always be a region of bad performance around DC.
4. the $I$ and $Q$ components of the analytic signal are not aligned due to the delay of the FIR Hilbert transformer

Let's address the different points above one by one. First of all, the coefficients of an ideal Hilbert transformer are given by

$$h[n]=\begin{cases}\frac{2}{n\pi},& n\textrm{ odd}\\0,&n\textrm{ even}\end{cases}\tag{1}$$

From $(1)$ you see that the coefficients must be anti-symmetric, which yours aren't; and also the scaling is wrong. The scaling is important because you combine the Hilbert transformed signal with the original signal to obtain the analytic signal.

A good approximation to the ideal Hilbert transform can be achieved by (symmetrically) truncating the coefficients $(1)$, and windowing them to minimize the Gibbs phenomenon. The necessary number of coefficients must be determined by trial and error for the given application.

For the given input signal you would need a very large number of coefficients due to its low frequency. Increasing the (average) frequency for a given number of coefficients will drastically improve the performance.

The $I$ and $Q$ components of the analytic signal can be aligned by throwing away the first $(N-1)/2$ samples of $Q$, where $N$ is the number of coefficients (including the zeros), assuming that $N$ is odd, which is the case if the coefficients $(1)$ are truncated symmetrically.

The plots below were created using Octave. I assume that the implementations of the command hilbert are similar in Octave and scipy. The red curves are the results from using hilbert, the blue ones from using the FIR Hilbert transformer. The top plot shows the amplitude estimates, and the bottom plot shows the frequency estimates.

The differences with your example are:

1. I used $21$ coefficients according to $(1)$, multiplied by a Hamming window.
2. I increased the frequency of the input signal by a factor $20$.
3. I correctly aligned the $I$ and $Q$ components, taking into account the delay of the FIR Hilbert transformer.

Here's the code:

idx = (0:511)';
w0 = 2*pi ./ (80 + idx/8);
w0 = 20*w0;
x0 = sin( w0 .* idx);
h = hilbert(x0);
N = 10;      % filter length of FIR Hilbert transformer is 2*N+1
nn = (1:N)';
g = 2*(sin(pi*nn/2)).^2./ (pi*nn);
g = [-g(N:-1:1);0;g];
hw = hamming(2*N+1);
g = g.*hw(:);    % windowed coefficients
I = x0;
Q = filter(g,1,x0);
h2 = I(1:end-N) + 1i*Q(N+1:end);
A1 = abs(h);
A2 = abs(h2);
ph1 = angle(h);
ph2 = angle(h2);
w1 = diff(unwrap(ph1));
w2 = diff(unwrap(ph2));

Note that the instantaneous frequency can also be computed without unwrapping the phase (which can be problematic for noisy signals), as mentioned in this answer.