turbo_rec.m

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

function [out_bits,berrs,ferrs,raw_berrs,raw_ferrs,varargout] = ...
    turbo_rec(Data_received, fading_coeff, N0_int, alpha_channel, alpha_enc,...
	      alpha_rsc, iterations, Trellis, Lut, ModMap, BitsLut, ...
	      int_input, PunctureF, LcOFlag, MaxFlag);

%------------------------------------------------------------------------------
% Receiver with soft demodulation and turbo-decoding.
% Demodulate and decode binary 1-D data vector with iterative MAP-algorithm.
%
% $RCSfile: turbo_rec.m,v $
% $Revision: 2.11 $
% $Date: 2004/07/29 14:33:15 $
% $Author: coffin $
%
% Format:
% -------
% [out_bits, berrs, ferrs, raw_berrs, [LcO1], [LcO2]] = ...
%   turbo_rec(Data_received, fading_coeff, N0, alpha_channel, alpha_enc, ...
%             alpha_rsc, iterations, Trellis, Lut, ModMap, BitsLut, ...
%             int_input, PunctureF, LcOFlag);
%
% Inputs:
% -------
% Data_received - Received data vector
% fading_coeff  - Channel path gains
% N0_int        - Total noise variance
% alpha_channel - Symbol-level interleaver after modulation
% alpha_enc     - Bit-level interleaver after turbo encoder
% alpha_rsc     - Bit-level interleaver between RSCs
% iterations     - Number of maximum decoding iterations 
% Trellis       - Structure holding most of the parameters related to code
% Lut           - Some look-up tables
% BitsLut       - Bit sequences related to constellation points
% ModMap        - Complex constellation points
% int_input     - Original information bits + tail for error calculation
% PunctureF     - Flag 0/1, whether code is punctured or not
% LcOFlag       - Flag 0/1, calculate soft estimates of coded bits
% 
% Outputs:
% --------
% out_bits      - Decoded information bits
% berrs         - Bit errors found
% ferrs         - Frame errors found
% raw_berrs     - Raw bit errors found
%
% Author: MVe, [e-mail]
% Date:   09.10.2002
%------------------------------------------------------------------------------

Infty = 1e100;				% for coherence with C++ implementation
Md = length(ModMap);			% # of constellation points
Nbits = log2(Md);			% # of bits / constellation points
L_Received = length(Data_received);

% De-interleaving
r(1,alpha_channel) = Data_received;
a(1,alpha_channel) = fading_coeff(1,:);	% fading coefficients
N0(:,alpha_channel) = N0_int;

%% -- Soft demodulator -- %%
% 
% Reference (for example):
%  Stefanov, A. and Duman, T. M., "Turbo Coded Modulation for
%  Systems with Transmission and Receive Antenna Diversity...".
%

%% Calculate log-likelihoods for bits b(l), l = 1,...,Nbits*L_data
%% conditioned to corresponding received symbol r(k), k = 1,...,L_data:
%% LL(b(l)|r(k)) = log(P(b(l)=1|r(k))) - log(P(b(l)=1|r(k)))
%
LL_tmp1 = soft_demodulator_vect_CC(r, a, ModMap, BitsLut, N0);

% $$$ LL_tmp1 = zeros(Nbits,L_Received);
% $$$ LL_tmp2 = zeros(1,Nbits*L_Received);
% $$$ 
% $$$ for k = 1:L_Received		% loop over received symbols
% $$$   
% $$$   for l = 1:Nbits		% loop over bits within one
% $$$                                 % received symbol
% $$$     
% $$$     sum_1 = -Infty;
% $$$     sum_0 = -Infty;
% $$$     
% $$$     for M = 1:Md		% loop over constellation points
% $$$       
% $$$       b_vect = BitsLut(M,:);	% c = f(b_vect), f() is modulation
% $$$       c = ModMap(M);		% c - constellation point
% $$$ 
% $$$       % Scaled Euclidian distance
% $$$       log_metric = -(abs( r(1,k) - a(1,k)*c )^2)/N0(k); 
% $$$ 
% $$$       % Symbol-level to bit-level a priori log-likelihood extraction
% $$$       if b_vect(l)==1
% $$$ 	sum_1 = maxx(sum_1,log_metric);	% maxx(x,y) = log(exp(x)+exp(y))
% $$$       else
% $$$ 	sum_0 = maxx(sum_0,log_metric);
% $$$       end % if b_vect(l)==1	
% $$$       
% $$$     end % for M = Md
% $$$ 
% $$$     % Log-likelihood ratio for the lth bit in the kth symbol
% $$$     LL_tmp1(l,k) = sum_1 - sum_0;
% $$$     
% $$$   end % for l = 1:Nbits
% $$$ 
% $$$ end % for k = 1:L_Received


