pabr-leansdr/src/leansdr/viterbi.h

299 wiersze
9.4 KiB
C++

// This file is part of LeanSDR Copyright (C) 2016-2018 <pabr@pabr.org>.
// See the toplevel README for more information.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
#ifndef LEANSDR_VITERBI_H
#define LEANSDR_VITERBI_H
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
// This is a generic implementation of Viterbi with explicit
// representation of the trellis. There is special support for
// convolutional coding, but the code can handle other schemes.
// TBD This is very inefficient. For a specific trellis all loops
// can be be unrolled.
namespace leansdr {
// TS is an integer type for a least NSTATES+1 states.
// NSTATES is the number of states (e.g. 2^(K-1)).
// TUS is an integer type for uncoded symbols (branch identifiers).
// NUS is the number of uncoded symbols.
// TCS is an integer type for coded symbols (branch labels).
// NCS is the number of coded symbols.
// TP is a type for representing paths.
// TPM, TBM are unsigned integer types for path/branch metrics.
// TPM is at least as wide as TBM.
template<typename TS, int NSTATES, typename TUS, int NUS, int NCS>
struct trellis {
static const int NOSTATE = NSTATES+1;
struct state {
struct branch {
TS pred; // Predecessor state or NOSTATE
TUS us; // Uncoded symbol
} branches[NCS]; // Incoming branches indexed by coded symbol
} states[NSTATES];
trellis() {
for ( TS s=0; s<NSTATES; ++s )
for ( int cs=0; cs<NCS; ++cs )
states[s].branches[cs].pred = NOSTATE;
}
// TBD Polynomial width should be a template parameter ?
void init_convolutional(const uint16_t G[]) {
if ( NCS & (NCS-1) ) {
fprintf(stderr, "NCS must be a power of 2\n");
exit(1);
}
// Derive number of polynomials from NCS.
int nG = log2i(NCS);
for ( TS s=0; s<NSTATES; ++s ) {
for ( TUS us=0; us<NUS; ++us ) {
// Run the convolutional encoder from state s with input us
uint64_t shiftreg = s; // TBD type
// Reverse bits
TUS us_rev = 0;
for ( int b=1; b<NUS; b*=2 ) if ( us & b ) us_rev |= (NUS/2/b);
shiftreg |= us_rev * NSTATES;
uint32_t cs = 0; // TBD type
for ( int g=0; g<nG; ++g )
cs = (cs<<1) | parity(shiftreg&G[g]);
shiftreg /= NUS; // Shift bits for 1 uncoded symbol
// [us] at state [s] emits [cs] and leads to state [shiftreg].
typename state::branch *b = &states[shiftreg].branches[cs];
if ( b->pred != NOSTATE ) {
fprintf(stderr, "Invalid convolutional code\n");
exit(1);
}
b->pred = s;
b->us = us;
}
}
}
void dump() {
for ( int s=0; s<NSTATES; ++s ) {
fprintf(stderr, "State %02x:", s);
for ( int cs=0; cs<NCS; ++cs ) {
typename state::branch *b = &states[s].branches[cs];
if ( b->pred == NOSTATE )
fprintf(stderr, " - ");
else
fprintf(stderr, " %02x+%x", b->pred, b->us);
}
fprintf(stderr, "\n");
}
}
};
// Interface that hides the templated internals.
template<typename TUS,
typename TCS,
typename TBM,
typename TPM>
struct viterbi_dec_interface {
virtual TUS update(TBM *costs, TPM *quality=NULL) = 0;
virtual TUS update(TCS s, TBM cost, TPM *quality=NULL) = 0;
};
template<typename TS, int NSTATES,
typename TUS, int NUS,
typename TCS, int NCS,
typename TBM, typename TPM,
typename TP>
struct viterbi_dec : viterbi_dec_interface<TUS,TCS,TBM,TPM> {
trellis<TS, NSTATES, TUS, NUS, NCS> *trell;
struct state {
TPM cost; // Metric of best path leading to this state
TP path; // Best path leading to this state
};
typedef state statebank[NSTATES];
state statebanks[2][NSTATES];
statebank *states, *newstates; // Alternate between banks
viterbi_dec(trellis<TS, NSTATES, TUS, NUS, NCS> *_trellis) :
trell(_trellis)
{
states = &statebanks[0];
newstates = &statebanks[1];
for ( TS s=0; s<NSTATES; ++s ) (*states)[s].cost = 0;
// Determine max value that can fit in TPM
max_tpm = (TPM)0 - 1;
if ( max_tpm < 0 ) {
// TPM is signed
for ( max_tpm=0; max_tpm*2+1>max_tpm; max_tpm=max_tpm*2+1 ) ;
}
}
// Update with full metric
TUS update(TBM costs[NCS], TPM *quality=NULL) {
TPM best_tpm = max_tpm, best2_tpm = max_tpm;
TS best_state = 0;
// Update all states
for ( int s=0; s<NSTATES; ++s ) {
TPM best_m = max_tpm;
typename trellis<TS,NSTATES,TUS,NUS,NCS>::state::branch *best_b = NULL;
// Select best branch
for ( int cs=0; cs<NCS; ++cs ) {
typename trellis<TS,NSTATES,TUS,NUS,NCS>::state::branch *b =
&trell->states[s].branches[cs];
if ( b->pred == trell->NOSTATE ) continue;
TPM m = (*states)[b->pred].cost + costs[cs];
if ( m <= best_m ) { // <= guarantees one match
best_m = m;
best_b = b;
}
}
(*newstates)[s].path = (*states)[best_b->pred].path;
(*newstates)[s].path.append(best_b->us);
(*newstates)[s].cost = best_m;
// Select best and second-best states
if ( best_m < best_tpm ) {
best_state = s;
best2_tpm = best_tpm;
best_tpm = best_m;
} else if ( best_m < best2_tpm )
best2_tpm = best_m;
}
// Swap banks
{ statebank *tmp=states; states=newstates; newstates=tmp; }
// Prevent overflow of path metrics
for ( TS s=0; s<NSTATES; ++s ) (*states)[s].cost -= best_tpm;
#if 0
// Observe that the min-max range remains bounded
fprintf(stderr,"-%2d = [", best_tpm);
for ( TS s=0; s<NSTATES; ++s ) fprintf(stderr," %d", (*states)[s].cost);
fprintf(stderr," ]\n");
#endif
// Return difference between best and second-best as quality metric.
if ( quality ) *quality = best2_tpm - best_tpm;
// Return uncoded symbol of best path
return (*states)[best_state].path.read();
}
// Update with partial metrics.
// The costs provided must be negative.
// The other symbols will be assigned a cost of 0.
TUS update(int nm, TCS cs[], TBM costs[], TPM *quality=NULL) {
TPM best_tpm = max_tpm, best2_tpm = max_tpm;
TS best_state = 0;
// Update all states
for ( int s=0; s<NSTATES; ++s ) {
// Select best branch among those for with metrics are provided
TPM best_m = max_tpm;
typename trellis<TS,NSTATES,TUS,NUS,NCS>::state::branch *best_b = NULL;
for ( int im=0; im<nm; ++im ) {
typename trellis<TS,NSTATES,TUS,NUS,NCS>::state::branch *b =
&trell->states[s].branches[cs[im]];
if ( b->pred == trell->NOSTATE ) continue;
TPM m = (*states)[b->pred].cost + costs[im];
if ( m <= best_m ) { // <= guarantees one match
best_m = m;
best_b = b;
}
}
if ( nm != NCS ) {
// Also scan the other branches.
// We actually rescan the branches with metrics.
// This works because costs are negative.
for ( int cs=0; cs<NCS; ++cs ) {
typename trellis<TS,NSTATES,TUS,NUS,NCS>::state::branch *b =
&trell->states[s].branches[cs];
if ( b->pred == trell->NOSTATE ) continue;
TPM m = (*states)[b->pred].cost;
if ( m <= best_m ) {
best_m = m;
best_b = b;
}
}
}
(*newstates)[s].path = (*states)[best_b->pred].path;
(*newstates)[s].path.append(best_b->us);
(*newstates)[s].cost = best_m;
// Select best states
if ( best_m < best_tpm ) {
best_state = s;
best2_tpm = best_tpm;
best_tpm = best_m;
} else if ( best_m < best2_tpm )
best2_tpm = best_m;
}
// Swap banks
{ statebank *tmp=states; states=newstates; newstates=tmp; }
// Prevent overflow of path metrics
for ( TS s=0; s<NSTATES; ++s ) (*states)[s].cost -= best_tpm;
#if 0
// Observe that the min-max range remains bounded
fprintf(stderr,"-%2d = [", best_tpm);
for ( TS s=0; s<NSTATES; ++s ) fprintf(stderr," %d", (*states)[s].cost);
fprintf(stderr," ]\n");
#endif
// Return difference between best and second-best as quality metric.
if ( quality ) *quality = best2_tpm - best_tpm;
// Return uncoded symbol of best path
return (*states)[best_state].path.read();
}
// Update with single-symbol metric.
// cost must be negative.
TUS update(TCS cs, TBM cost, TPM *quality=NULL) {
return update(1, &cs, &cost, quality);
}
void dump() {
fprintf(stderr, "[");
for ( TS s=0; s<NSTATES; ++s )
if ( states[s].cost )
fprintf(stderr, " %02x:%d", s, states[s].cost);
fprintf(stderr, "\n");
}
private:
TPM max_tpm;
};
// Paths (sequences of uncoded symbols) represented as bitstreams.
// NBITS is the number of bits per symbol.
// DEPTH is the number of symbols stored in the path.
// T is an unsigned integer type wider than NBITS*DEPTH.
template<typename T, typename TUS, int NBITS, int DEPTH>
struct bitpath {
T val;
bitpath() : val(0) { }
void append(TUS us) { val = (val<<NBITS) | us; }
TUS read() { return (val>>(DEPTH-1)*NBITS) & ((1<<NBITS)-1); }
};
} // namespace
#endif // LEANSDR_VITERBI_H