基于深度优先搜索算法的球形译码，具体文献可以参考03年information thoery上的那篇经典文章

function [zHat,wHat,nv,nf] = SEAL_det(mPhy,y,G,u,D,H,L); % % [zHat,wHat,nv,nf] = SEAL_det(mPhy,y,G,u,D,H,L) % % Given inputs - % % mPhy : the physical/spatial problem dimension, i.e., the number of transmit % antennas, the effective problem dimension is then m = 2*mPhy*L, % y : a real target signal (column) vector, % G : a real lattice generator matrix, \ LAST % u : a real lattice translation vector, > code % D : a real spherical shaping region squared radius, / % H : a real linear transform matrix, e.g., a block diagonal matrix of L % converted/stacked real MIMO channel matrices, and % L : the temporal problem dimension, i.e., the number of transmit periods % spanned by a codeword, % % SEAL_det returns - % % zHat : a possibly suboptimal solution to argmin |y - H*G*z|^2, % z \in Z_*^m % where Z_*^m is the set of integer vectors such that G*z is a valid % spherical LAST codeword, % wHat : the squared distance |y - H*G*zHat|^2, i.e., the solution node weight, % nv : the number of nodes expanded by the detector, % nf : the (approximate) number of floating points operations performed % (excluding pre-processing, e.g., QR factorization, and sorting), % % This real-valued depth-first stack-based sequential decoding algorithm uses a % decoding tree of m+1 levels where each node has an arbitrary and unknown % number of children. At each stage, the node under consideration is expanded % if its weight is less than the squared distance to nearest currently known % lattice point. This distance threshold is initially set to infinity. % % Because it is a depth-first traversal, we expand a node by computing its % first child. If it is a leaf node, clearly it cannot be further expanded. % In this case, we will have found a closer lattice point than that previously % known. Therefore we can adaptively reduce the distance threshold to reflect % this new discovery. % % If the weight of the node under consideration is larger than this distance % threshold, then the current search path is terminated because it cannot % possibly lead to a closest lattice point. Upon path termination, the next % node to be considered is the next sibling of its parent. % % Note that it is possible for a node to have no children, and also for the % next sibling to be invalidated by means of boundary control. % % Complex to real conversion by stacking should be done prior to calling % SEAL_det. In addition, if either lattice reduction assistance or MMSE % pre-processing are desired, these operations should also be applied in % advance. % % 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. % - The solution zHat is an integer vector; be sure to apply any necessary % scaling prior to calling SEAL_det so that this solution is appropriate. % % Example: % % mPhy = 2; % 2x2 MIMO system % L = 2; % code spans 2 time periods % m = 2*mPhy*L; % code matrix size is governed by the % % number of real space-time dimensions % y = randn(m,1); % generate random target vector % HPhyC = randn(mPhy,mPhy)+i*randn(mPhy,mPhy); % HPhyR = [real(HPhyC) -imag(HPhyC); imag(HPhyC) real(HPhyC)]; % H = kron(eye(L),HPhyR); % construct space-time x space-time channel % % load optimizedLASTcode; % get LAST code specification: G,u,D % [zHat,wHat,nv,nf] = SEAL_det(mPhy,y,G,u,D,H,L) % % Version 1.0, copyright 2006 by Karen Su. % % Please send comments and bug reports to [email protected]. % global w; global z; global cn; global zp; global d; global lb; global ub; global nf; zHat = []; wHat = Inf; nv = 0; nf = 0; % % A small number % EPS = 0.001; % % Basic error checking % [n,mChk] = size(H); m = 2*mPhy*L; % size of the channel matrix is 2*mPhy*L x 2*nPhy*L 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 % % Composite (effective) linear transform matrix for LAST decoding is % formed by taking the product of the MIMO channel matrix and the LAST % code generator matrix % F = H*G; % % Reduce problem to square if H is overdetermined; also factorize code % generator % [Fq,Fr] = qr(F); [Gq,Gr] = qr(G); wRoot = 0; yR = Fq'*(y-H*u); if m < n wRoot = norm(yR(m+1:n))^2; yR = yR(1:m); end % % A basic node data structure consists of % w : node weight % pyR : parent's residual target vector, elements 1:dim == py(1:dim,dim) % zA : accumulated symbol decision vector, length m-dim == z(dim+1:dim) % dim : node/residual dimension, root (m+1) to leaf (1) % cn : sibling ordinal % pw : parent node weight == w(dim+1) % ut : (scalar) target for siblings at this level == py(dim+1,dim) % zp : symbol decision of previous sibling % % In addition the following components are needed to perform boundary control % for spherical LAST decoding % paR : parent's shaping region residual offset vector, elements 1:dim == pa(1:dim,dim) % pd : parent's shaping region residual squared radius == d(dim+1) % lb : lower bound of admissible range for siblings at this level == lb(dim) % ub : upper bound of admissible range for siblings at this level == ub(dim) % uo : (scalar) offset for siblings at this level == pa(dim+1,dim) % % The algorithm requires only a fixed-size block of memory; we initialize it % with the root node. % w = zeros(1,m+1); py = zeros(m,m); z = zeros(1,m); dim = m; cn = zeros(1,m); zp = zeros(1,m); pa = zeros(m,m); d = zeros(1,m+1); lb = zeros(1,m); ub = zeros(1,m); dimp = dim+1; w(dimp) = wRoot; d(dimp) = D+EPS; % % Compute first child of root node % aR = -Gq'*u; pa(1:dim,dim) = aR; % Store parent (root) residual offset vector; tmp1 = pa(dim,dim)/Gr(dim,dim); % compute and store lower and upper bounds tmp2 = sqrt(d(dimp))/abs(Gr(dim,dim)); % for the range of admissible symbol decisions lb(dim) = ceil(tmp1-tmp2); % for siblings at this level. ub(dim) = floor(tmp1+tmp2); nf = nf + 9; py(1:dim,dim) = yR; % Store parent (root) residual target. compFirstChild(dim,dimp,lb(dim),ub(dim),Fr(dim,dim),Gr(dim,dim),py(dim,dim),pa(dim,dim)); backtrack = 0; while dim <= m if w(dim) < wHat % If the node under consideration has a smaller weight if dim == 1 % and it is a leaf node, then update zHat, wHat wHat = w(dim); zHat = z; backtrack = 1; % To backtrack, we have to compute the next sibling dim = dimp; % of the current node's parent dimp = dimp+1; if cn(dim) <= ub(dim)-lb(dim) % Check that there is at least one more sibling backtrack = 0; compNextSibling(dim,dimp,lb(dim),ub(dim),Fr(dim,dim),Gr(dim,dim),py(dim,dim),pa(dim,dim)); end nf = nf + 1; % % If the node under consideration has a sufficiently small weight and is % not a leaf node: % - then either we are backtracking (backtrack == 1, and so compute the % next sibling of its parent) % elseif backtrack & cn(dim) > ub(dim)-lb(dim) % No more siblings %backtrack = 1; dim = dimp; if dim <= m % Check that we have not backtracked to root dimp = dimp+1; if cn(dim) <= ub(dim)-lb(dim) % Check that there is at least one more sibling backtrack = 0;