% Total number of transmitted bits
L_tot = Nbits*L_Received;

LL_tmp2 = reshape(LL_tmp1,1,L_tot);

% Check if modulator has appended all-zeros tail
if Trellis.flag_3gpp
  if isfield(Trellis,'Nzeropad')
    LL_tmp2 = LL_tmp2(1,1:end-Trellis.Nzeropad);
    L_tot = L_tot-Trellis.Nzeropad;
  end % if isfield(Trellis,'Nzeropad')
else
  tmp = mod(length(int_input)*(3-PunctureF),Nbits);
  if tmp>0
    Ltail = Nbits-tmp;
    LL_tmp2 = LL_tmp2(1,1:end-Ltail);
    L_tot = L_tot-Ltail;
  end % if mod(L_total,k)
end % if Trellis.flag_3gpp

% De-interleave information and parity bits
LL = zeros(1,length(alpha_enc));
LL(1,alpha_enc) = LL_tmp2(1,1:end);

% Separate tail bits from coded sequence
if Trellis.flag_3gpp
  NEncoders = 2;
  LL_tail = reshape(LL(1,end-NEncoders*Trellis.Ntailbits_3gpp+1:end),...
		    Trellis.Ntailbits_3gpp,NEncoders).';
  LL = LL(1,1:end-NEncoders*Trellis.Ntailbits_3gpp);
  L_tot = L_tot-NEncoders*Trellis.Ntailbits_3gpp;
end % if Trellis.flag_3gpp


%% -- Turbo decoder (two component SISO-MAP algorithms) -- %%
% 
% Restrictions: Two 1/2-rate RSCs as CCs assumed.
%
% Reference: 
%   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.
%

% Separate a priori information to information and coded bit a
% priori LLRs for different decoders.
%
% Note. For AWGN channel P_k(c;I) = P_k(y(k)|x(c))*P(c)/P(y),
%       where x(c) is transmitted symbol. 

LL_I1 = zeros(2,L_tot/(2-PunctureF));	% bit reliabilities for 1st CC
LL_I2 = zeros(2,L_tot/(2-PunctureF));	% bit reliabilities for 2nd CC

if PunctureF
  error('DOES NOT WORK WITH PUNCTURING!!!!')
  LL_I1(1,1:end) = LL(1,1:2:end);	% LLRs of inf.bits for 1st MAP
else
  LL_I1(1,1:end) = LL(1,1:2:end);
end % if PunctureF

% $$$ LL_I2(1,1:end) = LL_I1(1,alpha_rsc);	% LLRs of inf.bits for 2nd MAP

if PunctureF
  error('DOES NOT WORK WITH PUNCTURING!!!!')
  LL_I1(2,1:2:end) = LL(1,2:4:end);	% LLRs of code bits for 1st MAP
  LL_I2(2,2:2:end) = LL(1,4:4:end);	% LLRs of code bits for 2ns MAP
else
  LL_I1(2,1:end) = LL(1,2:2:end);
% $$$   LL_I2(2,1:end) = LL(1,3:3:end);
end % if PunctureF



% Append tail bits in 3GPP compliant mode
if Trellis.flag_3gpp
  LL_I1 = [LL_I1,reshape(LL_tail(1,:),2,Trellis.Ntailbits_3gpp/2)];
  LL_I2 = [LL_I2,reshape(LL_tail(2,:),2,Trellis.Ntailbits_3gpp/2)];
%  L_tot = L_tot + NEncoders*Trellis.Ntailbits_3gpp;
end % if Trellis.flag_3gpp

% Prepare variables for output LLRs (information and coded bits)
LuO1 = zeros(1,L_tot/(2-PunctureF));
LcO1 = zeros(2,L_tot/(2-PunctureF));
% $$$ LuO2 = zeros(1,L_tot/(2-PunctureF));
% $$$ LcO2 = zeros(2,L_tot/(2-PunctureF));

