inagersh.m

控制系统计算机辅助设计——MATLAB语言与应用(源代码）

function inagersh(H,nij)
t=[0:.1:2*pi,2*pi]'; [nr,nc]=size(H); nw=nr/nc; ii0=1:nc; 
if nargin==1, ii=1:nc; jj=1:nc; 
else, ii=nij(1); jj=nij(2); end
for i=1:nc, circles{i}=[]; end
for k=1:nw  % 对各个频率获取逆 Nyquist 阵列
   Ginv=inv(H((k-1)*nc+1:k*nc,:)); nyq(:,:,k)=Ginv;
   for j=1:nc
      ij=find(ii0~=j);
      v=min([sum(abs(Ginv(ij,j))),sum(abs(Ginv(j,ij)))]);
      x0=real(Ginv(j,j)); y0=imag(Ginv(j,j));
      r=sum(abs(v)); % 计算 Gershgorin 圆的半径
      circles{j}=[circles{j}, x0+r*cos(t)+sqrt(-1)*(y0+r*sin(t))];
end,end
for i=ii, for j=jj   % 绘图
    if nargin==1, subplot(nc,nc,(i-1)*nc+j); end
    for k=1:nw, NN(k)=nyq(i,j,k); end
    if i==j,  % 对角图，带有 Gershgorin 带
       plot(real(NN),imag(NN),real(circles{i}),imag(circles{i}));
   else      % 非对角元素
       plot(real(NN),imag(NN))
   end, 
end,end