2

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.

robert bristow-johnson
  • 20,661
  • 4
  • 38
  • 76
Aleksandar Simonović
  • 93
  • 1
  • 7

2 Answers2

2

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/

Hilmar
  • 44,604
  • 1
  • 32
  • 63
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 with the least overall error 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 essentially a windowed Sinc function where the coefficients are all the peaks of the Sinc. The magnitude of the peaks of a Sinc is inversely proportional to sample index $n$. The impulse response of a perfect discrete Hilbert as given b elow also has a magnitude that is inversely proportional to $n$:

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

This is a non-causal impulse response with infinite frequency support, so to be realizable it is delayed and truncated. With the windowing approach to filter design we would multiply this impulse response by a window to minimize frequency domain aliasing and achieve the desired Hilbert filter coefficients. What the approach outlined above provides is an alternate path to doing this using the least squares optimum filter design algorithm. With the Discrete Hilbert for $n$ positive the impulse response is all positive and for n negative it would all be negative. So we get that directly through the modification of the optimally derived Half Band filter coefficients.

Note the similarity and differences with the continuous time Hilbert Transform which has an impulse response given as:

$$h(t) = \frac{1}{\pi t}$$

An intuition as to why every other coefficient is zero in the Discrete Hilbert (and as RBJ further details mathematically in the comments), consider how a continuous Hilbert in the frequency domain is a step function rotated to the $j$ axis, while the discrete Hilbert in the frequency domain is a square wave similarly rotated to the $j$ axis (when we extend the periodic frequency out beyond the first Nyquist zone). A square wave in one domain has every odd sample as zero in the other domain (for example, a square wave in time has non-zero odd harmonics only in frequency - same thing!)

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:

Hilbert Magnitude Lower

Hilbert Magnitude Upper

