splane3d.m

很多matlab的源代码

function splane3d(ss,ww,f,p)
% SPLANE3D Mesh 3-D plot of s-plane magnitude.
%
%	SPLANE3D(SIG,W,F) 3-D mesh plot of s-plane magnitude of STRING F
%	An example of a string variable is F = 's./(s.*s+0.16*s+25)' 
%
%       SPLANE3D(SIG,W,N,D) 3-D plot with N and D as arrays of H(s)=N(s)/D(s) 
%	NOTE: In this form, N and D must be in descending powers of s 
%
%	WITH s = sigma + j*w, SIG = limits on sigma = [SIGlow, SIGhigh]
%	W represents the limits on frequency w with the form [Wlo, Whigh]
%
%    	SPLANE3D (with no arguments) invokes the following example:
%
%	% Plot the 3-D magnitude of H(s) = s/(s*s+0.16s+25) 
%	% over the range -5<sigma<5 and -10<w<10
%	  >>splane3d([-5 5],[-10 10],[1 0],[1 .16 25])


% 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 splane3d,disp('Strike a key to see results of the example')
pause,splane3d([-5 5],[-10 10],[1 0],[1 .16 25]),return,end


if exist('version')==5,ml=4;else,ml=3;end
sc1=50;
if ml==3,sc=40;vr=version;if vr(1)=='S',sc=25;sc1=sc;end,else,sc=100;end,
ds=(ss(2)-ss(1))/sc1;dw=(ww(2)-ww(1))/sc;
ss=ss(1):ds:ss(2);w0=ww(1):dw:ww(2);
vx=matverch;
if vx<5,eval('[s1,w1]=meshdom(ss,w0);');
else, eval('[s1,w1]=meshgrid(ss,w0);');end

j=sqrt(-1);s=s1+j*w1;
if ml==4,w=fliplr(w0);else,w=ww(1):sc*dw/200:ww(2);end
if isstr(f),g=eval(f);else,g=polyval(f,s)./polyval(p,s);end,
g=abs(g);
kk=finite(g);i=find(kk==0);  
m=0;n=0;if ~isempty(i),[m,n]=size(i);g(i)=zeros(m,n);end,
gmx=.25*max(max(g));
g(i)=1e15*ones(m,n);
%end,
flag=0;i=find(ss>=0);
if ~isempty(i),flag=1;l=i(1);end
s=j*w;
if isstr(f),h=eval(f);else,h=polyval(f,s)./polyval(p,s);end,
h=abs(h);hm=max(h);

g1=g;g1=g1.*(g1<=hm)+hm*(g1>hm);
if isstr(f),ff=[' = ',f];else,ff='';end
disp('STRIKE A KEY BETWEEN PLOTS'), 
if ml==3,pause,end

if vx < 4, eval('clg');else,eval('clf');end

box='box';on='on';
if ml==3,
mesh(g,[35 45]),else
eval('mesh(s1,w1,g1),set(gca,box,on),view(35,45)')
if flag==1,
h0=h.*(h<=gmx)+gmx*(h>gmx);
eval('hold on,plot3(zeros(1,length(h)),w,h0),hold off'),
end,
end
title(['magnitude of F(s)' ff])
if flag==0,return,else,pause,end
g=g.*(g<=hm)+hm*(g>hm);
if ml==3,
mesh(g(:,1:l),[35 45]),else
[mm,nn]=size(g(:,l+1:50));g(:,l+1:50)=nan*ones(mm,nn);
ax=axis;
if ax(5)>0,ax(5)=0;end,
eval('mesh(s1,w1,g),axis(ax),set(gca,box,on)'),
eval('view(35,45),hold on,plot3(zeros(1,length(h)),w,h),hold off'),
end
title('LHP magnitude cutaway'),pause
if vx < 4, eval('clg');else,eval('clf');end


plot(w,h),title(['magnitude along jw-axis for F(s)' ff]),
hold off