Support non-standard constellations with DVB-S --viterbi.

master
pabr 2017-12-09 20:58:48 +01:00
rodzic 8f8ee3c02d
commit 838942f195
9 zmienionych plików z 655 dodań i 258 usunięć

Wyświetl plik

@ -1,3 +1,9 @@
HEAD
* leandvb, leandvbtx: DVB-S with non-standard constellations (with --viterbi).
* leandvb: Support all DVB-S code rates with --viterbi.
* leandvbtx: Support all DVB-S code rates.
* New viterbi_sync with simplified metrics (may degrade QPSK sensitivity).
2017-07-23 v1.2.0
* Fixed --hs mode on x86_64.
* Added --s8, --s16, --u16.

Wyświetl plik

@ -417,15 +417,14 @@ int run(config &cfg) {
if ( cfg.standard == config::DVB_S ) {
if ( cfg.constellation != cstln_lut<256>::QPSK &&
cfg.constellation != cstln_lut<256>::BPSK )
fail("Invalid constellation for DVB-S");
demod.cstln = new cstln_lut<256>(cfg.constellation);
fprintf(stderr, "Warning: non-standard constellation for DVB-S\n");
}
if ( cfg.standard == config::DVB_S2 ) {
// For DVB-S2 testing only.
// Constellation should be determined from PL signalling.
fprintf(stderr, "DVB-S2: Testing symbol sampler only.\n");
demod.cstln = make_dvbs2_constellation(cfg.constellation, cfg.fec);
}
demod.cstln = make_dvbs2_constellation(cfg.constellation, cfg.fec);
if ( cfg.hard_metric ) {
if ( cfg.verbose )
fprintf(stderr, "Using hard metric.\n");
@ -482,15 +481,14 @@ int run(config &cfg) {
deconvol_sync_simple *r_deconv = NULL;
if ( cfg.viterbi ) {
if ( cfg.fec != FEC12 )
fail("Viterbi only for code rate 1/2");
if ( cfg.constellation == cstln_lut<256>::QPSK ) {
viterbi_sync *r = new viterbi_sync(&sch, p_symbols, p_bytes);
if ( cfg.fastlock ) r->resync_period = 1;
} else {
viterbi_sync_bpsk *r = new viterbi_sync_bpsk(&sch, p_symbols, p_bytes);
if ( cfg.fastlock ) r->resync_period = 1;
}
if ( cfg.fec==FEC23 && (demod.cstln->nsymbols==4 ||
demod.cstln->nsymbols==64) ) {
if ( cfg.verbose ) fprintf(stderr, "Handling rate 2/3 as 4/6\n");
cfg.fec = FEC46;
}
viterbi_sync *r = new viterbi_sync(&sch, p_symbols, p_bytes,
demod.cstln, cfg.fec);
if ( cfg.fastlock ) r->resync_period = 1;
} else {
r_deconv = make_deconvol_sync_simple(&sch, p_symbols, p_bytes, cfg.fec);
r_deconv->fastlock = cfg.fastlock;
@ -921,7 +919,10 @@ void usage(const char *name, FILE *f, int c) {
" --tune HZ Bias frequency for demodulation\n"
" --drift Track frequency drift beyond safe limits\n"
" --standard S DVB-S (default), DVB-S2 (not implemented)\n"
" --const C QPSK (default), BPSK .. 32APSK (DVB-S2 only)\n"
" --const C QPSK (default),\n"
" BPSK .. 32APSK (DVB-S2),\n"
" 64APSKe (DVB-S2X),\n"
" 16QAM .. 256QAM (experimental)\n"
" --cr N/D Code rate 1/2 (default) .. 7/8 .. 9/10\n"
" --fastlock Synchronize more aggressively (CPU-intensive)\n"
" --sampler nearest, linear, rrc\n"
@ -1000,6 +1001,14 @@ int main(int argc, const char *argv[]) {
cfg.constellation = cstln_lut<256>::APSK16;
else if ( ! strcmp(argv[i], "32APSK" ) )
cfg.constellation = cstln_lut<256>::APSK32;
else if ( ! strcmp(argv[i], "64APSKe" ) )
cfg.constellation = cstln_lut<256>::APSK64E;
else if ( ! strcmp(argv[i], "16QAM" ) )
cfg.constellation = cstln_lut<256>::QAM16;
else if ( ! strcmp(argv[i], "64QAM" ) )
cfg.constellation = cstln_lut<256>::QAM64;
else if ( ! strcmp(argv[i], "256QAM" ) )
cfg.constellation = cstln_lut<256>::QAM256;
else usage(argv[0], stderr, 1);
}
else if ( ! strcmp(argv[i], "--cr") && i+1<argc ) {

Wyświetl plik

@ -41,6 +41,8 @@ private:
};
struct config {
cstln_lut<256>::predef constellation;
code_rate fec;
float power;
bool agc;
int interp;
@ -48,7 +50,8 @@ struct config {
float rolloff;
bool verbose, debug;
config()
: power(1.0), agc(false),
: constellation(cstln_lut<256>::QPSK),
power(1.0), agc(false),
interp(2), decim(1), rolloff(0.35),
verbose(false), debug(false)
{ }
@ -59,9 +62,10 @@ void run(config &cfg) {
sch.verbose = cfg.verbose;
sch.debug = cfg.debug;
unsigned long BUF_PACKETS = 12; // TBD Reduce copying
unsigned long BUF_BYTES = SIZE_RSPACKET*12;
unsigned long BUF_SYMBOLS = BUF_BYTES*8;
int buf_factor = 2;
unsigned long BUF_PACKETS = 12*buf_factor; // TBD Reduce copying
unsigned long BUF_BYTES = SIZE_RSPACKET*BUF_PACKETS;
unsigned long BUF_SYMBOLS = BUF_BYTES*8 * 2; // Worst case BPSK 1/2
unsigned long BUF_BASEBAND = 4096;
// TS PACKETS ON STDIN
@ -86,15 +90,22 @@ void run(config &cfg) {
// CONVOLUTIONAL CODER
cstln_lut<256> *cstln = make_dvbs2_constellation(cfg.constellation, cfg.fec);
int bits_per_symbol = log2(cstln->nsymbols);
if ( cfg.fec==FEC23 && (cstln->nsymbols==4 ||
cstln->nsymbols==64) ) {
if ( cfg.verbose ) fprintf(stderr, "Handling rate 2/3 as 4/6\n");
cfg.fec = FEC46;
}
pipebuf<u8> p_symbols(&sch, "symbols", BUF_SYMBOLS);
dvb_convol r_convol(&sch, p_mpegbytes, p_symbols);
dvb_convol r_convol(&sch, p_mpegbytes, p_symbols, cfg.fec, bits_per_symbol);
// IQ MAPPER
pipebuf<cf32> p_iqsymbols(&sch, "IQ symbols", BUF_SYMBOLS);
cstln_transmitter<f32,0> r_mod(&sch, p_symbols, p_iqsymbols);
cstln_lut<256> qpsk(cstln_lut<256>::QPSK);
r_mod.cstln = &qpsk;
r_mod.cstln = cstln;
// RESAMPLER
@ -152,6 +163,11 @@ void usage(const char *name, FILE *f, int c) {
fprintf(f, "Output float complex samples\n");
fprintf
(f, "\nOptions:"
" --const STRING QPSK (default),\n"
" BPSK .. 32APSK (DVB-S2),\n"
" 64APSKe (DVB-S2X),\n"
" 16QAM .. 256QAM (experimental)\n"
" --cr STRING 1/2, 2/3, 3/4, 5/6, 7/8\n"
" -f INTERP[/DECIM] Samples per symbols (default: 2)\n"
" --roll-off R RRC roll-off (defalt: 0.35)\n"
" --power P Output power (dB)\n"
@ -171,6 +187,42 @@ int main(int argc, char *argv[]) {
cfg.verbose = true;
else if ( ! strcmp(argv[i], "-d") )
cfg.debug = true;
else if ( ! strcmp(argv[i], "--cr") && i+1<argc ) {
++i;
// DVB-S
if ( ! strcmp(argv[i], "1/2" ) ) cfg.fec = FEC12;
else if ( ! strcmp(argv[i], "2/3" ) ) cfg.fec = FEC23;
else if ( ! strcmp(argv[i], "3/4" ) ) cfg.fec = FEC34;
else if ( ! strcmp(argv[i], "5/6" ) ) cfg.fec = FEC56;
else if ( ! strcmp(argv[i], "7/8" ) ) cfg.fec = FEC78;
// DVB-S2
else if ( ! strcmp(argv[i], "4/5" ) ) cfg.fec = FEC45;
else if ( ! strcmp(argv[i], "8/9" ) ) cfg.fec = FEC89;
else if ( ! strcmp(argv[i], "9/10" ) ) cfg.fec = FEC910;
else usage(argv[0], stderr, 1);
}
else if ( ! strcmp(argv[i], "--const") && i+1<argc ) {
++i;
if ( ! strcmp(argv[i], "BPSK" ) )
cfg.constellation = cstln_lut<256>::BPSK;
else if ( ! strcmp(argv[i], "QPSK" ) )
cfg.constellation = cstln_lut<256>::QPSK;
else if ( ! strcmp(argv[i], "8PSK" ) )
cfg.constellation = cstln_lut<256>::PSK8;
else if ( ! strcmp(argv[i], "16APSK" ) )
cfg.constellation = cstln_lut<256>::APSK16;
else if ( ! strcmp(argv[i], "32APSK" ) )
cfg.constellation = cstln_lut<256>::APSK32;
else if ( ! strcmp(argv[i], "64APSKe" ) )
cfg.constellation = cstln_lut<256>::APSK64E;
else if ( ! strcmp(argv[i], "16QAM" ) )
cfg.constellation = cstln_lut<256>::QAM16;
else if ( ! strcmp(argv[i], "64QAM" ) )
cfg.constellation = cstln_lut<256>::QAM64;
else if ( ! strcmp(argv[i], "256QAM" ) )
cfg.constellation = cstln_lut<256>::QAM256;
else usage(argv[0], stderr, 1);
}
else if ( ! strcmp(argv[i], "-f") && i+1<argc ) {
++i;
cfg.decim = 1;

Wyświetl plik

@ -203,6 +203,55 @@ namespace leansdr {
Thist hist;
}; // convol_poly2
// Generic BPSK..256QAM and puncturing
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");
}
private:
Thist hist;
int nhist;
Thist sersymb;
int nsersymb;
}; // convol_multipoly
} // namespace
#endif // LEANSDR_CONVOLUTIONAL_H

Wyświetl plik

@ -94,7 +94,7 @@ namespace leansdr {
private:
int logn;
int *bitrev;
complex<float> *omega, *omega_rev;
complex<T> *omega, *omega_rev;
float invsqrtn;
};

Wyświetl plik

@ -18,20 +18,22 @@ namespace leansdr {
// Generic deconvolution
enum code_rate {
FEC12, FEC23, FEC34, FEC56, FEC78, // DVB-S
FEC45, FEC89, FEC910, // DVB-S2
FEC12, FEC23, FEC46, FEC34, FEC56, FEC78, // DVB-S
FEC45, FEC89, FEC910, // DVB-S2
FEC_MAX
};
// Customize APSK radii according to code rate
cstln_lut<256> *make_dvbs2_constellation(cstln_lut<256>::predef c,
code_rate r) {
float gamma1=1, gamma2=1;
cstln_lut<256> * make_dvbs2_constellation(cstln_lut<256>::predef c,
code_rate r) {
float gamma1=1, gamma2=1, gamma3=1;
switch ( c ) {
case cstln_lut<256>::APSK16:
// EN 302 307, section 5.4.3, Table 9
switch ( r ) {
case FEC23: gamma1 = 3.15; break;
case FEC23:
case FEC46: gamma1 = 3.15; break;
case FEC34: gamma1 = 2.85; break;
case FEC45: gamma1 = 2.75; break;
case FEC56: gamma1 = 2.70; break;
@ -51,10 +53,14 @@ namespace leansdr {
default: fail("Code rate not supported with APSK32");
}
break;
case cstln_lut<256>::APSK64E:
// EN 302 307-2, section 5.4.5, Table 13f
gamma1 = 2.4; gamma2 = 4.3; gamma3 = 7;
break;
default:
break;
}
return new cstln_lut<256>(c, gamma1, gamma2);
return new cstln_lut<256>(c, gamma1, gamma2, gamma3);
}
// EN 300 421, section 4.4.3, table 2 Punctured code, G1=0171, G2=0133
@ -160,7 +166,7 @@ namespace leansdr {
}
void next_sync() {
if ( fastlock ) fatal("Bug: next_sync() called with fastlock");
if ( fastlock ) fail("Bug: next_sync() called with fastlock");
++locked;
if ( locked == &syncs[NSYNCS] ) {
locked = &syncs[0];
@ -466,6 +472,7 @@ namespace leansdr {
pY = 0x1; // 1
break;
case FEC23:
case FEC46:
pX = 0xa; // 1010 (Handle as FEC4/6, no half-symbols)
pY = 0xf; // 1111
break;
@ -493,27 +500,84 @@ namespace leansdr {
// CONVOLUTIONAL ENCODER
static const uint16_t polys_fec12[] = {
DVBS_G1, DVBS_G2 // X1Y1
};
static const uint16_t polys_fec23[] = {
DVBS_G1, DVBS_G2, DVBS_G2<<1 // X1Y1Y2
};
// Same code rate as 2/3, usable with QPSK
static const uint16_t polys_fec46[] = {
DVBS_G1, DVBS_G2, DVBS_G2<<1, // X1Y1Y2
DVBS_G1<<2, DVBS_G2<<2, DVBS_G2<<3 // X3Y3Y4
};
static const uint16_t polys_fec34[] = {
DVBS_G1, DVBS_G2, // X1Y1
DVBS_G2<<1, DVBS_G1<<2 // Y2X3
};
static const uint16_t polys_fec56[] = {
DVBS_G1, DVBS_G2, // X1Y1
DVBS_G2<<1, DVBS_G1<<2, // Y2X3
DVBS_G2<<3, DVBS_G1<<4 // Y4X5
};
static const uint16_t polys_fec78[] = {
DVBS_G1, DVBS_G2, // X1Y1
DVBS_G2<<1, DVBS_G2<<2, // Y2Y3
DVBS_G2<<3, DVBS_G1<<4, // Y4X5
DVBS_G2<<5, DVBS_G1<<6 // Y6X7
};
// FEC parameters, for convolutional coding only (not S2).
static struct fec_spec {
int bits_in; // Entering the convolutional coder
int bits_out; // Exiting the convolutional coder
const uint16_t *polys; // [bits_out]
} fec_specs[FEC_MAX] = {
[FEC12] = { 1, 2, polys_fec12 },
[FEC23] = { 2, 3, polys_fec23 },
[FEC46] = { 4, 6, polys_fec46 },
[FEC34] = { 3, 4, polys_fec34 },
[FEC56] = { 5, 6, polys_fec56 },
[FEC78] = { 7, 8, polys_fec78 },
};
struct dvb_convol : runnable {
typedef u8 uncoded_byte;
typedef u8 hardsymbol;
dvb_convol(scheduler *sch,
pipebuf<uncoded_byte> &_in,
pipebuf<hardsymbol> &_out)
pipebuf<hardsymbol> &_out,
code_rate fec,
int bits_per_symbol)
: runnable(sch, "dvb_convol"),
in(_in), out(_out,8)
in(_in), out(_out,64) // BPSK 7/8: 7 bytes in, 64 symbols out
{
fec_spec *fs = &fec_specs[fec];
if ( ! fs->bits_in ) fail("Unexpected FEC");
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");
}
void run() {
int count = min(in.readable(), out.writable()/8);
static const u8 remap[] = { 0, 1, 2, 3 };
convol.run(in.rd(), remap, out.wr(), count);
int count = min(in.readable(), out.writable()*convol.bps/
convol.bits_out*convol.bits_in/8);
// Process in multiples of the puncturing period and of 8 bits.
int chunk = convol.bits_in;
count = (count/chunk) * chunk;
convol.encode(in.rd(), out.wr(), count);
in.read(count);
out.written(count*8);
int nout = count*8/convol.bits_in*convol.bits_out/convol.bps;
out.written(nout);
}
private:
pipereader<uncoded_byte> in;
pipewriter<hardsymbol> out;
convol_poly2<uint32_t, 0171, 0133> convol;
convol_multipoly<uint16_t, 16> convol;
}; // dvb_convol
@ -1077,22 +1141,61 @@ namespace leansdr {
// VITERBI DECODING
// QPSK 1/2 only.
// Supports all code rates and constellations
// Simplified metric to support large constellations.
// This version implements puncturing by expanding the trellis.
// TBD Compare performance vs skipping updates in a 1/2 trellis.
struct viterbi_sync : runnable {
static const int TRACEBACK = 32; // Suitable for QPSK 1/2
typedef uint8_t TS, TCS, TUS;
typedef uint16_t TBM;
typedef uint32_t TPM;
typedef bitpath<uint32_t,TUS,1,TRACEBACK> dvb_path;
typedef viterbi_dec<TS,64, TUS,2, TCS,4, TBM, TPM, dvb_path> dvb_dec;
typedef trellis<TS,64, TUS,2, 4> dvb_trellis;
typedef int32_t TBM; // Only 16 bits per IQ, but several IQ per Viterbi CS
typedef int32_t TPM;
typedef viterbi_dec_interface<TUS,TCS,TBM,TPM> dvb_dec_interface;
// 1/2: 6 bits of state, 1 bit in, 2 bits out
typedef bitpath<uint32_t,TUS,1,32> path_12;
typedef trellis<TS,64, TUS,2, 4> trellis_12;
typedef viterbi_dec<TS,64, TUS,2, TCS,4, TBM, TPM, path_12> dvb_dec_12;
// 2/3: 6 bits of state, 2 bits in, 3 bits out
typedef bitpath<uint64_t,TUS,3,21> path_23;
typedef trellis<TS,64, TUS,4, 8> trellis_23;
typedef viterbi_dec<TS,64, TUS,4, TCS,8, TBM, TPM, path_23> dvb_dec_23;
// 4/6: 6 bits of state, 4 bits in, 6 bits out
typedef bitpath<uint64_t,TUS,4,16> path_46;
typedef trellis<TS,64, TUS,16, 64> trellis_46;
typedef viterbi_dec<TS,64, TUS,16, TCS,64, TBM, TPM, path_46> dvb_dec_46;
// 3/4: 6 bits of state, 3 bits in, 4 bits out
typedef bitpath<uint64_t,TUS,3,21> path_34;
typedef trellis<TS,64, TUS,8, 16> trellis_34;
typedef viterbi_dec<TS,64, TUS,8, TCS,16, TBM, TPM, path_34> dvb_dec_34;
// 5/6: 6 bits of state, 5 bits in, 6 bits out
typedef bitpath<uint64_t,TUS,5,12> path_56;
typedef trellis<TS,64, TUS,32, 64> trellis_56;
typedef viterbi_dec<TS,64, TUS,32, TCS,64, TBM, TPM, path_56> dvb_dec_56;
// QPSK 7/8: 6 bits of state, 7 bits in, 8 bits out
typedef bitpath<uint64_t,TUS,7,9> path_78;
typedef trellis<TS,64, TUS,128, 256> trellis_78;
typedef viterbi_dec<TS,64, TUS,128, TCS,256, TBM, TPM, path_78> dvb_dec_78;
private:
pipereader<softsymbol> in;
pipewriter<unsigned char> out;
static const int NSYNCS = 4;
dvb_dec *syncs[NSYNCS];
cstln_lut<256> *cstln;
fec_spec *fec;
int bits_per_symbol; // Bits per IQ symbol (not per coded symbol)
int nsyncs;
int nshifts;
struct sync {
int shift;
dvb_dec_interface *dec;
TCS *map; // [nsymbols]
} *syncs; // [nsyncs]
int current_sync;
static const int chunk_size = 128;
int resync_phase;
@ -1100,78 +1203,173 @@ namespace leansdr {
int resync_period;
viterbi_sync(scheduler *sch,
pipebuf<softsymbol> &_in,
pipebuf<unsigned char> &_out)
pipebuf<softsymbol> &_in, pipebuf<unsigned char> &_out,
cstln_lut<256> *_cstln, code_rate cr)
: runnable(sch, "viterbi_sync"),
in(_in), out(_out, chunk_size),
cstln(_cstln),
current_sync(0),
resync_phase(0),
resync_period(32) // 1/32 = 9% synchronization overhead
resync_period(32) // 1/32 = 9% synchronization overhead TBD
{
dvb_trellis *trell = new dvb_trellis();
uint64_t dvb_polynomials[] = { DVBS_G1, DVBS_G2 };
trell->init_convolutional(dvb_polynomials);
for ( int s=0; s<NSYNCS; ++s ) syncs[s] = new dvb_dec(trell);
bits_per_symbol = log2(cstln->nsymbols);
fec = &fec_specs[cr];
{ // 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");
}
int nconj;
switch ( cstln->nsymbols ) {
case 2: nconj = 1; break; // Conjugation is not relevant for BPSK
default: nconj = 2; break;
}
int nrotations;
switch ( cstln->nsymbols ) {
case 2:
case 4:
// For BPSK and QPSK, 180° rotation is handled as
// polarity inversion in mpeg_sync.
nrotations = cstln->nrotations/2;
break;
default:
nrotations = cstln->nrotations;
break;
}
nshifts = fec->bits_out / bits_per_symbol;
nsyncs = nconj * nrotations * nshifts;
// TBD Many HOM constellations are labelled in such a way
// that certain rot/conj combinations are equivalent to
// polarity inversion. We could reduce nsyncs.
syncs = new sync[nsyncs];
for ( int s=0; s<nsyncs; ++s ) {
// Bit pattern [shift|conj|rot]
int rot = s % nrotations;
int conj = (s/nrotations) % nconj;
int shift = s / nrotations / nconj;
syncs[s].shift = shift;
if ( shift ) // Reuse identical map
syncs[s].map = syncs[conj*nrotations+rot].map;
else
syncs[s].map = init_map(conj, 2*M_PI*rot/cstln->nrotations);
#if 0
fprintf(stderr, "sync %3d: conj%d offs%d rot%d/%d map:",
s, conj, syncs[s].shift, rot, cstln->nrotations);
for ( int i=0; i<cstln->nsymbols; ++i )
fprintf(stderr, " %2d", syncs[s].map[i]);
fprintf(stderr, "\n");
#endif
}
if ( cr == FEC12 ) {
trellis_12 *trell = new trellis_12();
trell->init_convolutional(fec->polys);
for ( int s=0; s<nsyncs; ++s ) syncs[s].dec = new dvb_dec_12(trell);
}
else if ( cr == FEC23 ) {
trellis_23 *trell = new trellis_23();
trell->init_convolutional(fec->polys);
for ( int s=0; s<nsyncs; ++s ) syncs[s].dec = new dvb_dec_23(trell);
}
else if ( cr == FEC46 ) {
trellis_46 *trell = new trellis_46();
trell->init_convolutional(fec->polys);
for ( int s=0; s<nsyncs; ++s ) syncs[s].dec = new dvb_dec_46(trell);
}
else if ( cr == FEC34 ) {
trellis_34 *trell = new trellis_34();
trell->init_convolutional(fec->polys);
for ( int s=0; s<nsyncs; ++s ) syncs[s].dec = new dvb_dec_34(trell);
}
else if ( cr == FEC56 ) {
trellis_56 *trell = new trellis_56();
trell->init_convolutional(fec->polys);
for ( int s=0; s<nsyncs; ++s ) syncs[s].dec = new dvb_dec_56(trell);
}
else if ( cr == FEC78 ) {
trellis_78 *trell = new trellis_78();
trell->init_convolutional(fec->polys);
for ( int s=0; s<nsyncs; ++s ) syncs[s].dec = new dvb_dec_78(trell);
}
else
fail("CR not supported");
}
inline TUS update_sync(int s, TBM m[4], TPM *discr) {
// EN 300 421, section 4.5 Baseband shaping and modulation
// EN 302 307, section 5.4.1
//
// IQ=10=(2) | IQ=00=(0)
// ----------+----------
// IQ=11=(3) | IQ=01=(1)
//
TBM vm[4];
switch ( s ) {
case 0: // Mapping for 0°
vm[0]=m[0]; vm[1]=m[1]; vm[2]=m[2]; vm[3]=m[3]; break;
case 1: // Mapping for 90°
vm[0]=m[2]; vm[2]=m[3]; vm[3]=m[1]; vm[1]=m[0]; break;
case 2: // Mapping for 0° conjugated
vm[0]=m[1]; vm[1]=m[0]; vm[2]=m[3]; vm[3]=m[2]; break;
case 3: // Mapping for 90° conjugated
vm[0]=m[3]; vm[2]=m[2]; vm[3]=m[0]; vm[1]=m[1]; break;
default:
return 0; // Avoid compiler warning
TCS *init_map(bool conj, float angle) {
// Each constellation has its own pattern for labels.
// Here we simply tabulate systematically.
TCS *map = new TCS[cstln->nsymbols];
float ca=cosf(angle), sa=sinf(angle);
for ( int i=0; i<cstln->nsymbols; ++i ) {
int8_t I = cstln->symbols[i].re;
int8_t Q = cstln->symbols[i].im;
if ( conj ) Q = -Q;
int8_t RI = I*ca - Q*sa;
int8_t RQ = I*sa + Q*ca;
cstln_lut<256>::result *pr = cstln->lookup(RI, RQ);
map[i] = pr->ss.symbol;
}
return syncs[s]->update(vm, discr);
return map;
}
inline TUS update_sync(int s, softsymbol *pin, TPM *discr) {
// Read one FEC ouput block
pin += syncs[s].shift;
TCS cs = 0;
TBM cost = 0;
for ( int i=0; i<nshifts; ++i,++pin ) {
cs = (cs<<bits_per_symbol) | syncs[s].map[pin->symbol];
cost += pin->cost;
}
return syncs[s].dec->update(cs, cost, discr);
}
void run() {
while ( in.readable()>=8*chunk_size && out.writable()>=chunk_size ) {
unsigned long totaldiscr[NSYNCS];
for ( int s=0; s<NSYNCS; ++s ) totaldiscr[s] = 0;
for ( int bytenum=0; bytenum<chunk_size; ++bytenum ) {
// Decode one byte
unsigned char byte = 0;
softsymbol *pin = in.rd();
// Number of FEC blocks to fill the bitpath depth.
// Before that we cannot discriminate between synchronizers
int discr_delay = 64 / fec->bits_in;
// Process [chunk_size] FEC blocks at a time
while ( in.readable() >= nshifts*chunk_size + (nshifts-1) &&
out.writable()*8 >= fec->bits_in*chunk_size ) {
TPM totaldiscr[nsyncs];
for ( int s=0; s<nsyncs; ++s ) totaldiscr[s] = 0;
uint64_t outstream = 0;
int nout = 0;
softsymbol *pin = in.rd();
for ( int blocknum=0; blocknum<chunk_size; ++blocknum,pin+=nshifts ) {
TPM discr;
TUS result = update_sync(current_sync, pin, &discr);
outstream = (outstream<<fec->bits_in) | result;
nout += fec->bits_in;
if ( blocknum >= discr_delay ) totaldiscr[current_sync] += discr;
if ( ! resync_phase ) {
// Every one in [resync_period] chunks, run all decoders
// and compute average quality metrics
for ( int b=0; b<8; ++b,++pin ) {
TUS bits[NSYNCS];
for ( int s=0; s<NSYNCS; ++s ) {
TPM discr;
bits[s] = update_sync(s, pin->metrics4, &discr);
if ( bytenum*8 > TRACEBACK ) totaldiscr[s] += discr;
}
byte = (byte<<1) | bits[current_sync];
}
} else {
// Otherwise run only the selected decoder
for ( int b=0; b<8; ++b,++pin ) {
TUS bit = update_sync(current_sync, pin->metrics4, NULL);
byte = (byte<<1) | bit;
// Every [resync_period] chunks, also run the other decoders.
for ( int s=0; s<nsyncs; ++s ) {
if ( s == current_sync ) continue;
TPM discr;
(void)update_sync(s, pin, &discr);
if ( blocknum >= discr_delay ) totaldiscr[s] += discr;
}
}
in.read(8);
out.write(byte);
while ( nout >= 8 ) {
out.write(outstream>>(nout-8));
nout -= 8;
}
} // chunk_size
in.read(chunk_size*nshifts);
if ( nout ) fail("overlapping out");
if ( ! resync_phase ) {
// Switch to another decoder ?
int best = current_sync;
for ( int s=0; s<NSYNCS; ++s )
for ( int s=0; s<nsyncs; ++s )
if ( totaldiscr[s] > totaldiscr[best] ) best = s;
if ( best != current_sync ) {
if ( sch->debug ) fprintf(stderr, "{%d->%d}", current_sync, best);
@ -1181,116 +1379,10 @@ namespace leansdr {
if ( ++resync_phase >= resync_period ) resync_phase = 0;
}
}
}; // viterbi_sync
// VITERBI DECODING - BPSK
// ETSI TR 101 198: The BPSK variant of DVB-S transmits R=I,Q,I,Q,...
// Code rate 1/2 only.
struct viterbi_sync_bpsk : runnable {
static const int TRACEBACK = 32; // Suitable for BPSK 1/2
typedef uint8_t TS, TCS, TUS;
typedef uint16_t TBM;
typedef uint32_t TPM;
typedef bitpath<uint32_t,TUS,1,TRACEBACK> dvb_path;
typedef viterbi_dec<TS,64, TUS,2, TCS,4, TBM, TPM, dvb_path> dvb_dec;
typedef trellis<TS,64, TUS,2, 4> dvb_trellis;
private:
pipereader<softsymbol> in;
pipewriter<unsigned char> out;
static const int NSYNCS = 2; // Alignment of symbol pairs
dvb_dec *syncs[2];
int current_sync;
static const int chunk_size = 128;
int resync_phase;
public:
int resync_period;
viterbi_sync_bpsk(scheduler *sch,
pipebuf<softsymbol> &_in,
pipebuf<unsigned char> &_out)
: runnable(sch, "viterbi_sync_bpsk"),
in(_in), out(_out, chunk_size),
current_sync(0),
resync_phase(0),
resync_period(32)
{
dvb_trellis *trell = new dvb_trellis();
uint64_t dvb_polynomials[] = { DVBS_G1, DVBS_G2 };
trell->init_convolutional(dvb_polynomials);
for ( int s=0; s<NSYNCS; ++s ) syncs[s] = new dvb_dec(trell);
}
inline TUS update_sync(int s, TBM mX[4], TBM mY[4], TPM *discr) {
// Reconstruct QPSK metrics from pairs of BPSK symbols.
// EN 300 421, section 4.5 Baseband shaping and modulation
// EN 302 307, section 5.4.1
//
// IQ=10=(2) | IQ=00=(0)
// ----------+----------
// IQ=11=(3) | IQ=01=(1)
//
TBM vm[4];
vm[0] = mX[0] + mY[0];
vm[1] = mX[0] + mY[1];
vm[2] = mX[1] + mY[0];
vm[3] = mX[1] + mY[1];
return syncs[s]->update(vm, discr);
}
void run() {
while ( in.readable()>=2*8*chunk_size+1 && // +1 for pair alignment
out.writable()>=chunk_size ) {
unsigned long totaldiscr[NSYNCS];
for ( int s=0; s<NSYNCS; ++s ) totaldiscr[s] = 0;
for ( int bytenum=0; bytenum<chunk_size; ++bytenum ) {
// Decode one byte
unsigned char byte = 0;
softsymbol *pin = in.rd();
if ( ! resync_phase ) {
// Every one in [resync_period] chunks, run all decoders
// and compute average quality metrics
for ( int b=0; b<8; ++b,pin+=2 ) {
TUS bits[NSYNCS];
for ( int s=0; s<NSYNCS; ++s ) {
TPM discr;
bits[s] = update_sync(s,
(pin+s)->metrics4, (pin+s+1)->metrics4,
&discr);
if ( bytenum*8 > TRACEBACK ) totaldiscr[s] += discr;
}
byte = (byte<<1) | bits[current_sync];
}
} else {
// Otherwise run only the selected decoder
for ( int b=0; b<8; ++b,pin+=2 ) {
TUS bit = update_sync(current_sync,
(pin+current_sync)->metrics4,
(pin+current_sync+1)->metrics4, NULL);
byte = (byte<<1) | bit;
}
}
in.read(16);
out.write(byte);
} // chunk_size
if ( ! resync_phase ) {
// Switch to another decoder ?
int best = current_sync;
for ( int s=0; s<NSYNCS; ++s )
if ( totaldiscr[s] > totaldiscr[best] ) best = s;
if ( best != current_sync ) {
if ( sch->debug ) fprintf(stderr, "{%d->%d}", current_sync, best);
current_sync = best;
}
}
if ( ++resync_phase >= resync_period ) resync_phase = 0;
}
}
}; // viterbi_sync_bpsk
} // namespace

Wyświetl plik

@ -67,6 +67,11 @@ namespace leansdr {
return parity((uint32_t)(x^(x>>32)));
}
int log2(uint64_t x) {
int n = -1;
for ( ; x; ++n,x>>=1 ) ;
return n;
}
// Pre-computed sin/cos for 16-bit angles

Wyświetl plik

@ -257,8 +257,8 @@ namespace leansdr {
}; // simple_agc
typedef unsigned short u_angle; // [0,2PI[ in 65536 steps
typedef signed short s_angle; // [-PI,PI[ in 65536 steps
typedef uint16_t u_angle; // [0,2PI[ in 65536 steps
typedef int16_t s_angle; // [-PI,PI[ in 65536 steps
// GENERIC CONSTELLATION DECODING BY LOOK-UP TABLE.
@ -268,8 +268,8 @@ namespace leansdr {
// Up to 256 symbols.
struct softsymbol {
unsigned short metrics4[4]; // For Viterbi QPSK
unsigned char symbol; // 000000IQ for QPSK
int16_t cost; // For Viterbi with TBM=int16_t
uint8_t symbol;
};
// Target RMS amplitude for AGC
@ -283,10 +283,20 @@ namespace leansdr {
struct cstln_lut {
complex<signed char> *symbols;
int nsymbols;
enum predef { BPSK, QPSK, PSK8, APSK16, APSK32 };
cstln_lut(predef type, float gamma1=1, float gamma2=1) {
int nrotations;
enum predef {
BPSK, // DVB-S2 (and DVB-S variant)
QPSK, // DVB-S
PSK8, APSK16, APSK32, // DVB-S2
APSK64E, // DVB-S2X
QAM16, QAM64, QAM256 // For experimentation only
};
cstln_lut(predef type, float gamma1=1, float gamma2=1, float gamma3=1) {
switch ( type ) {
case BPSK:
nrotations = 2;
nsymbols = 2;
symbols = new complex<signed char>[nsymbols];
#if 0 // BPSK at 0°
@ -301,6 +311,7 @@ namespace leansdr {
case QPSK:
// EN 300 421, section 4.5 Baseband shaping and modulation
// EN 302 307, section 5.4.1
nrotations = 4;
nsymbols = 4;
symbols = new complex<signed char>[nsymbols];
symbols[0] = polar(1, 4, 0.5);
@ -311,10 +322,11 @@ namespace leansdr {
break;
case PSK8:
// EN 302 307, section 5.4.2
nrotations = 8;
nsymbols = 8;
symbols = new complex<signed char>[nsymbols];
symbols[0] = polar(1, 8, 0);
symbols[1] = polar(1, 8, 7);
symbols[0] = polar(1, 8, 1);
symbols[1] = polar(1, 8, 0);
symbols[2] = polar(1, 8, 4);
symbols[3] = polar(1, 8, 5);
symbols[4] = polar(1, 8, 2);
@ -327,6 +339,7 @@ namespace leansdr {
// EN 302 307, section 5.4.3
float r1 = sqrtf(4 / (1+3*gamma1*gamma1));
float r2 = gamma1 * r1;
nrotations = 4;
nsymbols = 16;
symbols = new complex<signed char>[nsymbols];
symbols[0] = polar(r2, 12, 1.5);
@ -353,6 +366,7 @@ namespace leansdr {
float r1 = sqrtf(8 / (1+3*gamma1*gamma1+4*gamma2*gamma2));
float r2 = gamma1 * r1;
float r3 = gamma2 * r1;
nrotations = 4;
nsymbols = 32;
symbols = new complex<signed char>[nsymbols];
symbols[0] = polar(r2, 12, 1.5);
@ -390,6 +404,44 @@ namespace leansdr {
make_lut_from_symbols();
break;
}
case APSK64E: {
// EN 302 307-2, section 5.4.5, Table 13e
float r1 =
sqrtf(64 / (4+12*gamma1*gamma1+20*gamma2*gamma2+28*gamma3*gamma3));
float r2 = gamma1 * r1;
float r3 = gamma2 * r1;
float r4 = gamma3 * r1;
nrotations = 4;
nsymbols = 64;
symbols = new complex<signed char>[nsymbols];
polar2( 0, r4, 1.0/ 4, 7.0/4, 3.0/ 4, 5.0/ 4);
polar2( 4, r4, 13.0/28, 43.0/28, 15.0/28, 41.0/28);
polar2( 8, r4, 1.0/28, 55.0/28, 27.0/28, 29.0/28);
polar2(12, r1, 1.0/ 4, 7.0/ 4, 3.0/ 4, 5.0/ 4);
polar2(16, r4, 9.0/28, 47.0/28, 19.0/28, 37.0/28);
polar2(20, r4, 11.0/28, 45.0/28, 17.0/28, 39.0/28);
polar2(24, r3, 1.0/20, 39.0/20, 19.0/20, 21.0/20);
polar2(28, r2, 1.0/12, 23.0/12, 11.0/12, 13.0/12);
polar2(32, r4, 5.0/28, 51.0/28, 23.0/28, 33.0/28);
polar2(36, r3, 9.0/20, 31.0/20, 11.0/20, 29.0/20);
polar2(40, r4, 3.0/28, 53.0/28, 25.0/28, 31.0/28);
polar2(44, r2, 5.0/12, 19.0/12, 7.0/12, 17.0/12);
polar2(48, r3, 1.0/ 4, 7.0/ 4, 3.0/ 4, 5.0/ 4);
polar2(52, r3, 7.0/20, 33.0/20, 13.0/20, 27.0/20);
polar2(56, r3, 3.0/20, 37.0/20, 17.0/20, 23.0/20);
polar2(60, r2, 1.0/ 4, 7.0/ 4, 3.0/ 4, 5.0/ 4);
make_lut_from_symbols();
break;
}
case QAM16:
make_qam(16);
break;
case QAM64:
make_qam(64);
break;
case QAM256:
make_qam(256);
break;
default:
fail("Constellation not implemented");
}
@ -401,7 +453,7 @@ namespace leansdr {
inline result *lookup(float I, float Q) {
// Handling of overflows beyond the lookup table:
// - For BPSK/QPSK/8PSK we only care about the phase,
// so the following is fine.
// so the following is harmless and improves locking at low SNR.
// - For amplitude modulations this is not appropriate.
// However, if there is enough noise to cause overflow,
// demodulation would probably fail anyway.
@ -424,28 +476,68 @@ namespace leansdr {
float a = i * 2*M_PI / n;
return complex<signed char>(r*cosf(a)*cstln_amp, r*sinf(a)*cstln_amp);
}
// Helper function for some constellation tables
void polar2(int i, float r, float a0, float a1, float a2, float a3) {
float a[] = { a0, a1, a2, a3 };
for ( int j=0; j<4; ++j ) {
float phi = a[j] * M_PI;
symbols[i+j] = complex<signed char>(r*cosf(phi)*cstln_amp,
r*sinf(phi)*cstln_amp);
}
}
void make_qam(int n) {
nrotations = 4;
nsymbols = n;
symbols = new complex<signed char>[nsymbols];
int m = sqrtl(n);
float scale;
{ // Average power in first quadrant with unit grid
int q = m / 2;
float avgpower = 2*(q*0.25+(q-1)*q/2+(q-1)*q*(2*q-1)/6) / q;
scale = 1.0 / sqrtf(avgpower);
}
// Arbitrary mapping
int s = 0;
for ( int x=0; x<m; ++x )
for ( int y=0; y<m; ++y ) {
float I = x - (float)(m-1)/2;
float Q = y - (float)(m-1)/2;
symbols[s].re = I * scale * cstln_amp;
symbols[s].im = Q * scale * cstln_amp;
++s;
}
make_lut_from_symbols();
}
result lut[R][R];
void make_lut_from_symbols() {
for ( int I=-R/2; I<R/2; ++I )
for ( int Q=-R/2; Q<R/2; ++Q ) {
result *pr = &lut[I&(R-1)][Q&(R-1)];
unsigned int dmin = R*2;
unsigned char smin = 0;
// Simplified metric:
// Distance to nearest minus distance to second-nearest.
// Null at edge of decision regions
// => Suitable for Viterbi with partial metrics.
uint8_t nearest = 0;
int32_t cost=R*R*2, cost2=R*R*2;
for ( int s=0; s<nsymbols; ++s ) {
unsigned int d2 =
int32_t d2 =
(I-symbols[s].re)*(I-symbols[s].re) +
(Q-symbols[s].im)*(Q-symbols[s].im);
if ( d2 > 65535 ) fail("Unexpected constellation");
unsigned int d = sqrtf(d2);
if ( d < dmin ) { dmin=d; smin=s; }
if ( nsymbols <= 4 ) {
pr->ss.metrics4[s] = d2;
if ( d2 < cost ) {
cost2 = cost;
cost = d2;
nearest = s;
} else if ( d2 < cost2 ) {
cost2 = d2;
}
}
float ph_symbol = atan2f(symbols[smin].im,symbols[smin].re);
if ( cost > 32767 ) cost = 32767;
if ( cost2 > 32767 ) cost2 = 32767;
pr->ss.cost = cost - cost2;
pr->ss.symbol = nearest;
float ph_symbol = atan2f(symbols[pr->ss.symbol].im,
symbols[pr->ss.symbol].re);
float ph_err = atan2f(Q,I) - ph_symbol;
if ( dmin > 255 ) fail("dmin overflow");
pr->ss.symbol = smin;
pr->phase_error = (s32)(ph_err * 65536 / (2*M_PI)); // Mod 65536
}
}
@ -455,22 +547,11 @@ namespace leansdr {
void harden() {
for ( int i=0; i<R; ++i )
for ( int q=0; q<R; ++q ) {
unsigned short *m = lut[i][q].ss.metrics4;
int best;
if ( nsymbols == 2 ) {
best = (m[0]<m[1]) ? 0 : 1;
m[0] = (best==0) ? 0 : 1;
m[1] = (best==1) ? 0 : 1;
} else {
if ( m[0]<=m[1] && m[0]<=m[2] && m[0]<=m[3] ) best = 0;
else if ( m[1]<=m[2] && m[1]<=m[3] ) best = 1;
else if ( m[2]<=m[3] ) best = 2;
else best = 3;
for ( int s=0; s<4; ++s )
m[s] = hamming_weight((uint8_t)(s^best));
}
}
}
softsymbol *ss = &lut[i][q].ss;
if ( ss->cost < 0 ) ss->cost = -1;
if ( ss->cost > 0 ) ss->cost = 1;
} // for I,Q
}
}; // cstln_lut
@ -667,7 +748,7 @@ namespace leansdr {
float freq_alpha = 0.04;
float freq_beta = 0.0012 / omega * pll_adjustment;
float gain_mu = 0.02 / (cstln_amp*cstln_amp) * 2;
int max_meas = chunk_size/meas_decimation + 1;
// Large margin on output_size because mu adjustments
// can lead to more than chunk_size/min_omega symbols.
@ -677,7 +758,7 @@ namespace leansdr {
( !ss_out || ss_out ->writable()>=max_meas ) &&
( !mer_out || mer_out ->writable()>=max_meas ) &&
( !cstln_out || cstln_out->writable()>=max_meas ) ) {
sampler->update_freq(freqw);
complex<T> *pin=in.rd(), *pin0=pin, *pend=pin+chunk_size;
@ -687,7 +768,7 @@ namespace leansdr {
complex<float> sg; // Symbol before AGC;
complex<float> s; // For MER estimation and constellation viewer
complex<signed char> *cstln_point = NULL;
while ( pin < pend ) {
// Here mu is the time of the next symbol counted from 0 at pin.
if ( mu < 1 ) {
@ -894,7 +975,7 @@ namespace leansdr {
if ( ! freq_beta ) fail("Excessive oversampling");
float gain_mu = 0.02 / (cstln_amp*cstln_amp) * 2;
int max_meas = chunk_size/meas_decimation + 1;
// Largin margin on output_size because mu adjustments
// can lead to more than chunk_size/min_omega symbols.
@ -908,7 +989,7 @@ namespace leansdr {
cu8 s;
u_angle symbol_arg = 0; // Exported for constellation viewer
while ( pin < pend ) {
// Here mu is the time of the next symbol counted from 0 at pin.
if ( mu < 1 ) {

Wyświetl plik

@ -5,6 +5,12 @@
#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.
@ -34,20 +40,24 @@ namespace leansdr {
states[s].branches[cs].pred = NOSTATE;
}
void init_convolutional(uint64_t G[]) {
// 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;
for ( nG=1; (1<<nG)<NCS; ++nG ) ;
int nG = log2(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 | (us*NSTATES);
uint32_t cs = 0;
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
@ -65,7 +75,7 @@ namespace leansdr {
}
void dump() {
for ( TS s=0; s<NSTATES; ++s ) {
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];
@ -80,12 +90,23 @@ namespace leansdr {
};
// 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);
virtual TUS update(TCS s, TBM cost, TPM *quality=NULL);
};
template<typename TS, int NSTATES,
typename TUS, int NUS,
typename TCS, int NCS,
typename TBM, typename TPM,
typename TP>
struct viterbi_dec {
struct viterbi_dec : viterbi_dec_interface<TUS,TCS,TBM,TPM> {
trellis<TS, NSTATES, TUS, NUS, NCS> *trell;
struct state {
@ -102,17 +123,25 @@ namespace leansdr {
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 = -1LL, best2_tpm = -1LL;
TPM best_tpm = max_tpm, best2_tpm = max_tpm;
TS best_state = 0;
// Update all states
for ( TS s=0; s<NSTATES; ++s ) {
TPM best_m = -1LL;
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 ( TCS cs=0; cs<NCS; ++cs ) {
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;
@ -149,6 +178,77 @@ namespace leansdr {
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 )
@ -157,6 +257,9 @@ namespace leansdr {
fprintf(stderr, "\n");
}
private:
TPM max_tpm;
};
// Paths (sequences of uncoded symbols) represented as bitstreams.