1+ % BCJR algorithm for a half-rate systematic recursive convolutional code
2+ % having 1 memory element, a generator polynomial of [1,0] and a feedback
3+ % polynomial of [1,1]. For more information, see Section 1.3.2.2 of Rob's
4+ % thesis (http://eprints.ecs.soton.ac.uk/14980) or the BCJR paper
5+ % (http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=1055186).
6+ % Copyright (C) 2008 Robert G. Maunder
7+
8+ % This program is free software: you can redistribute it and/or modify it
9+ % under the terms of the GNU General Public License as published by the
10+ % Free Software Foundation, either version 3 of the License, or (at your
11+ % option) any later version.
12+
13+ % This program is distributed in the hope that it will be useful, but
14+ % WITHOUT ANY WARRANTY; without even the implied warranty of
15+ % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
16+ % Public License for more details.
17+
18+ % The GNU General Public License can be seen at http://www.gnu.org/licenses/.
19+
20+
21+
22+ function [aposteriori_uncoded_llrs , aposteriori_encoded1_llrs , aposteriori_encoded2_llrs ] = bcjr_decoder(apriori_uncoded_llrs , apriori_encoded1_llrs , apriori_encoded2_llrs )
23+
24+ if (length(apriori_uncoded_llrs ) ~= length(apriori_encoded1_llrs ) || length(apriori_encoded1_llrs ) ~= length(apriori_encoded2_llrs ))
25+ error(' LLR sequences must have the same length' );
26+ end
27+
28+
29+ % All calculations are performed in the logarithmic domain in order to
30+ % avoid numerical issues. These occur in the normal domain, because some of
31+ % the confidences can get smaller than the smallest number the computer can
32+ % store. See Section 1.3.2.4 of Rob's thesis for more information on this.
33+ %
34+ % A multiplication of two confidences is achieved using the addition of the
35+ % corresponding log-confidences. If A = log(a) and B = log(b), then
36+ % log(a*b) = A+B (Equation 1.17 in Rob's thesis).
37+ %
38+ % An addition of two confidences is achieved using the Jacobian logarithm
39+ % of the corresponding log-confidences. The Jacobian logarithm is defined
40+ % in the jac.m file. If A = log(a) and B = log(b), then
41+ % log(a+b) = max(A,B) + log(1+exp(-abs(A-B))) (Equation 1.19 in Rob's
42+ % thesis).
43+
44+ % Matrix to describe the trellis
45+ % Each row describes one transition in the trellis
46+ % Each state is allocated an index 1,2,3,... Note that this list starts
47+ % from 1 rather than 0.
48+ % FromState, ToState, UncodedBit, Encoded1Bit, Encoded2Bit
49+ transitions = [1 , 1 , 0 , 0 , 0 ;
50+ 1 , 2 , 1 , 1 , 1 ;
51+ 2 , 1 , 1 , 1 , 0 ;
52+ 2 , 2 , 0 , 0 , 1 ];
53+
54+ % Find the largest state index in the transitions matrix
55+ % In this example, we have two states since the code has one memory element
56+ state_count = max(max(transitions(: ,1 )),max(transitions(: ,2 )));
57+
58+ % Calculate the a priori transition log-confidences by adding the
59+ % log-confidences associated with each corresponding bit value. This is
60+ % similar to Equation 1.12 in Rob's thesis or Equation 9 in the BCJR paper.
61+ gammas= zeros(size(transitions ,1 ),length(apriori_uncoded_llrs ));
62+ for bit_index = 1 : length(apriori_uncoded_llrs )
63+ for transition_index = 1 : size(transitions ,1 )
64+ if transitions(transition_index , 3 )==0
65+ gammas(transition_index , bit_index ) = gammas(transition_index , bit_index ) + apriori_uncoded_llrs(bit_index )/2 ;
66+ % Dividing the LLR by 2 gives a value that is log-proportional to
67+ % the actual log-confidence it instills. Log-proportional is fine
68+ % for us though, because we're after the log-ratio of
69+ % confidences. Don't worry too much about this confusing
70+ % feature.
71+ else
72+ gammas(transition_index , bit_index ) = gammas(transition_index , bit_index ) - apriori_uncoded_llrs(bit_index )/2 ;
73+ end
74+
75+ if transitions(transition_index , 4 )==0
76+ gammas(transition_index , bit_index ) = gammas(transition_index , bit_index ) + apriori_encoded1_llrs(bit_index )/2 ;
77+ else
78+ gammas(transition_index , bit_index ) = gammas(transition_index , bit_index ) - apriori_encoded1_llrs(bit_index )/2 ;
79+ end
80+
81+ if transitions(transition_index , 5 )==0
82+ gammas(transition_index , bit_index ) = gammas(transition_index , bit_index ) + apriori_encoded2_llrs(bit_index )/2 ;
83+ else
84+ gammas(transition_index , bit_index ) = gammas(transition_index , bit_index ) - apriori_encoded2_llrs(bit_index )/2 ;
85+ end
86+ end
87+ end
88+
89+ % Forward recursion to calculate state log-confidences. This is similar to
90+ % Equation 1.13 in Rob's thesis or Equations 5 and 6 in the BCJR paper.
91+ alphas= zeros(state_count ,length(apriori_uncoded_llrs ));
92+ alphas= alphas - inf ;
93+ alphas(1 ,1 )=0 ; % We know that this is the first state
94+ for state_index = 2 : state_count
95+ alphas(state_index ,1 )=-inf ; % We know that this is *not* the first state (a log-confidence of minus infinity is equivalent to a confidence of 0)
96+ end
97+ for bit_index = 2 : length(apriori_uncoded_llrs )
98+ for transition_index = 1 : size(transitions ,1 )
99+ alphas(transitions(transition_index ,2 ),bit_index ) = jac(alphas(transitions(transition_index ,2 ),bit_index ),alphas(transitions(transition_index ,1 ),bit_index - 1 ) + gammas(transition_index , bit_index - 1 ));
100+ end
101+ end
102+
103+ % Backwards recursion to calculate state log-confidences. This is similar
104+ % to Equation 1.14 in Rob's thesis or Equations 7 and 8 in the BCJR paper.
105+ betas= zeros(state_count ,length(apriori_uncoded_llrs ));
106+ betas= betas - inf ;
107+ for state_index = 1 : state_count
108+ betas(state_index ,length(apriori_uncoded_llrs ))=0 ; % The final state could be any one of these
109+ end
110+ for bit_index = length(apriori_uncoded_llrs )-1 : -1 : 1
111+ for transition_index = 1 : size(transitions ,1 )
112+ betas(transitions(transition_index ,1 ),bit_index ) = jac(betas(transitions(transition_index ,1 ),bit_index ),betas(transitions(transition_index ,2 ),bit_index + 1 ) + gammas(transition_index , bit_index + 1 ));
113+ end
114+ end
115+
116+ % Calculate a posteriori transition log-confidences. This is similar to
117+ % Equation 1.15 in Rob's thesis or Equation 4 in the BCJR paper.
118+ deltas= zeros(size(transitions ,1 ),length(apriori_uncoded_llrs ));
119+ for bit_index = 1 : length(apriori_uncoded_llrs )
120+ for transition_index = 1 : size(transitions ,1 )
121+ deltas(transition_index , bit_index ) = alphas(transitions(transition_index ,1 ),bit_index ) + gammas(transition_index , bit_index ) + betas(transitions(transition_index ,2 ),bit_index );
122+ end
123+ end
124+
125+ % Calculate the a posteriori LLRs. This is similar to Equation 1.16 in
126+ % Rob's thesis.
127+ aposteriori_uncoded_llrs = zeros(1 ,length(apriori_uncoded_llrs ));
128+ for bit_index = 1 : length(apriori_uncoded_llrs )
129+ prob0= -inf ;
130+ prob1= -inf ;
131+ for transition_index = 1 : size(transitions ,1 )
132+ if transitions(transition_index ,3 )==0
133+ prob0 = jac(prob0 , deltas(transition_index ,bit_index ));
134+ else
135+ prob1 = jac(prob1 , deltas(transition_index ,bit_index ));
136+ end
137+ end
138+ aposteriori_uncoded_llrs(bit_index ) = prob0 - prob1 ;
139+ end
140+
141+ aposteriori_encoded1_llrs = zeros(1 ,length(apriori_uncoded_llrs ));
142+ for bit_index = 1 : length(apriori_uncoded_llrs )
143+ prob0= -inf ;
144+ prob1= -inf ;
145+ for transition_index = 1 : size(transitions ,1 )
146+ if transitions(transition_index ,4 )==0
147+ prob0 = jac(prob0 , deltas(transition_index ,bit_index ));
148+ else
149+ prob1 = jac(prob1 , deltas(transition_index ,bit_index ));
150+ end
151+ end
152+ aposteriori_encoded1_llrs(bit_index ) = prob0 - prob1 ;
153+ end
154+
155+ aposteriori_encoded2_llrs = zeros(1 ,length(apriori_uncoded_llrs ));
156+ for bit_index = 1 : length(apriori_uncoded_llrs )
157+ prob0= -inf ;
158+ prob1= -inf ;
159+ for transition_index = 1 : size(transitions ,1 )
160+ if transitions(transition_index ,5 )==0
161+ prob0 = jac(prob0 , deltas(transition_index ,bit_index ));
162+ else
163+ prob1 = jac(prob1 , deltas(transition_index ,bit_index ));
164+ end
165+ end
166+ aposteriori_encoded2_llrs(bit_index ) = prob0 - prob1 ;
167+ end
168+
169+ end
0 commit comments