iki.fi/o

Hilbert transform

Created by Olli Niemitalo on 2003-07-03, last modified 2014-08-28

This is probably the most efficient structure for implementing a Hilbert transform. Actually, it's not a Hilbert transform, but two all-pass IIR filters whose phase difference is approximately 90 degrees over a range of frequencies symmetric around Nyquist/2. Laurent de Soras uses these kind of filters in his HIIR resampling library, which also has functions for coefficient calculation. This note is available in Maple format: hilbert.mws

2003/07

90 degree phase difference IIR allpass pair

The basic building block is:

> H_sect := (z, a) -> (a^2 - z^(-2)) / (1 - a^2 * z^(-2));

[Maple Math]

Equation 1. Allpass section taking one multiplication when a^2 is a known coefficient.

Each such allpass section has poles located on the real axis at a and - a , and zeros at 1/ a and -1/ a .

We construct two allpass filters from such sections:

> H_1 := z -> H_sect(z, 0.6923878)*H_sect(z, 0.9360654322959)*H_sect(z, 0.9882295226860)*H_sect(z, 0.9987488452737)*z^(-1);

[Maple Math]

Equation 2. First allpass filter, delayed by one sample.

> H_2 := z -> H_sect(z, 0.4021921162426)*H_sect(z, 0.8561710882420)*H_sect(z, 0.9722909545651)*H_sect(z, 0.9952884791278);

[Maple Math]

Equation 3. Second allpass filter, has 90 (+- 0.7) degrees relative phase to first filter over a long range of frequencies.

> plot({argument(H_1(exp(I*freq))), argument(H_2(exp(I*freq)))}, freq=0..Pi, phase = -Pi..Pi);

Figure 1. Phases of H_1 (red) and H_2 (green).

> plot(argument(H_2(exp(I*freq))/H_1(exp(I*freq))), freq=0..Pi, phase = -Pi..Pi);

Figure 2. Phase difference of H_2 and H_1.

> plot({argument(H_2(exp(I*freq))/H_1(exp(I*freq))), Pi/2}, freq=0..Pi, phase = Pi/2*0.95..Pi/2*1.05);

Figure 3. Phase difference of H_2 and H_1 (green), detail near 90 degrees (red).

> plot(argument(H_2(exp(I*freq))/H_1(exp(I*freq))), freq=-0.005..0.005, phase = -Pi..Pi);

Figure 4. Phase difference of H_2 and H_1, detail near 0 Hz.

You can use the 90 degree phase difference of the filters as such, or you can construct a complex filter and analyze its properties:

> H := z -> 0.5*(H_2(z)+I*H_1(z));

[Maple Math]

Equation 4. Combined complex filter that will remove negative frequencies.

> plot(20*log10(abs(H(exp(I*freq)))), freq=-Pi..Pi, dB=-60..1, axes=boxed);

Figure 5. Frequency response of filter in Equation 4.

> plot(20*log10(abs(H(exp(I*freq), 0.4))), freq=-0.01..0.01, dB=-60..1);

Figure 6. Frequency response of filter in Equation 4, transition band detail

> plot(20*log10(abs(H(exp(I*freq), 0.4))), freq=0..0.01, dB=-0.0005..0.0005);

Figure 7. Frequency response of filter in Equation 4, transition band - passband detail.

> plot(20*log10(abs(H(exp(I*freq)))), freq=0..Pi, magnitude=-0.0005..0.0005);

Figure 8. Frequency response of filter in Equation 4, passband detail

Transition bandwidth is 0.002 times the width of passband, stopband is attenuated down to -44 dB and passband ripple is 0.0002 dB.

Plenty cheap for a total of 8 multiplications (plus final scaling by 0.5)!

The coefficients were found using a generic evolutionary algorithm. I believe that it would be possible to design coefficients for this filter structure using the software by Artur Krukowski, which finds coefficients for halfband filters: http://www.cmsa.wmin.ac.uk/~artur/Poly.html

2004/01/13

Ben Saucer has succesfully implemented this filter pair to be used in his audio effect. Here are some useful (edited) excerpts from our e-mail correspondence.

Here's a quick diagram of the allpass pair:

        ................ filter 1 .................
   +--> allpass --> allpass --> allpass --> allpass --> delay --> out1
   |
  in
   |    ................ filter 2 .................
   +--> allpass --> allpass --> allpass --> allpass ------------> out2 (+90 deg)

