jml.m

张贤达老师 通信信号处理4－9章的程序 这本书是张老师的经典之作。

function d=jml(y,R,A,alg,BITS)
% function d=jml(y,R,A,alg)
%
% finds the solution to the optimization problem
% d = arg min b'*A*R*A*b-2*y'*A*b for b constrained to binary vectors
%
% y = vector matched filter output
% R = signature correlation matrix
% A = diagonal amplitude matrix
% alg = 1 selects the standard brute force solution 
%         Important: you must define the global variable BITS which
%         is the K by 2^K matrix of all possible bit combinations 
%         (see makebinary.m). The brute force method may be faster for
%         smaller values of K, e.g. K<12.
% alg = 2 selects the depth-first branch and bound method. Much faster for
%         larger values of K.

% global BITS			% must be defined globally!
K=length(y);			% number of users
y=y(:);				% make sure we have a column
H=A*R*A;

if alg==1,			% brute force
  stat=2*BITS'*A*y-dot(BITS,H*BITS)';
  [junk,index]=max(stat);
  d=BITS(:,index);
else				% depth-first branch and bound
  U=chol(H(K:-1:1,K:-1:1));	% upper triangular
  L=U(K:-1:1,K:-1:1);		% lower triangular such that H=L'*L

  d=zeros(K,1);			% current best decision vector
  b=zeros(K,1);			% working decision vector
  y_tilde=L'\(A*y);		% whiten the MF outputs (fast inverse)
  V=zeros(K,K+1);		% residual whitened MF output vectors
  E=zeros(1,K+1);		% error

  k=0;				% initialize
  V(:,k+1)=-y_tilde;
  UPPER=inf;
  OPEN.bits=[];
  OPEN.level=[];

  done=0;			% flag to indicate when done
  skip3=0;			% flag to skip step 3
  % MAIN ALGORITHM BELOW
  while done==0,
    if skip3~=1,			% check to see if we skip this step
      k=k+1;				% step 3
      b=[b(1:k-1) ; -sign(V(k,k)); zeros(K-k,1)];
      if k<K,
        OPEN.bits=[OPEN.bits [b(1:k-1) ; sign(V(k,k)); zeros(K-k,1)]];
        OPEN.level=[OPEN.level k];
      end
    end
    V(:,k+1)=V(:,k)+b(k)*L(:,k);		% step 4
    E(k+1)=E(k)+V(k,k+1)^2;			% step 5
    if (E(k+1)>=UPPER)&(isempty(OPEN.bits)==0), % step 6
      % drop working decision vector and get new one from OPEN list
      b=OPEN.bits(:,end);		% new working bit vector
      k=OPEN.level(:,end);		% associated level (-1?)
      OPEN.bits=OPEN.bits(:,1:end-1);	% remove from list
      OPEN.level=OPEN.level(:,1:end-1);	% remove from list
      skip3=1;				% start at step 4
    elseif (E(k+1)<UPPER)&(k==K),	% step 7
      % found a solution better than UPPER but more tree to search
      d=b;				% update current best
      UPPER=E(k+1);			% update UPPER
      b=OPEN.bits(:,end);		% new working bit vector
      k=OPEN.level(:,end);		% associated level (-1?)
      OPEN.bits=OPEN.bits(:,1:end-1);	% remove from list
      OPEN.level=OPEN.level(:,1:end-1);	% remove from list
      skip3=1;				% start at step 4
    elseif (E(k+1)<UPPER)&(k~=K),	% step 8
      % just keep heading down the tree
      done=0;				% not done (probably not necessary)
      skip3=0;				% start at step 3
    elseif (k==K)&(E(k+1)<UPPER)&isempty(OPEN.bits), % step 9
      % tree has been exhausted - update current best and exit
      d=b;				% update current best solution
      done=1;				% flag to end loop
    else				% step 10 
      % tree has been exhausted - don't update current best and exit
      done=1;				% flag to end loop
    end
  end % while
  
end