logmapo_all.m

迭代多用户检测的MATLAB详细实现,非常全面,主要基于Vincent Poor和王晓东提出的算法实现,绝对的好东西.

function L_all = logmapo_all(rec_s,g,L_a,ind_dec)
% Copyright Nov 1998, Yufei Wu
% MPRG lab, Virginia Tech.
% for academic use only

% Log_MAP algorithm using straightforward method to compute branch metrics
% no approximation is used.
% Can be simplified to Max-Log-MAP by using approximation ln(e^x+e^y) = max(x,y).
% Input: rec_s: scaled received bits. 
%               rec_s = 0.5 * L_c * yk = ( 2 * a * rate * Eb/N0 ) * yk
%        g: code generator for the component RSC code, in binary matrix form.
%        L_a: a priori info. for the current decoder, 
%               scrambled version of extrinsic Inftyo. of the previous decoder.
%        ind_dec: index of decoder. Either 1 or 2. 
%               Encoder 1 is assumed to be terminated, while encoder 2 is open.
%
% Output: L_all: log-likelihood ratio of the symbols. Complete information.


% Total number of bits: Inftyo. + tail
L_total = length(rec_s)/2;
[n,K] = size(g); 
m = K - 1;
nstates = 2^m;          % number of states in the trellis

% Set up the trellis
[next_out, next_state, last_out, last_state] = trellis(g);

for state=1:nstates
    temp2=last_out(state,2:2:4);
  for bit=0:1
    b2_last_state(state,bit+1)=last_state(state,find(temp2==2*bit-1));
  end
   if (last_out(state,2)==-1)
       b2_last_out(state,:)=last_out(state,:);
   else
       b2_last_out(state,1:2)=last_out(state,3:4);
       b2_last_out(state,3:4)=last_out(state,1:2);
   end
end


Infty = 1e100;

% Initialization of Alpha
Alpha(1,1) = 1; 
Alpha(1,2:nstates) = zeros(1,nstates-1);

% Initialization of Beta
if ind_dec==1
   Beta(L_total,1) = 1;
   Beta(L_total,2:nstates) = zeros(1,nstates-1); 
elseif ind_dec==2
   Beta(L_total,1:nstates) = zeros(1,nstates);
else
   fprintf('ind_dec is limited to 1 and 2!\n');
end

% Trace forward, compute Alpha
for k = 2:L_total+1
    for state2 = 1:nstates
      gamma = zeros(1,nstates);
      gamma(last_state(state2,1)) = (0.5*(1-tanh(rec_s(2*k-3)/2))*0.5*(1+tanh(rec_s(2*k-2)/2)*last_out(state2,2)));%-log(1+exp(L_a(2*k-3)));
      gamma(last_state(state2,2)) = (0.5*(1+tanh(rec_s(2*k-3)/2))*0.5*(1+tanh(rec_s(2*k-2)/2)*last_out(state2,4)));% + L_a(2*k-3)-log(1+exp(L_a(2*k-3)));

 %     if(sum(exp(gamma+Alpha(k-1,:)))<1e-300)
 %        Alpha(k,state2)=-Infty;
 %     else
         Alpha(k,state2) = sum(  gamma.*Alpha(k-1,:) ) ;  
         %     end   
     end
    tempmax(k) = sum(Alpha(k,:));
    Alpha(k,:) = Alpha(k,:) / tempmax(k);
end     

% Trace backward, compute Beta
for k = L_total-1:-1:1
  for state1 = 1:nstates
     gamma = zeros(1,nstates);
     gamma(next_state(state1,1)) = (0.5*(1-tanh(rec_s(2*k+1)/2))*0.5*(1+tanh(rec_s(2*k+2)/2)*next_out(state1,2)));%-log(1+exp(L_a(2*k+1)));
     gamma(next_state(state1,2)) = (0.5*(1+tanh(rec_s(2*k+1)/2))*0.5*(1+tanh(rec_s(2*k+2)/2)*next_out(state1,4)));% +L_a(2*k+1)-log(1+exp(L_a(2*k+1)));
%     if(sum(exp(gamma+Beta(k+1,:)))<1e-300)
%        Beta(k,state1)=-Infty;
%     else
        Beta(k,state1) = sum(gamma.*Beta(k+1,:));
        %     end   
  end
  Beta(k,:) = Beta(k,:) / tempmax(k+1);
end

% Compute the soft output, log-likelihood ratio of symbols in the frame
for k = 1:L_total
  for state2 = 1:nstates
     gamma0 = (0.5*(1-tanh(rec_s(2*k-1)/2))*0.5*(1+tanh(rec_s(2*k)/2)*last_out(state2,2)))*L_a(2*k-1);%-log(1+exp(L_a(2*k-1)));
     gamma1 = (0.5*(1+tanh(rec_s(2*k-1)/2))*0.5*(1+tanh(rec_s(2*k)/2)*last_out(state2,4)))*L_a(2*k-1);% +L_a(2*k-1)-log(1+exp(L_a(2*k-1)));
     temp0(state2) = gamma0 * Alpha(k,last_state(state2,1)) * Beta(k,state2);
     temp1(state2) = gamma1 * Alpha(k,last_state(state2,2)) * Beta(k,state2);
     
     gamma00 = (0.5*(1+tanh(rec_s(2*k-1)/2)*b2_last_out(state2,1))*0.5*tanh(1+tanh(rec_s(2*k)/2)*b2_last_out(state2,2)))*L_a(2*k);% -log(1+exp(L_a(2*k)));
     gamma11 = (0.5*(1+tanh(rec_s(2*k-1)/2)*b2_last_out(state2,3))*0.5*tanh(1+tanh(rec_s(2*k)/2)*b2_last_out(state2,4)))*L_a(2*k);% +L_a(2*k)-log(1+exp(L_a(2*k)));
     
     temp00(state2)= gamma00 * Alpha(k,b2_last_state(state2,1)) * Beta(k,state2);
     temp01(state2)= gamma11 * Alpha(k,b2_last_state(state2,2)) * Beta(k,state2);
  end
  L_all(2*k-1) = log(sum(temp1)) - log(sum(temp0));
  L_all(2*k) = log(sum(temp01)) - log(sum(temp00));
end