WSJT-X/libm65/vit216.c

220 lines
5.4 KiB
C
Raw Normal View History

/* Viterbi decoder for arbitrary convolutional code
* viterbi27 and viterbi37 for the r=1/2 and r=1/3 K=7 codes are faster
* Copyright 1999 Phil Karn, KA9Q
* May be used under the terms of the GNU Public License
*/
/* Select code here */
#define V216
#ifdef V216
#define K 16 /* Constraint length */
#define N 2 /* Number of symbols per data bit */
#define Polys Poly216 /* Select polynomials here */
#endif
/* Rate 1/2 codes */
unsigned int Poly216[] = {0126723, 0152711}; /* k = 16 */
#include <memory.h>
#define NULL ((void *)0)
#define LONGBITS 32
#define LOGLONGBITS 5
#undef max
#define max(x,y) ((x) > (y) ? (x) : (y))
#define D (1 << max(0,K-LOGLONGBITS-1))
#define MAXNBITS 200 /* Maximum frame size (user bits) */
extern unsigned char Partab[]; /* Parity lookup table */
int Syms[1 << K];
int VDInit = 0;
int parity(int x)
{
x ^= (x >> 16);
x ^= (x >> 8);
return Partab[x & 0xff];
}
// Wrapper for calling "encode" from Fortran:
//void __stdcall ENCODE(
void enc216_(
unsigned char data[], // User data, 8 bits per byte
int *nbits, // Number of user bits
unsigned char symbols[], // Encoded one-bit symbols, 8 per byte
int *nsymbols, // Number of symbols
int *kk, // K
int *nn) // N
{
int nbytes;
nbytes=(*nbits+7)/8; // Always encode multiple of 8 information bits
enc216(symbols,data,nbytes,0,0); // Do the encoding
*nsymbols=(*nbits+K-1)*N; // Return number of encoded symbols
*kk=K;
*nn=N;
}
/* Convolutionally encode data into binary symbols */
enc216(unsigned char symbols[], unsigned char data[],
unsigned int nbytes, unsigned int startstate,
unsigned int endstate)
{
int i,j,k,n=-1;
unsigned int encstate = startstate;
for(k=0; k<nbytes; k++) {
for(i=7;i>=0;i--){
encstate = (encstate + encstate) + ((data[k] >> i) & 1);
for(j=0;j<N;j++) {
n=n+1;
symbols[n] = parity(encstate & Polys[j]);
}
}
}
// Flush out with zero tail. (No need, if tail-biting code.)
for(i=0; i<K-1;i++){
encstate = (encstate << 1) | ((endstate >> i) & 1);
for(j=0;j<N;j++) {
n=n+1;
symbols[n] = parity(encstate & Polys[j]);
}
}
return 0;
}
// Wrapper for calling "viterbi" from Fortran:
//void __stdcall VITERBI(
void vit216_(
unsigned char symbols[], /* Raw deinterleaved input symbols */
unsigned int *Nbits, /* Number of decoded information bits */
int mettab[2][256], /* Metric table, [sent sym][rx symbol] */
unsigned char ddec[], /* Decoded output data */
long *Metric /* Final path metric (bigger is better) */
){
long metric;
vit216(&metric,ddec,symbols,*Nbits,mettab,0,0);
*Metric=metric;
}
/* Viterbi decoder */
int vit216(
long *metric, /* Final path metric (returned value) */
unsigned char *data, /* Decoded output data */
unsigned char *symbols, /* Raw deinterleaved input symbols */
unsigned int nbits, /* Number of output bits */
int mettab[2][256], /* Metric table, [sent sym][rx symbol] */
unsigned int startstate, /* Encoder starting state */
unsigned int endstate /* Encoder ending state */
){
int bitcnt = -(K-1);
long m0,m1;
int i,j,sym,ipp;
int mets[1 << N];
unsigned long paths[(MAXNBITS+K-1)*D];
unsigned long *pp,mask;
long cmetric[1 << (K-1)],nmetric[1 << (K-1)];
memset(paths,0,sizeof(paths));
// Initialize on first time through:
if(!VDInit){
for(i=0;i<(1<<K);i++){
sym = 0;
for(j=0;j<N;j++)
sym = (sym << 1) + parity(i & Polys[j]);
Syms[i] = sym;
}
VDInit++;
}
// Keep only lower K-1 bits of specified startstate and endstate
startstate &= ~((1<<(K-1)) - 1);
endstate &= ~((1<<(K-1)) - 1);
/* Initialize starting metrics */
for(i=0;i< 1<<(K-1);i++)
cmetric[i] = -999999;
cmetric[startstate] = 0;
pp = paths;
ipp=0;
for(;;){ /* For each data bit */
/* Read input symbols and compute branch metrics */
for(i=0;i< 1<<N;i++){
mets[i] = 0;
for(j=0;j<N;j++){
mets[i] += mettab[(i >> (N-j-1)) & 1][symbols[j]];
}
}
symbols += N;
/* Run the add-compare-select operations */
mask = 1;
for(i=0;i< 1 << (K-1);i+=2){
int b1,b2;
b1 = mets[Syms[i]];
nmetric[i] = m0 = cmetric[i/2] + b1;
b2 = mets[Syms[i+1]];
b1 -= b2;
m1 = cmetric[(i/2) + (1<<(K-2))] + b2;
if(m1 > m0){
nmetric[i] = m1;
*pp |= mask;
}
m0 -= b1;
nmetric[i+1] = m0;
m1 += b1;
if(m1 > m0){
nmetric[i+1] = m1;
*pp |= mask << 1;
}
mask <<= 2;
if(mask == 0){
mask = 1;
pp++;
ipp++;
}
}
if(mask != 1){
pp++;
ipp++;
}
if(++bitcnt == nbits){
*metric = nmetric[endstate];
break;
}
memcpy(cmetric,nmetric,sizeof(cmetric));
}
/* Chain back from terminal state to produce decoded data */
if(data == NULL)
return 0;/* Discard output */
memset(data,0,(nbits+7)/8); /* round up in case nbits % 8 != 0 */
for(i=nbits-1;i >= 0;i--){
// int a0,a1;
pp -= D;
ipp -= D;
m0=endstate >> LOGLONGBITS;
m1=1L << (endstate & (LONGBITS-1));
if(pp[m0] & m1) {
// a0=nmetric[endstate];
endstate |= (1 << (K-1));
// a1=nmetric[endstate];
data[i>>3] |= 0x80 >> (i&7);
// printf("B %d %d %d %d\n",*metric,i,a0,a1);
}
endstate >>= 1;
}
return 0;
}