tfplot.m

很多matlab的源代码

function [hm,hp,hf] = tfplot(ty,b,a,w,pl1,pl2)
% TFPLOT Frequency response of CT or DT systems
%
%	[M,P,F]=TFPLOT(TY,N,D,F,PL1,PL2) Generates and plots frequency response 
%	TY = 's', for CT systems,  TY = 'z', for DT systems 
%	N, D are the NUM and DEN of H(s) or H(z) in descending powers 
%	F=[F1 F2 M] = [frequency range (HERTZ) and M points (DEF: M=201)
%	If TY='z', F defaults to F = [0 0.5 201]
%	PL1 = 0 (for linear), PL1 = 1 (for log) frequency [default: PL1 = 0]
%	PL2 = 0 (no plots), PL2 = 1 (mag), PL2 = 2 (phase, deg), PL2 > 2 (both)
%	[Default: PL2=3 if no input arguments given, otherwise PL2=0]
%
%	TFPLOT('z',N,D) plots the resp of H(z) over F=[0 0.5] (201 points)
%	TFPLOT('s',N,D) plots freq resp of H(s) by estimating the freq. range 
%       Output Arguments:
%	M = magnitude,  P = phase (radians),  F = frequency (Hertz).
%
%	TFPLOT (with no input arguments) invokes the following example:
%
%	% Plot the spectra of H(s)=s/(s*s+2s+2) 
%	  >>tfplot('s',[1 0],[1 2 2])  %Use default estimated range


% 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 tfplot,disp('Strike a key to see results of the example')
pause,
 vx=matverch;
 if vx < 4, eval('clg');else,eval('clf');end
tfplot('s',[1 0],[1 2 2]),return,end

if nargin<6,pl2=0;end,
if nargin<5,pl1=0;end
if nargout==0,if pl2==0,pl2=3;end,end
if nargin<4,
if ty=='z',w=[0 0.5 201];
else,rt=roots(a);ri=0;i=find(imag(rt)~=0);
if ~isempty(i),ri=rt(i);rt(i)=[];end
rf=max([min(abs(rt));min(sqrt(abs(ri)))]);w=[0 rf 201];
end,end
if length(w)==3,n=w(3);else,n=201;end,
ww=w(1:2);w1=min(ww);w2=max(ww);
if ty=='s',if pl1==1,if w1==0
error('0 frequency is illegal on log scale'),return,end,end,end

if ty=='s',fq=' frequency in Hertz';else,fq=' digital frequency F';end
pd='phase in degrees vs';pdw=['unwrapped ' pd];
mm='magnitude';
if pl1==1,mm=[mm ' in dB vs'];else,mm=[mm ' vs'];end
if pl1==0 | ty=='z',w=(linspace(w1,w2,n))';end
if ty=='s',if pl1==1
w1=floor(log10(w1));w2=ceil(log10(w2));
w=(logspace(w1,w2,n))';
end,end
if ty=='s',s=sqrt(-1)*2*pi*w;else,s=exp(sqrt(-1)*2*pi*w);end
h=polyval(b,s)./polyval(a,s);
mag=abs(h);ph=angle(h);phh=ph;h=[mag ph];
ph1=unwrap(ph)*180/pi;ph=ph*180/pi;
if pl1==1,mag=20*log10(mag);end,
if ty=='z',pl1=0;end
if pl2>2,disp('STRIKE A KEY FOR & BETWEEN PLOTS'),pause(2),end
if pl2==1 | pl2>2
if pl1~=1,plot(w,mag),else,semilogx(w,mag),end,
grid,title([mm fq]),
if pl2>2,pause,end
end
if pl2>1
if pl1~=1,subplot(211),plot(w,ph),grid,title([pd fq]),
subplot(212),plot(w,ph1),grid,title([pdw fq]),
else
subplot(211),semilogx(w,ph),grid,title([pd fq])
subplot(212),semilogx(w,ph1),grid,title([pdw fq])
end
end,subplot,hold off
if nargout==0,return,else,hm=mag;hp=phh;hf=w;end