function [v,xn,xmax,y] = simulateDSM(u,arg2,nlev,x0)%[v,xn,xmax,y] = simulateDSM(u,ABCD,nlev=2,x0=0)% or%[v,xn,xmax,y] = simulateDSM(u,ntf,nlev=2,x0=0)%%Compute the output of a general delta-sigma modulator with input u,%a structure described by ABCD, an initial state x0 (default zero) and %a quantizer with a number of levels specified by nlev.%Multiple quantizers are implied by making nlev an array,%and multiple inputs are implied by the number of rows in u.%%Alternatively, the modulator may be described by an NTF.%The NTF is zpk object. (The STF is assumed to be 1.)%The structure that is simulated is the block-diagional structure used by%zp2ss.m.fprintf(1,'Warning: You are running the non-mex version of simulateDSM.\n');fprintf(1,'Please compile the mex version with "mex simulateDSM.c"\n');if nargin<2    fprintf(1,'Error. simulateDSM needs at least two arguments.\n');    returnend% Handle the input argumentsparameters = {'u','arg2','nlev','x0'};defaults = [ NaN NaN 2 NaN ];for i=1:length(defaults)    parameter = char(parameters(i));    if i>nargin | ( eval(['isnumeric(' parameter ') '])  &  ...     eval(['any(isnan(' parameter ')) | isempty(' parameter ') ']) )	eval([parameter '=defaults(i);'])      endendnu = size(u,1);nq = length(nlev);if isobject(arg2) & strcmp(class(arg2),'zpk')    ntf.k = arg2.k;    ntf.zeros = arg2.z{:};    ntf.poles = arg2.p{:};    form = 2;    order = length(ntf.zeros);elseif isstruct(arg2)    if any(strcmp(fieldnames(arg2),'zeros'))	ntf.zeros = arg2.zeros;    else	error('No zeros field in the NTF.')    end    if any(strcmp(fieldnames(arg2),'poles'))	ntf.poles = arg2.poles;    else	error('No poles field in the NTF.')    end    form = 2;    order = length(ntf.zeros);elseif isnumeric(arg2)    if size(arg2,2) > 2 & size(arg2,2)==nq+size(arg2,1)	% ABCD dimesions OK	form = 1;	ABCD = arg2;	order = size(ABCD,1)-nq;    else	fprintf(1,'The ABCD argument does not have proper dimensions.\n');	if size(arg2,2) == 2		% Probably old (ver. 2) ntf form	    fprintf(1,'You appear to be using the old-style form of NTF specification.\n Automatic converstion to the new form will be done for this release only.\n');	    ntf.zeros = arg2(:,1);	    ntf.poles = arg2(:,2);	    form = 2;	    order = length(ntf.zeros);	else	    error('Exiting simulateDSM.')	end    endelse    error('The second argument is neither an ABCD matrix nor an NTF.\n');endif isnan(x0)    x0 = zeros(order,1);endif form==1    A = ABCD(1:order, 1:order);    B = ABCD(1:order, order+1:order+nu+nq);    C = ABCD(order+1:order+nq, 1:order);    D1= ABCD(order+1:order+nq, order+1:order+nu);else    [A,B2,C,D2] = zp2ss(ntf.poles,ntf.zeros,-1);	% A realization of 1/H    % Transform the realization so that C = [1 0 0 ...]    Sinv = orth([C' eye(order)])/norm(C); S = inv(Sinv);    C = C*Sinv;    if C(1)<0	S = -S;	Sinv = -Sinv;    end    A = S*A*Sinv; B2 = S*B2; C = [1 zeros(1,order-1)]; % C=C*Sinv;     D2 = 0;    % !!!! Assume stf=1    B1 = -B2;    D1 = 1;    B = [B1 B2];endN = length(u);v = zeros(nq,N);y = zeros(nq,N);if nargout > 1	% Need to store the state information    xn = zeros(order,N); endif nargout > 2	% Need to keep track of the state maxima    xmax = abs(x0);endfor i=1:N    y(:,i) = C*x0 + D1*u(:,i);    v(:,i) = ds_quantize(y(:,i),nlev);    x0 = A * x0 + B * [u(:,i);v(:,i)];    if nargout > 1	% Save the next state	xn(:,i) = x0;    end    if nargout > 2	% Keep track of the state maxima	xmax = max(abs(x0),xmax);    endendreturnfunction v = ds_quantize(y,n)%v = ds_quantize(y,n)%Quantize y to % an odd integer in [-n+1, n-1], if n is even, or% an even integer in [-n, n], if n is odd.%%This definition gives the same step height for both mid-rise%and mid-tread quantizers.if rem(n,2)==0	% mid-rise quantizer    v = 2*floor(0.5*y)+1;else 		% mid-tread quantizer    v = 2*floor(0.5*(y+1));end% Limit the outputfor qi=1:length(n)	% Loop for multiple quantizers    L = n(qi)-1;    i = v(qi,:)>L;     if any(i)	v(qi,i) = L;    end    i = v(qi,:)<-L;    if any(i)	v(qi,i) = -L;    endend