sfunid.m

本书是电子通信类的本科、研究生辅助教材

function  [sys, x0]  = sfunid(t,x,u,flag,order,npts,HowOften,offset,ts,method)
%SFUNID an S-function which performs system identification.
%
%	This M-file is designed to be used in a SIMULINK S-function block.
%	It stores up a buffer of input and output points of the system
%	and then uses the System Identification Toolbox to to identify a linear 
%	model.
%	
%	The input arguments are
%	order:          the orders of the model (usually a vector)
%	npts:           number of points to use in the fft (e.g. 128)
%	HowOften:       how often to plot the ffts (e.g. 64)
%	offset:         sample time offset (usually zeros)
%	ts:             how often to sample points (secs)
%	method:         the method to use which may one of:
% 	               'ar', 'arx', 'oe', 'armax', 'bj', 'pem'.
%
%	The system id. block displays two plots: the time history of actual
%	data versus predicted data and the associated error terms.
%
%	See also: SIMIDENT, MASKIDENT, SFUNC.

%	Copyright (c) 1990-94 by The MathWorks, Inc.
%	Andrew Grace 5-30-91.


if abs(flag) == 2   % Flag equals 2 on a real time hit 
        % Is it a sample hit ?
	sys = zeros(npts+1,3);
	sys(:) = x;
        if rem(t - offset + 1000*eps, ts) < 10000*eps
		x(1,1) = x(1,1) + 1;
		sys(x(1,1),:) = u(:)';
		div = x(1,1)/HowOften;
		if (div == round(div))
			u = [sys(x(1,1)+1:npts+1,1);sys(2:x(1,1),1)];
			y = [sys(x(1,1)+1:npts+1,2);sys(2:x(1,1),2)];
			e = [sys(x(1,1)+1:npts+1,3);sys(2:x(1,1),3)];
			tvec = (t - ts * npts + ts * (1:npts));

			if length(method)==2 & all(method(1:2) == 'ar')
				u = [];
			end

			th = feval(method, [y u], order);

			% Conversion to transfer function form:
			[A,B,C,D,F]=polyform(th);
	
			disp('')	
			disp('Transfer function:')
			d = conv(A,F);
			bl = length(B); dl = length(d);
			nl = max(bl,dl);
			num = [B zeros(1,nl-bl)];
			den = [d zeros(1,nl-dl)];
			if exist('printsys')	% If have control toolbox
				eval('printsys(num,den,''z'')');
			else
				num, den
			end
			if length(den)
				yest = filter(num,den,u);
			else 
				yest = zeros(npts,1);
			end

			disp('')
			disp('Noise model:')
			d = conv(A,D);
			al = length(A);
			cl = length(C);
			dl = length(d);
			nl = max(al,dl);
			num = [C zeros(1,nl-cl)];
			den = [d zeros(1,nl-dl)];
			if exist('printsys')	% If have control toolbox
				eval('printsys(num,den,''z'')');
			else
				num, den
			end

			if length(den)
				yest = yest + filter(num,den,e);
			end
				
			clf
			subplot(2,1,1)
			plot(tvec,y,'-',tvec,yest,'--');
			xlabel('Time (secs)')
			title('Actual output ''-'' vs. predicted model output ''--''')

			err = yest-y;
			subplot(2,1,2)
			plot(tvec,err);
			xlabel('Time (secs)')
			title('Error in predicted model')
			text(0.16,0.13,['Norm of error: ',num2str(norm(err))], 'sc');
			drawnow
			
		end
		if sys(1,1) == npts
			x(1,1) = 1;
		end
		sys(1,1) = x(1,1);
        end
elseif flag == 4
% Return next sample hit
        ns = (t - offset)/ts;  % Number of samples
        sys = offset + (1 + floor(ns + 1e-13*(1+ns)))*ts;
elseif flag  == 0, 
    sys = [0;3*npts+3;0;3;0;0]; 
    x0 = [[1 0 0];zeros(npts,3)];
	if HowOften > npts
	error('The number of points in the buffer must be more than the plot frequency')
	end
	% Initialize graph
	hold off
	clf
	title('Working - please wait');
	axis off
	drawnow
else 
        sys = [];
end