# IIR Hilbert Transformer

I'm beginning to explore discrete Hilbert transformers - ways to achieve 90°. phase shift across a band of perhaps 6 kHz at a 44.1 kHz sampling rate. I'm trying to stick with IIR filters in order to save computation time - in addition to the fact that I only need to phase shift just a band of frequencies.

I saw this post where the author Ross Wilkinson alludes to creating an IIR based Hilbert transformer by setting the transfer function's zeros to:

The RHS poles inside the unit circle (0 < x < 1) are the inverse of the zeros: \begin{equation} {(p_n)}^{k}_{n=0},\quad p_n = \frac{1}{q_n} \end{equation}

It would be great if somebody could point me towards another resource that would help explain where that comes from in greater detail.

FYI, I'm trying to make a single side band modulator for audio frequencies. It should be able to process in real time on a laptop, with a sampling rate of 44.1 kHz.

• it's a good reference, but you will notice that the IIR All-pass filters (APF) are pretty high order. lotsa poles and zeros. that might not be much savings over an FIR Hilbert Transformer which is a lot easier to understand and design. Feb 4 '17 at 23:14
• i am interested in this, too, and have posted to the music-dsp mailing list. i'll bet someone there has a good reference. Feb 5 '17 at 0:39
• I also saw this crazy implementation: webpages.charter.net/wa1sov/technical/allpass/allpass.html, where the author does some kind of elliptic polynomial. Feb 5 '17 at 0:51
• so, with the help of the music-dsp mailing list, i just rediscovered this page from Olli. it looks pretty good, but i don't know how he gets his coefficients. Feb 5 '17 at 3:01
• Thanks for that link from Olli. Comment #5 refers to a Hilbert transform API in csound - I have no idea what csound is, it seems like a programming language... The API then refers to a paper titled A Hilbert-Transformer Frequency Shifter for Audio, which seemed promising. That paper then refers to a chapter titled Special Filter Designs in Handbook for Digital Signal Processing, which I just purchased for $22. Feb 5 '17 at 21:03 ## 4 Answers This is achievable with two parallel all pass filters. The two all pass filters synthesize an odd ordered low pass filter whose pass band extends from -90º to +90º in the z-domain. (I will discuss this below). $$G_{lowpass}(z) = \frac{A_0(z)+z^{-1}A_1(z)}{2}$$ The low pass filter is then rotated by +90º so that its pass band extends from 0º to 180º, which approximates the Hilbert transform. Rotation mathematically is: $$H_{Hilbert}(z) =G_{lowpass}(-jz)$$ As a consequence, one of the all pass filters becomes completely imaginary - this is the Hilbert transformed signal path. By necessity, since the filter is approximating this Hilbert transform, then the imaginary branch must be applying a 90º phase shift relative to the output of $$A_0(z)$$. $$H_{Hillbert}(z)=G_{lowpass}(-jz)=\frac{A_0(z)+jz^{-1}A_1(z)}{2}$$ Example with Elliptical Filter The following example was taken from Handbook for Digital Signal Processing pg. 920 You can also use this program that I interpreted to design a new half-band elliptic filter with a higher order. Starting with a real, half-band elliptical low pass filter $$G(z)$$ whose frequency response satisfies $$1-\delta_1\leq|G(e^{j\omega})|\leq1\>\>\>\>\>for\>\>0\leq\omega\leq\omega_p$$ $$\>\>\>|G(e^{j\omega})|\leq\delta_2\>\>\>\>\>for\>\>\omega_s\leq\omega\leq\pi$$ To obtain half band we choose: $$\omega_s+\omega_p=\pi$$ $$(1-\delta_1)^{2}+\delta_2^{2}=1$$ But I just used the examples given seventh order poles or an elliptical low pass half-band filter: $$z=0,\>\>\>z=\pm j0.436688,\>\>\>z=\pm j0.743707,\>\>\>z=\pm j0.927758$$ $$G_{lowpass}(z) = \frac{A_0(z)+z^{-1}A_1(z)}{2}$$ Distribute the poles between the two all pass filters like so: $$A_0(z)=\frac{z^{-2}+0.190696}{1+0.190696z^{-2}}\cdot\frac{z^{-2}+0.860735}{1+0.860735z^{-2}}$$ $$A_1(z)=\frac{z^{-2}+0.553100}{1+0.553100^{-2}}$$ Rotate by 90º $$H_{Hilbert}(z)=G_{lowpass}(-jz)=\frac{A_0(-jz)+jz^{-1}A_1(-jz)}{2}$$ $$A_0(-jz)=\frac{0.190696-z^{-2}}{1-0.190696z^{-2}}\cdot\frac{0.860735-z^{-2}}{1-0.860735z^{-2}}$$ $$A_1(z)=\frac{0.553100-z^{-2}}{1-0.553100^{-2}}$$ Phases of $$A_0(z)$$ and $$A_1(z)$$: Magnitude Response of $$H_{Hilbert}(z)$$: Phase Difference Between $$A_0(z)$$ and $$A_1(z)$$ Branches: • I see - the poles come from the projection onto a straight line of a set of points along an ellipse. Interesting find! To replicate that phase difference curve, I had a quick trial using ellip() and ncauer() in Octave. The internal mathematics of the latter is ... painful, but I can confirm that it works Feb 20 '17 at 15:22 • By the way, I think the figure headings "Real Branch" and "Imag Branch" are misleading. As far as I can tell from my tests, all the poles (i.e. both branches) must be rotated by 90 degrees onto the real axis. Feb 21 '17 at 16:59 • Hi Ross, sorry for the late reply. Recall that the Hilbert transformer comes from a transformed half-band elliptic filter:$$G_{lowpass}(z)=\frac{A_0(z) +z^{-1} A_1(z)}{2}$$ The output of $$G_{lowpass}(z)$$ is completely real - its filter coefficients are also real. The output of $$H_{Hilbert}(z)=G_{lowpass}(-jz)$$ however, is complex: $$\frac{A_0(-jz)+jz^{-1}A_1(-jz)}{2}$$ While the filter coefficients of$A_1(-jz)$are real, because it is multiplied by$jz^{-1}$, the output data is imaginary. So, when I refer to Real and Imag branches, I'm referring to$A_0$and$A_1\$. Feb 25 '17 at 23:06

