球形译码的MATLAB仿真

1.问题描述：

球形译码的基本思想是在以一个矢量x 为中心的半径为d的多维球内搜索格点，通过限制或者减少搜索半径从而减少搜索的点数,进而使得计算时间减少。球形译码算法带来的优点在于它不需要象传统的最大似然译码算法那样需要在整个格内对所有的格点进行搜索,而只需要在一个事先设定的有限球形区域进行搜索,如果该区域所包含的点数相对于整个格内的总点数是相当小的,搜索时间就会大大减少。

影响球形译码的关键问题有：(1) 怎样选择搜索半径d。如果d太大,则球内会包含太多的点,复杂度就会接近或者达到最大似然译码的指数级复杂度。如果d 太小,则球内可能一个格点都不包含,那么球形译码算法将得不到合理的解。(2) 怎样才能判断一个点是否在球内。如果这种判断需要借助每一个格点和矢量之间的距离来判断的话,那么这种方法就不太理想,因为我们需要考察所有的点,所产生的计算量也是指数级的。

球形译码解决了第2个问题，此处均考虑信号为实数，因为复数可以通过增加一倍的维数，将实部和虚部分开，要判断一个点是否在半径为d的m维球内比较困难。若将m变为1，则从球退化为一个间距，这个点就相当于某根天线发送信号的实部或虚部，这样操作就简单很多，可以知道这个点是否在这个距离内。多根发送天线上的信号的实部和虚部分成很多维，每一维上有可能取值。球形译码算法相当于构建了一棵树，树的第k层节点对应的是落在半径为d，维数为k的球内的格点。

 

2.部分程序：

 

function y1 = spheredecodetoML(rev, C, H, codebook)
[m,n]=size(H);      %m为Nr，n为Nt
T = zeros(1,n);   

[Q,R] = qr(H);    %QR分解
R = R(1:n, :);    %R= R（1：Nt，：）
ro = Q'*rev;  %y'=Q1'y
for ii = 1:n  %R为上三角矩阵
    if R(ii,ii) == 0
        R(ii,ii) = 0.001;
    end
    for jj = ii+1:n
        R(ii,jj) = R(ii,jj)/R(ii,ii); %为后面的运算做铺垫
    end
    ro(ii) = ro(ii)/R(ii,ii);%相当于ym/Rmm，ym-1/Rm-1...
    R(ii,ii) = R(ii,ii)^2;
end
d = C;    %初始半径
T(n) = C; %T相当于是dm，m-1，m-2...的值，每维的半径
S = ro;   %相当于ym/Rmm
cn = length(codebook); %码本长度，4qam的话，就为2
len = zeros(1,n);      %每维可能值的个数 
Y = zeros(n, cn);      %存储
ind = zeros(1,n);      %作为指示标志
y1 = zeros(n,1);       %存储最后的数据输出
yt = [];               %只是一个中间存储变量，存储每层的可能结果 
while 1
    k = n;                      %从第n维开始
       while k >= 1
             U = sqrt(T(k)/R(k,k))+S(k);%上限值
             ii = cn;
             while ii >= 1 && U < codebook(ii)-10^(-11)
                   ii = ii - 1;
             end
             L = -sqrt(T(k)/R(k,k))+S(k);%下限值
             jj = 1;
             while jj <= cn && L > codebook(jj)+10^(-11)
                  jj = jj + 1;
             end
             len(k) = ii-jj+1;    %第K维的可能值的个数
             Y(k,:) = zeros(1, cn);
             yt = [];
             for kk = 1:len(k)
                 yt = [yt codebook(jj+kk-1)]; %第K维的可能值
             end
            [yt1, idex] = sort(abs(yt - S(k)));%按行升序排序
            yt = yt(idex);       %已经排序后的结果
            for kk = 1:len(k)
                Y(k,kk) = yt(kk);%从偏差最小的点开始存储，按序，存储信号的可能取值
            end
        
            ind(k)=0;%作为指示标志
            d1 = d+1;
            while d1 > d - 10^(-11)
                  ind(k) = ind(k) + 1;
                  while ind(k) > len(k)
                        if k == n
                           if y1~=zeros(n,1)
                              return       %结束，返回y1    
                           else
                              C=2*C;
                              y1=spheredecodetoML(rev, C, H, codebook);
                              return  
                           end   
                        end
                        
                        k = k+1;   %表示在第k维中没有找到值，则返回到k+1维中
                        ind(k) = ind(k) + 1; %在第k+1维中从下一个点开始搜
                  end
                  if k > 1
                     T(k-1) = T(k) - R(k,k)*(S(k)-Y(k,ind(k)))^2;%求dm-1的值，Y(k,ind(k))相当于第k维的取值   采用第二范数
                     %T(k-1)=max(T(k),R(k,k)*(S(k)-Y(k,ind(k)))^2);                                    % 采用无穷范数
                     tmp = 0;
                     for ii = k:n
                         tmp = tmp + R(k-1, ii)*Y(ii,ind(ii));
                     end
                     S(k-1) = ro(k-1) - tmp;  %
                     k = k - 1;
                     break;     
                  end
                  d1 = T(n) - T(1) + R(1,1)*(S(1) - Y(1, ind(1)))^2;%k=1时的d1值
              end %while d1 > d
              if d1 < d-10^(-11)
                 for ii=1:n
                     y1(ii) = Y(ii, ind(ii));
                 end
                 d = d1;
                 T(n)=d1;
                 break;
              end
          end %while k>=1
%     end   
end %while 1 k=n
 

3.仿真结论：

D104