3
$\begingroup$

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.

warping vs. no warping

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  
$\endgroup$
  • $\begingroup$ What is your problem specifically? You haven't made clear what you've tried and where it isn't working for you. $\endgroup$ – Jason R Apr 17 '13 at 12:45
  • $\begingroup$ I am having difficulty getting the 90 degree phase difference close to dc. The paper claims that the frequency range can be extended using the frequency warping described in the program. However I fail to get the phase difference below about 1kHz. $\endgroup$ – melinnde Apr 17 '13 at 12:56
  • $\begingroup$ I just tried to run your code. One thing that seems wrong is that your halfband design routine an order of N=5 is calculated whereas the paper use order 13. Do you know why? $\endgroup$ – niaren May 22 '13 at 8:21
3
$\begingroup$

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:

Figure 1 - Hilbert generator all-pass filter

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:

enter image description here

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:

enter image description here

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:

enter image description here

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.

$\endgroup$
  • $\begingroup$ i hadn't seen this post before. it's pretty good. $\endgroup$ – robert bristow-johnson Feb 4 '17 at 23:16
1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Thank you for the reply, I think I follow what you are saying and may try to attempt it. The method used above however, attempts to approach dc by sacrificing higher frequency content, pre-warping the half-band specification and shifting the poles of the filter to meet the specification. I've noticed one mistake in the code so far which I will rectify but what's troubling me most I can't get the results shown in the paper regardless of the fix $\endgroup$ – melinnde Apr 17 '13 at 14:27
  • $\begingroup$ You mention windowing the hilbert filter. You mean, simply construct the hilbert FIR filter, (lets say it was an FIR) and then literally apply a window to it? What benefit does that give over just no windowing, and convolving your input with the hilbert FIR? Thanks. $\endgroup$ – TheGrapeBeyond Aug 24 '13 at 3:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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