3
$\begingroup$

As an exercise with IIR allpass filters, I am trying to compensate the phase of a lowpass Linkwitz-Riley filter using the Matlab integrated function iirgrpdelay.

The code I am using is the following:

clear;
close all;
clc;

fs = 48000;
NFFT_half = 4097;

% Design Linkwitz-Riley 24
fcut = 300;
[b,a] = butter(2,fcut/(fs/2),'low');
b = conv(b,b);
a = conv(a,a);
[H,f] = freqz(b,a,NFFT_half,fs);

% Compute iirgrpdelay
N = 40;
gd = grpdelay(b,a,f,fs);
gd_inv = max(gd)-gd;
fn = f./(fs/2);
[num,den,tau] = iirgrpdelay(N,fn,[fn(1) fn(end)],gd_inv);
gd_eq = grpdelay(num,den,f,fs);

% Plot
figure;
semilogx(f,gd); hold on;
semilogx(f,gd_inv);
semilogx(f,gd_eq-tau); grid on;
legend("orig","des","approx");

If I set the cut frequency to fcut = 3000 Hz the approximation with IIR allpass filters is perfect, as you can see here:

enter image description here

But when I try to reduce the fcut to for example 300 Hz, Matlab gives the following error:

Error using iirgrpdelaymex
Poorly conditioned Hessian matrix. Cannot accurately compute the optimization because either the approximation error is extremely small (try reducing the number of
poles or zeros) or the filter specifications yield huge magnitude variations, such as mag=[1 1e9 0 0].

I tried to reduce/increase the order N of the allpass and to play with the additional input parameters of iirgrpdelay, but I can't manage to make it work with low cutoff frequencies, I always have the error reported above.

Does anyone have experience with this function? Are there any alternative techniques to compensate the phase/group delay using IIR allpass filters?

$\endgroup$
1
  • 2
    $\begingroup$ It's not an answer to your question, but I thought I'd check if you've seen Vicanek's "reverse IIR filtering". It's a fun idea, and they demonstrate its use making a linear-phase L-R crossover $\endgroup$
    – cloudfeet
    May 16, 2023 at 18:03

3 Answers 3

2
$\begingroup$

You mentioned that this is an exercise, so I won't comment on the practical value of what you're trying to do. For comments on this, see RBJ's answer. As noted in Dan's answer, if the ratio of cut-off frequency and sampling frequency is very low (as is the case for fcut=300), you're bound to run into problems. The group delay of the lowpass filter increases to a large value close to the cut-off frequency and then falls off quickly towards higher frequencies. In order to compensate for such large variation in group delay in a very small frequency range, you need a high order filter. IIR filters of high order ($20$ or more) are very difficult to design and often even more difficult to implement. The design process is difficult due to the highly non-linear nature of the design problem. Another difficulty is the specification of the desired phase. We have to use some bulk delay to get a desired phase or group delay response that can be approximated by a causal filter with reasonable accuracy, but there's no way of knowing that bulk delay in advance, and it will be different for each chosen filter order. In sum, I don't think that for this problem there is a method that can design an IIR allpass group delay equalizer that performs reasonably well (and that can be implemented).

Of course, the practical solution would be to use downsampling, but since this is meant to be an exercise in phase/delay equalization, I'll show you another approach which doesn't involve downsampling, which performs reasonably well, and which can be implemented. The solution I suggest is to use an FIR filter. The price you have to pay is a relatively large total delay and an increased computational cost. Another disadvantage is that, unlike with IIR filters, there exist no ideal FIR allpass filters. I.e., the magnitude of the frequency response can only be approximately constant, whereas an IIR allpass filter has a perfectly constant magnitude.

Let's see what a design with an FIR filter could look like. I chose a filter length of $201$ taps for the FIR equalizer, and a total desired delay of $200$ samples for the cascade of low pass filter and equalizer. For computing the filter coefficients I used the function lslevin.m, which computes least squares optimal FIR filters with prescribed magnitude and phase responses. The figure below shows the result. Clearly, the magnitude of the equalizer is not ideal but it's probably good enough for most applications. The maximum group delay error is less than $2$ samples, and it decreases with increasing frequency. Of course, the magnitude error as well as the group delay error could be further reduced by increasing the filter length and increasing the total desired delay accordingly.

