www.gusucode.com > MPSK,误码率性能仿真源码程序 > MPSK,误码率性能仿真源码程序/IET_MATLAB/functions/MSDSD.m

    % fucntion Multiple symbol differential sphere detection
% this code is based on Pseudocode for MSDSD, from paper [1]
% [1] Multiple symbol differential sphere decoding, L. Lampe, IEEE 2005

% inputs
% U : N*N upper triangular matrix 
% M : constellation size, example M=4 
% R : intial radious, example R=100 
%output
% shat : Maximum-likelihood decision 

% note that here (i and j) are used as a variable but (1i)=sqrt(-1)
function shat=MSDSD(U,M,R)

N=length(U);
shat=[];
m=zeros(1,N);
step=zeros(1,N);
n=zeros(1,N);
k=zeros(1,N);
d=zeros(1,N+1);

d(N)=abs(U(N,N));% initialize length
s(N)=1;%fix last component of s
i=N-1;%start with component i=N-1
k(i)=U(N-1,N);%sum of last N-i components

% find best cndidate point for ith component
[m(i),step(i),n(i)]=findBest(k(i),U(i,i),M);

% loop
while (1)  
    % update squared length Eq.7
    d(i)=sqrt(abs(k(i)+U(i,i)*exp(1i*2*pi*m(i)/M))^2+d(i+1)^2);
    
    % check sphere radius Eq.7 and constellation size
    if d(i)<R && n(i)<=M
        s(i)=exp(1i*2*pi*m(i)/M);%store candidate component
        if i~=1 % component 1 not reached yet
            i=i-1; %move down
            % add last N-i components,Eq.6 
            j=i+1:N;
            k(i)=sum(U(i,j).*s(j));
            [m(i),step(i),n(i)]=findBest(k(i),U(i,i),M);
        else  % first component reached
            shat=s; % best point so far
            R=d(i); % update sphere radius
            i=i+1;  % move up
            % next point examined for component i
            [m(i),step(i),n(i)]=findNext(m(i),step(i),n(i));
        end
    else    
        condition = true;
        while condition
            if i==N-1,return,end              
            i=i+1;
            condition=(n(i)==M);
        end
        [m(i),step(i),n(i)]=findNext(m(i),step(i),n(i));
    end
end    

end
% subfunction
    function [mi,stepi,ni]=findBest(ki,uii,M)
        p=M*(angle(-ki/uii))/2/pi;
        mi=round(p);
        ni=1;
        stepi=sign(p-mi);
        
    end
    

% subfunction
    function [mi,stepi,ni]=findNext(mi,stepi,ni)
        mi=mi+stepi;
        stepi=-stepi-sign(stepi);
        ni=ni+1;
    end