sfuncorr.m

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

function  [sys, x0]  = sfuncorr(t,x,u,flag,npts,HowOften,offset,ts,cross,biased)
%SFUNCORR an S-function which performs auto- and cross-correlation.
%
%	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
%	then plots the power spectral density of the input signal.
%	
%	The input arguments are:
%		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)
%		cross:     set to 1 for cross-correlation 0 for auto
%		biased:    set to 'biased' or 'unbiased'
%	
%	The cross correlator gives two plots: the time history, 
%	and the auto- or cross-correlation.
%
%	See also, SFUNC, XCORR, SFUNPSD.

%	Copyright (c) 1990-94 by The MathWorks, Inc.
%	Andrew Grace 5-30-91.
%	Revised Wes Wang 4-28-93, 8-17-93.

if abs(flag) == 2   % Flag equals 2 on a real time hit 
	% Is it a sample hit ?
	sys = zeros(npts+1+1,1+cross);
	sys(:) = x;
	h_fig = sys((1 + cross ) * (npts + 1 + 1));
	if h_fig <= 0
		% Graph initialization
		[sl_name,blocks] = get_param;
		[n_b, m_b]= size(blocks);
		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 exist
		ind = find(sl_name == setstr(10));
		for i = ind
			sl_name(i) = '_';
		end;
		Figures = get(0,'Chil');
		new_figure = 1;
		for i = 1:length(Figures)
			if strcmp(get(Figures(i), 'Type'), 'figure')
				if strcmp(get(Figures(i), 'Name'), sl_name)
					h_fig = Figures(i);
					handels = get(h_fig,'UserData');
					if (length(handels) == 7) & ~isempty(get(h_fig,'Child'))
						%handels = [h_axis1, h_plot1, h_titl1, h_axis2, h_plot2, h_titl2]
						for j=1:2
							set(handels(j*3-2),'Visible','off');
							set(handels(j*3-1),'XData',0,'YData',0,'Erasemode','none');
							set(handels(j*3),  'String','');
							xl = get(handels(j*3-2),'Xlabel');
							set(xl, 'String','');
						end;
						set(handels(7),'XData',0,'YData',0,'Erasemode','none');
						set(handels(3),'String','Working - please wait');
						new_figure = 0;
					end;
				end;
			end;
		end;
		if new_figure
			h_fig = figure('Unit','pixel','Pos',[100 100 500 400], ...
			'Name',sl_name);
			set(0, 'CurrentF', h_fig);
			handels(1) = subplot(211);
			zr = zeros(1+cross*(2-j),2);
			handels(2) = plot(0,0,'EraseMode','None');
			if cross
				hold on
				handels(7) = plot(0,0,'EraseMode','None');
				hold off
			else
				handels(7) = handels(2);
			end
			set(handels(1),'Visible','off');
			handels(3) = get(handels(1),'Title');
			set(handels(3),'String','Working - please wait');
			handels(4)=subplot(212);
			handels(5) = plot(0,0,'EraseMode','None');
			set(handels(4),'Visible','off');
			handels(6) = get(handels(4),'Title');
		end;
		nstates = (1 + cross ) * (npts + 1 + 1);
		sys(nstates) = h_fig;
		set(h_fig,'Nextplot','new');
		set(h_fig,'UserData',handels);
	end;
	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))
			handels = get(h_fig, 'UserData');
			lasta = [sys(x(1,1)+1:npts+1,1);sys(2:x(1,1),1)];
			if cross
				lastb = [sys(x(1,1)+1:npts+1,2);sys(2:x(1,1),2)];
				% Take real part to avoid small residual complex part
				y = real(xcorr(lasta,lastb,biased));
			else
				y = real(xcorr(lasta,biased));
			end
			tvec = (t - ts * sys(x(1,1)) + ts * (1:npts))';
			if cross
				set(handels(1),'Visible','on','Xlim',[min(tvec),max(tvec)],'Ylim', [min(min([lasta;lastb])), max(max([lasta;lastb]))]);
				set(handels(2),'XData',tvec,'YData',lasta)
				set(handels(7),'XData',tvec,'YData',lastb)
			else
				set(handels(1),'Visible','on','Xlim',[min(tvec),max(tvec)],'Ylim',[min(lasta) max(lasta)]);
				set(handels(2),'XData',tvec,'YData',lasta)
			end
			set(handels(3),'String','Time history')
			xl = get(handels(1), 'Xlabel');
			set(xl,'String','Time (secs)')
			set(handels(5), 'XData', tvec, 'YData', y(1:npts))
			set(handels(4),'Visible','on','Xlim',[min(tvec),max(tvec)],'Ylim',[min(y(1:npts)) max(y(1:npts))]);
			if cross
				set(handels(6),'String',['Cross correlation (', biased, ')'])
			else
				set(handels(6),'String',['Auto correlation (',biased, ')'])
			end
			xl = get(handels(4), 'Xlabel');
			set(xl,'String','Time (secs)')
			set(h_fig,'Color',get(h_fig,'Color'));
		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 		% Initialization
	%add one for figure handel
	nstates = (1 + cross ) * (npts + 1 + 1); 	% No. of discrete states
	sys = [0;nstates;0;(cross + 1);0;0];  		% Sizes of system (see SFUNC)
	if HowOften > npts
		error('The number of points in the buffer must be more than the plot frequency')
	end
	x0 = zeros(nstates,1);			% Initial conditions
	x0(1,1) = 1;				% Set counter to 1
	x0(nstates,1) = 0;
else 				% Other flag options ignored
	sys = [];
end

% end of sfuncorr