Hilbert transform filter for audio applications: Using IIR half-band parallel all pass structure

Question

Does anyone have any experience designing wideband IIR Hilbert transform filters with audio applications in mind? I am using the filter for single side band modulation audio effects discussed in this paper :-http://www.mikrocontroller.net/attachment/33905/Audio_Hilbert_WAR19.pdf.

I tried to implement the design technique discussed in this paper:- http://www.aes.org/e-lib/browse.cfm?elib=15680 but have not been able to achieve the results in the paper using frequency warping.

Any help on the topic is appreciated; please find MATLAB code below.

Update: I now achieve better results without warping with some loss in accuracy in the magnitude response.

Frequency Warping:

Half band frequencies:-

$\omega _{p}= \tan^{-1}\left [ \left ( \frac{\tan \frac{\omega _{1}}{2}}{\tan \frac{\omega _{2}}{2}} \right )^{\frac{1}{2}} \right ]$

$\omega _{s}= \tan^{-1}\left [ \left ( \frac{\tan \frac{\omega _{2}}{2}}{\tan \frac{\omega _{1}}{2}} \right )^{\frac{1}{2}} \right ]$

$ \text{Where } \omega _{1} \text{ and }\omega _{2} \text{ is the hilbert frequency range}$

Warping:-

$z^{-1}= \frac{k_{b}+z^{-1}}{1+k_{b}z^{-1}}$

$k_{b} = \frac{\beta -1}{\beta +1}$

$\beta = \left [ \tan \frac{\omega _{1}}{2}\tan \frac{\omega _{2}}{2} \right ]^{_{\frac{1}{2}}}$

$c_{new} = \frac{-p+k_{b}}{1-pk_{b}} $

$ \text{Where } c_{new} \text{ is the new filter coefficient and } p \text{ is the filter pole} $

Code

%% IIR Hilbert transformer pair
%specifacation
fs = 40000;
fsinv = 1/fs;
w1 = 200*fsinv*pi*2;
w2 = fsinv*15000*pi*2;

% prewarp spec
wp_hil = 2*atan(sqrt((tan(w1/2))/(tan(w2/2))))
ws_hil = 2*atan(sqrt((tan(w2/2))/(tan(w1/2))))

% halfband spec
%wp = (pi/2)-wp_hil; 
A_s = 80;

% Calculate poles for A0(z) and A1(z)
[alpha_dash] = IIR_Halfband(wp_hil ,ws_hil ,A_s); 

% Calculate poles for A0(z^2) and A1(z^2)
alpha = sort([1j*sqrt(alpha_dash) -1j*sqrt(alpha_dash)]);

% Convert to hilbert transform filter & interlace poles
var = 0;
idx_a = 1;
idx_b = 1;
A1_p = zeros((length(alpha))/2,1);
A0_p = zeros((length(alpha))/2,1);
for idx = length(alpha):-1:1
    if (var <2)
        A1_p(idx_b) = 1j*alpha(idx);
        var = var + 1;
        idx_b = idx_b +1;
    else
        A0_p(idx_a) = 1j*alpha(idx);
        var = var + 1;
        idx_a = idx_a +1;
        var = mod(var,4);
    end
end
A1_p = [0;A1_p]; % additional delay element for imajinary path

% Frequency warping
A1_p2 = zeros((length(A1_p)),1);
A0_p2 = zeros((length(A0_p)),1);

beta = sqrt((tan(w1/2))*(tan(w2/2)));
kb = (beta-1)/(beta+1);
for idx = 1:length(A1_p2)
    A1_p2(idx) = (-A1_p(idx)+ kb)/(-A1_p(idx)*kb +1);
end
for idx = 1:length(A0_p2)
    A0_p2(idx) = (-A0_p(idx)+ kb)/(-A0_p(idx)*kb +1);
end

% generate transfer functions
sys_i = tf([A1_p2(1) 1],[1 A1_p2(1)],fs,'Variable','z^-1');
sys_r =1;
for idx = 1:length(A0_p)
    sys_i = sys_i*tf([A1_p2(idx+1) 1],[1 A1_p2(idx+1)],fs,'Variable','z^-1');
    sys_r = sys_r*tf([A0_p2(idx) 1],[1 A0_p2(idx)],fs,'Variable','z^-1');
end

% plot filter atributes
%fvtool(cell2mat(sys_r.num),cell2mat(sys_r.den))
%fvtool(cell2mat(sys_i.num),cell2mat(sys_i.den))

% phase difference
[phi,w] = phasez(cell2mat(sys_r.num),cell2mat(sys_r.den)) ;
[phi2,w] = phasez(cell2mat(sys_i.num),cell2mat(sys_i.den)) ;
plot(w/(2*pi),180*((phi-phi2)/pi))

Half-band filter design function based on Program 13_9 in Digital Signal Processing: A Computer-Based Approach 3ed http://highered.mcgraw-hill.com/sites/0072865466/

% IIR Lowpass Half-band Filter Design
function[alpha] = IIR_Halfband(wp,ws,A_s) 
% Enter filter stopband specifications

%ws = pi - wp;
delta_s = 10^(-(A_s/20));
% Order estimation
r = tan(0.5*wp)/tan(0.5*ws);
r_hat = sqrt(1-(r^2));
q0 = 0.5*(1-sqrt(r_hat))/(1+sqrt(r_hat));
q = q0+2*(q0^5)+15*(q0^9)+150*(q0^13);
D = ((1-(delta_s^2))/(delta_s^2))^2;
N = ceil((log10(16*D))/(log10(1/q)))
if (mod(N,2) == 0)
    N = N+1
