ltr.m

很多matlab的源代码

function [N,D] = ltr(x)
% LTR Laplace Transform.
%
%	[N,D] = LTR(X) computes the Laplace transform. 
%     	For x(t) = Kexp(-at)(t^p)cos(wt+r), X has the form 
%	X = [K a p w r].  This is a 5 element ROW vector.
%	For N terms, X should be an Nx5 matrix with N such rows.
%
%	[N,D] =[num,den] of X(s).
%
%	LTR (with no input arguments) invokes the following example:
%
%  	% Find the LT of 4exp(-3t)+2texp(-3t) in rational function form
% 	  >>[NN,DD] = ltr([4 3 0;2 3 1])


% 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 ltr,disp('Strike a key to see results of the example')
pause,[NN,DD]=ltr([4 3 0;2 3 1]),return,end

% Remove rows whose first value is 0
xxa=x(:,1);i=find(xxa==0);
if ~isempty(i),x(i,:)=[];end
if isempty(x),N=[];D=[];return,end

[m,n]=size(x);z=zeros(m,5);z(1:m,1:n)=x;x=z;
xa=x(:,3);
if any(xa<0),error('powers of t cannot be negative'),return,end
xa=x(:,4);i=find(xa<0);
if ~isempty(i),x(i,4)=-x(i,4);x(i,5)=-x(i,5);end
%Sort by roots
z=x;
if m>1,xa=x(:,2);[mm,i]=sort(xa);for k=1:m,x(k,:)=z(i(k),:);end,end
while x(1,1)==0,x(1,:)=[];m=m-1;end
%Sort by frequency
z=x;za=z(:,2);x=[];
while ~isempty(za),
i=find(za==za(1));l=length(i);
zc=z(1:l,4);[ss,s]=sort(zc);
for k=1:l,x=[x;z(s(k),:)];end
za(1:l)=[];z(1:l,:)=[];
end
%Sort by power
z=x;za=pi*z(:,2)+sqrt(-1)*z(:,4);x=[];
while ~isempty(za),i=find(za==za(1));
l=length(i);zc=z(1:l,3);[ss,s]=sort(zc);
for k=1:l,x=[x;z(s(k),:)];end
za(1:l)=[];z(1:l,:)=[];
end
tf=ltrn(x(1,:));t2=tf;
if m>1,
for j=2:m,
if x(j,1)~=0
t1=t2(2,:);t2=ltrn(x(j,:));
den=conv(tf(2,:),t2(2,:));
n1=conv(tf(1,:),t2(2,:));
n2=conv(tf(2,:),t2(1,:));
if x(j,2)==x(j-1,2) & x(j,4)==x(j-1,4)
if x(j,3)<x(j-1,3),dc=t2(2,:);else,dc=t1;end
while dc(1)==0,dc(1)=[];end
n1=deconv(n1,dc);n2=deconv(n2,dc);den=deconv(den,dc);
end
l1=length(n1);l2=length(n2);ld=length(den);
n1=[zeros(1,ld-l1) n1];n2=[zeros(1,ld-l2) n2];tf=[n1+n2;den];
end,
end,
end
N=tf(1,:);D=tf(2,:);