siso_app.m

通信中常用的卷积码信道译码源码程序

function [LuO,LcO] = siso_app(LuI, LcI, Trellis, BitsLut, dec, LcOFlag);

%------------------------------------------------------------------------------
% Maximum a posteriori decoder with soft-inputs and soft-outputs. 
%
% Calculates Log-likelihood ratios of received bits. A priori information 
% of information and coded bits from the other decoder and soft democulator 
% accepted. Soft inputs, soft outputs. Log-domain calculations for numerical 
% stability. No approximations used in the calculation of logarithm
% of sum of exponents.
%
% Note: PCCC-schemes shouldn't use a priori LLRs: log(P(c;O)). See references.
% 
% Algorithm reference:
%   Benedetto, S., Divsalar, D., Montorsi, G. and Pollara, F.,
%   "A Soft-Input Soft-Output APP Module for Iterative Decoding of 
%    Concatenated Codes". IEEE Communications Letters, Vol.1, No.1, 1997.
%   Benedetto, S., Divsalar, D., Montorsi, G. and Pollara, F.,
%   "A Soft-Input Soft-Output Maximum A Posterioiri (MAP) Module to Decode
%    Parallel and Serial Concatenated Codes". TDA Progress Report 42-127.
% 
% Restrictions: Two 1/2-rate RSCs as CCs assumed.
%
% 
% Format:
% -------
% 
% [LuO,LcO] = siso_app(LuI, LcI, Trellis, BitsLut, dec, LcOFlag);
%
% Inputs:
% -------
%  LuI - log(P(u;I)), a priori log-likelihoods of the information bits 
%                     from other component decoder
%  LcI - log(P(c;I)), a priori log-likelihoods of the coded bits
%                     from demodulator (PCCC) or other component decoder(SCCC)
%  Trellis - Structure containing necessary information of CC trellises
%  BitsLut - Bit look-up table for edge output symbols
%  dec - Decoder 1 or 2
%  LcOFlag - Calculate likelihoods log(P(c;O)) or not
%
% Outputs:
% --------
%  LuO - log(P(u;O)), extrinsic bit information of the information bits
%  LcO - log(P(c;O)), extrinsic bit information of the coded bits
% 
% Author: MVe, [e-mail]
% Date:   14.10.2002
%------------------------------------------------------------------------------


% Format of 'edge' (transition connecting two states) for 1/2-rate CCs:
%  row = edge number
%                e_S         e_E          e_U           e_C1, e_C2
%  column = | start edge | end edge | input bit | output binary sequence |
%
[e_S, e_E, e_U, e_C1, e_C2] = deal(1,2,3,4,5);
edge(:,e_S:e_E) = Trellis.Edge(:,1:2);
edge(:,e_U) = Trellis.Edge(:,3)-1;		% binary input
edge(:,e_C1:e_C2) = BitsLut(Trellis.Edge(:,4),:);

NEncoders = 2;				% two CCs (see restrictions)

% Extract most used variables from struct 'Trellis'
number_of_states = Trellis.StNr;
k0 = Trellis.k;				% # of bits / input symbol
n0 = Trellis.n;				% # of bits / output symbol

% Define "infinity" (to keep things coherent with C++)
Infty = 1e100;

% Length of the input stream
L_input = size(LcI,2);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% -- alpha (forward recursion) -- %% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Note:
%
% L_e(u(n)) = log(P(u(n)==1) / P(u(n)==0)) = L(u(n)==1) - L(u(n)==0)
% 
%  ==> L(u(n)==1) = L_e(u(n)==1) - log(1 + exp(L_e(u(n)==1)))
%  ==> L(u(n)==0) = - log(1 + exp(L_e(u(n)==0)))
%       ==> L(u(n)) = u(n)*L_e(u(n)) - log(1 + exp(L_e(u(n))))  | u(n) = {0,1}
%
% L_e(u(n)) - extrinsic information from other component decoder ('siso_app')
%

% Format of 'alpha':
%  row = state
%  column = time
%
alpha = ones(number_of_states,L_input+1)*(-Infty);
alpha(1,1) = 0;