end;
% Readjusting ripple size
D = (1/16)*(10^(N*log10(1/q)));
delta_s = sqrt(1/(1+sqrt(D)));
% Computing filter coefficients
alpha = zeros(length((N-1)/2),1);
for k = 1:(N-1)/2,
    sum1 = 0;sum2 = 0;
    for i = 0:6,
        sum1 = sum1+((-1)^i)*(q^(i*(i+1)))*sin((2*i+1)*k*pi/N);
        if (i ~= 0)
        sum2 = sum2+((-1)^i)*(q^(i^2))*cos(2*pi*k*i/N);
        end;
    end;
    lambda(k) = (2*q^(0.25)*sum1)/(1+2*sum2);
    b(k) = sqrt((1-r*(lambda(k)^2))*(1-(lambda(k)^2)/r));
    c(k) = (2*b(k))/(1+(lambda(k)^2));
    alpha(k) = (2-c(k))/(2+c(k));
end;
alpha = sort(alpha);
end  

Answer 1

I realize this question is 2 years old now, and it looks like OP eventually achieved a pretty decent phase response - but at the cost of magnitude response accuracy (according to the update in the post).

However, since OP asked for help on the topic (not just the specific implementation), I believe I can add something useful.

It turns out there is actually quite an easy way to implement a Hilbert Transformer with an IIR parallel all-pass structure: A method based on a particular geometric sequence of poles and zeros in the complex plane.

Here it is:

The entire structure in Figure 1 is all-pass, pure real (jy = 0), with the sequence of x-coordinates for the zeros on the RHS of the z-plane (x>0) being generated with the formula:

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

For odd k.

(The example shows k=19; zeros at +/- 23.1407 and +/- 4.8105 have been omitted for clarity).

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}{z_n} \end{equation}

On the (x<0) LHS of the z-plane, the poles and zeros are simply the mirror image of those on the RHS.

A pair of all-pass filter sections can be obtained from the above layout by separating alternate zeros (by odd and even k), along with their corresponding p=1/z poles, to give the following structures:

Effectively, we have split the filter generator function into 2 parts:

\begin{equation} Filter A:\quad {(z_n)}^{(k-1)/2}_{n=0}, ‎‎z_n = e^{\frac{\pi}{2^{(2n)}}}\qquad Filter B:\quad {(z_n)}^{(k-1)/2}_{n=0},\quad ‎‎z_n = e^{\frac{\pi}{2^{(2n+1)}}} \end{equation}

The phase difference between the filter pair A and B is:

We can straighten those lines to a constant phase difference by adding a unit delay (${z^{-1}}$: itself an allpass filter with a pole at 0 and zero at ∞) in series with Filter A:

As you can see this gives you a ${\frac{\pi}{2}}$ = 90° phase difference between DC and the Nyquist frequency (at ${\pi}$ Hz), with sharp slopes at those 2 theoretical discontinuities.

Furthermore, the amplitude response is perfectly flat since the structure is all-pass.

The higher the k in the filter generator function, the closer to the unit circle the resulting pole/zero pairs, and the tighter the phase response near DC (with symmetrical behaviour near Nyquist) - at the expense of filter settling time.

As a bonus, you can rotate the positive and negative 'arms' of the combined filter in figure 1 into a conjugate '<' shaped pair mirrored in the x-axis. Dividing the filter in the same fashion previously explained, and inserting a unit delay in series on one branch, will give you a pair of outputs exhibiting a ${\pi}$ = 180° phase difference with a transition frequency which is dependent on the angle of the '<'.

By adding and subtracting the outputs of such a filter you get a low pass and high pass response, all from the same basic design.

Answer 2

This is a hard problem for a variety of reasons:

Firstly the ideal Hilbert Transformer is infinitely non causal, which is a fancy way of saying that the impulse response starts at $t = -\infty $. That's of course impossible to implement so some sort of an approximation is required. You need to find out over what frequency range and to what precision you need the approximation. For example for audio you could say "40Hz-16 kHz, phase = 90deg +- 1 deg, level = +- 0.1 dB".

Your IIR filter is causal by definition. In order to have any chance on matching the Hilbert transformer needs to be cascaded with a pre-delay, i.e. your target will need to be ideal Hilbert transformer cascaded with some bulk delay. The bulk delay will have to be applied to all parallel paths in your system as well and it will add to the overall latency. If you need a high precision Hilbert Transformer with very low latency, you are out of luck: Can't be done.

A quick & dirty way to determine the required pre-delay is to take a VERY long FIR version of the ideal HT and start truncating or windowing and determine at what window length you can meet the actual requirement. Then retry your IIR fit with a target of Hilbert + bulk delay.

FWIW: the reason what you see trouble at lower frequencies means that your pre-delay (if you have) is too short. The ideal Hilbert transformer requires an infinitely sharp transition of the phase from 0deg at DC to 90 deg at anthing above DC. That's off course impossible so there must be some sort of finite transition at low frequencies.

Answer 3

I created a version in python/signal to create filter:

def hilbert(order):
    order = math.floor(order/2)*2
    # See instructions for generating poles/zeros at https://dsp.stackexchange.com/questions/8692/hilbert-transform-filter-for-audio-applications-using-iir-half-band-parallel-al
    Ik = np.e ** (np.pi / (2 ** np.arange(1,order+1)))
    Ik = Ik * ((-1) ** np.arange(1,order+1))
    Isos = signal.zpk2sos(Ik, 1.0 / Ik, 1.0)
    h1 = signal.sosfreqz(Isos)
    k = 1/np.max(np.abs(h1))
    Isos = signal.zpk2sos(Ik, 1.0 / Ik, k)
    Qsos = signal.zpk2sos(-Ik, -1.0 / Ik, k)