asymp.m

来自「基于ＭＡＴＬＡＢ的辅助设计 简述了矩阵分析的重要性」· M 代码 · 共 100 行

	
function md = asymp(num,den,mx,wrang)
%
%ASYMP  ASYMP(NUM,DEN,MX,WRANG) is the function which is used to
%       calculate the asymptotic information for a given system
%       expressed by transfer function  [NUM,DEN].  Argument
%       MX is the log magnitude of the system over the frequency
%       range WRANG.
%

%    Author:  Ole Barup Sorensen, Rapid Data Ltd 

%    Copyright (c) 1989-94 by Rapid Data Ltd
%    Revision 10:36  07/02/94
	
[zer,pol,gain] = tf2zp(num,den);
zchk = exist('zer');
[nr,n0] = size(pol);
if zchk ~= 0, [nc,n0] = size(zer); end
[n0,nw] = size(wrang);
i1=1;
k=0;
w0 = -99;
for i = 1:nr
   nx = imag(pol(i));
   if nx > 1d-12
      w0 = abs(pol(i));
      k1 = -2;
   end
   if abs(nx) <= 1d-12
      w0 = abs(pol(i));
      k1 = -1;
   end
   if w0 ~= -99
      if w0 <= wrang(1)
         k = k + k1;
      elseif w0 < wrang(nw)
         i1 = i1 + 1;
         mata(i1,1) = w0;
        mata(i1,2) = k1;
      end
   end
   w0 = -99;
end
if zchk ~= 0
   for i=1:nc
      ny = imag(zer(i));
      if ny > 1d-12
         w0 = abs(zer(i));
        k1 = 2;
      end
      if abs(ny) <= 1d-12
         w0 = abs(zer(i));
         k1 = 1;
      end
      if w0 ~= -99
         if w0 <= wrang(1)
            k1 = k1 + 1;
         elseif w0 <= wrang(nw)
            i1 = i1 + 1;
            mata(i1,1) = w0;
            mata(i1,2) = k1;
         end
      end
      w0 = -99;
     end
   end
end
mata(1,1) = wrang(1);
mata(1,2) = k;
i1 =i1 + 1;
mata(i1,1) = wrang(nw);
mata(i1,2) = 0;
for i = 1 : i1
   for j = i+1 : i1
      if mata(i,1) > mata(j,1);
         mtemp = mata(i,1);
         mata(i,1) = mata(j,1);
         mata(j,1) = mtemp;
         mtemp = mata(i,2);
         mata(i,2) = mata(j,2);
         mata(j,2) = mtemp;
      end
   end
end
mata(1,3) = 2 * mx(1) - mx(2);
mata(1,2) = 20 * k;
for i = 2 : i1-1
   mata(i,2) = mata(i-1,2) + 20*mata(i,2);
   mata(i,3) = mata(i-1,3)+mata(i-1,2)*log10(mata(i,1)/mata(i-1,1));
end
k0 = 1;
for j = 1 : nw
   if wrang(j) > mata(k0+1,1)
      k0 = k0 + 1;
   end
   md(j) = mata(k0,3) + mata(k0,2) * log10(wrang(j)/mata(k0,1));
end