trl = poly2trellis( 3, [7 5], 7 ); win = 5; len = 100; inp = round(rand( 1, len ) ); enc = convenc( inp, trl ); msg = 2*enc - 1; msg = msg + 0.1*randn( 1, length(msg) ); llr = zeros( 1, len ); out = sovadec( msg, llr, trl, win ); hrd = vitdec( (msg>0), trl, win, 'trunc', 'hard' ); disp( [ inp', out', hrd' ] );
function out = sovadec( msg, llr, trl, win )% SOVADEC is an implementation of the soft input soft output Viterbi % algorithm. The algorithm can be called using% DEC = SOVADEC( MSG, LLR, TRELLIS, WIN )% where MSG is the soft input (codeword), LLR is a priori information% per bit about the bits (log likelihood ratios), TRELLIS is the % trellis structure describing the convolutional encoder used to % encode the original information message.% % The output of the function is the vector containing the soft % estimates of the originally encoded information. The implementation% is able to perform decoding using any trellis structure compatibile% with the standard Matlab POLY2TRELLIS function.%% WIN describes the size of the trellis buffer used to perform the% Viterbi algorithm. Thus, the decoder estimates the best path through% the trellis and searches back within this buffer at every decoding% time instant.% If WIN is omitted, the trellis of length N+1 (where N is the size% of the decoded message) is used.%% For the estimation of the reliability information only the second% best path (competitor) is used, even if there are more than two% paths merging at a particular state.%% The output of the decoding algorithm is of the following form% out = sign(inp) * log( P(inp=1|out) ) / log( P(inp=-1|out) )%% See also: POLY2TRELLIS, CONVENCO, VITDEC%% (c) 2003 Adrian Bohdanowicz%% $Id: sovadec.m,v 1.16 2003/06/15 14:57:29 adrian Exp $ enc = trellis2enc( trl ); if( enc.k == 1 ) out = sovadec_1N( msg, llr, trl, win ); % use 1/N optimized code (2x faster) else out = sovadec_KN( msg, llr, trl, win ); % uset K/N generic code endreturn; function enc = trellis2enc( trl ),% put the trellis structure into a more user friendly manner enc.k = log2( trl.numInputSymbols ); % number of inputs enc.n = log2( trl.numOutputSymbols ); % numbor of outputs enc.r = enc.k / enc.n; % code rate enc.ksym = trl.numInputSymbols; % number of possible input combinations enc.nsym = trl.numOutputSymbols; % number of possible output combinations enc.stat = trl.numStates; % number of encoder states % forward transitions: enc.next.states = trl.nextStates + 1; % NEXT states enc.next.output = trl.outputs; % NEXT outputs for i = 1:enc.ksym, % NEXT (binary) outputs enc.next.binout( :,:,i ) = 2*de2bi( oct2dec( trl.outputs(:,i) ), enc.n )-1; end % store possible binary outputs and inputs: enc.inp = de2bi( oct2dec( [0:enc.ksym-1] ), enc.k, 'left-msb' ); % all possible binary inputs enc.out = de2bi( oct2dec( [0:enc.nsym-1] ), enc.n, 'left-msb' ); % all possible binary outputs enc.bininp = 2*enc.inp-1;return;function out = sovadec_1N( msg, llr, trl, win ),% SOVADEC optimized for 1/N encoders (faster!) % error checking: if( ~istrellis( trl ) ), error( 'Incorrect input trellis!' ); end; if( nargin <= 2 ), error( 'Incorrect number of input args!' ); end; if( nargin == 3 ), win = length( llr )+1; end; % some parameters: INF = 9e9; % infinity enc = trellis2enc( trl ); % encoder parameters len = length( llr ) / enc.k; % number of decoding steps (total) win = min( [ win, len ] ); % trim the buffer if msg is short old = NaN; % to remember the last survivor % allocate memory for the trellis: metr = zeros( enc.stat, win+1 ) -INF; % path metric buffer metr( 1,1 ) = 0; % initial state => (0,0) surv = zeros( enc.stat, win+1 ); % survivor state buffer inpt = zeros( enc.stat, win+1 ); % survivor input buffer (dec. output) diff = zeros( enc.stat, win+1 ); % path metric difference comp = zeros( enc.stat, win+1 ); % competitor state buffer inpc = zeros( enc.stat, win+1 ); % competitor input buffer out = zeros( size(llr) ) + NaN; % hard output (bits) sft = zeros( size(llr) ) + INF; % soft output (sign with reliability) % decode all the bits: for i = 1:len, % indices + precalcuations: Cur = mod( i-1, win+1 ) +1; % curr trellis (cycl. buf) position Nxt = mod( i, win+1 ) +1; % next trellis (cycl. buf) position buf = msg( i*enc.n:-1:(i-1)*enc.n+1 ); % msg portion to be processed (reversed) llb = llr( (i-1)+1:i ); % SOVA: llr portion to be processed metr( :,Nxt ) = -INF -INF; % (2*) helps in initial stages (!!!) %% forward recursion: for s = 1:enc.stat, for j = 1:enc.ksym, nxt = enc.next.states( s, j ); % state after transition bin = enc.next.binout( s,:,j ); % transition output (encoder) mtr = bin*buf' + metr( s,Cur ); % transition metric mtr = mtr ... + enc.bininp( j )*(llb*enc.r)'; % SOVA if( metr( nxt,Nxt ) < mtr ), diff( nxt,Nxt ) = mtr - metr( nxt,Nxt ); % SOVA comp( nxt,Nxt ) = surv( nxt,Nxt ); % SOVA inpc( nxt,Nxt ) = inpt( nxt,Nxt ); % SOVA metr( nxt,Nxt ) = mtr; % store the metric surv( nxt,Nxt ) = s; % store the survival state inpt( nxt,Nxt ) = j-1; % store the survival input else dif = metr( nxt,Nxt ) - mtr; if( dif <= diff( nxt,Nxt ) ) diff( nxt,Nxt ) = dif; % SOVA comp( nxt,Nxt ) = s; % SOVA inpc( nxt,Nxt ) = j-1; % SOVA end end end end %% trace backwards: if( i < win ), continue; end; % proceed if the buffer has been filled; [ mtr, sur ] = max( metr( :,Nxt ) ); % find the intitial state (max metric) b = i; % temporary bit index clc = mod( Nxt-[1:win], win+1 ) +1; % indices in a 'cyclic buffer' operation for j = 1:win, % for all the bits in the buffer inp = inpt( sur, clc(j) ); % current bit-decoder output (encoder input) out( b ) = inp; % store the hard output tmp = clc( j ); cmp = comp( sur, tmp ); % SOVA: competitor state (previous) inc = inpc( sur, tmp ); % SOVA: competitor bit output dif = diff( sur, tmp ); % SOVA: corresp. path metric difference srv = surv( sur, tmp ); % SOVA: temporary survivor path state for k = j+1:win+1, % check all buffer bits srv and cmp paths if( inp ~= inc ), tmp = dif; idx = b - ( (k-1)-j ); % calculate index: [enc.k*(b-(k-1)+j-1)+1:enc.k*(b-(k-1)+j)] sft( idx ) = min( sft(idx), tmp ); % update LLRs for bits that are different end if( srv == cmp ), break; end; % stop if surv and comp merge (no need to continue) if( k == win+1 ), break; end; % stop if the end (otherwise: error) tmp = clc( k ); inp = inpt( srv, tmp ); % previous surv bit inc = inpt( cmp, tmp ); % previous comp bit srv = surv( srv, tmp ); % previous surv state cmp = surv( cmp, tmp ); % previous comp state end sur = surv( sur, clc(j) ); % state for the previous surv bit b = b - 1; % update bit index end end % provide soft output with +/- sign: out = (2*out-1) .* sft;return;function out = sovadec_KN( msg, llr, trl, win ) % error checking: if( ~istrellis( trl ) ), error( 'Incorrect input trellis!' ); end; if( nargin <= 2 ), error( 'Incorrect number of input args!' ); end; if( nargin == 3 ), win = length( llr )+1; end; % some parameters: INF = 9e9; % infinity enc = trellis2enc( trl ); % encoder parameters len = length( llr ) / enc.k; % number of decoding steps (total) win = min( [ win, len ] ); % trim the buffer if msg is short old = NaN; % to remember the last survivor % allocate memory for the trellis: metr = zeros( enc.stat, win+1 ) -INF; % path metric buffer metr( 1,1 ) = 0; % initial state => (0,0) surv = zeros( enc.stat, win+1 ); % survivor state buffer inpt = zeros( enc.stat, win+1 ); % survivor input buffer (dec. output) diff = zeros( enc.stat, win+1 ); % path metric difference comp = zeros( enc.stat, win+1 ); % competitor state buffer inpc = zeros( enc.stat, win+1 ); % competitor input buffer out = zeros( size(llr) ) + NaN; % hard output (bits) sft = zeros( size(llr) ) + INF; % soft output (sign with reliability) % decode all the bits: for i = 1:len, % indices + precalcuations: Cur = mod( i-1, win+1 ) +1; % curr trellis (cycl. buf) position Nxt = mod( i, win+1 ) +1; % next trellis (cycl. buf) position buf = msg( i*enc.n:-1:(i-1)*enc.n+1 ); % msg portion to be processed (reversed) llb = llr( (i-1)*enc.k+1:i*enc.k ); % SOVA: llr portion to be processed metr( :,Nxt ) = -INF -INF; % (2*) helps in initial stages (!!!) %% forward recursion: for s = 1:enc.stat, for j = 1:enc.ksym, nxt = enc.next.states( s, j ); % state after transition bin = enc.next.binout( s,:,j ); % transition output (encoder) mtr = bin*buf' + metr( s,Cur ); % transition metric mtr = mtr ... + enc.bininp(j,:)*(llb*enc.r)'; % SOVA if( metr( nxt,Nxt ) < mtr ), diff( nxt,Nxt ) = mtr - metr( nxt,Nxt ); % SOVA comp( nxt,Nxt ) = surv( nxt,Nxt ); % SOVA inpc( nxt,Nxt ) = inpt( nxt,Nxt ); % SOVA metr( nxt,Nxt ) = mtr; % store the metric surv( nxt,Nxt ) = s; % store the survival state inpt( nxt,Nxt ) = j-1; % store the survival input else dif = metr( nxt,Nxt ) - mtr; if( dif <= diff( nxt,Nxt ) ) diff( nxt,Nxt ) = dif; % SOVA comp( nxt,Nxt ) = s; % SOVA inpc( nxt,Nxt ) = j-1; % SOVA end end end end %% trace backwards: if( i < win ), continue; end; % proceed if the buffer has been filled; [ mtr, sur ] = max( metr( :,Nxt ) ); % find the intitial state (max metric) b = i; % temporary bit index clc = mod( Nxt-[1:win], win+1 ) +1; % indices in a 'cyclic buffer' operation for j = 1:win, % for all the bits in the buffer inp = inpt( sur, clc(j) ); % current bit-decoder output (encoder input) t = [ enc.k*(b-1)+1:enc.k*b ]; % compute the index out( t ) = enc.inp( inp+1,: ); % store the hard output cmp = comp( sur, clc(j) ); % SOVA: competitor state (previous) inc = inpc( sur, clc(j) ); % SOVA: competitor bit output dif = diff( sur, clc(j) ); % SOVA: corresp. path metric difference srv = surv( sur, clc(j) ); % SOVA: temporary survivor path state for k = j+1:win+1, % check all buffer bits srv and cmp paths inp = enc.inp( inp+1, : ); % convert to binary form inc = enc.inp( inc+1, : ); % convert to binary form tmp = ( inp == inc )*INF + dif; % for each different bit store the new dif idx = t - enc.k*( (k-1)-j ); % calculate index: [enc.k*(b-(k-1)+j-1)+1:enc.k*(b-(k-1)+j)] sft( idx ) = min( sft(idx), tmp ); % update LLRs for bits that are different if( srv == cmp ), break; end; % stop if surv and comp merge (no need to continue) if( k == win+1 ), break; end; % stop if the end (otherwise: error) inp = inpt( srv, clc(k) ); % previous surv bit srv = surv( srv, clc(k) ); % previous surv state inc = inpt( cmp, clc(k) ); % previous comp bit cmp = surv( cmp, clc(k) ); % previous comp state end sur = surv( sur, clc(j) ); % state for the previous surv bit b = b - 1; % update bit index end end % provide soft output with +/- sign: out = (2*out-1) .* sft;return;
