9
$\begingroup$

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:

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

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.

$\endgroup$
8

5 Answers 5

7
$\begingroup$

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: Pole distribution diagram $$ 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)$: All pass filter phase response

Magnitude Response of $H_{Hilbert}(z)$: Filter response

Phase Difference Between $A_0(z)$ and $A_1(z)$ Branches: Phase difference

$\endgroup$
3
  • $\begingroup$ 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 $\endgroup$ Commented Feb 20, 2017 at 15:22
  • $\begingroup$ 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. $\endgroup$ Commented Feb 21, 2017 at 16:59
  • $\begingroup$ 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$. $\endgroup$ Commented Feb 25, 2017 at 23:06
6
$\begingroup$

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.

$\endgroup$
3
  • $\begingroup$ well, i hope i can boost you by 10 points. (everyone should be able to comment, in my opinion.) $\endgroup$ Commented Feb 6, 2017 at 17:52
  • $\begingroup$ 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. $\endgroup$ Commented Feb 7, 2017 at 19:17
  • $\begingroup$ 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. $\endgroup$ Commented Feb 9, 2017 at 13:26
3
$\begingroup$

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

There are two ways to obtain a true Hilbert transformer by forward-backward (ping-pong) IIR filtering, here expressed using the frequency responses $H_{\text{ref}}(\omega)$ and $H_{\text{ref}+90^\circ}(\omega)$ of the all-pass branches, with the subscript denoting the approximate phase shift. The time reversal of a conjugate-symmetric frequency response $H(\omega)$ results in a frequency response $H(-\omega)$, which is the signature of a backwards pass on the left sides of:

$$H_{\text{ref}+90^\circ}(\omega)H_{\text{ref}}(-\omega) = H_{\text{ref}+90^\circ}(\omega)H_{-\text{ref}}(\omega) = H_{90^\circ}(\omega),$$

and

$$-H_{\text{ref}}(\omega)H_{\text{ref}+90^\circ}(-\omega) = -H_{\text{ref}}(\omega)H_{-\text{ref}-90^\circ}(\omega) = -H_{-90^\circ}(\omega) = H_{90^\circ}(\omega).$$

Furthermore, Martin Vicanek shows in a 2015—2022 article A New Reverse IIR Filtering Algorithm how to obtain a causal, truncated approximation of each the forward IIR filter impulse response and a causal delayed version of the backward IIR filter impulse response, using low-complexity FIR filters, to obtain an approximate delayed Hilbert transformer.

$\endgroup$
1
$\begingroup$

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.

$\endgroup$
0
$\begingroup$

I don't have the "reputation" yet to add reply comments unfortunately so I have to post here. First thanks to all the amazing people who have posted techniques and insights here. I've recently built three methods mentioned on this page. So I figured I should say thanks by sharing some observations, in case they are useful for others.

  1. Allpass pairs. I did 8, 12, and 16 stages, using coefficients generated by Laurent de Soras' excellent HIIR library as described in @Olli Niemitalo's post. Adding more stages improves behaviour nearer DC and Nyquist, at the price of more "latency". Latency is in quotes here because the response is not linear phase -- delay is frequency dependent. That also means that neither path will reconstruct the shape of the original waveform. But it's cheap and overall the delay is on the order of a few samples. Implementation was pretty straightforward following Olli's post.

  2. Reverse-IIR allpass technique from Martin Vicanek's paper referenced in Olli's post (https://www.vicanek.de/articles/ReverseIIR.pdf). The main point here is to get a linear phase output (closer to a true Hilbert), where one path is an unfiltered delayed copy of the input, and the other approximates two allpass chains, with one "reversed". The reversal converts 8 allpasses in series to 8 truncated-IIR sections in parallel. The total latency was 4095 samples (not 4096 as stated in the paper): the sum of the delays (2048 + 1024 + 512 + etc.) in the reversed IIR pass.

I think I got it right. The waveforms look good, the phase plot looks really great down to about 30Hz, but I'm definitely seeing more audible leakage between the two paths once I try to frequency shift them. I presume this is a side-effect of the truncated IIR.

Another thing maybe not clear in the paper -- I suspected that the one component that isn't filtered (the one scaled by C) should also be a 4095-sample delayed copy of the input to align it with the rest of the reversed allpass sections. I'm not 100% sure, but it does improve the SNR for me.

The latency of 4095 samples is a lot though. It being so large leads to a lot more "post-ringing" images with dynamic input. This method also has a higher CPU load and much bigger memory footprint for the delay sections in the reversed IIR (36855 floats total).

  1. The heterodyne techinque described by @hotpaw.

I used a windowed-sinc FIR of odd-numbered length with cutoff at samplerate/4 (half-band) for the filtering, so that I could get linear phase output. It worked, and reconstructed the shape of the original wave pretty well, though as with other methods the response starts to fall apart as you get to lower frequencies.

Using a longer FIR buys better performance at low frequencies at the cost of the latency. (The latency ends up being the length of the FIR minus one sample.)

E.g. a 193 samples FIR could reproduce 420Hz with an excellent 25% phase separation (at SR=48kHz) but starts to lose phase separation down below that. With a 577 sample FIR the threshold was down to 120Hz. Not as good as the reverse-IIR, but it's also 15% the latency. Increasing latency also increases CPU load for the FIR convolution, so at a certain point it won't make sense (and the low freq response has diminishing returns).

Speaking of performance, since I set the frequency at SR/4 there's a few optimizations: the FIR is zero at odd samples, except for the centre sample. So the filtering convolution code can be about half as expensive. It also means that the complex oscillator used for heterodyning outputs a sequence like [1 0], [0 1], [-1 0], [0 -1] so you can do it with no trigonometry calls and half the computations can be cut.

At first I also noticed that you lose half the spectrum with this method, but I found a way to resolve it -- I'm not sure if this was known or implied, or a new finding:

In the first pass the (real) input is filtered before heterodyning. This is the stage that effectively kills half your spectrum. I presumed that the imaginary component was zero. But what I found is that if instead you subtract the filtered part from (a delay compensated copy of) the input to get a HPF, you can use that HPF'd as your imaginary signal for the heterodyne. At the output this preserves the full spectrum of your input with a proper 25% separation. The latency compensation for that part is just half the [length of the FIR minus 1].

The output of the 0% phase rotated path is close to the original waveform, mostly linear phase, and the SNR looks better than the reverse-IIR method. It's also a lot cheaper and lower latency.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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