enter image description here

Below is the Matlab/Octave code to design the FIR equalizer shown above. The coefficient vectors a and b are the ones computed with the OP's code.

nf = 2048;
[H,w] = freqz(b,a,nf);
phi = angle(H);
gd = grpdelay(b,a,nf);

del = 200; phd = - phi - del*w;    % desired phase
h = lslevin(201,w,exp(1i*phd),ones(nf,1));
H = freqz(h,1,nf);
gdh = grpdelay(h,1,nf);
$\endgroup$
3
  • $\begingroup$ Could you please add to the answer the Matlab code you used to get the plots? $\endgroup$
    – skateskate
    May 13, 2023 at 16:22
  • $\begingroup$ @skateskate it's in the link lslevin.m $\endgroup$ May 14, 2023 at 3:00
  • 1
    $\begingroup$ @skateskate: I've added the code for designing the FIR equalizer. $\endgroup$
    – Matt L.
    May 14, 2023 at 10:59
4
$\begingroup$

First of all, there is no good reason to do this phase compensation with the low-pass half of a Linkwitz-Riley filter pair, unless you compensate it somehow with the high-pass counterpart. The LPF and HPF add (not just the magnitude but the complex transfer functions) to a flat-magnitude APF. If you phase compensate just the LPF, they will not add to an APF.

Now, in other contexts, if you have a, say, Butterworth LPF that has a little bump or lip in the group delay or phase delay response that you want to reduce with an APF, the best way to do that, in my opinion, is to manually adjust the resonant frequency and Q of the APF so that its group delay downslope (on the right) is matched to the upslope of the bump or lip. It brings the prominence of that bump down and makes for a virtually maximally flat group (or phase) delay.

$\endgroup$
2
  • $\begingroup$ That's a great tip! Makes a lot of sense, practical. $\endgroup$ May 12, 2023 at 22:49
  • $\begingroup$ If I am not wrong, the phase responses of the complementary LP and HP 24 dB/oct Linkwitz-Riley filters are the same, so if I can compute an "ideal" allpass filter that compensate all the phase for the LP filter, I can apply the same APF to the HP to get an ideal phase response also for the HP part, preserving the flat magnitude response when complex summed. Is it right? $\endgroup$
    – skateskate
    May 13, 2023 at 16:20
2
$\begingroup$

I believe the issue is the having poles and zeros near z=1 which results in very large magnitude variations as the tool is reporting. This comes from the very large ratio of sampling rate to cut frequency. One solution is to reduce the sampling rate with a band equalizer adjusting just the low frequency portion of your spectrum, with the high end of this band matched to a flat delay equalizer for the upper frequencies prior to combining the two. Also reducing the number of poles and zeros will help as the tool suggests, which means reducing the order, not increasing it. Another solution to reduce filter order while still getting the performance of a higher order solution is to implement the compensator as a cascade of lower order filters; second order sections are commonly done. For this approach the target group delays can be distributed amongst the sections such that the overall desired compensation is met when the filters are cascaded.

$\endgroup$
3
  • $\begingroup$ Could you please clarify why reducing the ratio of sampling rate to cut frequency could solve the problem? $\endgroup$
    – skateskate
    May 13, 2023 at 16:32
  • 1
    $\begingroup$ @skateskate I had provided some references in earlier comments regarding that but what Matt provided is more to the point: such ratios drive filter complexity and specifically the order of the filters. To see the issue consider how filter transfer functions are described with polynomials (the order is the order of the polynomial) and with you should see how very small errors can become significant when we raise them to large powers. Reducing the order reduces this issue significantly, either by factoring into smaller sections, reducing the sampling rate, or both. $\endgroup$ May 14, 2023 at 2:50
  • 1
    $\begingroup$ But I concur enthusiastically with Matt's other suggestion to simply use an FIR filter - it's what I tell everyone first getting exposed to filters--- only consider an IIR if you know you can't get it done with an FIR. Although I'm still in this camp, I have also learned a lot here from Robert Bistrow-Johnson, Hilmar (especially related to both of their experience in audio) and MattL on useful merits of IIR filters, starting with this related post of FIR vs IIR: dsp.stackexchange.com/questions/79400/… $\endgroup$ May 14, 2023 at 2:55

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.