zplane3d.m

很多matlab的源代码

function zplane3d(r,f,p)
% ZPLANE3D Mesh 3D-plot of z-plane magnitude
%
%	ZPLANE3D(R,F) Surface 3D-plot of z-plane magnitude of string function F.
%	An example of a string function would be F = ['z./(z-0.9)./(z+1.1)'];
%
%       ZPLANE3D(R,N,D) Uses the num and den of H(z) in descending order
%       R represents the limits [RLow, RHigh] on R where z = R*exp(j*theta).
%
%	ZPLANE3D (with no input arguments) invokes the following example:
%
%	%Plot the 3-D z-plane magnitude of H(z) = z/(z*z+0.2z-0.99)
%	>>zplane3d(2,[1 0],[1 0.2 -0.99])


% 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 zplane3d,disp('Strike a key to see results of the example')
pause,zplane3d(2,[1 0],[1 0.2 -0.99]),return,end

if exist('version')==5,ml=4;dd=50;else,ml=3;dd=20;
vr=version;if vr(1)=='S',dd=15;end,end

vx=matverch;

r=abs(r);ss=-r:r/dd:r;
if ml==4,ll=length(ss);ss(ll)=[];ss(1)=[];end
ww=ss;vw=[135 45];%50,45
%[s1,w1]=meshdom(ss,ww);
if vx<5,eval('[s1,w1]=meshdom(ss,ww);');
else, eval('[s1,w1]=meshgrid(ss,ww);');end
j=sqrt(-1);
z=s1+j*w1;
if isstr(f),g=eval(f);else,g=polyval(f,z)./polyval(p,z);end
g=abs(g);
%%%i=find(isnan(g)); %i=find(g==nan);
%%%if ~isempty(i),[m,n]=size(i);g(i)=1e15*ones(m,n);end
%%%i=find(finite(g));%i=find(g==inf);
%%%if ~isempty(i),[m,n]=size(i);g(i)=1e6*ones(m,n);end
th=2*pi*(-180:180)/360;z=exp(j*th);rez=real(z);imz=imag(z);


if isstr(f),h=eval(f);else,h=polyval(f,z)./polyval(p,z);end,
h=abs(h);hm=max(h);
if isstr(f),ff=[' = ',f];else,ff='';end,
w=s1.^2+w1.^2;
sc=1;g=g.*(g<=sc*hm)+sc*hm*(g>sc*hm);
h1=g;

i=find(w>r*r);if ~isempty(i),[m,n]=size(i);
if ml==3,g(i)=-0.2*hm*ones(m,n);else,g(i)=nan*ones(m,n);end,
end

disp('STRIKE A KEY BETWEEN PLOTS'),if ml==3,pause,end

if vx < 4, eval('clg');else,eval('clf');end
axis('square'),
if ml==3,mesh(g,[135 45]),else
eval('mesh(s1,w1,g)'),box='box';on='on';
%%%eval('axis([-r r -r r 0 ceil(hm)])')
eval('set(gca,box,on),view(vw),hold on,plot3(rez,imz,h)'),end
title(['z-plane magnitude of F(z)' ff]),pause,hold off %%Was [50 45]

if ml==3,hh=g;else,hh=h1;end
i=find(w>1);if isempty(i)
sc=0.5;g=g.*(g<=sc*hm)+sc*hm*(g>sc*hm);
if ml==3,mesh(g,[135 45]),else
eval('mesh(s1,w1,g),set(gca,box,on),view(vw)'),end
title('Blowup'),pause
else
[m,n]=size(i);if ml==3,hh(i)=-0.005*hm*ones(m,n);mesh(hh,[135 45]),else
hh(i)=nan*ones(m,n);
eval('mesh(s1,w1,hh),axis([-1 1 -1 1 0 ceil(hm)]);set(gca,box,on),view(vw)')
eval('hold on,plot3(rez,imz,h)'),
end,title('Inside unit circle'),pause
end

if ml==4,
if vx < 4, eval('clg');else,eval('clf');end
eval('axis([-1 1 -1 1 0 ceil(hm)]),set(gca,box,on),view(vw)')
eval('hold on,plot3(rez,imz,h),hold off'),title('around unit circle'),pause,end
plot(th/2/pi,h),title('magnitude along unit circle vs digital freq F')
axis('normal')