sfunyst.m

数字通信第四版原书的例程

function  [sys, x0, str, ts]  = sfunyst(t, x, u, flag, ax, ltype, npts, dt)
%SFUNYST S-function storage scope using graph window.
%
%	This M-file is designed to be used in a SIMULINK S-function block.
%	It stores the previous input points using discrete states
%	and then plots this against time.
%
% 	S-function Syntax (see SFUNC):
%		[SYS, X0]  = SFUNYST(T,X,U,FLAG,AX,LTYPE,NPTS)
%	where:
%		AX - initial graph axis
%		LTYPE - line type (e.g. 'r','r-.','x') (See PLOT)
%		NPTS  - number of storage points
%
%	See also, PLOT, SFUNY, SFUNXY, LORENZ2.

%	Copyright (c) 1990-94 by The MathWorks, Inc.
%	Andrew Grace 1-27-92.
%	Revised Wes Wang 4-28-93, 8-17-93, 12-6-93.

if abs(flag) == 2 | flag == 9 	% A real time hit or end of simulation (flag 9)
	npts = x(7);
  
	len_u = length(u);
	% add handle here
	h_fig = x(8);
	if (h_fig <= 0)
		if length(ax) < 4
			disp('Axis not defined correctly; using current axis');
			ax = [0 1 0 1];
		end

		% Initialize graph
		[sl_name,block] = get_param;
		[n_b, m_b]= size(block);
		if n_b < 1
			error('Cannot delete block while simulating')
		elseif n_b > 1
			error('Something wrong in get_param, You don''t have the current SIMULINK')
		end;
		
		% find out if the graphics window exists
		ind = find(sl_name == setstr(10));
		dash = '_';
		sl_name(ind)=dash(ones(size(ind)));

		% find out if the graphics window exists
		Figures = get(0,'Chil');
		new_figure = 1;
		i = 1;
		while ((new_figure) & (i <= length(Figures)))
			if strcmp(get(Figures(i), 'Type'), 'figure')
				if strcmp(get(Figures(i), 'Name'), sl_name)
					h_fig = Figures(i);
					handles = get(h_fig,'UserData');
					new_figure = 0;
					h_axis = handles(1);
					set(0,'CurrentF',h_fig);
					set(h_axis,'Xlim',ax(1:2),'Ylim',ax(3:4),...
						'UserData',[]);

					%input number has been changed
					delete(get(h_axis,'Child'));
					h_plot = [];
					h_plot2 = [];
					set(h_axis,'NextPLot','add');
					set(h_fig,'NextPlot','add');
					coll = ltype;
					for k=1:len_u
						[col, coll]=strtok(coll,'/');
						if length(col) <= 1
							col = 'w-';
						end;
						h_plot(k) = plot(0,0,col,'Erasemode','none','Color',...
							col(1),'LineStyle',col(2:length(col)));
						h_plot2(k) = plot(0,0,col,'Erasemode','none','Color',...
							col(1),'LineStyle',col(2:length(col)));
						set(h_plot(k),'UserData',[]);
					end;
					set(h_fig,'UserData',[h_axis, h_plot, h_plot2 ax(2:4), abs(ltype)]);
					set(h_axis,'NextPlot','new');
					set(h_fig,'NextPlot','new');
				end;
			end;
			i = i + 1;
		end;

		% the figure does not exist, and this is an end of simulation (flag 9) call
		% exit early
		if new_figure & (flag == 9)
			sys = [];
			return;
		end;

		% the figure does not exist, create a new figure
		if new_figure
			h_fig = figure('Unit','pixel','Pos',[100 100 400 300], ...
				'Name',sl_name, 'Number','off');;
			set(0, 'CurrentF', h_fig);
			h_axis = axes('Xlim', ax(1:2), 'Ylim', ax(3:4),... %'XTickLabels',[],'YTickLabels',[],...
				'UserData',[],'Box','on');
			set(h_axis,'NextPlot','add');
			coll = ltype;
			for k=1:len_u
				[col,coll] = strtok(coll,'/');
				if length(col) <= 1
					col = 'w-';
				end;
				h_plot(k) = plot(0,0,col,'Erasemode','none','Color',...
					col(1),'LineStyle',col(2:length(col)));
				h_plot2(k) = plot(0,0,col,'Erasemode','none','Color',...
					col(1),'LineStyle',col(2:length(col)));
				set(h_plot(k),'UserData',[]);
			end;
			set(h_fig,'UserData',[h_axis, h_plot, h_plot2, ax(2:4), abs(ltype)]);
			set(h_axis,'NextPlot','new');
			set(h_fig,'NextPlot','new');
			xl = get(h_axis,'XLabel');
			set(xl,'String','Time (second)');
		end;
		set(h_axis,'XGrid','on','YGrid','on');
		set(h_fig,'Color',get(h_fig,'Color'));
		set(h_fig,'NextPlot','new');
		% Initial conditions
		% Infinity flags used to signal first point.
		x = [0, -Inf, ax(1:4), npts, h_fig, 0];
		set(h_fig, 'UserData',[h_axis, h_plot, h_plot2, ax(2:4), abs(ltype)]);
	end;

	if ((max(get(0,'Child') == h_fig) > 0) | flag == 9)
		handles = get(h_fig,'UserData');
		len_c = length(ltype);
		len_h = length(handles);
		axis_color_reset = 0;

		if (t>x(2)) | ((len_u*2 + len_c) ~= (len_h - 4))
			axis_color_reset = 1;
		elseif max([ax(2:4) abs(ltype)] ~= handles(2*len_u+2:len_h))
			axis_color_reset = 1;
		end;

		h_axis = handles(1);
		h_plot = handles(2:len_u+1);
		h_plot2 = handles(len_u+2:len_u+len_u+1);
		axss = [handles(len_u+len_u+3:len_u+len_u+4) npts];

		if (length(ax) == 4) & (length(npts) ==1)
			if min(axss==[ax(3:4), npts]) <1
				axis_color_reset = 1;
