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:

$${(q_n)}^{k}_{n=0},\quad ‎‎q_n = e^{\frac{\pi}{2^n}}\quad$$

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

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.

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 wrote 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 – Ross Wilkinson 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. – Ross Wilkinson 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$. – Robby Wasabi 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.) – robert bristow-johnson 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. – Robby Wasabi 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. – Ross Wilkinson 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