multlinr.m

经典通信系统仿真书籍《通信系统与 MATLAB (Proakis)》源代码

function [sys,x0] = multlinr(t,x,u,flag,input,inpidx,table)
%MULTLINR is a multi-dimensional lookup table in S-function form.
%	This M-file is designed to be used in a SIMULINK S-function
%	block.  It performs a multi-dimensional lookup on the table
%	passed in through the parameter list. It makes use of the
%	variables INPUT, INPIDX, and TABLE, where 
%	INPUT(INPIDX(i-1):INPIDX(i)) is the i-th input vector.
%	The content in TABLE is arranged as FORTRAN storage for 
%	INPIDX(1) x INPIDX(2) x ... x INPIDX(N) variables, where
%	N = length(INPIDX).
%
%	The section of code under the FLAG == 3 call returns to SYS
%	a linearly interpolated intersection from the array TABLE.
%
%	All INPUT(INPIDX(i-1):INPIDX(i)) must increase monotonically.
%	Extrapolation is used for values outside the limits of the
%	input index based on the four closest points.

% Wes Wang 1/11/93
% Copyright (c) 1990-93 by The MathWorks, Inc.
% $Revision: 1.8 $  $Date: 1993/06/09 22:09:40 $

if flag == 3,
	% multilinear interpolation
        % Use "find" to indicate what xindex values x0 lies between
        N = length(inpidx);
	% k1: the beginning of i-th input
	k1 = 1;
	% k2: the end of i-th input
	k2 = inpidx(1);
	for i=1:N-1
		k1 = [k1 k1(i)+inpidx(i)];
		k2 = [k2 k2(i)+inpidx(i+1)];
	end;
	% find the index where u(i) is between xptrlo(i) and xptrhi(i)
	% except the case where u(i) is out of the range
	% In the case of out of range, xptrlo(i) = 1 or xptrhi(i) = inpidx(i)
	for i=1:N
		tmp = min(find(input(k1(i):k2(i)) > u(i)));
		if length(tmp) == 0 % when the input is larger than index
			tmp = inpidx(i);
		elseif tmp <= 1     % when the input is smaller than index
			tmp = 2;
		end;
		xptrhi(i) = tmp;
        end;
	xptrlo = xptrhi-1;
	% xlo and xhi are input at index related to xptrlo and xptrhi
	xlo = input(xptrlo+k1-1);
	xhi = input(xptrhi+k1-1);
	% xratiof ratio for lo
	xratiof = (u' - xlo)./(xhi - xlo);
	% xratiof ratio for hi
	xratiob = 1 - xratiof;
        %computing
	% sys = sum  (prod xratio?i(i))*f(x?(1),x?2(2),...,x?N(N))
	% where the combination of ?1,?2,....,?N exhausts all of possible 
	% combination binary, with 0-->lo, 1-->hi.
	sys = 0;
	tmp = 0;
	NN = 2^N;
	k3 = 2.^[0:N];
	for i = 1:NN
		for j = 1:N
			% using tmp as a binary test, if i-th bit is
			% 1, use low part, is 0, use high part.
			if rem(tmp,k3(j+1)) >= k3(j)
				xratio(j) = xratiob(j);
				xptr(j)   = xptrlo(j);
			else
				xratio(j) = xratiof(j);
				xptr(j)   = xptrhi(j);
			end;
		end;
		% the index for the table
		k4 = xptr(N) - 1;
		for j = 1:N-1
			k4 = k4 * inpidx(N - j) + xptr(N - j) - 1;
		end;
		k4 = table(k4 + 1);
		for j = 1:N
			k4 = k4 * xratio(j);
		end;
		sys = sys + k4;
		tmp = tmp+1;
	end;
	% ***************************************
	% Here ends the multilinear interpolation
	% ***************************************

elseif flag  == 0,
	% This part takes care of all initialization; it is used only once.
	x0 = [];
	% The system has no states, one output, and two inputs.
	sys = [0 0 1 length(inpidx) 0 1]';

else 
	% Flags not considered here are treated as unimportant.
	% Output is set to [].
	sys = [];

end