Magnitude response of allpass filter - group-delay equalization of an IIR filter
3 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
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?
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1540662/image.png)
0 comentarios
Respuesta aceptada
Paul
el 16 de Nov. de 2023
Editada: Paul
el 16 de Nov. de 2023
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)
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)))))
5 comentarios
Paul
el 17 de Nov. de 2023
[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.
Más respuestas (0)
Ver también
Categorías
Más información sobre Filter Design en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!