%				sys(9) = 4;                 %flag for replot
				set(h_fig, 'UserData',[h_axis, h_plot, h_plot2, ax(2:4), abs(ltype)]);
			end;
		end;

		discrete = any(get_param(gcb, 'sample time'));
		
		% Loop counter
		x_m = x(1);
		x(1) = x(1) + 1;

		if x(1) > npts
			x(1) = x(1) - npts;
		end;

		% h_axis UserData saves the time
		% h_plot UserData saves the previous input points

		sys = x;
		% Store last time and input point in circular buffer.
		buffer_t = get(h_axis,'UserData');
		len_b = length(buffer_t);
		if len_b < npts
			buffer_t = [buffer_t; t];
		else
			buffer_t(sys(1)) = t;
		end;
		set(h_axis,'UserData', buffer_t);
		buffer_y = [];
		for k=1:len_u
			tmp = get(h_plot(k),'UserData');
			if len_b < npts
				tmp = [tmp; u(k)];
			else
				tmp(sys(1)) = u(k);
			end;
			set(h_plot(k),'UserData',tmp);
			buffer_y = [buffer_y tmp];
		end;

		% If the input significantly exceeds the Y-axis
		% then count how many times this happens.
		if (1.1*u < sys(5)) | (1.1*u > sys(6))
			sys(9) = sys(9) + 1;
		end

		% If the time is greater than the X-axis or
		% the input has exceeded the Y-axis more than 3 times
		% or it is the end of simulation then rescale and replot.
		if (t > sys(2)) | (sys(9) > 3) | (flag == 9) | (axis_color_reset == 1)

			sys(9) = 0; % Set 'exceed Y-axis' counter to zero.

			if (sys(2) == -Inf) % Initial point
				if (flag == 9)
					sys = [];
				else
					sys(2) = ax(2);
				end
				return
			end

			newax = ax;
			if flag == 9
				factor = 1.2;
				sys(2) = t;
			elseif t > sys(2)
				% If time has exceeded axis limit then
				% be conservative in rescaling factor for Y axis:
				factor = 1.2;
				sys(2) = sys(2) + ax(2);
			else
				% otherwise double scaling:
				factor = 2;
			end
			newax(1,2) = sys(2);

			if ~isempty(buffer_t)
				newax(1,1) = min(buffer_t);
			end;

			% Work out new Y-axis ranges:
			minb = min(min(buffer_y));
			maxb = max(max(buffer_y));
			if ~isempty(minb)
				if minb < 0
					newax(1,3) = factor*minb;
				else
					newax(1,3) = (1/factor)*minb;
				end
				if maxb > 0
					newax(1,4) = factor*max(max(buffer_y));
				else
					newax(1,4) = (1/factor)*max(max(buffer_y));
				end
			else
				newax(1,3) = ax(3);
				newax(1,4) = ax(4);
			end;
			sys(5) = newax(1,3);
			sys(6) = newax(1,4);

			if newax(4) <= newax(3)
				if (newax(3) ~= 0)
					newax(4) = newax(3) * (1 + eps);
				else
					newax(4) = eps;
				end
			end;
			if newax(2) <= buffer_t(1)
				if (buffer_t(1) ~= 0)
					newax(2) = buffer_t(1) * (1 + eps);
				else
					newax(2) = eps;
				end
			end;
			

			len_b = length(buffer_t);
			if (sys(1) ~= len_b)  & (len_b == npts)
				tmp_t = [buffer_t(sys(1)+1:len_b); buffer_t(1:x(1))];
			else
				tmp_t = [buffer_t(1:len_b)];
			end;
			if (discrete)
				indx = [ 2:len_b; 2:len_b ];
				tmp_t = tmp_t([1; indx(:)]);
			end

			% Plot stored points and the y=0 line:
			for k=1:len_u
				if (sys(1) ~= len_b) & (len_b == npts)
					tmp = [buffer_y(sys(1)+1:len_b,k); buffer_y(1:x(1),k)];
				else
					tmp = [buffer_y(1:len_b,k)];
				end;
				if (discrete)
					indx = [ 1:(len_b - 1); 1:(len_b - 1) ];
					tmp = [tmp([indx(:); len_b])];
				end
				set(h_plot2(k), 'XData', tmp_t, 'YData',tmp);
			end;
			set(h_axis,'Xlim',newax(1:2),'Ylim',newax(3:4));
			set(h_fig,'Color',get(h_fig,'Color'));

		end		
		len_b = length(buffer_t);
		if (sys(1,1) == 1) & (len_b < npts) % Is it first point in cycle
			% Use none as Erasemode for faster plotting
			for k=1:len_u
				set(h_plot(k), 'XData',[t,t], 'YData',[u(k),u(k)]);
			end;
		else
			if x_m <= 0
				x_m = len_b;
			end;
			for k=1:len_u
				if (discrete)
					set(h_plot(k),'XData',[buffer_t(x_m),   t,               t   ],...
						      'YData',[buffer_y(x_m,k), buffer_y(x_m,k), u(k)]);
				else
					set(h_plot(k),'XData',[buffer_t(x_m),t],'YData',[buffer_y(x_m,k),u(k)]);
				end
			end;
			drawnow
		end

		if (sys(1) == npts) | (sys(1) == (npts + 1))
			sys(1,1) = 0;
		end
		if flag == 9
			sys = [];
		end
	else
		sys = x;
	end;
	
elseif flag  == 0, 	% Initialization

	% Return system sizes:
	test = length(u);
	if test <= 0
		test = -1;
	end;

	sys(1) = 0;			% 0 continuous states
%	sys(2) = 6 + 2*npt s+ 4;	% discrete states
	sys(2) = 9;			% discrete states
	sys(3) = 0;			% 0 outputs
	sys(4) = test;			% 1 input
	sys(5) = 0;			% 0 roots
	sys(6) = 0;			% no direct feedthrough
	sys(7) = 1;			% 1 sample time

	%x ==> [flag_for_replot, time_record, axes(x_l,x_h,y_l,y_h), npts, figure_handle, plot_control]
	if length(ax) ~= 4
		error('axis limit has not set correctly')
	end;
	if npts <=0
		npts = 1;
	end;
	x0 = [0, -Inf, ax(1), ax(2), ax(3), ax(4), npts, 0, 0];

	if (nargin > 7)
		ts = [dt, 0];
	else
		ts = [-1, 0];
	end
	
elseif (flag == 4) & (nargin > 7)
	if (dt > 0)
		sys = t + dt;
	else
		sys = [];
	end;
else
	% Other flag options not defined.
	sys = [];
end