talg_det.m

Stack-based sequential decoder for M-QAM modulated MIMO-type problems, i.e., of fixed tree depth.

function [xHat,wHat,nv,nn,nf,nd] = Talg_det(T,m,xMax,y,H,cplx)
%
% [xHat,wHat,nv,nn,nf,nd] = Talg_det(T,m,xMax,y,H,cplx)
%
% Given inputs -
%
% T    : the maximum number of nodes to retain on the stack,
% m    : the problem dimension, i.e., the number of symbols to decode,
% xMax : a parameter specifying the admissible solution space, e.g.,
%        each element of xHat can take values in {-xMax+1,..,-1,0,1,..,xMax},
% y    : a real or complex target signal (column) vector,
% H    : a real or complex linear transform matrix, e.g., a MIMO channel, and
% cplx : flag specifying if xHat should be (0) real- or (1) complex-valued.
%
% Talg_det returns -
%
% xHat : a possibly suboptimal solution to   argmin    |y - H*x|^2, 
%                                          x \in Z_*^m
%        where Z_* is the set of integers such that -xMax < Z_* <= xMax, 
% wHat : the squared distance |y - H*xHat|^2, i.e., the solution node weight,
% nv   : the number of nodes expanded by the detector,
% nn   : the number of nodes left on the stack upon termination, 
% nf   : the (approximate) number of floating points operations performed
%        (excluding pre-processing, e.g., QR factorization, and sorting), and
% nd   : the number of nodes dropped during detector execution.
%
% The real-valued version of this priority first stack-based sequential
% decoding algorithm is based on a decoding tree of m+1 levels where each node
% has 2*xMax children.  At each stage, the node with the smallest weight is
% expanded by computing its first child and next sibling nodes.  Not all
% computed nodes are eventually expanded; in particular nv nodes are expanded
% and nn+nv+nd nodes are computed.  The stack maintains a list of those nodes
% that are known, i.e., have been computed, but that have not yet been
% expanded. 
%
% If T is exceeded during the replacement of a node by its first child and 
% next sibling, then the last node in the heap is dropped.  Note that this 
% node may not necessarily have the largest weight, since the stack is not 
% totally ordered.
%
% If T == 0, we use a heap having a reasonable maximum size of 32768 nodes.
% This limit is set due to practical computing considerations; at the caller's
% risk it can be exceeded by explicitly setting T to a value larger than 2^15.
%
% If xMax == 0, we do not apply (rectangular, or any) boundary control. In
% other words, Talg_det behaves as a lattice decoder.
%
% For more sophisticated operation, xMax may also be a vector of length m.
% Then each node at the beginning of stage j in the tree, where the root node
% is at the beginning of stage m and the leaf nodes are found at the end of
% stage 1, has 2*xMax(j) children.  Equivalently, symbol xHat(j) is drawn from
% {-xMax(j)+1,..,-1,0,1,..,xMax(j)}.
%
% If cplx == 1, we consider a tree of 2*m+1 levels with each node still having 
% 2*xMax children.  Alternatively, either xMax should be a complex-valued
% vector, or imag(xMax) will be taken to be equal to real(xMax).
%
% If either lattice reduction assistance or MMSE pre-processing are desired,
% these operations should be applied in advance of calling Talg_det.
%
% Notes:
%
% - Size(H,1) is expected to be equal to length(y).
% - Size(H,1) is expected to be greater than or equal to size(H,2).
% - Size(H,2) is expected to be equal to m.
% - Length(xMax) is expected to be equal to either 1 or m.
% - Because of Matlab memory issues that appear to be un-trappable (R14SP1), 
%   there is a practical upper bound of 2^15 encouraged for T.
%
% Example:
%
%   T    = 4;                        % max stack size
%   m    = 4;                        % 4x4 MIMO system
%   xMax = 2;                        % 16-QAM modulation
%   cplx = 1;
%   y    = randn(m,1)+i(randn(m,1);  % generate random complex target vector
%   H    = randn(m,m)+i*randn(m,m);  % generate random complex channel
%   [xHat,wHat,nv,nn,nf,nd] = Talg_det(T,m,xMax,y,H,cplx);
%
% Version 1.2, copyright 2006 by Karen Su.
%
% Please send comments and bug reports to [e-mail].
%

% Version 1.1
%   Minor cosmetic alterations.
% Version 1.2
%   Minor cosmetic alterations; scalar | changed to ||; added provision for
%   input xMax to be a vector (different order modulation for each symbol).

xHat = [];
wHat = Inf;
nv   = 0;
nn   = 0;
nf   = 0;
nd   = 0;

%
% Basic error checking
%
[n,mChk] = size(H);
if n ~= length(y) || mChk ~= m
  fprintf('Error: Decoding failed, invalid dimensions for H.\n');
  return;
end
if n < m
  fprintf('Error: Cannot solve under-determined problem, n (%i) < m (%i).\n', n, m);
  return;
end
if length(xMax) == 1
  xMax = xMax*ones(m,1);
elseif size(xMax,2) == m
  xMax = xMax.';
elseif size(xMax,1) ~= m
  fprintf('Error: Decoding failed, invalid dimensions for xMax.\n');
  return;
end

%
% Complex -> real conversion if necessary
%
if cplx == 1 
  m = 2*m;
  n = 2*n;
  y = [real(y);imag(y)];
  H = [[real(H),-imag(H)];[imag(H),real(H)]];
  if isreal(xMax)
    xMax = [xMax;xMax];
  else
    xMax = [real(xMax);imag(xMax)];
  end
end

%
% Reduce problem to square if H is overdetermined
%
[Hq,Hr] = qr(H);
wRoot   = 0;
yR      = Hq'*y;
if m < n
  wRoot = norm(yR(m+1:n))^2;
  yR    = yR(1:m);
end

%
% A node data structure consists of
%   w   : node weight
%   pR  : parent's residual target vector, elements 1:dim
%   xA  : accumulated symbol decision vector, length m-dim
%   dim : node/residual dimension, root (m) to leaf (0)
%   cn  : sibling ordinal
%   pw  : parent node weight
%   ut  : (scalar) target for siblings at this level
%   xp  : symbol decision of previous sibling
%
% We initialize the algorithm with the root node.  Note: wRoot and yR 
% defined previously
%
xA  = [];
dim = m;
cn  = 1;
pw  = [];
ut  = [];
xp  = [];

xMin = -xMax+1;
B    = 2*xMax;                          

%
% Compute first child node of root note
%
x0 = yR(dim)/Hr(dim,dim);                 % The first child is found by 
x1 = round(x0);                           % quantization and symbol-independent
if xMax(dim) > 0
  x1 = min([x1 xMax(dim)]);               % boundary control.
  x1 = max([x1 xMin(dim)]);
end

xA = x1;                                  % First child decision stored,
w  = wRoot+(yR(dim)-Hr(dim,dim)*x1)^2;    % its weight/cost is computed and

dim = dim-1;                              % its dimension and its
pR  = yR(1:dim);                          % parent's residual target stored.

cn  = 1;                                  % First child has ordinal 1.
pw  = wRoot;                              % Weight of parent node is stored.

ut  = yR(dim+1);               % (Scalar) target for siblings at child dim,
                               % i.e., costs of sibling symbol decisions are
                               % computed with respect to this reference.

xp  = x0;                      % For first child, store pre-quantization
                               % data as previous sibling "symbol" decision.
nf = nf + 6;
nv = nv + 1;

%
% Initialize storage
%
type = 1;                      % Heap to be initialized for lattice decoding.
if T == 0                      % Use maximum size heap if T not specified.
  T = 2^15;
end
hp = heap('init',m,T,type);

while dim > 0 
  nf = nf + 1;

  %
  % Compute first child node 
  %
  dimp = dim+1;                    % Derive dimension of parent node.

  yR = pR - Hr(1:dim,dimp)*xA(1);  % Compute current node residual target; 
                                   % note that this vector is only computed
                                   % upon expansion of a node.

  x0 = yR(dim)/Hr(dim,dim); 
  x1 = round(x0);                  % Quantization and boundary control
  if xMax(dim) > 0
    x1 = min([x1 xMax(dim)]);
    x1 = max([x1 xMin(dim)]);
  end
  xAc  = [x1; xA];                          % Store first child decision;
  wc   = w+(yR(dim)-Hr(dim,dim)*x1)^2;      % compute its weight and
  dimc = dim-1;                             % its dimension;
  yRc  = yR(1:dimc);                        % store parent residual target.
  cnc  = 1;                                 % First child ordinal is 1.
  pwc  = w;                                 % Store parent weight,
  utc  = yR(dim);                           % local target reference, and
  xpc  = x0;                                % pre-quantization data.
  nf   = nf + 8 + 2*dim;

  %
  % Replace root node with first child
  %
  ns = heap('replacemin',hp,wc,yRc,xAc,dimc,cnc,pwc,utc,xpc);
  %nswaps = nswaps + ns;
  nf = nf + 1;
  %heap('printcontents',hp);

  %
  % Compute next sibling if necessary
  %
  if xMax(dimp) == 0 || cn < B(dimp)
    xn = xA(1)-cn*sign(xA(1)-xp);  % The next sibling symbol decision can be
    if xMax(dimp) > 0
      if xn > xMax(dimp)           % derived from the current and previous ones,
        xn = xMax(dimp)-cn;        % followed by boundary control. Note that
      elseif xn < xMin(dimp)       % it is dependent on the child ordinal, but
        xn = xMin(dimp)+cn;        % is still symbol-independent.
      end
    end
    xAs  = [xn; xA(2:end)];                 % Store next sibling decision.
    ws   = pw+(ut-Hr(dimp,dimp)*xn)^2;      % Compute sibling node weight.
    dims = dim;                             % Siblings have same dimension,
    yRs  = pR;                              % same parent residual target,
    cns  = cn+1;                            % one larger child ordinal,
    pws  = pw;                              % same parent node weight, and
    uts  = ut;                              % same local target reference.
    xps  = xA(1);                           % Previous symbol decision is that
                                            % of current sibling.
    nf = nf + 10;

    %
    % Insert new node and check if T exceeded
    %
    nn = heap('getsize',hp);
    if nn == T
      nd = nd + 1;
    end
    ns = heap('insert',hp,ws,yRs,xAs,dims,cns,pws,uts,xps);
    %nswaps = nswaps + ns;
    nf = nf + 1;
  end
  nf = nf + 1;
  nv = nv + 1;
  %heap('printcontents',hp);

  %
  % Get next node from heap
  %
  [w,pR,xA,dim,cn,pw,ut,xp] = heap('getmin',hp);
  %nswaps = nswaps + ns;
  nf = nf + 1;
end
nf = nf + 1;

%
% Delete heap 
%
nn = heap('getsize',hp);
heap('delete',hp,'hp');

%
% Return smallest weight leaf node; an ML solution
%
wHat = w;
xHat = xA;

%
% Undo any complex -> real transformation that may have been done before
%
if cplx == 1
  m = m/2;
  xHat = xHat(1:m)+i*xHat(m+1:2*m);
end