ssresp.m

很多matlab的源代码

function ys = ssresp(ty,tfn,tfd,x,ty1)
% SSRESP Steady-state response to sinusoidal inputs in symbolic form.
%
%	YS = SSRESP(TY,NUM,DEN,X) Steady-state response in symbolic form 
%	NUM and DEN are coefficients of H=NUM/DEN in descending powers 
%	TY = 's'if NUM and DEN are for CT or H(s)=NUM/DEN
%	TY = 'z' if NUM and DEN are coefficients of DT or H(z)=NUM/DEN,
%	X = [K, w, r] where X(t) = Kcos(wt+r)  OR X(n) = Kcos(wn+r).
%	For N inputs, X = Nx3 matrix of the form [K1 w1 r1; K2 w2 r2; ...]
%	TY1=0 to NOT combine terms at identical frequencies (Def TY1=1)
%
%	YS = output in symbolic form (a string variable of t or n)
%	NOTE: YS may be evaluated using EVAL
%
%	SSRESP (with no input arguments) invokes the following example:
%
%	% Find and plot the ss response of H(s)=1/(s+2) to the input 
%	% x(t) = 2*cos(t - pi/6) - 4*sin(2*t + pi/4) 
%	  >>yst = ssresp('s',1,[1 2],[2 1 -pi/6;-4 2 -pi/4]) 
%	  >>t=0:0.05:10;yss=eval(yst);plot(t,yss)


% 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 ssresp,disp('Strike a key to see results of the example')
pause,yst=ssresp('s',1,[1 2],[2 1 -pi/6;-4 2 -pi/4]),
t=0:0.1:20;yss=eval(yst);
 vx=matverch;
 if vx < 4, eval('clg');else,eval('clf');end
plot(t,yss),axesn,grid,pause,return,end
 
%ADD STABILITY WARNINGS
if nargin<5,ty1=1;end
ln=length(tfn);ld=length(tfd);[m,n]=size(x);z=zeros(m,3);z(1:m,1:n)=x;x=z;n=3;
warn=0;if ty=='s',if tfd(ld)==0,warn=1;end,end
if ty=='s',b=real(roots(tfd));
if any(b>=0),warn=1;end,end
if ty=='z',if sum(tfd)==0,warn=1;end,end
if ty=='z',b=abs(abs(roots(tfd))-1);
if any(b<1e-4),warn=1;end,end
if warn==1,disp('WARNING: TF is not stable'),end

xa=x(:,2);i=find(xa<0);
if ~isempty(i),x(i,2)=-x(i,2);x(i,3)=-x(i,3);end
z=x;if m>1,[y,i]=sort(x(:,2));
for k=1:n,x(:,k)=z(i,k);end

%Combine inputs at identical frequencies if ty1=1
if ty1==1
m=0;j=sqrt(-1);z=[];
while ~isempty(y)
m=m+1;a1=0;w=y(1);i=find(y==w);
if w==0,if i==1,z=[z;x(1) 0 0];else,z=[z;sum(x(i,:))];end
else,a1=x(i,1).*exp(j*x(i,3));a1=sum(a1);z=[z;abs(a1) w angle(a1)];end
x(i,:)=[];y(i)=[];end,end
end

w=z(1,2);if w==0,i=2;
if ty=='s',p=[tfn(ln)/tfd(ld) 0 0];else,p=[sum(tfn)/sum(tfd) 0 0];end
else,i=1;p=[];end
j=sqrt(-1);
for k=i:m
w=z(k,2);if ty=='s',s=j*w;else,s=exp(j*w);end
h=polyval(tfn,s)./polyval(tfd,s);mag=abs(h);ph=angle(h);
if ~isempty(p),p=[p;mag w ph];else,p=[mag w ph];end,end

y=[z(:,1).*p(:,1) z(:,2) rem(z(:,3)+p(:,3),2*pi)];

%Now for the output
if ty=='s',nt='t';else,nt='n';end
a=num2str(y(1,1));if y(1,2)==0,ys=a;else
w=num2str(y(1,2));f=num2str(y(1,2)/2/pi);ps=num2str(y(1,3)/pi);py=eval(ps);
py=py.*(abs(py)>100*eps);
if py<0,s=[];elseif py>0,s='+';else,s=[];ps=[];end
ys=[a '*cos(2*pi*' f '*' nt s ps '*pi)'];end

for k=2:m
a=num2str(y(k,1));pa=eval(a);if pa~=0
if pa<0,s1=[];else,s1='+';end
w=num2str(y(k,2));f=num2str(y(k,2)/2/pi);ps=num2str(y(k,3)/pi);py=eval(ps);
py=py.*(abs(py)>100*eps);
if py<0,s=[];elseif py>0,s='+';else,s=[];ps=[];end
ys=[ys s1 a '*cos(2*pi*' f '*' nt s ps '*pi)'];end,end