Magnitude response of allpass filter - group-delay equalization of an IIR filter
3 ビュー (過去 30 日間)
古いコメントを表示
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 件のコメント
採用された回答
Paul
2023 年 11 月 16 日
編集済み: Paul
2023 年 11 月 16 日
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 件のコメント
Paul
2023 年 11 月 17 日
[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.
その他の回答 (0 件)
参考
カテゴリ
Help Center および File Exchange で Filter Design についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!