Magnitude response of allpass filter - group-delay equalization of an IIR filter

In iirgrpdelay help page, the following example was provided for a group delay equalization of an IIR filter:

[be,ae] = ellip(4,1,40,0.2);

f = 0:0.001:0.2;

g = grpdelay(be,ae,f,2); % Equalize only the passband.

a = max(g)-g;

[num,den]=iirgrpdelay(8, f, [0 0.2], a);

Since, iirgrpdelay returns an allpass IIR filter, I was expecting that magnitude of the filter is 0 dB. However, magnitude is not 0 dB for the entire frequency range but something noisy in the passband (attached screen capture). Am I missing something here?

Answer 1

[Paul]

Hi Jay,

Here is the filter

[be,ae] = ellip(4,1,40,0.2);
f = 0:0.001:0.2;
g = grpdelay(be,ae,f,2);   % Equalize only the passband.
a = max(g)-g;
[num,den]=iirgrpdelay(8, f, [0 0.2], a);

Here is its frequency response. Note that the deviations are around 1e-9 dB, so the magnitude is very close to 1

freqz(num,den)

Plote the magnnitude on an absolute scale to better visualize closeness to 1

[h,w] = freqz(num,den);
plot(w,abs(h))

Using the sos form cleans up the "noise", but still not a perfect allpass.

sos = tf2sos(num,den);
[h,w] = freqz(sos);
plot(w,abs(h))

Check the poles and zeros.

zplane(num,den)

The magnitudes of the poles and zeros should be reciprocals of each other. They are close, but not exact.

diff(sort([abs(roots(den)) 1./abs(roots(num))]),1,2)

ans = 8×1

1.0e-10 * -0.0564 -0.0564 0.1980 0.1980 -0.1374 -0.1374 0.0207 0.0207

I'm suprised that iirgrpdelay doesn't enforce that relationship. I checked the code and iirgrpdelay actually develops the filter in sos form, which should be far superior, and then converts to num,den as the last step. That's surprising to me, i.e., why not just stay in sos form? I checked in the debugger, and the low frequency deviation from unity of the sos form inside iirgrpdelay was no more than 1e-14 and usually smaller, so quite a bit better.

Anyway, we can verify with the Symbolic Math Toolbox that the very small deviation at low frequency is not just numerical error in freqz with the sos form.

syms z
numz = simplify(poly2sym(num,z)/z^8);
denz = simplify(poly2sym(den,z)/z^8);
syms h(z)
h(z) = simplify(numz/denz);

figure
plot(w,double(abs(h(exp(1j*w)))))

[Jay]

Hi Paul,

Thanks for the detail navigations and explanations! It clarifies many of questions I had in mind.

One more question, perhaps a little off topic though. Which way would you recommend for the combined magnitude response? In the example, only combined phase response was provided. What if am I supposed to do if I want to check magnitude response for the combined filters (e.g., the response fo both the original filter and the allpass filter)? I found cascade function but it requires filter objects as the arguments..

Thanks again for the time and effort you put.

[Paul]

I'm not sure what you mean by "combined filters" because I didn't see anything in the example in the question (or on iirgrpdelay) that I would call "combined phase." Do you mean the frequency response of be(z)*num(z)/ae(z)/den(z), i.e., the response of the original filter and the allpass filter connected in series?

[Jay]

Yes, I meant the response of the original filter and the allpass filter connected in seriese. You indicated the equation in z-domain and this is only cascade of two filters. I was wondering if there is any fuction so that the calculation is more compact way? What if I have 8 sos filters connected in seriese. What am I supposed to do?

[Paul]

[be,ae] = ellip(4,1,40,0.2);
sosel   = tf2sos(be,ae);

f = 0:0.001:0.2;
g = grpdelay(be,ae,f,2);   % Equalize only the passband.
a = max(g)-g;
[num,den] = iirgrpdelay(8, f, [0 0.2], a);
sosgrp    = tf2sos(num,den);

One way to get the frequency response of the series connection is to compute the individual responses and multiply them

[hel,w] = freqz(sosel);
hgrp    = freqz(sosgrp,w);
h1 = hel.*hgrp;

If you want to get a single SOS implementation, you can just append one to the other

sosboth = vertcat(sosel,sosgrp);
h2 = freqz(sosboth,w);
figure
subplot(211);
plot(w,abs(h1),w,abs(h2)),grid
subplot(212);
plot(w,180/pi*angle(h1),w,180/pi*angle(h2)),grid

However, the sos form constructed from scratch follows a set of rules (see Algorithms) that are likely not obeyed by that simple concatenation. If you want the "true" sos form, one alternative would be to use sos2zp to get the zpk form of each sos filter, concatenate the zeros, concatenate the poles, multiply the gains, and then go back to sos with zp2sos.

Also, make approprate changes if you use the two output form of tf2sos to keep the system gain of each sos filter separate from the dynamic sections.

[Jay]

Thanks Paul! It helps a lot. Thanks for your time and efforts in the conversations.