VE7UWU
/
sdrangel
mirror of https://github.com/f4exb/sdrangel.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

DATV demodulator: leansdr: removed calls to exit

This commit is contained in:
f4exb 2018-02-26 01:02:33 +01:00
parent b8147ffacc
commit a483b58028
8 changed files with 1205 additions and 877 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

473
plugins/channelrx/demoddatv/leansdr/convolutional.h
View File

@ -1,103 +1,131 @@
#ifndef LEANSDR_CONVOLUTIONAL_H
#define LEANSDR_CONVOLUTIONAL_H
namespace leansdr {
namespace leansdr
{
// ALGEBRAIC DECONVOLUTION
// ALGEBRAIC DECONVOLUTION
// QPSK 1/2 only.
// This is a straightforward implementation, provided for reference.
// deconvol_poly2 is functionally equivalent and much faster.
// QPSK 1/2 only.
// This is a straightforward implementation, provided for reference.
// deconvol_poly2 is functionally equivalent and much faster.
template<typename Tin, // Input IQ symbols
typename Thist, // Input shift register (IQIQIQ...)
Thist POLY_DECONVOL, // Taps (IQIQIQ...)
Thist POLY_ERRORS>
struct deconvol_poly {
template<typename Tin, // Input IQ symbols
typename Thist, // Input shift register (IQIQIQ...)
Thist POLY_DECONVOL, // Taps (IQIQIQ...)
Thist POLY_ERRORS>
struct deconvol_poly
{
typedef u8 hardsymbol;
// Support soft of float input
inline u8 SYMVAL(const hardsymbol *s) { return *s; }
inline u8 SYMVAL(const softsymbol *s) { return s->symbol; }
inline u8 SYMVAL(const hardsymbol *s)
{
return *s;
}
inline u8 SYMVAL(const softsymbol *s)
{
return s->symbol;
}
typedef u8 decoded_byte;
deconvol_poly() : hist(0) { }
deconvol_poly() :
hist(0)
{
}
// Remap and deconvolve [nb*8] symbols into [nb] bytes.
// Return estimated number of bit errors.
int run(const Tin *pin, const u8 remap[], decoded_byte *pout, int nb) {
int nerrors = 0;
int halfway = nb / 2;
for ( ; nb--; ++pout ) {
decoded_byte byte = 0;
for ( int bit=8; bit--; ++pin) {
hist = (hist<<2) | remap[SYMVAL(pin)];
byte = (byte<<1) | parity(hist&POLY_DECONVOL);
if ( nb < halfway )
nerrors += parity(hist&POLY_ERRORS);
}
*pout = byte;
}
return nerrors;
int run(const Tin *pin, const u8 remap[], decoded_byte *pout, int nb)
{
int nerrors = 0;
int halfway = nb / 2;
for (; nb--; ++pout)
{
decoded_byte byte = 0;
for (int bit = 8; bit--; ++pin)
{
hist = (hist << 2) | remap[SYMVAL(pin)];
byte = (byte << 1) | parity(hist & POLY_DECONVOL);
if (nb < halfway)
nerrors += parity(hist & POLY_ERRORS);
}
*pout = byte;
}
return nerrors;
}
private:
private:
Thist hist;
}; // deconvol_poly
};
// deconvol_poly
// ALGEBRAIC DECONVOLUTION, OPTIMIZED
// ALGEBRAIC DECONVOLUTION, OPTIMIZED
// QPSK 1/2 only.
// Functionally equivalent to deconvol_poly,
// but processing 32 bits in parallel.
// QPSK 1/2 only.
// Functionally equivalent to deconvol_poly,
// but processing 32 bits in parallel.
template<typename Tin, // Input IQ symbols
typename Thist, // Input shift registers (one for I, one for Q)
typename Tpoly, // Taps (interleaved IQIQIQ...)
Tpoly POLY_DECONVOL,
Tpoly POLY_ERRORS>
struct deconvol_poly2 {
template<typename Tin, // Input IQ symbols
typename Thist, // Input shift registers (one for I, one for Q)
typename Tpoly, // Taps (interleaved IQIQIQ...)
Tpoly POLY_DECONVOL, Tpoly POLY_ERRORS>
struct deconvol_poly2
{
typedef u8 hardsymbol;
// Support instanciation of template with soft of float input
inline u8 SYMVAL(const hardsymbol *s) { return *s; }
inline u8 SYMVAL(const softsymbol *s) { return s->symbol; }
inline u8 SYMVAL(const hardsymbol *s)
{
return *s;
}
inline u8 SYMVAL(const softsymbol *s)
{
return s->symbol;
}
typedef u8 decoded_byte;
deconvol_poly2() : inI(0), inQ(0) { }
// Remap and deconvolve [nb*8] symbols into [nb] bytes.
// Return estimated number of bit errors.
deconvol_poly2() :
inI(0), inQ(0)
{
}
int run(const Tin *pin, const u8 remap[], decoded_byte *pout, int nb) {
if ( nb & (sizeof(Thist)-1) )
fail("Must deconvolve sizeof(Thist) bytes at a time");
nb /= sizeof(Thist);
unsigned long nerrors = 0;
int halfway = nb / 2;
Thist histI=inI, histQ=inQ;
for ( ; nb--; ) {
// This is where we convolve bits in parallel.
Thist wd = 0; // decoded bits
Thist we = 0; // error bits (should be 0)
// Remap and deconvolve [nb*8] symbols into [nb] bytes.
// Return estimated number of bit errors or -1 if error.
int run(const Tin *pin, const u8 remap[], decoded_byte *pout, int nb)
{
if (nb & (sizeof(Thist) - 1)) {
fail("deconvol_poly2::run", "Must deconvolve sizeof(Thist) bytes at a time");
return -1;
}
nb /= sizeof(Thist);
unsigned long nerrors = 0;
int halfway = nb / 2;
Thist histI = inI, histQ = inQ;
for (; nb--;)
{
// This is where we convolve bits in parallel.
Thist wd = 0; // decoded bits
Thist we = 0; // error bits (should be 0)
#if 0
// Trust gcc to unroll and evaluate the bit tests at compile-time.
for ( int bit=sizeof(Thist)*8; bit--; ++pin ) {
u8 iq = remap[SYMVAL(pin)];
histI = (histI<<1) | (iq>>1);
histQ = (histQ<<1) | (iq&1);
if ( POLY_DECONVOL & ((Tpoly)2<<(2*bit)) ) wd ^= histI;
if ( POLY_DECONVOL & ((Tpoly)1<<(2*bit)) ) wd ^= histQ;
if ( POLY_ERRORS & ((Tpoly)2<<(2*bit)) ) we ^= histI;
if ( POLY_ERRORS & ((Tpoly)1<<(2*bit)) ) we ^= histQ;
}
// Trust gcc to unroll and evaluate the bit tests at compile-time.
for ( int bit=sizeof(Thist)*8; bit--; ++pin )
{
u8 iq = remap[SYMVAL(pin)];
histI = (histI<<1) | (iq>>1);
histQ = (histQ<<1) | (iq&1);
if ( POLY_DECONVOL & ((Tpoly)2<<(2*bit)) ) wd ^= histI;
if ( POLY_DECONVOL & ((Tpoly)1<<(2*bit)) ) wd ^= histQ;
if ( POLY_ERRORS & ((Tpoly)2<<(2*bit)) ) we ^= histI;
if ( POLY_ERRORS & ((Tpoly)1<<(2*bit)) ) we ^= histQ;
}
#else
// Unroll manually.
// Unroll manually.
#define LOOP(bit) { \
u8 iq = remap[SYMVAL(pin)]; \
histI = (histI<<1) | (iq>>1); \
@ -108,150 +136,231 @@ namespace leansdr {
if ( POLY_ERRORS & ((Tpoly)1<<(2*bit)) ) we ^= histQ; \
++pin; \
}
// Don't shift by more than the operand width
switch ( sizeof(Thist)*8 ) {
// Don't shift by more than the operand width
switch (sizeof(Thist) * 8)
{
#if 0 // Not needed yet - avoid compiler warnings
case 64:
LOOP(63); LOOP(62); LOOP(61); LOOP(60);
LOOP(59); LOOP(58); LOOP(57); LOOP(56);
LOOP(55); LOOP(54); LOOP(53); LOOP(52);
LOOP(51); LOOP(50); LOOP(49); LOOP(48);
LOOP(47); LOOP(46); LOOP(45); LOOP(44);
LOOP(43); LOOP(42); LOOP(41); LOOP(40);
LOOP(39); LOOP(38); LOOP(37); LOOP(36);
LOOP(35); LOOP(34); LOOP(33); LOOP(32);
// Fall-through
case 64:
LOOP(63); LOOP(62); LOOP(61); LOOP(60);
LOOP(59); LOOP(58); LOOP(57); LOOP(56);
LOOP(55); LOOP(54); LOOP(53); LOOP(52);
LOOP(51); LOOP(50); LOOP(49); LOOP(48);
LOOP(47); LOOP(46); LOOP(45); LOOP(44);
LOOP(43); LOOP(42); LOOP(41); LOOP(40);
LOOP(39); LOOP(38); LOOP(37); LOOP(36);
LOOP(35); LOOP(34); LOOP(33); LOOP(32);
// Fall-through
#endif
case 32:
LOOP(31); LOOP(30); LOOP(29); LOOP(28);
LOOP(27); LOOP(26); LOOP(25); LOOP(24);
LOOP(23); LOOP(22); LOOP(21); LOOP(20);
LOOP(19); LOOP(18); LOOP(17); LOOP(16);
// Fall-through
case 16:
LOOP(15); LOOP(14); LOOP(13); LOOP(12);
LOOP(11); LOOP(10); LOOP( 9); LOOP( 8);
// Fall-through
case 8:
LOOP( 7); LOOP( 6); LOOP( 5); LOOP( 4);
LOOP( 3); LOOP( 2); LOOP( 1); LOOP( 0);
break;
default:
fail("Thist not supported");
}
case 32:
LOOP(31)
;
LOOP(30)
;
LOOP(29)
;
LOOP(28)
;
LOOP(27)
;
LOOP(26)
;
LOOP(25)
;
LOOP(24)
;
LOOP(23)
;
LOOP(22)
;
LOOP(21)
;
LOOP(20)
;
LOOP(19)
;
LOOP(18)
;
LOOP(17)
;
LOOP(16)
;
// Fall-through
case 16:
LOOP(15)
;
LOOP(14)
;
LOOP(13)
;
LOOP(12)
;
LOOP(11)
;
LOOP(10)
;
LOOP(9)
;
LOOP(8)
;
// Fall-through
case 8:
LOOP(7)
;
LOOP(6)
;
LOOP(5)
;
LOOP(4)
;
LOOP(3)
;
LOOP(2)
;
LOOP(1)
;
LOOP(0)
;
break;
default:
fail("deconvol_poly2::run", "Thist not supported (1)");
return -1;
}
#undef LOOP
#endif
switch ( sizeof(Thist)*8 ) {
switch (sizeof(Thist) * 8)
{
#if 0 // Not needed yet - avoid compiler warnings
case 64:
*pout++ = wd >> 56;
*pout++ = wd >> 48;
*pout++ = wd >> 40;
*pout++ = wd >> 32;
// Fall-through
case 64:
*pout++ = wd >> 56;
*pout++ = wd >> 48;
*pout++ = wd >> 40;
*pout++ = wd >> 32;
// Fall-through
#endif
case 32:
*pout++ = wd >> 24;
*pout++ = wd >> 16;
// Fall-through
case 16:
*pout++ = wd >> 8;
// Fall-through
case 8:
*pout++ = wd;
break;
default:
fail("Thist not supported");
}
// Count errors when the shift registers are full
if ( nb < halfway ) nerrors += hamming_weight(we);
}
inI = histI;
inQ = histQ;
return nerrors;
case 32:
*pout++ = wd >> 24;
*pout++ = wd >> 16;
// Fall-through
case 16:
*pout++ = wd >> 8;
// Fall-through
case 8:
*pout++ = wd;
break;
default:
fail("deconvol_poly2::run", "Thist not supported (2)");
return -1;
}
// Count errors when the shift registers are full
if (nb < halfway)
nerrors += hamming_weight(we);
}
inI = histI;
inQ = histQ;
return nerrors;
}
private:
private:
Thist inI, inQ;
}; // deconvol_poly2
};
// deconvol_poly2
// CONVOLUTIONAL ENCODER
// CONVOLUTIONAL ENCODER
// QPSK 1/2 only.
// QPSK 1/2 only.
template<typename Thist, uint64_t POLY1, uint64_t POLY2>
struct convol_poly2 {
template<typename Thist, uint64_t POLY1, uint64_t POLY2>
struct convol_poly2
{
typedef u8 uncoded_byte;
typedef u8 hardsymbol;
convol_poly2() : hist(0) { }
convol_poly2() :
hist(0)
{
}
// Convolve [count] bytes into [count*8] symbols, and remap.
void run(const uncoded_byte *pin, const u8 remap[],
hardsymbol *pout, int count) {
for ( ; count--; ++pin ) {
uncoded_byte b = *pin;
for ( int bit=8; bit--; ++pout ) {
hist = (hist>>1) | (((b>>bit)&1)<<6);
u8 s = (parity(hist&POLY1)<<1) | parity(hist&POLY2);
*pout = remap[s];
}
}
void run(const uncoded_byte *pin, const u8 remap[], hardsymbol *pout, int count)
{
for (; count--; ++pin)
{
uncoded_byte b = *pin;
for (int bit = 8; bit--; ++pout)
{
hist = (hist >> 1) | (((b >> bit) & 1) << 6);
u8 s = (parity(hist & POLY1) << 1) | parity(hist & POLY2);
*pout = remap[s];
}
}
}
private:
private:
Thist hist;
}; // convol_poly2
};
// convol_poly2
// Generic BPSK..256QAM and puncturing
// Generic BPSK..256QAM and puncturing
template<typename Thist, int HISTSIZE>
struct convol_multipoly {
template<typename Thist, int HISTSIZE>
struct convol_multipoly
{
typedef u8 uncoded_byte;
typedef u8 hardsymbol;
int bits_in, bits_out, bps;
const Thist *polys; // [bits_out]
convol_multipoly()
: bits_in(0), bits_out(0), bps(0),
hist(0), nhist(0), sersymb(0), nsersymb(0)
{ }
void encode(const uncoded_byte *pin, hardsymbol *pout, int count) {
if ( !bits_in || !bits_out || !bps )
fatal("convol_multipoly not configured");
hardsymbol symbmask = (1<<bps) - 1;
for ( ; count--; ++pin ) {
uncoded_byte b = *pin;
for ( int bit=8; bit--; ) {
hist = (hist>>1) | ((Thist)((b>>bit)&1)<<(HISTSIZE-1));
++nhist;
if ( nhist == bits_in ) {
for ( int p=0; p<bits_out; ++p ) {
int b = parity((Thist)(hist & polys[p]));
sersymb = (sersymb<<1) | b;
}
nhist = 0;
nsersymb += bits_out;
while ( nsersymb >= bps ) {
hardsymbol s = (sersymb >> (nsersymb-bps)) & symbmask;
*pout++ = s;
nsersymb -= bps;
}
}
}
}
// Ensure deterministic output size
// TBD We can relax this
if ( nhist || nsersymb ) fatal("partial run");
convol_multipoly() :
bits_in(0), bits_out(0), bps(0), polys(0), hist(0), nhist(0), sersymb(0), nsersymb(0)
{
}
private:
void encode(const uncoded_byte *pin, hardsymbol *pout, int count)
{
if (!bits_in || !bits_out || !bps)
{
fatal("leansdr::convol_multipoly::encode: convol_multipoly not configured");
return;
}
hardsymbol symbmask = (1 << bps) - 1;
for (; count--; ++pin)
{
uncoded_byte b = *pin;
for (int bit = 8; bit--;)
{
hist = (hist >> 1) | ((Thist) ((b >> bit) & 1) << (HISTSIZE - 1));
++nhist;
if (nhist == bits_in)
{
for (int p = 0; p < bits_out; ++p)
{
int b = parity((Thist) (hist & polys[p]));
sersymb = (sersymb << 1) | b;
}
nhist = 0;
nsersymb += bits_out;
while (nsersymb >= bps)
{
hardsymbol s = (sersymb >> (nsersymb - bps)) & symbmask;
*pout++ = s;
nsersymb -= bps;
}
}
}
}
// Ensure deterministic output size
// TBD We can relax this
if (nhist || nsersymb) {
fatal("leansdr::convol_multipoly::encode: partial run");
}
}
private:
Thist hist;
int nhist;
Thist sersymb;
int nsersymb;
}; // convol_multipoly
};
// convol_multipoly
} // namespace
}// namespace
#endif // LEANSDR_CONVOLUTIONAL_H

6
plugins/channelrx/demoddatv/leansdr/dsp.h
View File

@ -362,8 +362,10 @@ struct fir_resampler: runnable
tap_multiplier(1),
freq_tol(0.1)
{
if (decim != 1)
fail("fir_resampler: decim not implemented"); // TBD
if (decim != 1) {
fail("fir_resampler::fir_resampler", "decim not implemented"); // TBD
return;
}
shifted_coeffs = new T[ncoeffs];
set_freq(0);
}

66
plugins/channelrx/demoddatv/leansdr/dvb.h
View File

@ -59,7 +59,8 @@ inline cstln_lut<256> * make_dvbs2_constellation(cstln_lut<256>::predef c,
gamma1 = 2.57;
break;
default:
fail("Code rate not supported with APSK16");
fail("cstln_lut<256>::make_dvbs2_constellation", "Code rate not supported with APSK16");
return 0;
}
break;
case cstln_lut<256>::APSK32:
@ -87,7 +88,8 @@ inline cstln_lut<256> * make_dvbs2_constellation(cstln_lut<256>::predef c,
gamma2 = 4.30;
break;
default:
fail("Code rate not supported with APSK32");
fail("cstln_lut<256>::make_dvbs2_constellation", "Code rate not supported with APSK32");
return 0;
}
break;
case cstln_lut<256>::APSK64E:
@ -213,7 +215,10 @@ struct deconvol_sync: runnable
void next_sync()
{
if (fastlock)
fail("Bug: next_sync() called with fastlock");
{
fail("deconvol_sync::next_sync", "Bug: next_sync() called with fastlock");
return;
}
++locked;
if (locked == &syncs[NSYNCS])
{
@ -325,7 +330,10 @@ private:
if (d == 0x00000000fb11d640LL)
d2 = 0x00fbea3c8679c980LL;
if (d2 == d)
fail("Alt polynomial not provided");
{
fail("deconvol_sync::inverse_convolution", "Alt polynomial not provided");
return;
}
deconv2[b] = d2;
}
@ -346,17 +354,27 @@ private:
int expect = (b == i) ? 1 : 0;
int d = parity(iq & deconv[b]);
if (d != expect)
fail("Failed to inverse convolutional coding");
{
fail("deconvol_sync::inverse_convolution", "Failed to inverse convolutional coding");
}
int d2 = parity(iq & deconv2[b]);
if (d2 != expect)
fail("Failed to inverse convolutional coding (alt)");
{
fail("deconvol_sync::inverse_convolution", "Failed to inverse convolutional coding (alt)");
}
}
if (traceback > sizeof(iq_t) * 8)
fail("Bug: traceback exceeds register size");
{
fail("deconvol_sync::inverse_convolution", "Bug: traceback exceeds register size");
}
if (log2(deconv[b]) + 1 > traceback)
fail("traceback insufficient for deconvolution");
{
fail("deconvol_sync::inverse_convolution", "traceback insufficient for deconvolution");
}
if (log2(deconv2[b]) + 1 > traceback)
fail("traceback insufficient for deconvolution (alt)");
{
fail("deconvol_sync::inverse_convolution", "traceback insufficient for deconvolution (alt)");
}
}
}
@ -668,14 +686,20 @@ struct dvb_convol: runnable
{
fec_spec *fs = &fec_specs[fec];
if (!fs->bits_in)
fail("Unexpected FEC");
{
fail("dvb_convol::dvb_convol", "Unexpected FEC");
return;
}
convol.bits_in = fs->bits_in;
convol.bits_out = fs->bits_out;
convol.polys = fs->polys;
convol.bps = bits_per_symbol;
// FEC must output a whole number of IQ symbols
if (convol.bits_out % convol.bps)
fail("Code rate not suitable for this constellation");
{
fail("dvb_convol::dvb_convol", "Code rate not suitable for this constellation");
return;
}
}
void run()
@ -1187,7 +1211,10 @@ struct rs_decoder: runnable
if (sizeof(Tbyte) == 1)
memcpy(pout, pin, SIZE_TSPACKET);
else
fail("Erasures not implemented");
{
fail("rs_decoder::run", "Erasures not implemented");
return;
}
u8 synd[16];
bool corrupted = rs.syndromes(pin, synd);
@ -1448,7 +1475,10 @@ public:
{ // Sanity check: FEC block size must be a multiple of label size.
int symbols_per_block = fec->bits_out / bits_per_symbol;
if (bits_per_symbol * symbols_per_block != fec->bits_out)
fail("Code rate not suitable for this constellation");
{
fail("viterbi_sync::viterbi_sync", "Code rate not suitable for this constellation");
return;
}
}
int nconj;
switch (cstln->nsymbols)
@ -1554,7 +1584,10 @@ public:
syncs[s].dec = new dvb_dec_78(trell);
}
else
fail("CR not supported");
{
fail("viterbi_sync::viterbi_sync", "CR not supported");
return;
}
}
@ -1640,7 +1673,10 @@ public:
} // chunk_size
in.read(chunk_size * nshifts);
if (nout)
fail("overlapping out");
{
fail("viterbi_sync::run", "overlapping out");
return;
}
if (!resync_phase)
{
// Switch to another decoder ?

35
plugins/channelrx/demoddatv/leansdr/framework.h
View File

@ -13,13 +13,11 @@ namespace leansdr
inline void fatal(const char *s)
{
perror(s);
exit(1);
}
inline void fail(const char *s)
inline void fail(const char *f, const char *s)
{
fprintf(stderr, "leansdr::fail: %s\n", s);
exit(1);
fprintf(stderr, "leansdr::%s: %s\n", f, s);
}
//////////////////////////////////////////////////////////////////////
@ -114,16 +112,19 @@ struct scheduler
void add_pipe(pipebuf_common *p)
{
if (npipes == MAX_PIPES) {
fail("MAX_PIPES");
if (npipes == MAX_PIPES)
{
fail("scheduler::add_pipe", "MAX_PIPES");
return;
}
pipes[npipes++] = p;
}
void add_runnable(runnable_common *r)
{
if (nrunnables == MAX_RUNNABLES) {
fail("MAX_RUNNABLES");
if (nrunnables == MAX_RUNNABLES)
{
fail("scheduler::add_runnable", "MAX_RUNNABLES");
}
runnables[nrunnables++] = r;
}
@ -211,8 +212,10 @@ struct pipebuf: pipebuf_common
int add_reader()
{
if (nrd == MAX_READERS) {
fail("too many readers");
if (nrd == MAX_READERS)
{
fail("pipebuf::add_reader", "too many readers");
return nrd;
}
rds[nrd] = wr;
return nrd++;
@ -295,8 +298,8 @@ struct pipewriter
{
if (buf.end < buf.wr)
{
fprintf(stderr, "leansdr::pipewriter::writable: Bug: overflow to %s\n", buf.name);
exit(1);
fprintf(stderr, "leansdr::pipewriter::writable: overflow in %s buffer\n", buf.name);
return 0;
}
unsigned long delta = buf.end - buf.wr;
@ -317,8 +320,8 @@ struct pipewriter
{
if (buf.wr + n > buf.end)
{
fprintf(stderr, "leansdr::pipewriter::written: Bug: overflow to %s\n", buf.name);
exit(1);
fprintf(stderr, "leansdr::pipewriter::written: overflow in %s buffer\n", buf.name);
return;
}
buf.wr += n;
buf.total_written += n;
@ -378,8 +381,8 @@ struct pipereader
{
if (buf.rds[id] + n > buf.wr)
{
fprintf(stderr, "leansdr::pipereader::read: Bug: underflow from %s\n", buf.name);
exit(1);
fprintf(stderr, "leansdr::pipereader::read: underflow in %s buffer\n", buf.name);
return;
}
buf.rds[id] += n;
buf.total_read += n;

555
plugins/channelrx/demoddatv/leansdr/generic.h
View File

@ -6,7 +6,8 @@
#include "leansdr/math.h"
namespace leansdr {
namespace leansdr
{
//////////////////////////////////////////////////////////////////////
// Simple blocks
@ -16,103 +17,144 @@ namespace leansdr {
// If the file descriptor is seekable, data can be looped.
template<typename T>
struct file_reader : runnable {
file_reader(scheduler *sch, int _fdin, pipebuf<T> &_out)
: runnable(sch, _out.name),
loop(false),
fdin(_fdin), out(_out)
{
}
void run() {
size_t size = out.writable() * sizeof(T);
if ( ! size ) return;
again:
ssize_t nr = read(fdin, out.wr(), size);
if ( nr < 0 ) fatal("read");
if ( ! nr ) {
if ( ! loop ) return;
if ( sch->debug ) fprintf(stderr, "%s looping\n", name);
off_t res = lseek(fdin, 0, SEEK_SET);
if ( res == (off_t)-1 ) fatal("lseek");
goto again;
struct file_reader: runnable
{
file_reader(scheduler *sch, int _fdin, pipebuf<T> &_out) :
runnable(sch, _out.name), loop(false), fdin(_fdin), out(_out)
{
}
void run()
{
size_t size = out.writable() * sizeof(T);
if (!size)
return;
// Always stop at element boundary (may block)
size_t partial = nr % sizeof(T);
size_t remain = partial ? sizeof(T)-partial : 0;
while ( remain ) {
if ( sch->debug ) fprintf(stderr, "+");
ssize_t nr2 = read(fdin, (char*)out.wr()+nr, remain);
if ( nr2 <= 0 ) fatal("partial read");
nr += nr2;
remain -= nr2;
again: ssize_t nr = read(fdin, out.wr(), size);
if (nr < 0)
{
fatal("leansdr::file_reader::run: read");
return;
}
if (!nr)
{
if (!loop)
return;
if (sch->debug)
fprintf(stderr, "leansdr::file_reader::run: %s looping\n", name);
off_t res = lseek(fdin, 0, SEEK_SET);
if (res == (off_t) -1)
{
fatal("leansdr::file_reader::run: lseek");
return;
}
goto again;
}
// Always stop at element boundary (may block)
size_t partial = nr % sizeof(T);
size_t remain = partial ? sizeof(T) - partial : 0;
while (remain)
{
if (sch->debug)
fprintf(stderr, "+");
ssize_t nr2 = read(fdin, (char*) out.wr() + nr, remain);
if (nr2 <= 0)
{
fatal("leansdr::file_reader::run: partial read");
return;
}
nr += nr2;
remain -= nr2;
}
out.written(nr / sizeof(T));
}
out.written(nr / sizeof(T));
}
bool loop;
bool loop;
private:
int fdin;
pipewriter<T> out;
int fdin;
pipewriter<T> out;
};
// [file_writer] writes raw data from a [pipebuf] to a file descriptor.
template<typename T>
struct file_writer : runnable {
file_writer(scheduler *sch, pipebuf<T> &_in, int _fdout) :
runnable(sch, _in.name),
in(_in), fdout(_fdout) {
}
void run() {
int size = in.readable() * sizeof(T);
if ( ! size ) return;
int nw = write(fdout, in.rd(), size);
if ( ! nw ) fatal("pipe");
if ( nw < 0 ) fatal("write");
if ( nw % sizeof(T) ) fatal("partial write");
in.read(nw/sizeof(T));
}
struct file_writer: runnable
{
file_writer(scheduler *sch, pipebuf<T> &_in, int _fdout) :
runnable(sch, _in.name), in(_in), fdout(_fdout)
{
}
void run()
{
int size = in.readable() * sizeof(T);
if (!size)
return;
int nw = write(fdout, in.rd(), size);
if (!nw)
{
fatal("leansdr::file_writer::run: pipe");
return;
}
if (nw < 0)
{
fatal("leansdr::file_writer::run: write");
return;
}
if (nw % sizeof(T))
{
fatal("leansdr::file_writer::run:partial write");
return;
}
in.read(nw / sizeof(T));
}
private:
pipereader<T> in;
int fdout;
pipereader<T> in;
int fdout;
};
// [file_printer] writes data from a [pipebuf] to a file descriptor,
// with printf-style formatting and optional scaling.
template<typename T>
struct file_printer : runnable {
file_printer(scheduler *sch, const char *_format,
pipebuf<T> &_in, int _fdout,
int _decimation=1) :
runnable(sch, _in.name),
scale(1), decimation(_decimation),
in(_in), format(_format), fdout(_fdout), phase(0) {
}
void run() {
int n = in.readable();
T *pin=in.rd(), *pend=pin+n;
for ( ; pin<pend; ++pin ) {
if ( ++phase >= decimation ) {
phase -= decimation;
char buf[256];
int len = snprintf(buf, sizeof(buf), format, (*pin)*scale);
if ( len < 0 ) fatal("obsolete glibc");
int nw = write(fdout, buf, len);
if ( nw != len ) fatal("partial write");
}
struct file_printer: runnable
{
file_printer(scheduler *sch, const char *_format, pipebuf<T> &_in, int _fdout, int _decimation = 1) :
runnable(sch, _in.name), scale(1), decimation(_decimation), in(_in), format(_format), fdout(_fdout), phase(0)
{
}
in.read(n);
}
T scale;
int decimation;
void run()
{
int n = in.readable();
T *pin = in.rd(), *pend = pin + n;
for (; pin < pend; ++pin)
{
if (++phase >= decimation)
{
phase -= decimation;
char buf[256];
int len = snprintf(buf, sizeof(buf), format, (*pin) * scale);
if (len < 0)
{
fatal("leansdr::file_printer::run: obsolete glibc");
return;
}
int nw = write(fdout, buf, len);
if (nw != len)
{
fatal("leansdr::file_printer::run: partial write");
return;
}
}
}
in.read(n);
}
T scale;
int decimation;
private:
pipereader<T> in;
const char *format;
int fdout;
int phase;
pipereader<T> in;
const char *format;
int fdout;
int phase;
};
// [file_carrayprinter] writes all data available from a [pipebuf]
@ -120,229 +162,246 @@ private:
// Special case for complex.
template<typename T>
struct file_carrayprinter : runnable {
file_carrayprinter(scheduler *sch,
const char *_head,
const char *_format,
const char *_sep,
const char *_tail,
pipebuf< complex<T> > &_in, int _fdout) :
runnable(sch, _in.name),
scale(1), fixed_size(0), in(_in),
head(_head), format(_format), sep(_sep), tail(_tail),
fout(fdopen(_fdout,"w")) {
}
void run() {
int n, nmin = fixed_size ? fixed_size : 1;
while ( (n=in.readable()) >= nmin ) {
if ( fixed_size ) n = fixed_size;
if ( fout ) {
fprintf(fout, head, n);
complex<T> *pin = in.rd();
for ( int i=0; i<n; ++i ) {
if ( i ) fprintf(fout, "%s", sep);
fprintf(fout, format, pin[i].re*scale, pin[i].im*scale);
}
fprintf(fout, "%s", tail);
}
fflush(fout);
in.read(n);
struct file_carrayprinter: runnable
{
file_carrayprinter(scheduler *sch, const char *_head, const char *_format, const char *_sep, const char *_tail, pipebuf<complex<T> > &_in, int _fdout) :
runnable(sch, _in.name), scale(1), fixed_size(0), in(_in), head(_head), format(_format), sep(_sep), tail(_tail), fout(fdopen(_fdout, "w"))
{
}
}
T scale;
int fixed_size; // Number of elements per batch, or 0.
void run()
{
int n, nmin = fixed_size ? fixed_size : 1;
while ((n = in.readable()) >= nmin)
{
if (fixed_size)
n = fixed_size;
if (fout)
{
fprintf(fout, head, n);
complex<T> *pin = in.rd();
for (int i = 0; i < n; ++i)
{
if (i)
fprintf(fout, "%s", sep);
fprintf(fout, format, pin[i].re * scale, pin[i].im * scale);
}
fprintf(fout, "%s", tail);
}
fflush(fout);
in.read(n);
}
}
T scale;
int fixed_size; // Number of elements per batch, or 0.
private:
pipereader< complex<T> > in;
const char *head, *format, *sep, *tail;
FILE *fout;
pipereader<complex<T> > in;
const char *head, *format, *sep, *tail;
FILE *fout;
};
template<typename T, int N>
struct file_vectorprinter : runnable {
file_vectorprinter(scheduler *sch,
const char *_head,
const char *_format,
const char *_sep,
const char *_tail,
pipebuf<T[N]> &_in, int _fdout) :
runnable(sch, _in.name), scale(1), in(_in),
head(_head), format(_format), sep(_sep), tail(_tail) {
fout = fdopen(_fdout,"w");
if ( ! fout ) fatal("fdopen");
}
void run() {
while ( in.readable() >= 1 ) {
fprintf(fout, head, N);
T (*pin)[N] = in.rd();
for ( int i=0; i<N; ++i ) {
if ( i ) fprintf(fout, "%s", sep);
fprintf(fout, format, (*pin)[i]*scale);
}
fprintf(fout, "%s", tail);
in.read(1);
struct file_vectorprinter: runnable
{
file_vectorprinter(scheduler *sch, const char *_head, const char *_format, const char *_sep, const char *_tail, pipebuf<T[N]> &_in, int _fdout) :
runnable(sch, _in.name), scale(1), in(_in), head(_head), format(_format), sep(_sep), tail(_tail)
{
fout = fdopen(_fdout, "w");
if (!fout)
{
fatal("leansdr::file_vectorprinter::file_vectorprinter: fdopen");
}
}
fflush(fout);
}
T scale;
void run()
{
while (in.readable() >= 1)
{
fprintf(fout, head, N);
T (*pin)[N] = in.rd();
for (int i = 0; i < N; ++i)
{
if (i)
fprintf(fout, "%s", sep);
fprintf(fout, format, (*pin)[i] * scale);
}
fprintf(fout, "%s", tail);
in.read(1);
}
fflush(fout);
}
T scale;
private:
pipereader<T[N]> in;
const char *head, *format, *sep, *tail;
FILE *fout;
pipereader<T[N]> in;
const char *head, *format, *sep, *tail;
FILE *fout;
};
// [itemcounter] writes the number of input items to the output [pipebuf].
// [Tout] must be a numeric type.
template<typename Tin, typename Tout>
struct itemcounter : runnable {
itemcounter(scheduler *sch, pipebuf<Tin> &_in, pipebuf<Tout> &_out)
: runnable(sch, "itemcounter"),
in(_in), out(_out) {
}
void run() {
if ( out.writable() < 1 ) return;
unsigned long count = in.readable();
if ( ! count ) return;
out.write(count);
in.read(count);
}
struct itemcounter: runnable
{
itemcounter(scheduler *sch, pipebuf<Tin> &_in, pipebuf<Tout> &_out) :
runnable(sch, "itemcounter"), in(_in), out(_out)
{
}
void run()
{
if (out.writable() < 1)
return;
unsigned long count = in.readable();
if (!count)
return;
out.write(count);
in.read(count);
}
private:
pipereader<Tin> in;
pipewriter<Tout> out;
pipereader<Tin> in;
pipewriter<Tout> out;
};
// [decimator] forwards 1 in N sample.
template<typename T>
struct decimator : runnable {
unsigned int d;
struct decimator: runnable
{
unsigned int d;
decimator(scheduler *sch, int _d, pipebuf<T> &_in, pipebuf<T> &_out)
: runnable(sch, "decimator"),
d(_d),
in(_in), out(_out) {
}
void run() {
unsigned long count = min(in.readable()/d, out.writable());
T *pin=in.rd(), *pend=pin+count*d, *pout=out.wr();
for ( ; pin<pend; pin+=d, ++pout )
*pout = *pin;
in.read(count*d);
out.written(count);
}
decimator(scheduler *sch, int _d, pipebuf<T> &_in, pipebuf<T> &_out) :
runnable(sch, "decimator"), d(_d), in(_in), out(_out)
{
}
void run()
{
unsigned long count = min(in.readable() / d, out.writable());
T *pin = in.rd(), *pend = pin + count * d, *pout = out.wr();
for (; pin < pend; pin += d, ++pout)
*pout = *pin;
in.read(count * d);
out.written(count);
}
private:
pipereader<T> in;
pipewriter<T> out;
pipereader<T> in;
pipewriter<T> out;
};
// [rate_estimator] accumulates counts of two quantities
// and periodically outputs their ratio.
// [rate_estimator] accumulates counts of two quantities
// and periodically outputs their ratio.
template<typename T>
struct rate_estimator : runnable {
template<typename T>
struct rate_estimator: runnable
{
int sample_size;
rate_estimator(scheduler *sch,
pipebuf<int> &_num, pipebuf<int> &_den,
pipebuf<float> &_rate)
: runnable(sch, "rate_estimator"),
sample_size(10000),
num(_num), den(_den), rate(_rate),
acc_num(0), acc_den(0) {
rate_estimator(scheduler *sch, pipebuf<int> &_num, pipebuf<int> &_den, pipebuf<float> &_rate) :
runnable(sch, "rate_estimator"), sample_size(10000), num(_num), den(_den), rate(_rate), acc_num(0), acc_den(0)
{
}
void run() {
if ( rate.writable() < 1 ) return;
int count = min(num.readable(), den.readable());
int *pnum=num.rd(), *pden=den.rd();
for ( int n=count; n--; ++pnum,++pden ) {
acc_num += *pnum;
acc_den += *pden;
}
num.read(count);
den.read(count);
if ( acc_den >= sample_size ) {
rate.write((float)acc_num / acc_den);
acc_num = acc_den = 0;
}
void run()
{
if (rate.writable() < 1)
return;
int count = min(num.readable(), den.readable());
int *pnum = num.rd(), *pden = den.rd();
for (int n = count; n--; ++pnum, ++pden)
{
acc_num += *pnum;
acc_den += *pden;
}
num.read(count);
den.read(count);
if (acc_den >= sample_size)
{
rate.write((float) acc_num / acc_den);
acc_num = acc_den = 0;
}
}
private:
private:
pipereader<int> num, den;
pipewriter<float> rate;
T acc_num, acc_den;
};
};
// SERIALIZER
// SERIALIZER
template<typename Tin, typename Tout>
struct serializer : runnable {
serializer(scheduler *sch, pipebuf<Tin> &_in, pipebuf<Tout> &_out)
: nin(max((size_t)1,sizeof(Tin)/sizeof(Tout))),
nout(max((size_t)1,sizeof(Tout)/sizeof(Tin))),
in(_in), out(_out,nout)
template<typename Tin, typename Tout>
struct serializer: runnable
{
serializer(scheduler *sch, pipebuf<Tin> &_in, pipebuf<Tout> &_out) :
nin(max((size_t) 1, sizeof(Tin) / sizeof(Tout))), nout(max((size_t) 1, sizeof(Tout) / sizeof(Tin))), in(_in), out(_out, nout)
{
if ( nin*sizeof(Tin) != nout*sizeof(Tout) )
fail("serializer: incompatible sizes");
if (nin * sizeof(Tin) != nout * sizeof(Tout))
{
fail("serializer::serializer", "incompatible sizes");
return;
}
}
void run() {
while ( in.readable()>=nin && out.writable()>=nout ) {
memcpy(out.wr(), in.rd(), nout*sizeof(Tout));
in.read(nin);
out.written(nout);
}
void run()
{
while (in.readable() >= nin && out.writable() >= nout)
{
memcpy(out.wr(), in.rd(), nout * sizeof(Tout));
in.read(nin);
out.written(nout);
}
}
private:
private:
int nin, nout;
pipereader<Tin> in;
pipewriter<Tout> out;
}; // serializer
};
// serializer
// [buffer_reader] reads from a user-supplied buffer.
// [buffer_reader] reads from a user-supplied buffer.
template<typename T>
struct buffer_reader : runnable {
buffer_reader(scheduler *sch, T *_data, int _count, pipebuf<T> &_out)
: runnable(sch, "buffer_reader"),
data(_data), count(_count), out(_out), pos(0) {
template<typename T>
struct buffer_reader: runnable
{
buffer_reader(scheduler *sch, T *_data, int _count, pipebuf<T> &_out) :
runnable(sch, "buffer_reader"), data(_data), count(_count), out(_out), pos(0)
{
}
void run() {
int n = min(out.writable(), (unsigned long)(count-pos));
memcpy(out.wr(), &data[pos], n*sizeof(T));
pos += n;
out.written(n);
void run()
{
int n = min(out.writable(), (unsigned long) (count - pos));
memcpy(out.wr(), &data[pos], n * sizeof(T));
pos += n;
out.written(n);
}
private:
private:
T *data;
int count;
pipewriter<T> out;
int pos;
}; // buffer_reader
};
// buffer_reader
// [buffer_writer] writes to a user-supplied buffer.
// [buffer_writer] writes to a user-supplied buffer.
template<typename T>
struct buffer_writer : runnable {
buffer_writer(scheduler *sch, pipebuf<T> &_in, T *_data, int _count)
: runnable(sch, "buffer_reader"),
in(_in), data(_data), count(_count), pos(0) {
template<typename T>
struct buffer_writer: runnable
{
buffer_writer(scheduler *sch, pipebuf<T> &_in, T *_data, int _count) :
runnable(sch, "buffer_reader"), in(_in), data(_data), count(_count), pos(0)
{
}
void run() {
int n = min(in.readable(), (unsigned long)(count-pos));
memcpy(&data[pos], in.rd(), n*sizeof(T));
in.read(n);
pos += n;
void run()
{
int n = min(in.readable(), (unsigned long) (count - pos));
memcpy(&data[pos], in.rd(), n * sizeof(T));
in.read(n);
pos += n;
}
private:
private:
pipereader<T> in;
T *data;
int count;
int pos;
}; // buffer_writer
};
// buffer_writer
} // namespace
}// namespace
#endif // LEANSDR_GENERIC_H

445
plugins/channelrx/demoddatv/leansdr/rs.h
View File

@ -5,88 +5,118 @@
#define DEBUG_RS 0
namespace leansdr {
namespace leansdr
{
// Finite group GF(2^N).
// Finite group GF(2^N).
// GF(2) is the ring ({0,1},+,*).
// GF(2)[X] is the ring of polynomials with coefficients in GF(2).
// P(X) is an irreducible polynomial of GF(2)[X].
// N is the degree of P(x).
// P is the bitfield representation of P(X), with degree 0 at LSB.
// GF(2)[X]/(P) is GF(2)[X] modulo P(X).
// (GF(2)[X]/(P), +) is a group with 2^N elements.
// (GF(2)[X]/(P)*, *) is a group with 2^N-1 elements.
// (GF(2)[X]/(P), +, *) is a field with 2^N elements, noted GF(2^N).
// Te is a C++ integer type for representing elements of GF(2^N).
// "0" is 0
// "1" is 1
// "2" is X
// "3" is X+1
// "4" is X^2
// Tp is a C++ integer type for representing P(X) (1 bit larger than Te).
// ALPHA is a primitive element of GF(2^N). Usually "2"=[X] is chosen.
// GF(2) is the ring ({0,1},+,*).
// GF(2)[X] is the ring of polynomials with coefficients in GF(2).
// P(X) is an irreducible polynomial of GF(2)[X].
// N is the degree of P(x).
// P is the bitfield representation of P(X), with degree 0 at LSB.
// GF(2)[X]/(P) is GF(2)[X] modulo P(X).
// (GF(2)[X]/(P), +) is a group with 2^N elements.
// (GF(2)[X]/(P)*, *) is a group with 2^N-1 elements.
// (GF(2)[X]/(P), +, *) is a field with 2^N elements, noted GF(2^N).
// Te is a C++ integer type for representing elements of GF(2^N).
// "0" is 0
// "1" is 1
// "2" is X
// "3" is X+1
// "4" is X^2
// Tp is a C++ integer type for representing P(X) (1 bit larger than Te).
// ALPHA is a primitive element of GF(2^N). Usually "2"=[X] is chosen.
template<typename Te, typename Tp, Tp P, int N, Te ALPHA>
struct gf2x_p {
gf2x_p() {
if ( ALPHA != 2 ) fail("alpha!=2 not implemented");
// Precompute log and exp tables.
Tp alpha_i = 1;
for ( Tp i=0; i<(1<<N); ++i ) {
lut_exp[i] = alpha_i;
lut_exp[((1<<N)-1)+i] = alpha_i;
lut_log[alpha_i] = i;
alpha_i <<= 1; // Multiply by alpha=[X] i.e. increase degrees
if ( alpha_i & (1<<N) ) alpha_i ^= P; // Modulo P iteratively
}
template<typename Te, typename Tp, Tp P, int N, Te ALPHA>
struct gf2x_p
{
gf2x_p()
{
if (ALPHA != 2)
{
fail("gf2x_p::gf2x_p", "alpha!=2 not implemented");
return;
}
// Precompute log and exp tables.
Tp alpha_i = 1;
for (Tp i = 0; i < (1 << N); ++i)
{
lut_exp[i] = alpha_i;
lut_exp[((1 << N) - 1) + i] = alpha_i;
lut_log[alpha_i] = i;
alpha_i <<= 1; // Multiply by alpha=[X] i.e. increase degrees
if (alpha_i & (1 << N))
alpha_i ^= P; // Modulo P iteratively
}
}
static const Te alpha = ALPHA;
inline Te add(Te x, Te y) { return x ^ y; } // Addition modulo 2
inline Te sub(Te x, Te y) { return x ^ y; } // Subtraction modulo 2
inline Te mul(Te x, Te y) {
if ( !x || !y ) return 0;
return lut_exp[lut_log[x] + lut_log[y]];
inline Te add(Te x, Te y)
{
return x ^ y;
} // Addition modulo 2
inline Te sub(Te x, Te y)
{
return x ^ y;
} // Subtraction modulo 2
inline Te mul(Te x, Te y)
{
if (!x || !y)
return 0;
return lut_exp[lut_log[x] + lut_log[y]];
}
inline Te div(Te x, Te y) {
//if ( ! y ) fail("div"); // TODO
if ( ! x ) return 0;
return lut_exp[lut_log[x] + ((1<<N)-1) - lut_log[y]];
inline Te div(Te x, Te y)
{
//if ( ! y ) fail("div"); // TODO
if (!x)
return 0;
return lut_exp[lut_log[x] + ((1 << N) - 1) - lut_log[y]];
}
inline Te inv(Te x) {
// if ( ! x ) fail("inv");
return lut_exp[((1<<N)-1) - lut_log[x]];
inline Te inv(Te x)
{
// if ( ! x ) fail("inv");
return lut_exp[((1 << N) - 1) - lut_log[x]];
}
inline Te exp(Te x) { return lut_exp[x]; }
inline Te log(Te x) { return lut_log[x]; }
private:
Te lut_exp[(1<<N)*2]; // Wrap to avoid indexing modulo 2^N-1
Te lut_log[1<<N];
};
inline Te exp(Te x)
{
return lut_exp[x];
}
inline Te log(Te x)
{
return lut_log[x];
}
private:
Te lut_exp[(1 << N) * 2]; // Wrap to avoid indexing modulo 2^N-1
Te lut_log[1 << N];
};
// Reed-Solomon for RS(204,188) shortened from RS(255,239).
// Reed-Solomon for RS(204,188) shortened from RS(255,239).
struct rs_engine {
struct rs_engine
{
// EN 300 421, section 4.4.2, Field Generator Polynomial
// p(X) = X^8 + X^4 + X^3 + X^2 + 1
gf2x_p<unsigned char, unsigned short, 0x11d, 8, 2> gf;
u8 G[17]; // { G_16, ..., G_0 }
rs_engine() {
// EN 300 421, section 4.4.2, Code Generator Polynomial
// G(X) = (X-alpha^0)*...*(X-alpha^15)
for ( int i=0; i<=16; ++i ) G[i] = (i==16) ? 1 : 0; // Init G=1
for ( int d=0; d<16; ++d ) {
// Multiply by (X-alpha^d)
// G := X*G - alpha^d*G
for ( int i=0; i<=16; ++i )
G[i] = gf.sub((i==16)?0:G[i+1], gf.mul(gf.exp(d),G[i]));
}
rs_engine()
{
// EN 300 421, section 4.4.2, Code Generator Polynomial
// G(X) = (X-alpha^0)*...*(X-alpha^15)
for (int i = 0; i <= 16; ++i)
G[i] = (i == 16) ? 1 : 0; // Init G=1
for (int d = 0; d < 16; ++d)
{
// Multiply by (X-alpha^d)
// G := X*G - alpha^d*G
for (int i = 0; i <= 16; ++i)
G[i] = gf.sub((i == 16) ? 0 : G[i + 1], gf.mul(gf.exp(d), G[i]));
}
#if DEBUG_RS
fprintf(stderr, "RS generator:");
for ( int i=0; i<=16; ++i ) fprintf(stderr, " %02x", G[i]);
fprintf(stderr, "\n");
fprintf(stderr, "RS generator:");
for ( int i=0; i<=16; ++i ) fprintf(stderr, " %02x", G[i]);
fprintf(stderr, "\n");
#endif
}
@ -96,163 +126,190 @@ namespace leansdr {
// By convention coefficients are listed by decreasing degree here,
// so we can evaluate syndromes of the shortened code without
// prepending with 51 zeroes.
bool syndromes(const u8 *poly, u8 *synd) {
bool corrupted = false;
for ( int i=0; i<16; ++i ) {
synd[i] = eval_poly_rev(poly, 204, gf.exp(i));
if ( synd[i] ) corrupted = true;
}
return corrupted;
bool syndromes(const u8 *poly, u8 *synd)
{
bool corrupted = false;
for (int i = 0; i < 16; ++i)
{
synd[i] = eval_poly_rev(poly, 204, gf.exp(i));
if (synd[i])
corrupted = true;
}
return corrupted;
}
u8 eval_poly_rev(const u8 *poly, int n, u8 x) {
// poly[0]*x^(n-1) + .. + poly[n-1]*x^0 with Hörner method.
u8 acc = 0;
for ( int i=0; i<n; ++i ) acc = gf.add(gf.mul(acc,x), poly[i]);
return acc;
u8 eval_poly_rev(const u8 *poly, int n, u8 x)
{
// poly[0]*x^(n-1) + .. + poly[n-1]*x^0 with Hörner method.
u8 acc = 0;
for (int i = 0; i < n; ++i)
acc = gf.add(gf.mul(acc, x), poly[i]);
return acc;
}
// Evaluation with coefficients listed by increasing degree.
u8 eval_poly(const u8 *poly, int deg, u8 x) {
// poly[0]*x^0 + .. + poly[deg]*x^deg with Hörner method.
u8 acc = 0;
for ( ; deg>=0; --deg ) acc = gf.add(gf.mul(acc,x), poly[deg]);
return acc;
u8 eval_poly(const u8 *poly, int deg, u8 x)
{
// poly[0]*x^0 + .. + poly[deg]*x^deg with Hörner method.
u8 acc = 0;
for (; deg >= 0; --deg)
acc = gf.add(gf.mul(acc, x), poly[deg]);
return acc;
}
// Append parity symbols
void encode(u8 msg[204]) {
// TBD Avoid copying
u8 p[204];
memcpy(p, msg, 188);
memset(p+188, 0, 16);
// p = msg*X^16
void encode(u8 msg[204])
{
// TBD Avoid copying
u8 p[204];
memcpy(p, msg, 188);
memset(p + 188, 0, 16);
// p = msg*X^16
#if DEBUG_RS
fprintf(stderr, "uncoded:");
for ( int i=0; i<204; ++i ) fprintf(stderr, " %d", p[i]);
fprintf(stderr, "\n");
fprintf(stderr, "uncoded:");
for ( int i=0; i<204; ++i ) fprintf(stderr, " %d", p[i]);
fprintf(stderr, "\n");
#endif
// Compute remainder modulo G
for ( int d=0; d<188; ++d ) {
// Clear monomial of degree d
if ( ! p[d] ) continue;
u8 k = gf.div(p[d], G[0]);
// p(X) := p(X) - k*G(X)*X^(188-d)
for ( int i=0; i<=16; ++i )
p[d+i] = gf.sub(p[d+i], gf.mul(k,G[i]));
}
// Compute remainder modulo G
for (int d = 0; d < 188; ++d)
{
// Clear monomial of degree d
if (!p[d])
continue;
u8 k = gf.div(p[d], G[0]);
// p(X) := p(X) - k*G(X)*X^(188-d)
for (int i = 0; i <= 16; ++i)
p[d + i] = gf.sub(p[d + i], gf.mul(k, G[i]));
}
#if DEBUG_RS
fprintf(stderr, "coded:");
for ( int i=0; i<204; ++i ) fprintf(stderr, " %d", p[i]);
fprintf(stderr, "\n");
fprintf(stderr, "coded:");
for ( int i=0; i<204; ++i ) fprintf(stderr, " %d", p[i]);
fprintf(stderr, "\n");
#endif
memcpy(msg+188, p+188, 16);
memcpy(msg + 188, p + 188, 16);
}
// Try to fix errors in pout[].
// If pin[] is provided, errors will be fixed in the original
// message too and syndromes will be updated.
bool correct(u8 synd[16], u8 pout[188],
u8 pin[204]=NULL, int *bits_corrected=NULL) {
// Berlekamp - Massey
// http://en.wikipedia.org/wiki/Berlekamp%E2%80%93Massey_algorithm#Code_sample
u8 C[16] = { 1,0,0,0, 0,0,0,0, 0,0,0,0, 0,0,0,0 }; // Max degree is L
u8 B[16] = { 1,0,0,0, 0,0,0,0, 0,0,0,0, 0,0,0,0 };
int L = 0;
int m = 1;
u8 b = 1;
for ( int n=0; n<16; ++n ) {
u8 d = synd[n];
for ( int i=1; i<=L; ++i ) d ^= gf.mul(C[i], synd[n-i]);
if ( ! d ) {
++m;
} else if ( 2*L <= n ) {
u8 T[16];
memcpy(T, C, sizeof(T));
for ( int i=0; i<16-m; ++i )
C[m+i] ^= gf.mul(d, gf.mul(gf.inv(b),B[i]));
L = n + 1 - L;
memcpy(B, T, sizeof(B));
b = d;
m = 1;
} else {
for ( int i=0; i<16-m; ++i )
C[m+i] ^= gf.mul(d, gf.mul(gf.inv(b),B[i]));
++m;
}
}
// L is the number of errors
// C of degree L is now the error locator polynomial (Lambda)
bool correct(u8 synd[16], u8 pout[188], u8 pin[204] = NULL, int *bits_corrected = NULL)
{
// Berlekamp - Massey
// http://en.wikipedia.org/wiki/Berlekamp%E2%80%93Massey_algorithm#Code_sample
u8 C[16] =
{ 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }; // Max degree is L
u8 B[16] =
{ 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
int L = 0;
int m = 1;
u8 b = 1;
for (int n = 0; n < 16; ++n)
{
u8 d = synd[n];
for (int i = 1; i <= L; ++i)
d ^= gf.mul(C[i], synd[n - i]);
if (!d)
{
++m;
}
else if (2 * L <= n)
{
u8 T[16];
memcpy(T, C, sizeof(T));
for (int i = 0; i < 16 - m; ++i)
C[m + i] ^= gf.mul(d, gf.mul(gf.inv(b), B[i]));
L = n + 1 - L;
memcpy(B, T, sizeof(B));
b = d;
m = 1;
}
else
{
for (int i = 0; i < 16 - m; ++i)
C[m + i] ^= gf.mul(d, gf.mul(gf.inv(b), B[i]));
++m;
}
}
// L is the number of errors
// C of degree L is now the error locator polynomial (Lambda)
#if DEBUG_RS
fprintf(stderr, "[L=%d C=",L);
for ( int i=0; i<16; ++i ) fprintf(stderr, " %d", C[i]);
fprintf(stderr, "]\n");
fprintf(stderr, "[S=");
for ( int i=0; i<16; ++i ) fprintf(stderr, " %d", synd[i]);
fprintf(stderr, "]\n");
fprintf(stderr, "[L=%d C=",L);
for ( int i=0; i<16; ++i ) fprintf(stderr, " %d", C[i]);
fprintf(stderr, "]\n");
fprintf(stderr, "[S=");
for ( int i=0; i<16; ++i ) fprintf(stderr, " %d", synd[i]);
fprintf(stderr, "]\n");
#endif
// Forney
// http://en.wikipedia.org/wiki/Forney_algorithm (2t=16)
// Compute Omega
u8 omega[16];
memset(omega, 0, sizeof(omega));
// TODO loops
for ( int i=0; i<16; ++i )
for ( int j=0; j<16; ++j )
if ( i+j < 16 ) omega[i+j] ^= gf.mul(synd[i], C[j]);
// Forney
// http://en.wikipedia.org/wiki/Forney_algorithm (2t=16)
// Compute Omega
u8 omega[16];
memset(omega, 0, sizeof(omega));
// TODO loops
for (int i = 0; i < 16; ++i)
for (int j = 0; j < 16; ++j)
if (i + j < 16)
omega[i + j] ^= gf.mul(synd[i], C[j]);
#if DEBUG_RS
fprintf(stderr, "omega=");
for ( int i=0; i<16; ++i ) fprintf(stderr, " %d", omega[i]);
fprintf(stderr, "\n");
fprintf(stderr, "omega=");
for ( int i=0; i<16; ++i ) fprintf(stderr, " %d", omega[i]);
fprintf(stderr, "\n");
#endif
// Compute Lambda'
u8 Cprime[15];
for ( int i=0; i<15; ++i )
Cprime[i] = (i&1) ? 0 : C[i+1];
// Compute Lambda'
u8 Cprime[15];
for (int i = 0; i < 15; ++i)
Cprime[i] = (i & 1) ? 0 : C[i + 1];
#if DEBUG_RS
fprintf(stderr, "Cprime=");
for ( int i=0; i<15; ++i ) fprintf(stderr, " %d", Cprime[i]);
fprintf(stderr, "\n");
fprintf(stderr, "Cprime=");
for ( int i=0; i<15; ++i ) fprintf(stderr, " %d", Cprime[i]);
fprintf(stderr, "\n");
#endif
// Find zeroes of C by exhaustive search?
// TODO Chien method
int roots_found = 0;
for ( int i=0; i<255; ++i ) {
u8 r = gf.exp(i); // Candidate root alpha^0..alpha^254
u8 v = eval_poly(C, L, r);
if ( ! v ) {
// r is a root X_k^-1 of the error locator polynomial.
u8 xk = gf.inv(r);
int loc = (255-i) % 255; // == log(xk)
// Find zeroes of C by exhaustive search?
// TODO Chien method
int roots_found = 0;
for (int i = 0; i < 255; ++i)
{
u8 r = gf.exp(i); // Candidate root alpha^0..alpha^254
u8 v = eval_poly(C, L, r);
if (!v)
{
// r is a root X_k^-1 of the error locator polynomial.
u8 xk = gf.inv(r);
int loc = (255 - i) % 255; // == log(xk)
#if DEBUG_RS
fprintf(stderr, "found root=%d, inv=%d, loc=%d\n", r, xk, loc);
fprintf(stderr, "found root=%d, inv=%d, loc=%d\n", r, xk, loc);
#endif
if ( loc < 204 ) {
// Evaluate e_k
u8 num = gf.mul(xk, eval_poly(omega, L, r));
u8 den = eval_poly(Cprime, 14, r);
u8 e = gf.div(num, den);
// Subtract e from coefficient of degree loc.
// Note: Coeffients listed by decreasing degree in pin[] and pout[].
if ( bits_corrected ) *bits_corrected += hamming_weight(e);
if ( loc >= 16 ) pout[203-loc] ^= e;
if ( pin ) pin[203-loc] ^= e;
}
if ( ++roots_found == L ) break;
}
}
if ( pin )
return syndromes(pin, synd);
else
return false;
if (loc < 204)
{
// Evaluate e_k
u8 num = gf.mul(xk, eval_poly(omega, L, r));
u8 den = eval_poly(Cprime, 14, r);
u8 e = gf.div(num, den);
// Subtract e from coefficient of degree loc.
// Note: Coeffients listed by decreasing degree in pin[] and pout[].
if (bits_corrected)
*bits_corrected += hamming_weight(e);
if (loc >= 16)
pout[203 - loc] ^= e;
if (pin)
pin[203 - loc] ^= e;
}
if (++roots_found == L)
break;
}
}
if (pin)
return syndromes(pin, synd);
else
return false;
}
};
};
} // namespace

23
plugins/channelrx/demoddatv/leansdr/sdr.h
View File

@ -531,7 +531,8 @@ struct cstln_lut
make_qam(256);
break;
default:
fail("Constellation not implemented");
fail("cstln_lut::cstln_lut", "Constellation not implemented");
return;
}
}
struct result
@ -936,7 +937,10 @@ struct cstln_receiver: runnable
void run()
{
if (!cstln)
fail("constellation not set");
{
fail("cstln_lut::run", "constellation not set");
return;
}
// Magic constants that work with the qa recordings.
float freq_alpha = 0.04;
@ -1190,7 +1194,10 @@ struct fast_qpsk_receiver: runnable
signed long freq_alpha = 0.04 * 65536;
signed long freq_beta = 0.0012 * 256 * 65536 / omega * pll_adjustment;
if (!freq_beta)
fail("Excessive oversampling");
{
fail("fast_qpsk_receiver::run", "Excessive oversampling");
return;
}
float gain_mu = 0.02 / (cstln_amp * cstln_amp) * 2;
@ -1431,7 +1438,10 @@ struct cstln_transmitter: runnable
void run()
{
if (!cstln)
fail("constellation not set");
{
fail("cstln_transmitter::run", "constellation not set");
return;
}
int count = min(in.readable(), out.writable());
u8 *pin = in.rd(), *pend = pin + count;
complex<Tout> *pout = out.wr();
@ -1524,8 +1534,9 @@ struct cnr_fft: runnable
avgpower(NULL),
phase(0)
{
if (bandwidth > 0.25)
fail("CNR estimator requires Fsampling > 4x Fsignal");
if (bandwidth > 0.25) {
fail("cnr_fft::cnr_fft", "CNR estimator requires Fsampling > 4x Fsignal");
}
}
float bandwidth;

479
plugins/channelrx/demoddatv/leansdr/viterbi.h
View File

@ -11,271 +11,322 @@
// TBD This is very inefficient. For a specific trellis all loops
// can be be unrolled.
namespace leansdr {
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.
// 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;
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
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;
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);
void init_convolutional(const uint16_t G[])
{
if (NCS & (NCS - 1))
{
fprintf(stderr, "leansdr::trellis::init_convolutional: NCS must be a power of 2\n");
return;
}
// 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;
}
}
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, "leansdr::trellis::init_convolutional: Invalid convolutional code\n");
return;
}
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");
}
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;
virtual TUS update(int nm, TCS cs[], TBM costs[], TPM *quality=NULL)=0;
};
// 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;
virtual TUS update(int nm, TCS cs[], TBM costs[], 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> {
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
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)
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 ) ;
}
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;
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");
// 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();
// 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;
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");
// 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();
// 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);
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");
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:
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.
// 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 {
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); }
};
bitpath() :
val(0)
{
}
void append(TUS us)
{
val = (val << NBITS) | us;
}
TUS read()
{
return (val >> (DEPTH - 1) * NBITS) & ((1 << NBITS) - 1);
}
};
} // namespace