I have insufficient reputation to answer in the comments, so here goes:

I believe Olli calculated his coefficients using some kind of genetic algorithm (I don't know the details).

All I did was plot (from Olli's coefficients) the resulting pole/zero positions in the z-plane, and then take the logarithm to transform into the s-plane.

Olli's poles and zeros almost formed a geometric pattern, so I played around with them in the s-plane until I made a nice pattern, then took the exponential to re-transform back into the z-plane.

Maybe it's already been done elsewhere - I don't know, so I don't have a reference for you. I was just curious to know whether there was an analytical way to achieve similar results to Olli.

By the way it turns out that Olli's 8th order effort beats mine :-), although with mine you can achieve arbitrarily wide bandwidth using ever higher orders, and rotate the filters in the z-plane to make low-pass / high-pass summation pairs (that doesn't work with Olli's).

You can also base the formula on multiples of numbers different from pi, which changes the response characteristics slightly - e.g you can achieve a higher bandwidth with a given filter order, at the expense of greater oscillation in the frequency response curve.

• well, i hope i can boost you by 10 points. (everyone should be able to comment, in my opinion.) Feb 6 '17 at 17:52
• Ross, my theory is that the coefficients are calculated using polynomials such as those used in the Chebyshev and elliptical filters - only instead of shaping magnitude, they are used to shape phase. I'm hoping that this textbook, specifically a chapter in it titled Special Filter Designs by Phillip A. Regalia, will give me some background on how to do that, and where it comes from. Feb 7 '17 at 19:17
• Hi Robby. Looking at the distribution of poles/zeros given by the genetic algorithm (see link in post), it feels like there ought to be some kind of generator function. It's certainly not pi/2^n, or any variant I tried. A related question: Given the generator, what happens to the x-values of the existing poles when you increase from order-8 to order-10? It's possible they will all move to accommodate the additional poles. With pi/2^n you can simply add more terms without affecting the others, so you can build a filter with multiple taps. Feb 9 '17 at 13:26

I did use Differential Evolution to calculate the coefficients. But you can re-design the filter pair easily using the HIIR library by Laurent de Soras (its source code will automatically unzip to a subdirectory hiir). You can use this C++ HilbertDesign.cpp source and compile with g++ using the compile-command quoted on the first line:

// -*- compile-command: "g++ HilbertDesign.cpp -std=c++11 -msse -I. -g3 -O0 hiir/PolyphaseIir2Designer.cpp -o HilbertDesign" -*-
#include <stdio.h>
#include "hiir/PolyphaseIir2Designer.h"

const int numCoefs = 8; // Number of coefficients, must be even
double transition = 2*20.0/44100; // Sampling frequency is 44.1 kHz. Approx. 90 deg phase difference band is from 20 Hz to 22050 Hz - 20 Hz. The transition bandwidth is twice 20 Hz.

double coefs[numCoefs];

int main() {
hiir::PolyphaseIir2Designer::compute_coefs_spec_order_tbw (coefs, numCoefs, transition);
printf("Phase reference path c coefficients:\n");
for (int i = 1; i < numCoefs; i += 2) {
printf("%.20f,", coefs[i]);
}
printf("\n+90 deg path c coefficients:\n");
for (int i = 0; i < numCoefs; i += 2) {
printf("%.20f,", coefs[i]);
}
printf("\n");
return 0;
}


Running HilbertDesign outputs:

Phase reference path c coefficients:
0.47944111608296202665,0.87624358989504858020,0.97660296916871658368,0.99749940412203375040,
+90 deg path c coefficients:
0.16177741706363166219,0.73306690130335572242,0.94536301966806279840,0.99060051416704042460,


You'd use a coefficient $$c$$ in a section with transfer function:

$$H(z) = \frac{c - z^{-2}}{1 - cz^{-2}},\tag{1}$$

implementable by:

$$y[k] = c\,\left(x[k] + y[k-2]\right) - x[k-2],\tag{2}$$

with input $$x$$ and output $$y$$.

Each of the two paths is a cascade of a number of such sections. The phase reference path should be delayed by an additional one-sample delay. The phase difference between the two paths will be approximately 90 degrees over a band centered about sampling frequency / 4. This symmetry gives the somewhat sparse transfer function (Eq. 1) with a simple computational implementation of the filter (Eq. 2). The coefficients are optimal in equiripple sense.

I used HIIR v1.20. The c coefficients are equivalent to my a^2 coefficients. The square roots of the coefficients from HIIR are quite close to my original pole locations, a, here interleaved for the two paths which gives them in ascending order:

a               sqrt(c)          a-sqrt(c)
0.4021921162    0.402215635     -2.35187886997168E-005
0.6923878       0.6924168658    -2.9065827921726E-005
0.8561710882    0.8561932617    -2.21734129581819E-005
0.9360654323    0.9360788374    -1.34051398409385E-005
0.9722909546    0.972297804     -6.84943738937793E-006
0.9882295227    0.9882322446    -2.72186419525422E-006
0.9952884791    0.9952891611    -0.000000682
0.9987488453    0.9987489195    -7.4186057630321E-008


Another possible option of only slightly more computational complexity than running something like 2 sets of 10 sections, might be to simply lowpass the real input, complex heterodyne that result up (or down) in frequency by the lowpass cutoff frequency to IQ, lowpass both the I and Q again, then complex heterodyne back down (or up) to complex baseband. That should chop off one sideband with the same transition shape as the low pass filter. Heterodyne by slightly more than the low-pass bandwidth (add half the transition width) to reject more of the carrier from the SSB result.

This will work with any cascade of biquad units that implement a your favorite IIR low-pass filter, plus two sets of complex vector multiplications to do the heterodyne. A recursion (or two) can generate the complex heterodyne mixer input or you can use a sinewave table in memory.