berrs = zeros(1,iterations);
ferrs = zeros(1,iterations);

LuI1 = zeros(1,L_tot/(2-PunctureF));
% $$$ LuI2 = zeros(1,L_tot/(2-PunctureF));

%%%%%%%%%%%%%%%%%%%%%%
%% -- Iterations -- %%
%%%%%%%%%%%%%%%%%%%%%%

% $$$ LuI1_2 = zeros(1,L_tot/(3-PunctureF));
% $$$ LuI2_2 = zeros(1,L_tot/(3-PunctureF));


for i1 = 1:iterations
  
  %% 1st component decoder %%
  %
  CDecoder = 1;
  [LuO1,LcO1] = ...
      siso_app_C(LuI1, LL_I1, Trellis, Lut.Outputs_bin, CDecoder, LcOFlag,...
		 MaxFlag);

% $$$   [LuO1_2,LcO1_2] = ...
% $$$       siso_app_norm_C(LuI1_2, LL_I1, Trellis, Lut.Outputs_bin, CDecoder, ...
% $$$ 		      LcOFlag, MaxFlag);

  
% $$$   % Extrinsic information for 2nd MAP decoder.
% $$$   if Trellis.flag_3gpp
% $$$     % Puncture tail bits away
% $$$     LuO1 = LuO1(1,1:end-Trellis.Ntailbits_3gpp/2);
% $$$   end % if Trellis.flag_3gpp
% $$$   LuI2(1,:) = LuO1(1,alpha_rsc);  
  
% $$$   LuI2_2(1,:) = LuO1_2(1,alpha_rsc);  

  %% 2nd component decoder %%
  %
% $$$   CDecoder = 2;
% $$$   [LuO2,LcO2] = ...
% $$$       siso_app_C(LuI2, LL_I2, Trellis, Lut.Outputs_bin, CDecoder, LcOFlag,...
% $$$ 		 MaxFlag);

% $$$   [LuO2_2,LcO2_2] = ...
% $$$       siso_app_norm_C(LuI2_2, LL_I2, Trellis, Lut.Outputs_bin, CDecoder, ...
% $$$ 		      LcOFlag, MaxFlag);


  LuO1_tot(i1,:) = abs(LuO1);  
  LcO1_1_tot(i1,:) = abs(LcO1(1,:));
  LcO1_2_tot(i1,:) = abs(LcO1(2,:));  
% $$$   LuO2_tot(i1,:) = abs(LuO2);
% $$$   LcO2_1_tot(i1,:) = abs(LcO2(1,:));
% $$$   LcO2_2_tot(i1,:) = abs(LcO2(2,:));

  
% $$$   % Extrinsic information for 1st MAP decoder.
% $$$   if Trellis.flag_3gpp
% $$$     % Puncture tail bits away
% $$$     LuO2 = LuO2(1,1:end-Trellis.Ntailbits_3gpp/2);
% $$$   end % if Trellis.flag_3gpp
% $$$   LuI1(1,alpha_rsc)=LuO2(1,:);

% $$$   LuI1_2(1,alpha_rsc)=LuO2_2(1,:);  

  % Full LLR information of the information bits.
  LL_O1 = LuO1+LuI1;

  decisions(1,:) = (LL_O1(1,1:end-Trellis.m)>0);	% hard decisions

  % Number of bit errors in current iteration
  berrs(i1) = length(find(decisions~=int_input));

  % Count frame errors for the current iteration
  if berrs(i1)>0
    ferrs(i1) = ferrs(i1)+1;
  else
    berrs(i1:iterations)=0;    
%    break;				% break iterations for-loop
  end % if berrs(j)>0

end % for i1 = 1:iterations

% Output LLRs of coded bits
if LcOFlag
  varargout{1} = LcO1;
% $$$   varargout{2} = LcO2;
end % if LcOFlag

if Trellis.flag_3gpp
  raw_berrs = length(find((LL_I1(1,1:end-Trellis.Ntailbits_3gpp/2)>0)~=...
			  int_input));  
else
  raw_berrs = length(find((LL_I1(1,1:end-Trellis.m)>0)~=int_input));
  if raw_berrs>0
    raw_ferrs = 1;
  else
    raw_ferrs = 0;
  end % if raw_berrs>0 
end % if Trellis.flag_3gpp

out_bits = decisions;