Implementing FIR Hilbert Transformer in Python

Question

I want to implement FIR Hilbert Transformer in Python using window method. I'm new to digital signal processing so please point out to me all mistakes you notice.

Ideal Hilbert transformer is described as: $H(e^{j\omega}) = \begin{cases} -j & \text{, $0 < \omega < \pi$} \\ j & \text{, $-\pi < \omega < 0$} \\ \end{cases} $

from which we conclude that ideal Hilbert Transformer is all pass filter $|H(e^{j\omega})| = 1$.

I know that FIR filter can be implemented using scipy.signal.firwin function.

What I don't know is which parameters to pass to function in order to get Hilbert filter aproximation.

Answer 1

You can't implement an ideal Hilbert transformer: its impulse response is non causal and infinite in time.

So you can only implement an "approximate" one and the best way to do this depends on the specific requirement of your application: what's your frequency range of interest, how much magnitude and/or phase deviation can you tolerate, are you constrained in latency or computational resources, do you need to maintain causality and transients, etc.

Python has a Hilbert transform function https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.hilbert.html but without understanding your requirements first, it will be difficult to control the parameters properly.

Very good article on Hilbert transform and analytical signals in the discrete domain: http://andrewduncan.net/air/

Answer 2

There are routines that will provide the Hilbert coefficients directly, but an approach I like to use given its simplicity and clarity in functionality is to transform a Half Band filter to a Hilbert as follows:

Step 1: Estimate the number of taps needed from the specifications using these commonly used estimators. The one Marcus Mueller provided would be well suited for the OP's case given there is a passband ripple and stopband requirement. Using that along with 45 dB rejection, the estimated number of coefficients (taps) needed came to 269 (also for the estimator, the passband edge of 0.99 is changed to 0.98 due to the Half Band to Hilbert transformation). Nearly half of these 269 coefficients will be zero, and the coefficients are symmetric such that it can be implemented with approximately 69 multiplications.)

Step 2: If an absolute min/max requirement is needed, use the Remez algorithm to design an optimum half band filter in minimax sense. (Typically I would prefer to use the least squares algorithm for an optimum solution in the least squared sense which would result in fewer taps but have higher peak error. A windowing design approach can be done as well in which case I would choose the Kaiser window and window the Sinc function that defines this HB filter; the results will be quite similar to least squares). In Python the Remez solution is done as follows:

coeff = sig.remez(N, [0, pb, 0.5- pb, .5], [1, 1, 0, 0])

Alternatively for least squares use:

coeff = sig.firls(N, [0, 2 * pb, 2* (0.5- pb), 1], [1, 1, 0, 0])

Where pb is the pass band edge normalized to $F_s = 1$, where $F_s$ is the sampling rate and $N$ is the total number of coefficients which must be an odd number.

Step 3: Transform the $N$ coefficients numbered $1:N$ for a Half Band filter (which should be odd) to the coefficients for a Hilbert Transform by setting the center coefficient $\lceil N/2 \rceil$ of the Half Band to 0, all coefficients that appear earlier than the center coefficient (coefficients $1: \lfloor N/2 \rfloor$) to be negative the absolute value of the half band, and all the coefficients that appear later than the center coefficient (coefficients $\lceil N/2 \rceil+1: N$) to be positive. Multiply all coefficients by two to normalize the scale.

The reason this works is the Half Band filter is (a windowed version of) a Sinc function where the coefficients are all the peaks of the Sinc. The peaks of the Sinc go down at rate $1/n$. The impulse response of a Hibert is a non-causal $1/n$ (with $\pi$ scaling but for clarity of the main point I leave that out). For for n positive it would all be positive and for n negative it would all be negative. So we get that directly through the modification of the Half Band filter coefficients as optimally derived (and as RBJ further details mathematically in the comments). We make it causal by truncating the Sinc and shifting in out to $n= \lceil N/2 \rceil$. For this reason we need to feed the input into a tracking filter which is just a simple delay by $n= \lceil N/2 \rceil$ samples, which will be in quadrature (as desired) from the output of the filter above.

Doing the above came up with the following result for amplitude balance with a perfect quadrature balance across the band: