function L_all = sova(rec_s, g, L_a, ind_dec)
% This function implememts Soft Output Viterbi Algorithm in trace
back mode
% Input:
% rec_s: scaled received bits. rec_s(k) = 0.5 * L_c(k) * y(k)
% L_c = 4 * a * Es/No, reliability value of the channel
% y: received bits
% g: encoder generator matrix in binary form,
g(1,:) for feedback, g(2,:) for feedforward
% L_a: a priori information about the info. bits. Extrinsic info.
from the previous
% component decoder
% ind_dec: index of the component decoder.
% =1: component decoder 1; The trellis is terminated to all zero
state
% =2: component decoder 2; The trellis is not perfectly
terminated.
% Output:
% L_all: log ( P(x=1|y) ) / ( P(x=-1|y) )
%
% Copyright: Yufei Wu, Nov. 1998
% MPRG lab, Virginia Tech
% for academic use only
% Frame size, info. + tail bits
L_total = length(L_a);
[n,K] = size(g);
m = K - 1;
nstates = 2^m;
Infty = 1e10;
% SOVA window size. Make decision after 'delta' delay. Decide
bit k when received bits
% for bit (k+delta) are processed. Trace back from (k+delta) to
k.
delta =
30;
% Set up the trellis defined by g.
[next_out, next_state, last_out, last_state] = trellis(g);
% Initialize path metrics to -Infty
for t=1:L_total+1
for state=1:nstates
path_metric(state,t) = -Infty;
end
end
% Trace forward to compute all the path metrics
path_metric(1,1) = 0;
for t=1:L_total
y = rec_s(2*t-1:2*t);
for state=1:nstates
sym0 = last_out(state,1:2);
sym1 = last_out(state,3:4);
state0 = last_state(state,1);
state1 = last_state(state,2);
Mk0 = y*sym0' - L_a(t)/2 + path_metric(state0,t);
Mk1 = y*sym1' + L_a(t)/2 + path_metric(state1,t);
if Mk0>Mk1
path_metric(state,t+1)=Mk0;
Mdiff(state,t+1) = Mk0 - Mk1;
prev_bit(state, t+1) = 0;
else
path_metric(state,t+1)=Mk1;
Mdiff(state,t+1) = Mk1 - Mk0;
prev_bit(state,t+1) = 1;
end
end
end
% For decoder 1, trace back from all zero state,
% for decoder two, trace back from the most likely state
if ind_dec == 1
mlstate(L_total+1) = 1;
else
mlstate(L_total+1) = find(
path_metric(:,L_total+1)==max(path_metric(:,L_total+1)) );
end
% Trace back to get the estimated bits, and the most likely
path
for t=L_total:-1:1
est(t) =
prev_bit(mlstate(t+1),t+1);
mlstate(t) =
last_state(mlstate(t+1), est(t)+1);
end
% Find the minimum delta that corresponds to a compitition
path with different info. bit
estimation. % Give the soft output
for t=1:L_total
llr = Infty;
for i=0:delta
if t+i
bit = 1-est(t+i);
temp_state = last_state(mlstate(t+i+1), bit+1);
for j=i-1:-1:0
bit = prev_bit(temp_state,t+j+1);
temp_state = last_state(temp_state, bit+1);
end
if bit~=est(t)
llr = min( llr,Mdiff(mlstate(t+i+1), t+i+1) );
end
end
end
L_all(t) = (2*est(t) - 1) *
llr;
end ----------------------------------------------------------
function [next_out, next_state, last_out, last_state] =
trellis(g)
% copyright Nov. 1998 Yufei Wu
% MPRG lab, Virginia Tech
% for academic use only
% set up the trellis given code generator g
% g given in binary matrix form. e.g. g = [ 1 1 1; 1 0 1 ];
% next_out(i,1:2): trellis next_out (systematic bit; parity
bit) when input = 0, state = i; next_out(i,j) = -1 or 1
% next_out(i,3:4): trellis next_out (systematic
bit; parity bit) when input = 1, state = i;
% next_state(i,1): next state when input = 0, state = i;
next_state(i,i) = 1,...2^m
% next_state(i,2): next state when input = 1, state = i;
% last_out(i,1:2): trellis last_out (systematic bit; parity bit)
when input = 0, state = i; last_out(i,j) = -1 or 1
% last_out(i,3:4): trellis last_out (systematic
bit; parity bit) when input = 1, state = i;
% last_state(i,1): previous state that comes to state i when info.
bit = 0;
% last_state(i,2): previous state that comes to state i when info.
bit = 1;
[n,K] = size(g);
m = K - 1;
max_state = 2^m;
% set up next_out and next_state matrices for systematic
code
for state=1:max_state
state_vector = bin_state(
state-1, m );
% when receive a 0
d_k = 0;
a_k = rem( g(1,:)*[0
state_vector]', 2 );
[out_0, state_0] =
encode_bit(g, a_k, state_vector);
out_0(1) = 0;
% when receive a 1
d_k = 1;
a_k = rem( g(1,:)*[1
state_vector]', 2 );
[out_1, state_1] =
encode_bit(g, a_k, state_vector);
out_1(1) = 1;
next_out(state,:) = 2*[out_0
out_1]-1;
next_state(state,:) =
[(int_state(state_0)+1) (int_state(state_1)+1)];
end
% find out which two previous states can come to present
state
last_state = zeros(max_state,2);
for bit=0:1
for state=1:max_state
last_state(next_state(state,bit+1), bit+1)=state;
last_out(next_state(state, bit+1), bit*2+1:bit*2+2) ...
= next_out(state, bit*2+1:bit*2+2);
end
end
-------------------------------------------------------------
% This script simulates the classical turbo encoding-decoding
system.
% It simulates parallel concatenated convolutional codes.
% Two component rate 1/2 RSC (Recursive Systematic Convolutional)
component encoders are assumed.
% First encoder is terminated with tails bits. (Info + tail) bits
are scrambled and passed to
% the second encoder, while second encoder is left open without
tail bits of itself.
%
% Random information bits are modulated into +1/-1, and transmitted
through a AWGN channel.
% Interleavers are randomly generated for each frame.
%
% Log-MAP algorithm without quantization or approximation is
used.
% By making use of ln(e^x+e^y) = max(x,y) +
ln(1+e^(-abs(x-y))),
% the Log-MAP can be simplified with a look-up table for the
correction function.
% If use approximation ln(e^x+e^y) = max(x,y), it becomes
MAX-Log-MAP.
%
% Copyright Nov 1998, Yufei Wu
% MPRG lab, Virginia Tech.
% for academic use only
clear all
% Write display messages to a text file
diary turbo_logmap.txt
% Choose decoding algorithm
dec_alg = input(' Please enter the decoding algorithm. (0:Log-MAP,
1:SOVA) default
0 ');
if isempty(dec_alg)
dec_alg = 0;
end
% Frame size
L_total = input(' Please enter the frame size (= info + tail,
default: 400) ');
if isempty(L_total)
L_total =
400; % infomation bits plus tail bits
end
% Code generator
g = input(' Please enter code generator: ( default: g = [1 1 1; 1 0
1 ]
) ');
if isempty(g)
g = [ 1 1 1;
1 0 1 ];
end
%g = [1 1 0 1; 1 1 1 1];
%g = [1 1 1 1 1; 1 0 0 0 1];
[n,K] = size(g);
m = K - 1;
nstates = 2^m;
%puncture = 0, puncturing into rate 1/2;
%puncture = 1, no puncturing
puncture = input(' Please choose punctured / unpunctured (0/1):
default
0 ');
if isempty(puncture)
puncture =
0;
end
% Code rate
rate = 1/(2+puncture);
% Fading amplitude; a=1 in AWGN channel
a = 1;
% Number of iterations
niter = input(' Please enter number of iterations for each frame:
default
5 ');
if isempty(niter)
niter = 5;
end % Number of frame errors to count as a stop criterior
ferrlim = input(' Please enter number of frame errors to terminate:
default
15 ');
if isempty(ferrlim)
ferrlim = 15;
end
EbN0db = input(' Please enter Eb/N0 in dB : default
[2.0] ');
if isempty(EbN0db)
EbN0db = [2.0];
end
fprintf('\n\n----------------------------------------------------\n');
if dec_alg == 0
fprintf(' === Log-MAP decoder
=== \n');
else
fprintf(' === SOVA decoder ===
\n');
end
fprintf(' Frame size = %6d\n',L_total);
fprintf(' code generator: \n');
for i = 1:n
for j =
1:K
fprintf( '%6d', g(i,j));
end
fprintf('\n');
end if puncture==0
fprintf(' Punctured, code rate
= 1/2 \n');
else
fprintf(' Unpunctured, code
rate = 1/3 \n');
end
fprintf(' iteration number = %6d\n',
niter);
fprintf(' terminate frame errors = %6d\n', ferrlim);
fprintf(' Eb / N0 (dB) = ');
for i = 1:length(EbN0db)
fprintf('%10.2f',EbN0db(i));
end
fprintf('\n----------------------------------------------------\n\n');
fprintf('+ + + + Please be patient. Wait a while to get the result.
+ + + +\n');
for nEN = 1:length(EbN0db)
en =
10^(EbN0db(nEN)/10); % convert Eb/N0 from unit db to normal numbers
L_c = 4*a*en*rate;
% reliability value of the channel
sigma = 1/sqrt(2*rate*en);
% standard deviation of AWGN noise
% Clear bit error counter and frame error counter
errs(nEN,1:niter) =
zeros(1,niter);
nferr(nEN,1:niter) =
zeros(1,niter);
nframe =
0; % clear
counter of transmitted frames
while nferr(nEN,
niter)
nframe = nframe +
1; x = round(rand(1,
L_total-m)); % info. bits
[temp, alpha] =
sort(rand(1,L_total)); % random interleaver mapping
en_output = encoderm( x, g, alpha, puncture ) ; % encoder output
(+1/-1)
r = en_output+sigma*randn(1,L_total*(2+puncture)); % received
bits
yk = demultiplex(r,alpha,puncture); % demultiplex to get input for
decoder 1 and 2
% Scale the received
bits rec_s = 0.5*L_c*yk;
% Initialize extrinsic
information L_e(1:L_total) = zeros(1,L_total);
for iter = 1:niter
% Decoder one
L_a(alpha) = L_e; % a priori info.
if dec_alg == 0
L_all = logmapo(rec_s(1,:), g, L_a, 1); %
complete info.
else L_all = sova0(rec_s(1,:), g, L_a, 1); % complete
info.
end L_e = L_all - 2*rec_s(1,1:2:2*L_total) - L_a; %
extrinsic info.
% Decoder
two L_a = L_e(alpha); % a priori info.
if dec_alg == 0
L_all = logmapo(rec_s(2,:), g, L_a, 2); %
complete info. else
L_all = sova0(rec_s(2,:), g, L_a, 2); % complete
info.
end
L_e = L_all - 2*rec_s(2,1:2:2*L_total) - L_a; %
extrinsic info.
% Estimate the info.
bits xhat(alpha) = (sign(L_all)+1)/2;
% Number of bit errors in current iteration
err(iter) = length(find(xhat(1:L_total-m)~=x));
% Count frame errors for the current iteration
if err(iter)>0
nferr(nEN,iter) = nferr(nEN,iter)+1;
end end %iter
% Total number of bit errors for all iterations
errs(nEN,1:niter) = errs(nEN,1:niter) + err(1:niter);
if rem(nframe,3)==0 | nferr(nEN, niter)==ferrlim
% Bit error rate
ber(nEN,1:niter) = errs(nEN,1:niter)/nframe/(L_total-m);
% Frame error rate
fer(nEN,1:niter) = nferr(nEN,1:niter)/nframe;
% Display intermediate results in
process fprintf('************** Eb/N0 = %5.2f db **************\n',
EbN0db(nEN));
fprintf('Frame size = %d, rate 1/%d. \n', L_total,
2+puncture);
fprintf('%d frames transmitted, %d frames in error.\n', nframe,
nferr(nEN, niter));
fprintf('Bit Error Rate (from iteration 1 to iteration %d):\n',
niter);
for i=1:niter
fprintf('%8.4e ', ber(nEN,i));
end
fprintf('\n');
fprintf('Frame Error Rate (from iteration 1 to iteration %d):\n',
niter);
for i=1:niter
fprintf('%8.4e ', fer(nEN,i));
end
fprintf('\n');
fprintf('***********************************************\n\n');
% Save intermediate results
save turbo_sys_demo EbN0db ber fer
end
end %while
end %nEN
diary off
-----------------------------------------------------