MATLAB / Octaveで，IIRフィルタの極・零点の位置を算出 -> z平面上にプロット

z平面上にプロットするのであれば， zplane関数で簡単にできてしまいます．
零点-極プロット - MATLAB zplane - MathWorks 日本
zplane関数から，極・零点の位置を抜き出す方法がわからなかった．
そこで，算出しプロットしようと思い立った次第です．

今回は，residue関数のお世話になりました．
部分分数展開と多項式係数間の変換 - MATLAB residue
<追記 20130823> root関数でも位置を求められそうです．
MATLAB 音声信号処理 極・零点配置

<ソースコード>

 

%% set filter coefficient
b=[1.0 0.0 0.0 -0.5];
a=[1.0 0.0 1.0];

%% calculate zero-pole position
[r_b, p_b, k_b]=residue(1,b);
[r_a, p_a, k_a]=residue(1,a);
origin_flag=0;
origin_num=0;
if(length(p_a) > length(p_b))
    origin_flag=-1;
    origin_num=length(p_a)-length(p_b);
    p_b=[p_b', zeros(origin_num,1)']';
elseif(length(p_a) < length(p_b))
    origin_flag=1;
    origin_num=length(p_b)-length(p_a);
    p_a=[p_a', zeros(origin_num,1)']';
end

%% plot
figure;plot( real(p_b), imag(p_b), 'o');
axis([-1.2 1.2 -1.2 1.2]);
axis equal
grid on
hold on
plot( real(p_a), imag(p_a), 'x' );
plot(cos(0:pi/30:2*pi), sin(0:pi/30:2*pi), 'Color', 'black');
legend('zero','pole');
legend('boxon');
xlabel('Real part');
ylabel('Imaginary part');
if(origin_flag ~= 0)
    text(0.05, 0.05, {origin_num});
end
hold off