for k = 2:L_input+1			% time / frequency
  for e = 1:number_of_states*2^k0	% loop over edges

      alpha(edge(e,e_E),k) = ...
	  maxx(alpha(edge(e,e_E),k), ...
	       alpha(edge(e,e_S),k-1)+ ...
	       edge(e,e_U)*LuI(k-1)-maxx(0,LuI(k-1))+ ...
	       edge(e,e_C1)*LcI(1,k-1)-maxx(0,LcI(1,k-1))+ ...
	       edge(e,e_C2)*LcI(2,k-1)-maxx(0,LcI(2,k-1)));

  end % for e = 1:number_of_states*2^k0
end % for k = 2:L_input+1



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% -- beta (backward recursion) -- %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Format of 'beta':
%  row = state
%  column = time
%
beta = ones(number_of_states,L_input+1)*(-Infty);

% First encoder is terminated, the others are not.
if dec == 1
  beta(1,L_input+1) = 0;
else
  beta(:,L_input+1) = alpha(:,L_input+1); 
end % if dec == 1


for k = L_input:(-1):1			% time / frequency
  for e = 1:number_of_states*2^k0	% loop over edges

      beta(edge(e,e_S),k) = ...
	  maxx(beta(edge(e,e_S),k), ...
	       beta(edge(e,e_E),k+1)+ ...
	       edge(e,e_U)*LuI(k)-maxx(0,LuI(k))+ ...
	       edge(e,e_C1)*LcI(1,k)-maxx(0,LcI(1,k))+ ...
	       edge(e,e_C2)*LcI(2,k)-maxx(0,LcI(2,k)));

  end % for e = 1:number_of_states*2^k0
end % for k = L_input:-1:1


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% -- LuO(k,l)=log(P(u(k,l);O)) -- %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Note: Works only for 1/2-rate CCs (extensions if needed).
%

LuO_tmp = ones(2,L_input+1)*(-Infty);
LuO = ones(k0,L_input)*(-Infty);

for k = 2:L_input+1
  for e = 1:number_of_states*2^k0

    LuO_tmp(edge(e,e_U)+1,k) = ...
	maxx(LuO_tmp(edge(e,e_U)+1,k), ...
	     alpha(edge(e,e_S),k-1)+ ...
	     edge(e,e_C1)*LcI(1,k-1)-maxx(0,LcI(1,k-1))+ ...
	     edge(e,e_C2)*LcI(2,k-1)-maxx(0,LcI(2,k-1))+ ...
	     beta(edge(e,e_E),k));
    
  end % for e = 1:number_of_states*2^k0
end % for k = 1:L_input

LuO(1,1:L_input) = LuO_tmp(2,2:end)-LuO_tmp(1,2:end);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% -- LcO(k,l)=log(P(c(k,l);O)) -- %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Note: Works only for 1/2-rate CCs  (extensions if needed).
%

LcO_tmp = ones(2,L_input+1,n0)*(-Infty);
LcO = ones(n0,L_input)*(-Infty);

if LcOFlag ~= 0

  for k = 2:L_input+1
    for e = 1:number_of_states*2^k0

      % first output bit
      LcO_tmp(edge(e,e_U)+1,k,1) = ...
	  maxx(LcO_tmp(edge(e,e_U)+1,k,1), ...
	       alpha(edge(e,e_S),k-1)+ ...
	       edge(e,e_U)*LuI(1,k-1)-maxx(0,LuI(1,k-1))+ ...
	       edge(e,e_C2)*LcI(2,k-1)-maxx(0,LcI(2,k-1))+ ...
	       beta(edge(e,e_E),k));
      
      % second output bit
      LcO_tmp(edge(e,e_U)+1,k,2) = ...
	  maxx(LcO_tmp(edge(e,e_U)+1,k,2), ...
	       alpha(edge(e,e_S),k-1)+ ...
	       edge(e,e_U)*LuI(1,k-1)-maxx(0,LuI(1,k-1))+ ...
	       edge(e,e_C1)*LcI(1,k-1)-maxx(0,LcI(1,k-1))+ ...
	       beta(edge(e,e_E),k));
      
    end % for e = 1:number_of_states*2^k0
  end % for k = 1:L_input
  
  for i1 = 1:n0
    LcO(i1,1:L_input) = LcO_tmp(2,2:end,i1)-LcO_tmp(1,2:end,i1);
  end % for i1 = 1:n0
  
end % if LcOFlag ~= 0