sphere_decoding.m

两发两收的空间多路复用信道的matlab仿真程序。

function [transmit_detected,erasure]=sphere_decoding(R,H,N,ModScheme)
% R: received signal vector(column vector)
% H: MIMO channel matrix with nR*nT
% N: signal number to be detected, i.e. Tx Antenna number
% transmit_detected: decoded column vector
% erasure: 1-decoding sucess, 0-decoding failure

switch ModScheme
      case 'BPSK'
          amplitude=1;
          NumBit = 2;
          constellation = [-1, 1];
      case 'QPSK'
          amplitude=sqrt(2);
          NumBit = 4;
          constellation = [-1, 1];
      case '8PSK'
          constellation = exp(j*[pi/8,3*pi/8,5*pi/8,7*pi/8,-pi/8,-3*pi/8,-5*pi/8,-7*pi/8]);
          NumBit = 8;
      case '16QAM'    
          amplitude=sqrt(10);
          NumBit = 16;
          constellation = [-3, -1, 1, 3];
      case '64QAM'
          amplitude=sqrt(42);
          NumBit = 64;
          constellation = [-7 -5 -3, -1, 1, 3, 5, 7];
      case '256QAM'
          amplitude=sqrt(170);
          NumBit = 256;
          constellation = [-15 -13 -11 -9 -7 -5 -3, -1, 1, 3, 5, 7 9 11 13 15];
      otherwise
          disp('Unknown  modualtion scheme.')
end

len=length(constellation);
R=R*amplitude;
n=2*N;
r=[real(R)' imag(R)'];h=[real(H)' imag(H)';-imag(H)' real(H)'];

h_chol=chol(h*h');var_temp1=diag(diag(h_chol));var_temp2=kron(diag(var_temp1),ones(1,n));
q=(h_chol.*var_temp1.*var_temp1+h_chol-var_temp1)./var_temp2;
for k1=1:n
    cc(k1)=(norm(2*h(k1,:))/2)^2;
end
c=min(cc);
step=sqrt(c);
u_sup=[-0.4326   -1.6656    0.1253    0.2877   -1.1465    1.1909    1.1892   -0.0376];
for k2=1:100000000      % for different radius c
    d_square=c;
    T=zeros(1,n);T(n)=c;
    rol=r*pinv(h);
    s=rol;
    L=zeros(1,n);u=zeros(1,n);
    epslon=zeros(1,n);
    
    ii=n;
    L(ii)=floor(sqrt(T(ii)/q(ii,ii))+s(ii));
    u(ii)=ceil_odd(-sqrt(T(ii)/q(ii,ii))+s(ii),constellation,L(ii),len);
    
    for k1=1:10000000   % for some radius c
        %     if u(ii)>L(ii)
        if u(ii)>constellation(len)|u(ii)>L(ii)
            if ii==n
                break
            else
                ii=ii+1;
            end
        elseif ii>1     % Iteratively Update Bounds
            epslon(ii)=rol(ii)-u(ii);
            T(ii-1)=T(ii)-q(ii,ii)*(s(ii)-u(ii))^2;
            s(ii-1)=rol(ii-1)+sum(q(ii-1,ii:n).*epslon(ii:n));
            ii=ii-1;
            L(ii)=floor(sqrt(T(ii)/q(ii,ii))+s(ii));
            u(ii)=ceil_odd(-sqrt(T(ii)/q(ii,ii))+s(ii),constellation,L(ii),len)-2;
        else        % ii==1
            d_sup_square=T(n)-T(1)+q(1,1)*(s(1)-u(1))^2;
            if d_sup_square<d_square
                u_sup=u;
                d_square=d_sup_square;
                T(n)=d_sup_square;
                ii=n;
                L(ii)=floor(sqrt(T(ii)/q(ii,ii))+s(ii));
                u(ii)=ceil_odd(-sqrt(T(ii)/q(ii,ii))+s(ii),constellation,L(ii),len)-2;
            else
%                 u_sup=u;
            end     % end if d_sup_square
        end     % end if u(ii)
        u(ii)=u(ii)+2;
    end     % end of k1
    % k1
    % u_sup=u;
    lim=IS_LATTICE(constellation,n,u_sup,len);
    if lim==1
        erasure=0;  % decoding sucess
        break        
    elseif c>200000
        erasure=1;  % decoding failure
        break
    else
        c=c+2*step*sqrt(c)+step^2;  % increase search radius
    end % end if lim
end     % end of k2
k1;
k2;
u_sup=u_sup';
transmit_detected=u_sup(1:N)+j*u_sup(N+1:n);
transmit_detected=transmit_detected/amplitude;