dtsim.m

ADSP TOOLBOX: Version 2.0 and gui m-files

function [yt,yzs,yzi] = dtsim(b,a,u,ic)
% DTSIM Recursive solution of difference equation.
%
%	[YT,YZS,YZI]= DTSIM(B,A,U,IC) Recursive solution of difference eqs 
%	B, A contain the coefficient arrays of the BACKWARD difference eq:
%           a0y[n]+a1y[n-1]+...+aNy[n-N] = b0x[n]+b1x[n-1]+...bNx[n-N]
%	IC = {y[-1], y[-2],...y[-N]} are initial conditions [Default: IC = 0].
%	U is the input sequence    Y is the output sequence. 
%       Y=Total, YZS=zero-state, YZI=zero-input response
%	Y = DTSIM(B,A,K) for SCALAR K returns impulse response h[n] to K terms.
%
%	DTSIM (with no input arguments) invokes the following example:
%
%	% Find and plot the response of y[n]+0.5y[n-1]=2x[n] with y[-1]=2 
%       % to the input x[n]=cos(0.4*n)u[n] for n=0 to n=40
%       >>nn = 0:40;uin = cos(0.4*nn);
%       >>yt = dtsim([2 0],[1 0.5],uin,2);
%       >>dtplot(nn,yt,'o')


% 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 dtsim,disp('Strike a key to see results of the example')
pause,nn=0:40;uin=cos(0.4*nn);ytt=dtsim([2 0],[1 0.5],uin,2);
dtplot(nn,ytt,'o'),return,end

%% HANDLES causal and noncausal systems

if nargin<4,ic=0;end,
while b(1)==0,b(1)=[];end
while a(1)==0,a(1)=[];end
if a(1)~=1,b=b/a(1);a=a/a(1);end
la=length(a);lb=length(b);lu=length(u);
b=[zeros(1,la-lb) b];
if lb>la,ls=int2str(lb-la);disp(' ')
disp(['WARNING: Noncausal system. Response index = Input index - ' ls ])
end

if nargin==3,
if lu==1,yt=filter(b,a,[1 zeros(1,u-1)]);
else, yt=filter(b,a,u);end
yzs=yt;yzi=[];return,end
%if lb>la,error('DTSIM requires causal system for nonzero IC'),return,end
ic=ic(:).';n=length(ic);m=la-1;
if n>m,error('Too many initial conditions'),return,end
ic=[ic zeros(1,m-n)];
if ~any(ic),yt=filter(b,a,u);yzs=yt;yzi=[];return,end
d=a;d(1)=[];d=-d;h=zeros(m);
for j=1:m,h(j,m-j+1:m)=ic(1:j)*d(j);end
nic=sum(h);nic=nic(m:-1:1);


if lb>la,ind=lb-la;else,ind=0;end
z=[zeros(1,ind)];nic=[z nic];
yzs=filter(b,a,u);
yzi=filter(b,a,0*u,nic);
if ind>0,add=[ic z];add=add(1:ind);add=fliplr(add);
yzi(1:ind)=add;end
yt=yzs+yzi;