bodelin.m

很多matlab的源代码

function [ht,ha] = bodelin(n,d,ty)
% BODELIN Asymptotic and true bode magnitude and phase plots.
%
%	[HT,HA]=BODELIN(N,D,TY) Returns true and asymptotic Bode plots of 
%	H(s)=N(s)/D(s) with coefficients N and D in descending order of s.  
%	TY is a plotting argument: TY=0 (magnitude), TY=1(phase), TY>1(both).
%	DEFAULT: TY=2 (plot both)
%
%	HT = [freq dB phase] in three columns (actual)
%	HA = [freq dB phase] in three columns (asymptotic)
%	USAGE: Use semilogx(HA(:,1),HA(:,2)) to plot asymptotic mag. etc.
%
%       BODELIN (with no input arguments) invokes the following example:
%
%         % Bode and true magnitude plot for H(s) = (s+1)/(s*s+10s+100)
%         >>bodelin([1 1],[1 10 100],0)


% ADSP Toolbox: Version 2.0 
% For use with "Analog and Digital Signal Processing", 2nd Ed.
% Published by PWS Publishing Co.
%
% Ashok Ambardar, EE Dept. MTU, Houghton, MI 49931, USA
% http://www.ee.mtu/faculty/akambard.html
% e-mail: [e-mail]
% Copyright (c) 1998


if nargin==0,help bodelin,disp('Strike a key to see results of the example')
pause,bodelin([1 1],[1 10 100],0);return,end

v=exist('version');if nargin<3,ty=2;end
while d(1)==0,d(1)=[];end,n=n/d(1);d=d/d(1);
rn=flipud(cplxpair(roots(n)));rd=flipud(cplxpair(roots(d)));
sn=1;sd=1;
%if any(any(n<0)),sn=-1;end,if any(any(d<0)),sd=-1;end,s=sn*sd;
s=1;if n(1)<0,s=-1;n=-n;end
zn=0;zd=0;i=find(rn==0);j=find(rd==0);if ~isempty(i);zn=length(i);rn(i)=[];end
if ~isempty(j);zd=length(j);rd(j)=[];end,t=real([rn;rd]);
if any(t>0),error('Cannot plot. TF has some roots in RHP'),return,end
add=zn-zd;xn=abs(rn);xd=abs(rd);ln=length(xn);ld=length(xd);
f=[xn;xd];i=find(f==0);if ~isempty(i);f(i)=[];end,l=length(f);
if l==0,w1=-1;w2=1;f=[];else
f=log10(f);f=sort([f;f+1;f-1;]);lf=3*l;w1=floor(f(1))-1;w2=ceil(f(lf))+1;
if w2<=f(l),w2=w2+1;end,if f(1)<=w1,w1=w1-1;end,end
f=[w1;f;w2];nn=length(f);m=zeros(nn,1);p=m;

%Asymptotic magnitude and phase
for i=1:ln,f0=log10(xn(i));
m=m+(f-f0).*(f>=f0);p1=(f-f0+1).*(f>=f0-1 & f<=f0+1)+2*(f>f0+1);p=p+p1;
if imag(xn(i))~=0,i=i+1;m=m+(f-f0).*(f>=f0);p=p+p1;end,end
for j=1:ld,f0=log10(xd(j));
m=m-(f-f0).*(f>=f0);p1=(f-f0+1).*(f>=f0-1 & f<=f0+1)+2*(f>f0+1);p=p-p1;
if imag(xd(j))~=0,j=j+1;m=m-(f-f0).*(f>=f0);p=p-p1;end,end
p=45*p;m=m+add*f;p=p+90*rem(add+2*(s==-1),4);
while p(1)<=-180;p=p+360;end,while p(1)>180;p=p-360;end %MAKE MAT COMP

if n(ln+1)>0;m=m+log10(abs(n(ln+1)));end %dc gain
if d(ld+1)>0;m=m-log10(abs(d(ld+1)));end,m=20*m;

%Plots
f1=logspace(f(1),f(nn),2);w=logspace(f(1),f(nn),200);
ss=sqrt(-1)*w;h=polyval(n*s,ss)./polyval(d,ss);
mag=20*log10(abs(h));ph=180*unwrap(angle(h))/pi;
ht=[w(:) mag(:) ph(:)];ha=[(10).^(f(:)) m(:) p(:)];
if v==5,f=(10).^f;end
if ty~=1
mn=min(m)-5;mm=max(m)+5;hm=5*ceil(max(mag)/5);hn=5*floor(min(mag)/5);
mmax=max(mm,hm);mmin=min(mn,hn);m1=[mmin mmin];m2=[mmax mmax];
%semilogx(f1,m1,'.',f1,m2,'.'),grid,hold on,
ax=[min(f) max(f) mmin mmax];
if v~=5,axis(ax),end,semilogx(w,mag,'-.'),if v==5,axis(ax),end
grid,hold on,plot(f,m,'-'),hold off
end
if ty>0,if ty>1,disp('Strike a key between plots'),pause,end
mn=min(p)-10;mm=max(p)+10;pm=5*ceil(max(ph)/5);pn=5*floor(min(ph)/5);
mmax=max(mm,pm);mmin=min(mn,pn);m1=[mmin mmin];m2=[mmax mmax];
%semilogx(f1,m1,'.',f1,m2,'.'),grid,hold on,
ax=[min(f) max(f) mmin mmax];
if v~=5,axis(ax),end,semilogx(w,ph,'-.'),if v==5,axis(ax),end,
grid,hold on,plot(f,p,'-'),hold off,end