求（2，1，7）维特比译码器的matlab程序

Posted 2023-05-05

非常紧急，请高手帮忙！！！

参考技术A 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
end
return;

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;

实时应用的维特比算法

【中文标题】实时应用的维特比算法

【英文标题】：Viterbi algorithm for real-time applications

【发布时间】：2011-12-22 07:34:31

【问题描述】：

我知道，给定一个 HMM 和一个观察，维特比算法可以猜测产生这个观察的隐藏状态序列。但是你想实时使用它的情况呢？我的意思是逐步找到隐藏状态。每次输入上的观察符号时，都会猜测隐藏状态，而无需知道接下来的整个观察序列。 我想将它用于实时运行的音频应用程序，因此观察将是每个时间帧的音频特征值的序列。

【问题讨论】：

你找到答案了吗？

【参考方案1】：

如果你有兴趣预测 T 时刻的隐藏状态，当你看到观察 O_T 时，你就有数据 O_1, ..., O_T-1, O_T。现在最可能的状态是向前向后找到，其中向后变量只是 1，因为我们看不到未来。总之，我们有 P(我们在时间 T 处于隐藏状态 i) = \alpha_T(i) / P(O_1, ..., O_T | \lambda)，其中 P(O_1, ..., O_T| \ λ) = \sum_i=1^n \alpha_T(i)。那么 P 的所有 i 的最大索引（我们在时间 T 处于隐藏状态 i）将是您的隐藏状态。

正式符号请参考http://courses.media.mit.edu/2010fall/mas622j/ProblemSets/ps4/tutorial.pdf。

如果这是您想要的，或者您有其他想法，请告诉我。如果您只想实时找到最佳状态序列，只需计算 alpha 变量，无需为此展望未来。

【讨论】：

在给定时间=t（现在）之前的状态和观察序列的情况下，如何预测时间=t+1（下一个时间步）的状态？

@Anil Vaitla 因为在时间T+1 的新观察可能会修改维特比路径，是否有可靠的方法可以找到从哪个时间t<T 无论新观察如何，维特比路径都不会改变?