We can use cookbook formulas to convert an allpass section into code. A general IIR recurrence relation:

  out(t) = a0*in(t) + a1*in(t-1) + a2*in(t-2) + ...
         + b1*out(t-1) + b2*out(t-2) + ...

results in the transfer function:

         a0 + a1*z^-1 + a2*z^-2 + ...
  H(z) = ----------------------------
          1 - b1*z^-1 - b2*z^-2 - ...

The allpass section in question has the following transfer function:

          a^2 - z^-2
  H(z) = ------------
         1 - a^2 z^-2

We want to convert this into the recurrence relation. According to the cookbook formulas and the above transfer function:

  a0 = a^2, a2 = -1, b2 = a^2, rest of coefficients zero

  =>  out(t) = a^2*in(t) - in(t-2) + a^2*out(t-2)

which simplifies to the one-multiplication allpass section:

  out(t) = a^2*(in(t) + out(t-2)) - in(t-2)

8 Comments »

  1. High quality stuff! I’ve used this filter for extracting audio envelopes as abs(out1)+abs(out2) and a peak-hold filter.

    It would be interesting to learn more about the theory behind it, and how the coefficients were calculated. Krukowskis site appears to have been taken offline permanently =(

    Comment by David Revelj — 2010-05-05 @ 13:30

  2. David, I don’t know well the theory behind this. Waybackmachine has some PDF’s at Artur Krukowski’s web page. The first of them cites harris, f., M. d’Oreye de Lantremange and A. G. Constantinides, “Digital signal processing with efficient polyphase recursive all-pass filters”, presented at International Conference on Signal Processing, Florence, Italy, September 4-6, 1991.

    I too have used this for envelope extraction, as sqrt(out1^2+out2^2) followed by a smoothing filter. This formulation gives a flat envelope for any audio frequency sinusoid. I recall that there were some problems from the large group delay at low frequencies. Unfortunately I did not put here a phase plot to highlight that.

    For calculation of the coefficients I used Differential Evolution optimization with a minmax cost function. It can be guided a bit by the fact that the a’s of the two filters are always interleaved. It is enough to look at one half of the spectrum, the other half will behave exactly the same. This symmetry is also reflected in the low amount of computations that the filters require, so it should not be broken even when requirements for the other half of the frequency response are more relaxed. It is better to keep the symmetrical structure.

    Comment by Olli Niemitalo — 2010-05-05 @ 14:42

  3. Very nice blog!
    The figure 1 doesn’t appear.

    Comment by Harry Lewis — 2011-08-30 @ 15:56

  4. Thanks Harry,
    figure 1 should be visible now

    Comment by Olli Niemitalo — 2011-09-01 @ 17:26

  5. Great presentation! It would be interesting to compare this to the Hutchins design using six 1st order all-passes. E.g., as implemented in Csound:

    http://www.csounds.com/manual/html/hilbert.html
    http://sourceforge.net/p/csound/csound5-git/ci/08a709e36bd9911e99ec16d1677fa9c62d3c19b6/tree/Opcodes/ugsc.c#l92

    Comment by ljp — 2013-02-17 @ 14:30

  6. hey Olli, first time i stumbled upon your blog. this pair of APFs to compute something like an analytic signal appears very useful (like for frequency shifting). i might put it to use someday.

    i noticed your page on Wikipedia (small world). anyway, i might disagree with you about the Wikipedia editor “Bob K”. i *don’t* think he’s been as useful as he is prolific. i think he has crapped up many, many articles related to signal processing and the sampling theorem. but i’m glad to see you at the English wikipedia.

    bestest,

    r b-j

    Comment by robert bristow-johnson — 2014-07-14 @ 21:39

  7. oh, BTW, were you (yehar) the person that had that old SHArC page on the web? does it still live somewhere?

    Comment by robert bristow-johnson — 2014-07-14 @ 21:41

  8. Hi RBJ! I don’t have a SHARC page, and don’t remember having had one. Maybe a pile of documents in one dir but nothing more. By the way by some extra modulation trickery it is possible to run this filter back and forth to get a linear-phase negative frequency removal filter. All the best and see you around!

    Comment by Olli Niemitalo — 2014-08-28 @ 00:19

RSS feed for comments on this post. TrackBack URL

Leave a comment

Powered by WordPress