0
$\begingroup$

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.

$\endgroup$
3
  • 1
    $\begingroup$ Being new to signal processing, I recommend starting with easier problems. Also, SciPy already includes a Hilbert transformer; why do you want to implement your own? $\endgroup$ – MBaz May 18 at 17:39
  • $\begingroup$ That's one of tasks I need to do, I understand theory, but I'm not sure how to implement. I would appreciate help. Thanks $\endgroup$ – Aleksandar Simonović May 18 at 17:41
  • 1
    $\begingroup$ Hilbert transformers are not, strictly speaking, all-pass filters (APF) because they don't pass DC. $\endgroup$ – robert bristow-johnson May 19 at 3:02
2
$\begingroup$

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/

$\endgroup$
4
  • $\begingroup$ These are the requirements that I have: Oscillations in bandwidth must be less than 0.05dB, Attenuation less than 3db in frequency range [0.01pi, 0.99pi] $\endgroup$ – Aleksandar Simonović May 18 at 19:55
  • $\begingroup$ This is what i wrote: N = 50 hh = signal.firwin(2*N+1, [0.01,0.99], window='hamming', pass_zero = False) $\endgroup$ – Aleksandar Simonović May 18 at 19:56
  • $\begingroup$ That's a very nice link (the 2nd). +1. @AleksandarSimonović That link should make you understand everything. $\endgroup$ – a concerned citizen May 18 at 22:03
  • $\begingroup$ I was happy to see the reference to Andrew Duncan's analytic impulse stuff. I was both impressed and jealous when that came out. $\endgroup$ – robert bristow-johnson May 20 at 21:25
0
$\begingroup$

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:

Hilbert Magnitude Lower

Hilbert Magnitude Upper

$\endgroup$
24
  • 1
    $\begingroup$ 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}$ $\endgroup$ – robert bristow-johnson May 20 at 17:02
  • 1
    $\begingroup$ 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). $\endgroup$ – robert bristow-johnson May 23 at 2:18
  • 1
    $\begingroup$ 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. $\endgroup$ – robert bristow-johnson May 23 at 2:21
  • 1
    $\begingroup$ 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. $\endgroup$ – robert bristow-johnson May 23 at 2:34
  • 1
    $\begingroup$ 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. $\endgroup$ – robert bristow-johnson May 23 at 2:38

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.