Dan Boschen
  • 50,942
  • 2
  • 57
  • 135
  • i think the (still acausal) halfband Hilbert transformer can always be thought of as a windowed ideal Hilbert transformer. $$ h[n] = \frac{1 - (-1)^n}{\pi n}w[n] $$ whatever is the width of the window $w[n]$ is twice the delay you have to put in to make the filter realizable. then delay the unhilberted signal by the same amount. $w[0]=1$ and $w[-n]=w[n]$. That will be a halfband filter because every even sample of $h[n]$ is zero. – robert bristow-johnson May 20 '21 at 04:26
  • @robertbristow-johnson My post was about the direct relationship of a regular non-Hilbert Halfband filter to the Hilbert Half-band; I just didn't create the Half-band using a standard "window" design approach even though the result can be similar (Kaiser window and least squares). – Dan Boschen May 20 '21 at 04:46
  • i think i agree with you. i just wanted to point out that this halfband can be looked at as a consequence that every other sample of $h[n]$ is zero. – robert bristow-johnson May 20 '21 at 05:42
  • @robertbristow-johnson I see, yes I follow you now. – Dan Boschen May 20 '21 at 14:54
  • and, Dan, in continuous time, it's not a $\operatorname{sinc}(\cdot)$ function for the Hilbert. the numerator is not $\sin(\pi t)$ but is $\frac{1}{2}\big(1-\cos(\pi t)\big)$. That's what goes over the $\pi t$ in the denominator. – robert bristow-johnson May 20 '21 at 16:36
  • @Robert I believe it works if you use a Sinc in the way I described? The Sinc is what is used to make the half-band, and if you just use the peaks of the Sinc they go down at 1/n as needed for the Hilbert. (With the sign change I described, which is what the cosine does directly I assume, but the point was how to get it directly from the Half Band solution, that other algorithms such as the least squares algorithm I demonstrated provide directly.) – Dan Boschen May 20 '21 at 16:53
  • it doesn't work. the sinc() has even symmetry. Hilberts have odd symmetry. i can't edit my comment anymore, but i was wrong about the $\frac12$ factor. otherwise, i am quite sure that i got it right. – robert bristow-johnson May 20 '21 at 16:58
  • as in general the Hilbert is simply $h[n] = \frac{1}{\pi n}$ and then modified for causal implementation through truncation, time offset and windowing. – Dan Boschen May 20 '21 at 16:59
  • @robertbristow-johnson my answer does have odd symmetry- I think you missed how I made all the samples to the left of center negative. – Dan Boschen May 20 '21 at 17:00
  • 1
    Bandlimited wire: $$ G(f) = \operatorname{rect}(f) $$ then $g(t) = \operatorname{sinc}(t)$ Then the bandlimited Hilbert: $$ H(f) = -j\operatorname{sgn}(f)\operatorname{rect}(f) $$ and then $h(t)=\frac{1-\cos(\pi t)}{\pi t}$ – robert bristow-johnson May 20 '21 at 17:02
  • Note too how the result for the impulse response derived from the non-bandlimited case (simply $h[n] = \frac{1}{\pi n}$) is equal to a decimated by two version of the impulse response from the bandlimited derivation- skipping the zeros basically. – Dan Boschen May 20 '21 at 17:06
  • Sampling the continuous-time Hilbert impulse response $$h(t)=\frac{1}{\pi t}$$ results in aliasing in the frequency domain. If this is a digital filter with a continuous delay parameter that is at least as long as half of the window width, then you need to start out with a simple bandlimited model. – robert bristow-johnson May 20 '21 at 20:53
  • It's well known that to get halfband symmetry, you need to somehow zero every other sample of $h[n]$ except possibly for the center tap which sometimes is used as a "pedestal" to raise or lower the whole thing. If you start out with a bandlimited Hilbert, which you must if you are windowing and sampling, there simply is no $\operatorname{sinc}(\cdot)$ in the math for the Hilbert part. – robert bristow-johnson May 21 '21 at 14:09
  • @robertbristow-johnson The peaks of the Sinc go down at 1/n and I make all the positive ones positive and all the negative ones negative. Wala! (The only reason I refer to a Sinc is because that is equivalent to what would be returned by the least squares algorithm after some windowing). – Dan Boschen May 21 '21 at 14:28
  • But the $\operatorname{sinc}(\cdot)$ isn't in there in the first place. There is no $\sin(\cdot)$ in the numerator to multiply $w[n]$ in the first place. We both agree that there is a $\frac{1}{n}$ factor. – robert bristow-johnson May 21 '21 at 16:40
  • @robertbristow-johnson The relationship I used was $\text{sgn n} \times \text{abs}(0.5 \text{sinc}(\pi n/2)) = 1/(\pi n)$ with n as odd integers excluding n=0 (which is set to 0). – Dan Boschen May 21 '21 at 18:42
  • well, $0$ ain't an odd integer in the first place. so what does that relation have with the Hilbert transform? – robert bristow-johnson May 22 '21 at 21:31
  • @robertbristow-johnson yes correct. It would either be n odd for the decimated solution or every n otherwise with n=0 set to 0. Stay tuned and I will update my answer on the mathematical relationship to go from Sinc to Hilbert and why it works out the way I show. – Dan Boschen May 22 '21 at 23:20
  • 1
    Look, Dan: If you're going to be talking about the sinc function (in the time domain), it's about a corresponding bandlimited counterpart in the frequency domain. But the Hilbert is not one single rect() function but two of them with opposite sign. Now we know that $$ \frac{1-\cos(\pi t)}{\pi t} = \frac{2\sin^2(\frac{\pi}{2}t)}{\pi t} = \sin(\tfrac{\pi}{2}t)\operatorname{sinc}(\tfrac{t}{2}) $$ but that is a consequence of the original analysis (and it's at half the frequency). – robert bristow-johnson May 23 '21 at 02:18
  • 1
    then notice that the $\sin(\tfrac{\pi}{2}t)$ puts in the zeros for every even integer and flips the sign for every odd integer, which is what accomplishes your abs() operation. but the abs() is not the correct math. it just happens to be correct for integer n, but for integer n, the sinc() function for the bandlimited wire is simply the Kronecker delta. – robert bristow-johnson May 23 '21 at 02:21
  • 1
    So here's the cool thing: when you window either (the wire and the Hilbert transformer) to get an FIR, that window spectrum is convolved with the edges of the rect functions, and at high frequencies that convolving does the same thing to both real (wire) and imaginary (Hilbert) parts, so the angle is unaffected. that means that the phase increment (which is instantaneous frequency) is also not affected near the bandedge. – robert bristow-johnson May 23 '21 at 02:34
  • 1
    It also means that multiplying by $e^{j \omega_0 n}$ will correctly offset each frequency component cleanly, even those close to Nyquist that have amplitude affected by that convolving by the window spectrum. So you get a nice frequency shifter that doesn't have aliasing. But down by DC the wire and the Hilbert work differently, so to prevent crap down there, you precede both the wire and the Hilbert transformer with a simple DC-blocking HPF so that there is no content down there by DC. – robert bristow-johnson May 23 '21 at 02:38
  • @robertbristow-johnson very nice, you did it for me. The Hilbert naturally blocks DC but the wire doesn't (which is what you meant by "work differently") and you need the phase to match so I see that you would precede both with the same filter. – Dan Boschen May 23 '21 at 03:26