DATV demod: switched to work branch copy of leansdr

This commit is contained in:
f4exb 2019-03-17 21:31:42 +01:00
parent 7b9cb0e9fe
commit d4fe404dd6
31 changed files with 8844 additions and 2451 deletions

View File

@ -7,7 +7,12 @@ set(datv_SOURCES
datvdemodgui.cpp
datvdemodplugin.cpp
datvideostream.cpp
datvideorender.cpp
datvideorender.cpp
leansdr/dvb.cpp
leansdr/filtergen.cpp
leansdr/framework.cpp
leansdr/math.cpp
leansdr/sdr.cpp
)
set(datv_HEADERS
@ -15,7 +20,12 @@ set(datv_HEADERS
datvdemodgui.h
datvdemodplugin.h
datvideostream.h
datvideorender.h
datvideorender.h
leansdr/dvb.h
leansdr/filtergen.h
leansdr/framework.h
leansdr/math.h
leansdr/sdr.h
)
set(datv_FORMS

View File

@ -28,12 +28,87 @@ namespace leansdr {
static const int DEFAULT_GUI_DECIMATION = 64;
static inline cstln_lut<eucl_ss, 256> * make_dvbs2_constellation(cstln_lut<eucl_ss, 256>::predef c,
code_rate r)
{
float gamma1 = 1, gamma2 = 1, gamma3 = 1;
switch (c)
{
case cstln_lut<eucl_ss, 256>::APSK16:
// EN 302 307, section 5.4.3, Table 9
switch (r)
{
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;
case FEC89:
gamma1 = 2.60;
break;
case FEC910:
gamma1 = 2.57;
break;
default:
fail("cstln_lut<256>::make_dvbs2_constellation: Code rate not supported with APSK16");
return 0;
}
break;
case cstln_lut<eucl_ss, 256>::APSK32:
// EN 302 307, section 5.4.4, Table 10
switch (r)
{
case FEC34:
gamma1 = 2.84;
gamma2 = 5.27;
break;
case FEC45:
gamma1 = 2.72;
gamma2 = 4.87;
break;
case FEC56:
gamma1 = 2.64;
gamma2 = 4.64;
break;
case FEC89:
gamma1 = 2.54;
gamma2 = 4.33;
break;
case FEC910:
gamma1 = 2.53;
gamma2 = 4.30;
break;
default:
fail("cstln_lut<eucl_ss, 256>::make_dvbs2_constellation: Code rate not supported with APSK32");
return 0;
}
break;
case cstln_lut<eucl_ss, 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<eucl_ss, 256>(c, gamma1, gamma2, gamma3);
}
template<typename T> struct datvconstellation: runnable
{
T xymin, xymax;
unsigned long decimation;
unsigned long pixels_per_frame;
cstln_lut<256> **cstln; // Optional ptr to optional constellation
long pixels_per_frame;
cstln_lut<eucl_ss, 256> **cstln; // Optional ptr to optional constellation
TVScreen *m_objDATVScreen;
pipereader<complex<T> > in;
unsigned long phase;

View File

@ -531,34 +531,34 @@ void DATVDemod::InitDATVFramework()
switch(m_objRunning.enmModulation)
{
case BPSK:
m_objCfg.constellation = leansdr::cstln_lut<256>::BPSK;
m_objCfg.constellation = leansdr::cstln_lut<leansdr::eucl_ss, 256>::BPSK;
break;
case QPSK:
m_objCfg.constellation = leansdr::cstln_lut<256>::QPSK;
m_objCfg.constellation = leansdr::cstln_lut<leansdr::eucl_ss, 256>::QPSK;
break;
case PSK8:
m_objCfg.constellation = leansdr::cstln_lut<256>::PSK8;
m_objCfg.constellation = leansdr::cstln_lut<leansdr::eucl_ss, 256>::PSK8;
break;
case APSK16:
m_objCfg.constellation = leansdr::cstln_lut<256>::APSK16;
m_objCfg.constellation = leansdr::cstln_lut<leansdr::eucl_ss, 256>::APSK16;
break;
case APSK32:
m_objCfg.constellation = leansdr::cstln_lut<256>::APSK32;
m_objCfg.constellation = leansdr::cstln_lut<leansdr::eucl_ss, 256>::APSK32;
break;
case APSK64E:
m_objCfg.constellation = leansdr::cstln_lut<256>::APSK64E;
m_objCfg.constellation = leansdr::cstln_lut<leansdr::eucl_ss, 256>::APSK64E;
break;
case QAM16:
m_objCfg.constellation = leansdr::cstln_lut<256>::QAM16;
m_objCfg.constellation = leansdr::cstln_lut<leansdr::eucl_ss, 256>::QAM16;
break;
case QAM64:
m_objCfg.constellation = leansdr::cstln_lut<256>::QAM64;
m_objCfg.constellation = leansdr::cstln_lut<leansdr::eucl_ss, 256>::QAM64;
break;
case QAM256:
m_objCfg.constellation = leansdr::cstln_lut<256>::QAM256;
m_objCfg.constellation = leansdr::cstln_lut<leansdr::eucl_ss, 256>::QAM256;
break;
default:
m_objCfg.constellation = leansdr::cstln_lut<256>::BPSK;
m_objCfg.constellation = leansdr::cstln_lut<leansdr::eucl_ss, 256>::BPSK;
break;
}
@ -648,7 +648,7 @@ void DATVDemod::InitDATVFramework()
// Generic constellation receiver
p_symbols = new leansdr::pipebuf<leansdr::softsymbol>(m_objScheduler, "PSK soft-symbols", BUF_SYMBOLS);
p_symbols = new leansdr::pipebuf<leansdr::eucl_ss>(m_objScheduler, "PSK soft-symbols", BUF_SYMBOLS);
p_freq = new leansdr::pipebuf<leansdr::f32> (m_objScheduler, "freq", BUF_SLOW);
p_ss = new leansdr::pipebuf<leansdr::f32> (m_objScheduler, "SS", BUF_SLOW);
p_mer = new leansdr::pipebuf<leansdr::f32> (m_objScheduler, "MER", BUF_SLOW);
@ -682,7 +682,7 @@ void DATVDemod::InitDATVFramework()
return;
}
m_objDemodulator = new leansdr::cstln_receiver<leansdr::f32>(
m_objDemodulator = new leansdr::cstln_receiver<leansdr::f32, leansdr::eucl_ss>(
m_objScheduler,
sampler,
*p_preprocessed,
@ -694,8 +694,8 @@ void DATVDemod::InitDATVFramework()
if (m_objCfg.standard == DVB_S)
{
if ( m_objCfg.constellation != leansdr::cstln_lut<256>::QPSK
&& m_objCfg.constellation != leansdr::cstln_lut<256>::BPSK )
if ( m_objCfg.constellation != leansdr::cstln_lut<leansdr::eucl_ss, 256>::QPSK
&& m_objCfg.constellation != leansdr::cstln_lut<leansdr::eucl_ss, 256>::BPSK )
{
qWarning("DATVDemod::InitDATVFramework: non-standard constellation for DVB-S");
}

View File

@ -25,11 +25,7 @@ class DownChannelizer;
#define rfFilterFftLength 1024
#ifndef LEANSDR_FRAMEWORK
#define LEANSDR_FRAMEWORK
//LeanSDR
#include "leansdr/framework.h"
#include "leansdr/generic.h"
#include "leansdr/dsp.h"
@ -41,10 +37,10 @@ class DownChannelizer;
#include "leansdr/hdlc.h"
#include "leansdr/iess.h"
#endif
#include "datvconstellation.h"
#include "datvvideoplayer.h"
#include "datvideostream.h"
#include "datvideorender.h"
#include "channel/channelsinkapi.h"
#include "dsp/basebandsamplesink.h"
@ -60,9 +56,6 @@ class DownChannelizer;
#include "util/message.h"
#include "util/movingaverage.h"
#include "datvideostream.h"
#include "datvideorender.h"
#include <QMutex>
enum DATVModulation { BPSK, QPSK, PSK8, APSK16, APSK32, APSK64E, QAM16, QAM64, QAM256 };
@ -83,7 +76,7 @@ struct config
bool cnr; // Measure CNR
unsigned int decim; // Decimation, 0=auto
float Fm; // QPSK symbol rate (Hz)
leansdr::cstln_lut<256>::predef constellation;
leansdr::cstln_lut<leansdr::eucl_ss, 256>::predef constellation;
leansdr::code_rate fec;
float Ftune; // Bias frequency for the QPSK demodulator (Hz)
bool allow_drift;
@ -109,7 +102,7 @@ struct config
cnr(false),
decim(0),
Fm(2e6),
constellation(leansdr::cstln_lut<256>::QPSK),
constellation(leansdr::cstln_lut<leansdr::eucl_ss, 256>::QPSK),
fec(leansdr::FEC12),
Ftune(0),
allow_drift(false),
@ -323,7 +316,7 @@ private:
};
unsigned long m_lngExpectedReadIQ;
unsigned long m_lngReadIQ;
long m_lngReadIQ;
//************** LEANDBV Parameters **************
@ -372,7 +365,7 @@ private:
float *coeffs_sampler;
int ncoeffs_sampler;
leansdr::pipebuf<leansdr::softsymbol> *p_symbols;
leansdr::pipebuf<leansdr::eucl_ss> *p_symbols;
leansdr::pipebuf<leansdr::f32> *p_freq;
leansdr::pipebuf<leansdr::f32> *p_ss;
leansdr::pipebuf<leansdr::f32> *p_mer;
@ -386,7 +379,7 @@ private:
leansdr::file_writer<leansdr::cf32> *r_ppout;
//GENERIC CONSTELLATION RECEIVER
leansdr::cstln_receiver<leansdr::f32> *m_objDemodulator;
leansdr::cstln_receiver<leansdr::f32, leansdr::eucl_ss> *m_objDemodulator;
// DECONVOLUTION AND SYNCHRONIZATION
leansdr::pipebuf<leansdr::u8> *p_bytes;

View File

@ -21,7 +21,7 @@
#include <QMediaMetaData>
#include "datvdemodgui.h"
#include "datvideostream.h"
//#include "datvideostream.h"
#include "device/devicesourceapi.h"
#include "device/deviceuiset.h"

View File

@ -0,0 +1,250 @@
// This file is part of LeanSDR Copyright (C) 2016-2018 <pabr@pabr.org>.
// See the toplevel README for more information.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
#ifndef LEANSDR_BCH_H
#define LEANSDR_BCH_H
#include "leansdr/discrmath.h"
namespace leansdr
{
// Interface to hide the template parameters
struct bch_interface
{
virtual void encode(const uint8_t *msg, size_t msgbytes, uint8_t *out) = 0;
virtual int decode(uint8_t *cw, size_t cwbytes) = 0;
}; // bch_interface
// BCH error correction.
// T: Unsigned type for packing binary polynomials.
// N: Number of parity bits.
// NP: Width of the polynomials supplied.
// DP: Actual degree of the minimum polynomials (all must be same).
// TGF: Unsigned type for syndromes (must be wider than DP).
// GFTRUNCGEN: Generator polynomial for GF(2^DP), with X^DP omitted.
template <typename T, int N, int NP, int DP, typename TGF, int GFTRUNCGEN>
struct bch_engine : bch_interface
{
bch_engine(const bitvect<T, NP> *polys, int _npolys)
: npolys(_npolys)
{
// Build the generator polynomial (product of polys[]).
g = 1;
for (int i = 0; i < npolys; ++i)
g = g * polys[i];
// Convert the polynomials to truncated representation
// (with X^DP omitted) for use with divmod().
truncpolys = new bitvect<T, DP>[npolys];
for (int i = 0; i < npolys; ++i)
truncpolys[i].copy(polys[i]);
// Check which polynomial contains each root.
// Note: The DVB-S2 polynomials are numbered so that
// syndpoly[2*i]==i, but we don't use that property.
syndpolys = new int[2 * npolys];
for (int i = 0; i < 2 * npolys; ++i)
{
int j;
for (j = 0; j < npolys; ++j)
if (!eval_poly(truncpolys[j], true, 1 + i))
break;
if (j == npolys)
fail("Bad polynomials/root");
syndpolys[i] = j;
}
}
// Generate BCH parity bits.
void encode(const uint8_t *msg, size_t msgbytes, uint8_t *out)
{
bitvect<T, N> parity = shiftdivmod(msg, msgbytes, g);
// Output as bytes, coefficient of highest degree first
for (int i = N / 8; i--; ++out)
*out = parity.v[i / sizeof(T)] >> ((i & (sizeof(T) - 1)) * 8);
}
// Decode BCH.
// Return number of bits corrected, or -1 on failure.
int decode(uint8_t *cw, size_t cwbytes)
{
//again:
bool corrupted = false;
// Divide by individual polynomials.
// TBD Maybe do in parallel, scanning cw only once.
bitvect<T, DP> rem[npolys];
for (int j = 0; j < npolys; ++j)
{
rem[j] = divmod(cw, cwbytes, truncpolys[j]);
}
// Compute syndromes.
TGF S[2 * npolys];
for (int i = 0; i < 2 * npolys; ++i)
{
// Compute R(alpha^(1+i)), exploiting the fact that
// R(x)=Q(x)g_j(X)+rem_j(X) and g_j(alpha^(1+i))=0
// for some j that we already determined.
// TBD Compute even exponents using conjugates.
S[i] = eval_poly(rem[syndpolys[i]], false, 1 + i);
if (S[i])
corrupted = true;
}
if (!corrupted)
return 0;
#if 0
fprintf(stderr, "synd:");
for ( int i=0; i<2*npolys; ++i ) fprintf(stderr, " %04x", S[i]);
fprintf(stderr, "\n");
#endif
// S_j = R(alpha_j) = 0+E(alpha_j) = sum(l=1..L)((alpha^j)^i_l)
// where i_1 .. i_L are the degrees of the non-zero coefficients of E.
// S_j = sum(l=1..L)((alpha^i_l)^j) = sum(l=1..L)(X_l^j)
// Berlekamp - Massey
// http://en.wikipedia.org/wiki/Berlekamp%E2%80%93Massey_algorithm
// TBD More efficient to work with logs of syndromes ?
int NN = 2 * npolys;
TGF C[NN] = {
1,
},
B[NN] = {
1,
};
int L = 0, m = 1;
TGF b = 1;
for (int n = 0; n < NN; ++n)
{
TGF d = S[n];
for (int i = 1; i <= L; ++i)
d = GF.add(d, GF.mul(C[i], S[n - i]));
if (d == 0)
++m;
else
{
TGF d_div_b = GF.mul(d, GF.inv(b));
if (2 * L <= n)
{
TGF tmp[NN];
memcpy(tmp, C, sizeof(tmp));
for (int i = 0; i < NN - m; ++i)
C[m + i] = GF.sub(C[m + i], GF.mul(d_div_b, B[i]));
L = n + 1 - L;
memcpy(B, tmp, sizeof(B));
b = d;
m = 1;
}
else
{
for (int i = 0; i < NN - m; ++i)
C[m + i] = GF.sub(C[m + i], GF.mul(d_div_b, B[i]));
++m;
}
}
}
// L is the number of errors.
// C of degree L is the error locator polynomial (Lambda).
// C(X) = sum(l=1..L)(1-X_l*X).
#if 0
fprintf(stderr, "C[%d]=", L);
for ( int i=0; i<NN; ++i ) fprintf(stderr, " %04x", C[i]);
fprintf(stderr, "\n");
#endif
// Forney
// http://en.wikipedia.org/wiki/Forney_algorithm
// Simplified because coefficients are in GF(2).
// Find zeroes of C by exhaustive search.
// TODO Chien method
int roots_found = 0;
for (int i = 0; i < (1 << DP) - 1; ++i)
{
// Candidate root ALPHA^i
TGF v = eval_poly(C, L, i);
if (!v)
{
// ALPHA^i is a root of C, i.e. the inverse of an X_l.
int loc = (i ? (1 << DP) - 1 - i : 0); // exponent of inverse
// Reverse because cw[0..cwbytes-1] is stored MSB first
int rloc = cwbytes * 8 - 1 - loc;
if (rloc < 0)
{
// This may happen if the code is used truncated.
return -1;
}
cw[rloc / 8] ^= 128 >> (rloc & 7);
++roots_found;
if (roots_found == L)
break;
}
}
if (roots_found != L)
return -1;
return L;
}
private:
// Eval a GF(2)[X] polynomial at a power of ALPHA.
TGF eval_poly(const bitvect<T, DP> &poly, bool is_trunc, int rootexp)
{
TGF acc = 0;
int re = 0;
for (int i = 0; i < DP; ++i)
{
if (poly[i])
acc = GF.add(acc, GF.exp(re));
re += rootexp;
if (re >= (1 << DP) - 1)
re -= (1 << DP) - 1; // mod 2^DP-1 incrementally
}
if (is_trunc)
acc = GF.add(acc, GF.exp(re));
return acc;
}
// Eval a GF(2^16)[X] polynomial at a power of ALPHA.
TGF eval_poly(const TGF *poly, int deg, int rootexp)
{
TGF acc = 0;
int re = 0;
for (int i = 0; i <= deg; ++i)
{
acc = GF.add(acc, GF.mul(poly[i], GF.exp(re)));
re += rootexp;
if (re >= (1 << DP) - 1)
re -= (1 << DP) - 1; // mod 2^DP-1 incrementally
}
return acc;
}
bitvect<T, DP> *truncpolys;
int npolys;
int *syndpolys;
bitvect<T, N> g;
// Finite group for syndrome calculations
gf2n<TGF, DP, 2, GFTRUNCGEN> GF;
}; // bch_engine
} // namespace leansdr
#endif // LEANSDR_BCH_H

View File

@ -1,6 +1,25 @@
// This file is part of LeanSDR Copyright (C) 2016-2018 <pabr@pabr.org>.
// See the toplevel README for more information.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
#ifndef LEANSDR_CONVOLUTIONAL_H
#define LEANSDR_CONVOLUTIONAL_H
#include "framework.h"
#include "sdr.h"
namespace leansdr
{
@ -10,30 +29,21 @@ namespace leansdr
// 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>
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.
@ -42,9 +52,11 @@ struct deconvol_poly
{
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)];
@ -52,16 +64,17 @@ struct deconvol_poly
if (nb < halfway)
nerrors += parity(hist & POLY_ERRORS);
}
*pout = byte;
}
return nerrors;
}
private:
private:
Thist hist;
};
// deconvol_poly
}; // deconvol_poly
// ALGEBRAIC DECONVOLUTION, OPTIMIZED
@ -69,223 +82,186 @@ private:
// 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>
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)
{
}
deconvol_poly2() : inI(0), inQ(0) {}
// Remap and deconvolve [nb*8] symbols into [nb] bytes.
// Return estimated number of bit errors or -1 if error.
// Return estimated number of bit errors.
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;
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)
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.
#define LOOP(bit) { \
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; \
++pin; \
}
#define LOOP(bit) \
{ \
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; \
++pin; \
}
// Don't shift by more than the operand width
switch (sizeof(Thist) /* 8*/)
switch (sizeof(Thist) * 8)
{
#if 0 // Not needed yet - avoid compiler warnings
case 8:
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
#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
#endif
case 4:
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)
;
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 2:
LOOP(15)
;
LOOP(14)
;
LOOP(13)
;
LOOP(12)
;
LOOP(11)
;
LOOP(10)
;
LOOP(9)
;
LOOP(8)
;
case 16:
LOOP(15);
LOOP(14);
LOOP(13);
LOOP(12);
LOOP(11);
LOOP(10);
LOOP(9);
LOOP(8);
// Fall-through
case 1:
LOOP(7)
;
LOOP(6)
;
LOOP(5)
;
LOOP(4)
;
LOOP(3)
;
LOOP(2)
;
LOOP(1)
;
LOOP(0)
;
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;
fail("Thist not supported");
}
#undef LOOP
#endif
switch (sizeof(Thist) /* 8*/)
switch (sizeof(Thist) * 8)
{
#if 0 // Not needed yet - avoid compiler warnings
case 8:
*pout++ = wd >> 56;
*pout++ = wd >> 48;
*pout++ = wd >> 40;
*pout++ = wd >> 32;
// Fall-through
#if 0 // Not needed yet - avoid compiler warnings
case 64:
*pout++ = wd >> 56;
*pout++ = wd >> 48;
*pout++ = wd >> 40;
*pout++ = wd >> 32;
// Fall-through
#endif
case 4:
case 32:
*pout++ = wd >> 24;
*pout++ = wd >> 16;
// Fall-through
case 2:
case 16:
*pout++ = wd >> 8;
// Fall-through
case 1:
case 8:
*pout++ = wd;
break;
default:
fail("deconvol_poly2::run", "Thist not supported (2)");
return -1;
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;
}
private:
private:
Thist inI, inQ;
};
// deconvol_poly2
}; // deconvol_poly2
// CONVOLUTIONAL ENCODER
// QPSK 1/2 only.
template<typename Thist, uint64_t POLY1, uint64_t 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)
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);
@ -294,50 +270,55 @@ struct convol_poly2
}
}
}
private:
private:
Thist hist;
};
// convol_poly2
}; // convol_poly2
// Generic BPSK..256QAM and puncturing
template<typename Thist, int HISTSIZE>
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]
const Thist *polys; // [bits_out]
convol_multipoly() :
bits_in(0), bits_out(0), bps(0), polys(0), hist(0), nhist(0), sersymb(0), nsersymb(0)
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("leansdr::convol_multipoly::encode: convol_multipoly not configured");
return;
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));
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]));
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;
@ -350,17 +331,17 @@ struct convol_multipoly
// Ensure deterministic output size
// TBD We can relax this
if (nhist || nsersymb) {
fatal("leansdr::convol_multipoly::encode: partial run");
fatal("partial run");
}
}
private:
private:
Thist hist;
int nhist;
Thist sersymb;
int nsersymb;
};
// convol_multipoly
}; // convol_multipoly
}// namespace
} // namespace leansdr
#endif // LEANSDR_CONVOLUTIONAL_H
#endif // LEANSDR_CONVOLUTIONAL_H

View File

@ -0,0 +1,73 @@
// This file is part of LeanSDR Copyright (C) 2016-2018 <pabr@pabr.org>.
// See the toplevel README for more information.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
#ifndef LEANSDR_CRC_H
#define LEANSDR_CRC_H
#include <stdint.h>
#include "leansdr/discrmath.h"
namespace leansdr
{
// EN 302 307-1 section 5.1.4 CRC-8 encoder
struct crc8_engine
{
crc8_engine()
{
// Precompute
// EN 302 307-1 5.1.4 Figure 2
bitvect<uint8_t, 8> g = POLY_DVBS2_CRC8;
for (int u = 0; u < 256; ++u)
{
uint8_t u8 = u;
bitvect<uint8_t, 8> c = shiftdivmod(&u8, 1, g);
table[u] = c.v[0];
}
}
uint8_t compute(const uint8_t *buf, int len)
{
uint8_t c = 0;
for (; len--; ++buf)
c = table[c ^ *buf];
return c;
}
private:
static const uint8_t POLY_DVBS2_CRC8 = 0xd5; // (1)11010101
uint8_t table[256];
};
// CRC-32 ITU V.42 for FINGERPRINT
uint32_t crc32(const uint8_t *buf, int len)
{
static const uint32_t poly = 0xedb88320;
uint32_t c = 0xffffffff;
for (int i = 0; i < len; ++i)
{
c ^= buf[i];
for (int bit = 8; bit--;)
c = (c & 1) ? (c >> 1) ^ poly : (c >> 1);
}
return c ^ 0xffffffff;
}
} // namespace leansdr
#endif // LEANSDR_CRC_H

View File

@ -0,0 +1,268 @@
// This file is part of LeanSDR Copyright (C) 2018 <pabr@pabr.org>.
// See the toplevel README for more information.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
#ifndef LEANSDR_DISCRMATH_H
#define LEANSDR_DISCRMATH_H
#include <cstddef>
namespace leansdr
{
// Polynomials with N binary coefficients.
// This is designed to compile into trivial inline code.
// Bits are packed into words of type T.
// Unused most-significant bits of the last word are 0.
// T must be an unsigned type.
template <typename T, int N>
struct bitvect
{
static const size_t WSIZE = sizeof(T) * 8;
static const size_t NW = (N + WSIZE - 1) / WSIZE;
T v[NW];
bitvect() {}
bitvect(T val)
{
v[0] = val;
for (int i = 1; i < NW; ++i)
v[i] = 0;
}
// Assign from another bitvect, with truncation or padding
template <int M>
bitvect<T, N> &copy(const bitvect<T, M> &a)
{
int nw = (a.NW < NW) ? a.NW : NW;
for (int i = 0; i < nw; ++i)
v[i] = a.v[i];
if (M < N)
for (int i = a.NW; i < NW; ++i)
v[i] = 0;
if (M > N)
truncate_to_N();
return *this;
}
bool operator[](unsigned int i) const
{
return v[i / WSIZE] & ((T)1 << (i & (WSIZE - 1)));
}
// Add (in GF(2))
template <int M>
bitvect<T, N> &operator+=(const bitvect<T, M> &a)
{
int nw = (a.NW < NW) ? a.NW : NW;
for (int i = 0; i < nw; ++i)
v[i] ^= a.v[i];
if (M > N)
truncate_to_N();
return *this;
}
// Multiply by X^n
bitvect<T, N> &operator<<=(unsigned int n)
{
if (n >= WSIZE)
fail("shift exceeds word width");
T mask = ~(((T)-1) << n);
for (int i = NW - 1; i > 0; --i)
v[i] = (v[i] << n) | ((v[i - 1] >> (WSIZE - n)) & mask);
v[0] <<= n;
if (N != NW * WSIZE)
truncate_to_N();
return *this;
}
private:
// Call this whenever bits N .. NW*WSIZE-1 may have been set.
inline void truncate_to_N()
{
v[NW - 1] &= ((T)-1) >> (NW * WSIZE - N);
}
}; // bitvect
// Return m modulo (X^N+p).
// p is typically a generator polynomial of degree N
// with the highest-coefficient monomial ommitted.
// m is packed into nm words of type Tm.
// The MSB of m[0] is the highest-order coefficient of m.
template <typename T, int N, typename Tm>
bitvect<T, N> divmod(const Tm *m, size_t nm, const bitvect<T, N> &p)
{
bitvect<T, N> res = 0;
const Tm bitmask = (Tm)1 << (sizeof(Tm) * 8 - 1);
for (; nm--; ++m)
{
Tm mi = *m;
for (int bit = sizeof(Tm) * 8; bit--; mi <<= 1)
{
// Multiply by X, save outgoing coeff of degree N
bool resN = res[N - 1];
res <<= 1;
// Add m[i]
if (mi & bitmask)
res.v[0] ^= 1;
// Modulo X^N+p
if (resN)
res += p;
}
}
return res;
}
// Return (m*X^N) modulo (X^N+p).
// Same as divmod(), slightly faster than appending N zeroes.
template <typename T, int N, typename Tm>
bitvect<T, N> shiftdivmod(const Tm *m, size_t nm, const bitvect<T, N> &p,
T init = 0)
{
bitvect<T, N> res;
for (int i = 0; i < res.NW; ++i)
res.v[i] = init;
const Tm bitmask = (Tm)1 << (sizeof(Tm) * 8 - 1);
for (; nm--; ++m)
{
Tm mi = *m;
for (int bit = sizeof(Tm) * 8; bit--; mi <<= 1)
{
// Multiply by X, save outgoing coeff of degree N
bool resN = res[N - 1];
res <<= 1;
// Add m[i]*X^N
resN ^= (bool)(mi & bitmask);
// Modulo X^N+p
if (resN)
res += p;
}
}
return res;
}
template <typename T, int N>
bool operator==(const bitvect<T, N> &a, const bitvect<T, N> &b)
{
for (int i = 0; i < a.NW; ++i)
if (a.v[i] != b.v[i])
return false;
return true;
}
// Add (in GF(2))
template <typename T, int N>
bitvect<T, N> operator+(const bitvect<T, N> &a, const bitvect<T, N> &b)
{
bitvect<T, N> res;
for (int i = 0; i < a.NW; ++i)
res.v[i] = a.v[i] ^ b.v[i];
return res;
}
// Polynomial multiplication.
template <typename T, int N, int NB>
bitvect<T, N> operator*(bitvect<T, N> a, const bitvect<T, NB> &b)
{
bitvect<T, N> res = 0;
for (int i = 0; i < NB; ++i, a <<= 1)
if (b[i])
res += a;
// TBD If truncation is needed, do it only once at the end.
return res;
}
// Finite group GF(2^N), for small 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).
// 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 encoding elements of GF(2^N).
// Binary coefficients are packed, with degree 0 at LSB.
// (Te)0 is 0
// (Te)1 is 1
// (Te)2 is X
// (Te)3 is X+1
// (Te)4 is X^2
// TRUNCP is the encoding of the generator, with highest monomial omitted.
// ALPHA is a primitive element of GF(2^N). Usually "2"=[X] is chosen.
template <typename Te, int N, Te ALPHA, Te TRUNCP>
struct gf2n
{
typedef Te element;
static const Te alpha = ALPHA;
gf2n()
{
if (ALPHA != 2)
fail("alpha!=2 not implemented");
// Precompute log and exp tables.
Te alpha_i = 1; // ALPHA^0
for (int i = 0; i < (1 << N); ++i)
{
lut_exp[i] = alpha_i; // ALPHA^i
lut_exp[((1 << N) - 1) + i] = alpha_i; // Wrap to avoid modulo 2^N-1
lut_log[alpha_i] = i;
bool overflow = alpha_i & (1 << (N - 1));
alpha_i <<= 1; // Multiply by alpha=[X] i.e. increase degrees
alpha_i &= ~((~(Te)0) << N); // In case Te is wider than N bits
if (overflow)
alpha_i ^= TRUNCP; // Modulo P iteratively
}
}
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 (!x)
return 0;
return lut_exp[lut_log[x] + ((1 << N) - 1) - lut_log[y]];
}
inline Te inv(Te x)
{
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]; // Extra room for wrapping modulo 2^N-1
Te lut_log[1 << N];
}; // gf2n
} // namespace leansdr
#endif // LEANSDR_DISCRMATH_H

View File

@ -1,9 +1,25 @@
// This file is part of LeanSDR Copyright (C) 2016-2018 <pabr@pabr.org>.
// See the toplevel README for more information.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
#ifndef LEANSDR_DSP_H
#define LEANSDR_DSP_H
#include <math.h>
#include "leansdr/framework.h"
#include "leansdr/math.h"
#include <math.h>
namespace leansdr
{
@ -14,15 +30,15 @@ namespace leansdr
// [cconverter] converts complex streams between numric types,
// with optional ofsetting and rational scaling.
template<typename Tin, int Zin, typename Tout, int Zout, int Gn, int Gd>
struct cconverter: runnable
template <typename Tin, int Zin, typename Tout, int Zout, int Gn, int Gd>
struct cconverter : runnable
{
cconverter(scheduler *sch, pipebuf<complex<Tin> > &_in,
pipebuf<complex<Tout> > &_out) :
runnable(sch, "cconverter"), in(_in), out(_out)
cconverter(scheduler *sch, pipebuf<complex<Tin>> &_in,
pipebuf<complex<Tout>> &_out)
: runnable(sch, "cconverter"),
in(_in), out(_out)
{
}
void run()
{
unsigned long count = min(in.readable(), out.writable());
@ -30,25 +46,23 @@ struct cconverter: runnable
complex<Tout> *pout = out.wr();
for (; pin < pend; ++pin, ++pout)
{
pout->re = Zout + (pin->re - (Tin) Zin) * Gn / Gd;
pout->im = Zout + (pin->im - (Tin) Zin) * Gn / Gd;
pout->re = Zout + (pin->re - (Tin)Zin) * Gn / Gd;
pout->im = Zout + (pin->im - (Tin)Zin) * Gn / Gd;
}
in.read(count);
out.written(count);
}
private:
pipereader<complex<Tin> > in;
pipewriter<complex<Tout> > out;
private:
pipereader<complex<Tin>> in;
pipewriter<complex<Tout>> out;
};
template<typename T>
template <typename T>
struct cfft_engine
{
const unsigned int n;
cfft_engine(unsigned int _n) :
n(_n), invsqrtn(1.0 / sqrt(n))
const int n;
cfft_engine(int _n) : n(_n), invsqrtn(1.0 / sqrt(n))
{
// Compute log2(n)
logn = 0;
@ -56,29 +70,28 @@ struct cfft_engine
++logn;
// Bit reversal
bitrev = new int[n];
for (unsigned int i = 0; i < n; ++i)
for (int i = 0; i < n; ++i)
{
bitrev[i] = 0;
for (int b = 0; b < logn; ++b)
bitrev[i] = (bitrev[i] << 1) | ((i >> b) & 1);
}
// Float constants
omega = new complex<T> [n];
omega_rev = new complex<T> [n];
for (unsigned int i = 0; i < n; ++i)
omega = new complex<T>[n];
omega_rev = new complex<T>[n];
for (int i = 0; i < n; ++i)
{
float a = 2.0 * M_PI * i / n;
omega_rev[i].re = (omega[i].re = cosf(a));
omega_rev[i].im = -(omega[i].im = sinf(a));
}
}
void inplace(complex<T> *data, bool reverse = false)
{
// Bit-reversal permutation
for (unsigned int i = 0; i < n; ++i)
for (int i = 0; i < n; ++i)
{
unsigned int r = bitrev[i];
int r = bitrev[i];
if (r < i)
{
complex<T> tmp = data[i];
@ -100,7 +113,7 @@ struct cfft_engine
complex<T> &w = om[k * dom];
complex<T> &dqk = data[q + k];
complex<T> x(w.re * dqk.re - w.im * dqk.im,
w.re * dqk.im + w.im * dqk.re);
w.re * dqk.im + w.im * dqk.re);
data[q + k].re = data[p + k].re - x.re;
data[q + k].im = data[p + k].im - x.im;
data[p + k].re = data[p + k].re + x.re;
@ -111,7 +124,7 @@ struct cfft_engine
if (reverse)
{
float invn = 1.0 / n;
for (unsigned int i = 0; i < n; ++i)
for (int i = 0; i < n; ++i)
{
data[i].re *= invn;
data[i].im *= invn;
@ -119,24 +132,22 @@ struct cfft_engine
}
}
private:
private:
int logn;
int *bitrev;
complex<T> *omega, *omega_rev;
float invsqrtn;
};
template<typename T>
struct adder: runnable
template <typename T>
struct adder : runnable
{
adder(scheduler *sch, pipebuf<T> &_in1, pipebuf<T> &_in2, pipebuf<T> &_out) :
runnable(sch, "adder"),
in1(_in1),
in2(_in2),
out(_out)
adder(scheduler *sch,
pipebuf<T> &_in1, pipebuf<T> &_in2, pipebuf<T> &_out)
: runnable(sch, "adder"),
in1(_in1), in2(_in2), out(_out)
{
}
void run()
{
int n = out.writable();
@ -144,8 +155,7 @@ struct adder: runnable
n = in1.readable();
if (in2.readable() < n)
n = in2.readable();
T *pin1 = in1.rd(), *pin2 = in2.rd(), *pout = out.wr(), *pend = pout
+ n;
T *pin1 = in1.rd(), *pin2 = in2.rd(), *pout = out.wr(), *pend = pout + n;
while (pout < pend)
*pout++ = *pin1++ + *pin2++;
in1.read(n);
@ -153,24 +163,22 @@ struct adder: runnable
out.written(n);
}
private:
private:
pipereader<T> in1, in2;
pipewriter<T> out;
};
template<typename Tscale, typename Tin, typename Tout>
struct scaler: runnable
template <typename Tscale, typename Tin, typename Tout>
struct scaler : runnable
{
Tscale scale;
scaler(scheduler *sch, Tscale _scale, pipebuf<Tin> &_in, pipebuf<Tout> &_out) :
runnable(sch, "scaler"),
scale(_scale),
in(_in),
out(_out)
scaler(scheduler *sch, Tscale _scale,
pipebuf<Tin> &_in, pipebuf<Tout> &_out)
: runnable(sch, "scaler"),
scale(_scale),
in(_in), out(_out)
{
}
void run()
{
unsigned long count = min(in.readable(), out.writable());
@ -182,23 +190,20 @@ struct scaler: runnable
out.written(count);
}
private:
private:
pipereader<Tin> in;
pipewriter<Tout> out;
};
// [awgb_c] generates complex white gaussian noise.
template<typename T>
struct wgn_c: runnable
template <typename T>
struct wgn_c : runnable
{
wgn_c(scheduler *sch, pipebuf<complex<T> > &_out) :
runnable(sch, "awgn"),
stddev(1.0),
out(_out)
wgn_c(scheduler *sch, pipebuf<complex<T>> &_out)
: runnable(sch, "awgn"), stddev(1.0), out(_out)
{
}
void run()
{
int n = out.writable();
@ -220,21 +225,17 @@ struct wgn_c: runnable
}
out.written(n);
}
float stddev;
private:
pipewriter<complex<T> > out;
private:
pipewriter<complex<T>> out;
};
template<typename T>
struct naive_lowpass: runnable
template <typename T>
struct naive_lowpass : runnable
{
naive_lowpass(scheduler *sch, pipebuf<T> &_in, pipebuf<T> &_out, int _w) :
runnable(sch, "lowpass"),
in(_in),
out(_out),
w(_w)
naive_lowpass(scheduler *sch, pipebuf<T> &_in, pipebuf<T> &_out, int _w)
: runnable(sch, "lowpass"), in(_in), out(_out), w(_w)
{
}
@ -257,25 +258,23 @@ struct naive_lowpass: runnable
out.written(count);
}
private:
private:
pipereader<T> in;
pipewriter<T> out;
int w;
};
template<typename T, typename Tc>
struct fir_filter: runnable
template <typename T, typename Tc>
struct fir_filter : runnable
{
fir_filter(scheduler *sch, int _ncoeffs, Tc *_coeffs, pipebuf<T> &_in, pipebuf<T> &_out, unsigned int _decim = 1) :
runnable(sch, "fir_filter"),
freq_tap(NULL),
tap_multiplier(1),
freq_tol(0.1),
ncoeffs(_ncoeffs),
coeffs(_coeffs),
in(_in),
out(_out),
decim(_decim)
fir_filter(scheduler *sch, int _ncoeffs, Tc *_coeffs,
pipebuf<T> &_in, pipebuf<T> &_out,
unsigned int _decim = 1)
: runnable(sch, "fir_filter"),
ncoeffs(_ncoeffs), coeffs(_coeffs),
in(_in), out(_out),
decim(_decim),
freq_tap(NULL), tap_multiplier(1), freq_tol(0.1)
{
shifted_coeffs = new T[ncoeffs];
set_freq(0);
@ -292,16 +291,15 @@ struct fir_filter: runnable
if (fabs(current_freq - new_freq) > freq_tol)
{
if (sch->verbose)
fprintf(stderr, "Shifting filter %f -> %f\n", current_freq,
new_freq);
fprintf(stderr, "Shifting filter %f -> %f\n",
current_freq, new_freq);
set_freq(new_freq);
}
}
unsigned long count = min((in.readable() - ncoeffs) / decim,
out.writable());
T *pin = in.rd() + ncoeffs, *pend = pin + count * decim, *pout =
out.wr();
out.writable());
T *pin = in.rd() + ncoeffs, *pend = pin + count * decim, *pout = out.wr();
// TBD use coeffs when current_freq=0 (fewer mults if float)
for (; pin < pend; pin += decim, ++pout)
{
@ -316,25 +314,20 @@ struct fir_filter: runnable
out.written(count);
}
public:
float *freq_tap;
float tap_multiplier;
float freq_tol;
private:
private:
unsigned int ncoeffs;
Tc *coeffs;
pipereader<T> in;
pipewriter<T> out;
unsigned int decim;
T *shifted_coeffs;
float current_freq;
void set_freq(float f)
{
for (unsigned int i = 0; i < ncoeffs; ++i)
{
float a = 2 * M_PI * f * (i - (ncoeffs / 2.0f));
float a = 2 * M_PI * f * (i - ncoeffs / 2);
float c = cosf(a), s = sinf(a);
// TBD Support T=complex
shifted_coeffs[i].re = coeffs[i] * c;
@ -342,30 +335,29 @@ private:
}
current_freq = f;
}
};
// fir_filter
public:
float *freq_tap;
float tap_multiplier;
float freq_tol;
}; // fir_filter
// FIR FILTER WITH INTERPOLATION AND DECIMATION
template<typename T, typename Tc>
struct fir_resampler: runnable
template <typename T, typename Tc>
struct fir_resampler : runnable
{
fir_resampler(scheduler *sch, int _ncoeffs, Tc *_coeffs, pipebuf<T> &_in, pipebuf<T> &_out, int _interp = 1, int _decim = 1) :
runnable(sch, "fir_resampler"),
ncoeffs(_ncoeffs),
coeffs(_coeffs),
interp(_interp),
decim(_decim),
in(_in),
out(_out, interp),
freq_tap(NULL),
tap_multiplier(1),
freq_tol(0.1)
fir_resampler(scheduler *sch, int _ncoeffs, Tc *_coeffs,
pipebuf<T> &_in, pipebuf<T> &_out,
int _interp = 1, int _decim = 1)
: runnable(sch, "fir_resampler"),
ncoeffs(_ncoeffs), coeffs(_coeffs),
interp(_interp), decim(_decim),
in(_in), out(_out, interp),
freq_tap(NULL), tap_multiplier(1), freq_tol(0.1)
{
if (decim != 1) {
fail("fir_resampler::fir_resampler", "decim not implemented"); // TBD
return;
}
if (decim != 1)
fail("fir_resampler: decim not implemented"); // TBD
shifted_coeffs = new T[ncoeffs];
set_freq(0);
}
@ -381,8 +373,8 @@ struct fir_resampler: runnable
if (fabs(current_freq - new_freq) > freq_tol)
{
if (sch->verbose)
fprintf(stderr, "Shifting filter %f -> %f\n", current_freq,
new_freq);
fprintf(stderr, "Shifting filter %f -> %f\n",
current_freq, new_freq);
set_freq(new_freq);
}
}
@ -390,7 +382,7 @@ struct fir_resampler: runnable
if (in.readable() * interp < ncoeffs)
return;
unsigned long count = min((in.readable() * interp - ncoeffs) / interp,
out.writable() / interp);
out.writable() / interp);
int latency = (ncoeffs + interp) / interp;
T *pin = in.rd() + latency, *pend = pin + count, *pout = out.wr();
// TBD use coeffs when current_freq=0 (fewer mults if float)
@ -410,20 +402,21 @@ struct fir_resampler: runnable
out.written(count * interp);
}
public:
float *freq_tap;
float tap_multiplier;
float freq_tol;
private:
private:
unsigned int ncoeffs;
Tc *coeffs;
int interp, decim;
pipereader<T> in;
pipewriter<T> out;
public:
float *freq_tap;
float tap_multiplier;
float freq_tol;
private:
T *shifted_coeffs;
float current_freq;
void set_freq(float f)
{
for (int i = 0; i < ncoeffs; ++i)
@ -436,9 +429,8 @@ private:
}
current_freq = f;
}
};
// fir_resampler
}; // fir_resampler
}// namespace
} // namespace leansdr
#endif // LEANSDR_DSP_H
#endif // LEANSDR_DSP_H

View File

@ -0,0 +1,46 @@
#include "dvb.h"
namespace leansdr
{
deconvol_sync_simple *make_deconvol_sync_simple(scheduler *sch,
pipebuf<eucl_ss> &_in,
pipebuf<u8> &_out,
enum code_rate rate)
{
// EN 300 421, section 4.4.3 Inner coding
uint32_t pX, pY;
switch (rate)
{
case FEC12:
pX = 0x1; // 1
pY = 0x1; // 1
break;
case FEC23:
case FEC46:
pX = 0xa; // 1010 (Handle as FEC4/6, no half-symbols)
pY = 0xf; // 1111
break;
case FEC34:
pX = 0x5; // 101
pY = 0x6; // 110
break;
case FEC56:
pX = 0x15; // 10101
pY = 0x1a; // 11010
break;
case FEC78:
pX = 0x45; // 1000101
pY = 0x7a; // 1111010
break;
default:
//fail("Code rate not implemented");
// For testing DVB-S2 constellations.
fprintf(stderr, "Code rate not implemented; proceeding anyway\n");
pX = pY = 1;
}
return new deconvol_sync_simple(sch, _in, _out, DVBS_G1, DVBS_G2, pX, pY);
}
} // leansdr

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,20 @@
#include <stdio.h>
#include "filtergen.h"
namespace leansdr
{
namespace filtergen
{
void dump_filter(const char *name, int ncoeffs, float *coeffs)
{
fprintf(stderr, "%s = [", name);
for (int i = 0; i < ncoeffs; ++i)
fprintf(stderr, "%s %f", (i ? "," : ""), coeffs[i]);
fprintf(stderr, " ];\n");
}
} // filtergen
} // leansdr

View File

@ -1,90 +1,123 @@
// This file is part of LeanSDR Copyright (C) 2016-2018 <pabr@pabr.org>.
// See the toplevel README for more information.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
#ifndef LEANSDR_FILTERGEN_H
#define LEANSDR_FILTERGEN_H
#include <math.h>
namespace leansdr {
namespace filtergen {
template<typename T>
void normalize_power(int n, T *coeffs, float gain=1) {
float s2 = 0;
for ( int i=0; i<n; ++i ) s2 = s2 + coeffs[i]*coeffs[i]; // TBD complex
if ( s2 ) gain /= gen_sqrt(s2);
for ( int i=0; i<n; ++i ) coeffs[i] = coeffs[i] * gain;
}
#include "framework.h"
template<typename T>
void normalize_dcgain(int n, T *coeffs, float gain=1) {
float s = 0;
for ( int i=0; i<n; ++i ) s = s + coeffs[i];
if ( s ) gain /= s;
for ( int i=0; i<n; ++i ) coeffs[i] = coeffs[i] * gain;
}
namespace leansdr
{
namespace filtergen
{
// Generate coefficients for a sinc filter.
// https://en.wikipedia.org/wiki/Sinc_filter
template <typename T>
void normalize_power(int n, T *coeffs, float gain = 1)
{
float s2 = 0;
for (int i = 0; i < n; ++i)
s2 = s2 + coeffs[i] * coeffs[i]; // TBD complex
if (s2)
gain /= gen_sqrt(s2);
for (int i = 0; i < n; ++i)
coeffs[i] = coeffs[i] * gain;
}
template<typename T>
int lowpass(int order, float Fcut, T **coeffs, float gain=1) {
int ncoeffs = order + 1;
*coeffs = new T[ncoeffs];
for ( int i=0; i<ncoeffs; ++i ) {
float t = i - (ncoeffs-1)*0.5;
float sinc = 2*Fcut * (t ? sin(2*M_PI*Fcut*t)/(2*M_PI*Fcut*t) : 1);
#if 0 // Hamming
template <typename T>
void normalize_dcgain(int n, T *coeffs, float gain = 1)
{
float s = 0;
for (int i = 0; i < n; ++i)
s = s + coeffs[i];
if (s)
gain /= s;
for (int i = 0; i < n; ++i)
coeffs[i] = coeffs[i] * gain;
}
template <typename T>
void cancel_dcgain(int n, T *coeffs)
{
float s = 0;
for (int i = 0; i < n; ++i)
s = s + coeffs[i];
for (int i = 0; i < n; ++i)
coeffs[i] -= s / n;
}
// Generate coefficients for a sinc filter.
// https://en.wikipedia.org/wiki/Sinc_filter
template <typename T>
int lowpass(int order, float Fcut, T **coeffs, float gain = 1)
{
int ncoeffs = order + 1;
*coeffs = new T[ncoeffs];
for (int i = 0; i < ncoeffs; ++i)
{
float t = i - (ncoeffs - 1) * 0.5;
float sinc = 2 * Fcut * (t ? sin(2 * M_PI * Fcut * t) / (2 * M_PI * Fcut * t) : 1);
#if 0 // Hamming
float alpha = 25.0/46, beta = 21.0/46;
float window = alpha - beta*cos(2*M_PI*i/order);
#else
float window = 1;
float window = 1;
#endif
(*coeffs)[i] = sinc * window;
}
normalize_dcgain(ncoeffs, *coeffs, gain);
return ncoeffs;
(*coeffs)[i] = sinc * window;
}
normalize_dcgain(ncoeffs, *coeffs, gain);
return ncoeffs;
}
// Generate coefficients for a RRC filter.
// https://en.wikipedia.org/wiki/Root-raised-cosine_filter
// Generate coefficients for a RRC filter.
// https://en.wikipedia.org/wiki/Root-raised-cosine_filter
template<typename T>
int root_raised_cosine(int order, float Fs, float rolloff, T **coeffs) {
float B = rolloff, pi = M_PI;
int ncoeffs = (order+1) | 1;
*coeffs = new T[ncoeffs];
for ( int i=0; i<ncoeffs; ++i ) {
int t = i - ncoeffs/2;
float c;
if ( t == 0 )
c = sqrt(Fs) * (1-B+4*B/pi);
else {
float tT = t * Fs;
float den = pi*tT*(1-(4*B*tT)*(4*B*tT));
if ( ! den )
c = B*sqrt(Fs/2) * ( (1+2/pi)*sin(pi/(4*B)) +
(1-2/pi)*cos(pi/(4*B)) );
else
c = sqrt(Fs) * ( sin(pi*tT*(1-B)) +
4*B*tT*cos(pi*tT*(1+B)) ) / den;
}
(*coeffs)[i] = c;
}
normalize_dcgain(ncoeffs, *coeffs);
return ncoeffs;
template <typename T>
int root_raised_cosine(int order, float Fs, float rolloff, T **coeffs)
{
float B = rolloff, pi = M_PI;
int ncoeffs = (order + 1) | 1;
*coeffs = new T[ncoeffs];
for (int i = 0; i < ncoeffs; ++i)
{
int t = i - ncoeffs / 2;
float c;
if (t == 0)
c = sqrt(Fs) * (1 - B + 4 * B / pi);
else
{
float tT = t * Fs;
float den = pi * tT * (1 - (4 * B * tT) * (4 * B * tT));
if (!den)
c = B * sqrt(Fs / 2) * ((1 + 2 / pi) * sin(pi / (4 * B)) + (1 - 2 / pi) * cos(pi / (4 * B)));
else
c = sqrt(Fs) * (sin(pi * tT * (1 - B)) + 4 * B * tT * cos(pi * tT * (1 + B))) / den;
}
(*coeffs)[i] = c;
}
normalize_dcgain(ncoeffs, *coeffs);
return ncoeffs;
}
// Dump filter coefficients for matlab/octave
void dump_filter(const char *name, int ncoeffs, float *coeffs);
// Dump filter coefficients for matlab/octave
} // namespace filtergen
} // namespace leansdr
inline void dump_filter(const char *name, int ncoeffs, float *coeffs) {
fprintf(stderr, "%s = [", name);
for ( int i=0; i<ncoeffs; ++i )
fprintf(stderr, "%s %f", (i?",":""), coeffs[i]);
fprintf(stderr, " ];\n");
}
} // namespace
} // namespace
#endif // LEANSDR_FILTERGEN_H
#endif // LEANSDR_FILTERGEN_H

View File

@ -0,0 +1,16 @@
#include "framework.h"
namespace leansdr
{
void fatal(const char *s)
{
perror(s);
}
void fail(const char *s)
{
fprintf(stderr, "** %s\n", s);
}
} // leansdr

View File

@ -1,26 +1,37 @@
// This file is part of LeanSDR Copyright (C) 2016-2018 <pabr@pabr.org>.
// See the toplevel README for more information.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
#ifndef LEANSDR_FRAMEWORK_H
#define LEANSDR_FRAMEWORK_H
#include <cstddef>
#include <algorithm>
#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#ifndef VERSION
#define VERSION "undefined"
#endif
namespace leansdr
{
inline void fatal(const char *s)
{
perror(s);
}
inline void fail(const char *f, const char *s)
{
fprintf(stderr, "leansdr::%s: %s\n", f, s);
}
void fatal(const char *s);
void fail(const char *s);
//////////////////////////////////////////////////////////////////////
// DSP framework
@ -50,25 +61,25 @@ struct pipebuf_common
virtual void dump(std::size_t *total_bufs)
{
(void) total_bufs;
(void)total_bufs;
}
pipebuf_common(const char *_name) :
name(_name)
const char *name;
pipebuf_common(const char *_name) : name(_name)
{
}
virtual ~pipebuf_common()
{
}
const char *name;
};
struct runnable_common
{
runnable_common(const char *_name) :
name(_name)
const char *name;
runnable_common(const char *_name) : name(_name)
{
}
@ -83,13 +94,12 @@ struct runnable_common
virtual void shutdown()
{
}
#ifdef DEBUG
~runnable_common()
{ fprintf(stderr, "Deallocating %s !\n", name);}
{
fprintf(stderr, "Deallocating %s !\n", name);
}
#endif
const char *name;
};
struct window_placement
@ -105,39 +115,35 @@ struct scheduler
runnable_common *runnables[MAX_RUNNABLES];
int nrunnables;
window_placement *windows;
bool verbose, debug;
bool verbose, debug, debug2;
scheduler() :
npipes(0), nrunnables(0), windows(nullptr), verbose(false), debug(false)
scheduler() : npipes(0),
nrunnables(0),
windows(NULL),
verbose(false),
debug(false),
debug2(false)
{
std::fill(pipes, pipes + MAX_PIPES, nullptr);
std::fill(runnables, runnables + MAX_RUNNABLES, nullptr);
}
void add_pipe(pipebuf_common *p)
{
if (npipes == MAX_PIPES)
{
fail("scheduler::add_pipe", "MAX_PIPES");
return;
}
fail("MAX_PIPES");
pipes[npipes++] = p;
}
void add_runnable(runnable_common *r)
{
if (nrunnables == MAX_RUNNABLES)
{
fail("scheduler::add_runnable", "MAX_RUNNABLES");
}
fail("MAX_RUNNABLES");
runnables[nrunnables++] = r;
}
void step()
{
for (int i = 0; i < nrunnables; ++i) {
for (int i = 0; i < nrunnables; ++i)
runnables[i]->run();
}
}
void run()
@ -148,26 +154,23 @@ struct scheduler
{
step();
unsigned long long h = hash();
if (h == prev_hash) {
if (h == prev_hash)
break;
}
prev_hash = h;
}
}
void shutdown()
{
for (int i = 0; i < nrunnables; ++i) {
for (int i = 0; i < nrunnables; ++i)
runnables[i]->shutdown();
}
}
unsigned long long hash()
{
unsigned long long h = 0;
for (int i = 0; i < npipes; ++i) {
for (int i = 0; i < npipes; ++i)
h += (1 + i) * pipes[i]->hash();
}
return h;
}
@ -175,26 +178,26 @@ struct scheduler
{
fprintf(stderr, "\n");
std::size_t total_bufs = 0;
for (int i = 0; i < npipes; ++i) {
for (int i = 0; i < npipes; ++i)
pipes[i]->dump(&total_bufs);
}
fprintf(stderr, "leansdr::scheduler::dump Total buffer memory: %ld KiB\n", (unsigned long) total_bufs / 1024);
fprintf(stderr, "Total buffer memory: %ld KiB\n",
(unsigned long)total_bufs / 1024);
}
};
struct runnable: runnable_common
struct runnable : runnable_common
{
runnable(scheduler *_sch, const char *name) :
runnable_common(name), sch(_sch)
runnable(scheduler *_sch, const char *name) : runnable_common(name), sch(_sch)
{
sch->add_runnable(this);
}
protected:
protected:
scheduler *sch;
};
template<typename T>
struct pipebuf: pipebuf_common
template <typename T>
struct pipebuf : pipebuf_common
{
T *buf;
T *rds[MAX_READERS];
@ -207,9 +210,13 @@ struct pipebuf: pipebuf_common
return sizeof(T);
}
pipebuf(scheduler *sch, const char *name, unsigned long size) :
pipebuf_common(name), buf(new T[size]), nrd(0), wr(buf), end(
buf + size), min_write(1), total_written(0), total_read(0)
pipebuf(scheduler *sch, const char *name, unsigned long size) : pipebuf_common(name),
buf(new T[size]),
nrd(0), wr(buf),
end(buf + size),
min_write(1),
total_written(0),
total_read(0)
{
sch->add_pipe(this);
}
@ -217,10 +224,7 @@ struct pipebuf: pipebuf_common
int add_reader()
{
if (nrd == MAX_READERS)
{
fail("pipebuf::add_reader", "too many readers");
return nrd;
}
fail("too many readers");
rds[nrd] = wr;
return nrd++;
}
@ -229,16 +233,12 @@ struct pipebuf: pipebuf_common
{
T *rd = wr;
for (int i = 0; i < nrd; ++i)
{
if (rds[i] < rd) {
if (rds[i] < rd)
rd = rds[i];
}
}
memmove(buf, rd, (wr - rd) * sizeof(T));
wr -= rd - buf;
for (int i = 0; i < nrd; ++i) {
for (int i = 0; i < nrd; ++i)
rds[i] -= rd - buf;
}
}
long long hash()
@ -248,71 +248,53 @@ struct pipebuf: pipebuf_common
void dump(std::size_t *total_bufs)
{
if (total_written < 10000) {
fprintf(stderr, "leansdr::pipebuf::dump: .%-16s : %4ld/%4ld", name, total_read, total_written);
} else if (total_written < 1000000) {
fprintf(stderr, "leansdr::pipebuf::dump: .%-16s : %3ldk/%3ldk", name, total_read / 1000, total_written / 1000);
} else {
fprintf(stderr, "leansdr::pipebuf::dump: .%-16s : %3ldM/%3ldM", name, total_read / 1000000, total_written / 1000000);
}
if (total_written < 10000)
fprintf(stderr, ".%-16s : %4ld/%4ld", name, total_read,
total_written);
else if (total_written < 1000000)
fprintf(stderr, ".%-16s : %3ldk/%3ldk", name, total_read / 1000,
total_written / 1000);
else
fprintf(stderr, ".%-16s : %3ldM/%3ldM", name, total_read / 1000000,
total_written / 1000000);
*total_bufs += (end - buf) * sizeof(T);
unsigned long nw = end - wr;
fprintf(stderr, "leansdr::pipebuf: %6ld writable %c,", nw, (nw < min_write) ? '!' : ' ');
fprintf(stderr, " %6ld writable %c,", nw, (nw < min_write) ? '!' : ' ');
T *rd = wr;
for (int j = 0; j < nrd; ++j)
{
if (rds[j] < rd) {
if (rds[j] < rd)
rd = rds[j];
}
}
fprintf(stderr, "leansdr::pipebuf::dump: %6d unread (", (int) (wr - rd));
for (int j = 0; j < nrd; ++j) {
fprintf(stderr, "leansdr::pipebuf: %d", (int) (wr - rds[j]));
}
fprintf(stderr, "leansdr::pipebuf::dump: )\n");
fprintf(stderr, " %6d unread (", (int)(wr - rd));
for (int j = 0; j < nrd; ++j)
fprintf(stderr, " %d", (int)(wr - rds[j]));
fprintf(stderr, " )\n");
}
unsigned long min_write;
unsigned long total_written, total_read;
#ifdef DEBUG
~pipebuf()
{ fprintf(stderr, "Deallocating %s !\n", name);}
{
fprintf(stderr, "Deallocating %s !\n", name);
}
#endif
};
template<typename T>
template <typename T>
struct pipewriter
{
pipebuf<T> &buf;
pipewriter(pipebuf<T> &_buf, unsigned long min_write = 1) :
buf(_buf)
pipewriter(pipebuf<T> &_buf, unsigned long min_write = 1) : buf(_buf)
{
if (min_write > buf.min_write) {
if (min_write > buf.min_write)
buf.min_write = min_write;
}
}
/** Return number of items writable at this->wr, 0 if full. */
unsigned long writable()
// Return number of items writable at this->wr, 0 if full.
long writable()
{
if (buf.end < buf.wr)
{
fprintf(stderr, "leansdr::pipewriter::writable: overflow in %s buffer\n", buf.name);
return 0;
}
unsigned long delta = buf.end - buf.wr;
if (delta < buf.min_write) {
if (buf.end < buf.min_write + buf.wr)
buf.pack();
}
return delta;
return buf.end - buf.wr;
}
T *wr()
@ -324,9 +306,9 @@ struct pipewriter
{
if (buf.wr + n > buf.end)
{
fprintf(stderr, "leansdr::pipewriter::written: overflow in %s buffer\n", buf.name);
return;
fprintf(stderr, "Bug: overflow to %s\n", buf.name);
}
buf.wr += n;
buf.total_written += n;
}
@ -340,38 +322,36 @@ struct pipewriter
// Convenience functions for working with optional pipes
template<typename T>
pipewriter<T> *opt_writer(pipebuf<T> *buf)
template <typename T>
pipewriter<T> *opt_writer(pipebuf<T> *buf, unsigned long min_write = 1)
{
return buf ? new pipewriter<T>(*buf) : NULL;
return buf ? new pipewriter<T>(*buf, min_write) : NULL;
}
template<typename T>
bool opt_writable(pipewriter<T> *p, unsigned int n = 1)
template <typename T>
bool opt_writable(pipewriter<T> *p, int n = 1)
{
return (p == NULL) || p->writable() >= n;
}
template<typename T>
template <typename T>
void opt_write(pipewriter<T> *p, T val)
{
if (p) {
if (p)
p->write(val);
}
}
template<typename T>
template <typename T>
struct pipereader
{
pipebuf<T> &buf;
int id;
pipereader(pipebuf<T> &_buf) :
buf(_buf), id(_buf.add_reader())
pipereader(pipebuf<T> &_buf) : buf(_buf), id(_buf.add_reader())
{
}
unsigned long readable()
long readable()
{
return buf.wr - buf.rds[id];
}
@ -385,9 +365,9 @@ struct pipereader
{
if (buf.rds[id] + n > buf.wr)
{
fprintf(stderr, "leansdr::pipereader::read: underflow in %s buffer\n", buf.name);
return;
fprintf(stderr, "Bug: underflow from %s\n", buf.name);
}
buf.rds[id] += n;
buf.total_read += n;
}
@ -395,7 +375,8 @@ struct pipereader
// Math functions for templates
template<typename T> T gen_sqrt(T x);
template <typename T>
T gen_sqrt(T x);
inline float gen_sqrt(float x)
{
return sqrtf(x);
@ -411,7 +392,8 @@ inline long double gen_sqrt(long double x)
return sqrtl(x);
}
template<typename T> T gen_abs(T x);
template <typename T>
T gen_abs(T x);
inline float gen_abs(float x)
{
return fabsf(x);
@ -427,7 +409,8 @@ inline long int gen_abs(long int x)
return labs(x);
}
template<typename T> T gen_hypot(T x, T y);
template <typename T>
T gen_hypot(T x, T y);
inline float gen_hypot(float x, float y)
{
return hypotf(x, y);
@ -438,7 +421,8 @@ inline long double gen_hypot(long double x, long double y)
return hypotl(x, y);
}
template<typename T> T gen_atan2(T y, T x);
template <typename T>
T gen_atan2(T y, T x);
inline float gen_atan2(float y, float x)
{
return atan2f(y, x);
@ -449,13 +433,13 @@ inline long double gen_atan2(long double y, long double x)
return atan2l(y, x);
}
template<typename T>
template <typename T>
T min(const T &x, const T &y)
{
return (x < y) ? x : y;
}
template<typename T>
template <typename T>
T max(const T &x, const T &y)
{
return (x < y) ? y : x;
@ -470,6 +454,6 @@ typedef signed char s8;
typedef signed short s16;
typedef signed long s32;
} // namespace
} // namespace leansdr
#endif // LEANSDR_FRAMEWORK_H
#endif // LEANSDR_FRAMEWORK_H

View File

@ -1,14 +1,26 @@
// This file is part of LeanSDR Copyright (C) 2016-2018 <pabr@pabr.org>.
// See the toplevel README for more information.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
#ifndef LEANSDR_GENERIC_H
#define LEANSDR_GENERIC_H
#include <errno.h>
#include <fcntl.h>
#include <sys/types.h>
#ifdef _MSC_VER
#include <stdlib.h>
#include <io.h>
#else
#include <unistd.h>
#endif
#include "leansdr/math.h"
@ -22,11 +34,14 @@ namespace leansdr
// [file_reader] reads raw data from a file descriptor into a [pipebuf].
// If the file descriptor is seekable, data can be looped.
template<typename T>
struct file_reader: runnable
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)
file_reader(scheduler *sch, int _fdin, pipebuf<T> &_out)
: runnable(sch, _out.name),
loop(false),
filler(NULL),
fdin(_fdin), out(_out)
{
}
void run()
@ -35,32 +50,29 @@ struct file_reader: runnable
if (!size)
return;
#ifdef _MSC_VER
again: ssize_t nr = _read(fdin, out.wr(), size);
#else
again: ssize_t nr = read(fdin, out.wr(), size);
#endif
if (nr < 0)
again:
ssize_t nr = read(fdin, out.wr(), size);
if (nr < 0 && errno == EWOULDBLOCK)
{
fatal("leansdr::file_reader::run: read");
if (filler)
{
if (sch->debug)
fprintf(stderr, "U");
out.write(*filler);
}
return;
}
if (nr < 0)
fatal("read(file_reader)");
if (!nr)
{
if (!loop)
return;
if (sch->debug)
fprintf(stderr, "leansdr::file_reader::run: %s looping\n", name);
#ifdef _MSC_VER
off_t res = _lseek(fdin, 0, SEEK_SET);
#else
fprintf(stderr, "%s looping\n", name);
off_t res = lseek(fdin, 0, SEEK_SET);
#endif
if (res == (off_t) -1)
{
fatal("leansdr::file_reader::run: lseek");
return;
}
if (res == (off_t)-1)
fatal("lseek");
goto again;
}
@ -71,16 +83,9 @@ struct file_reader: runnable
{
if (sch->debug)
fprintf(stderr, "+");
#ifdef _MSC_VER
ssize_t nr2 = _read(fdin, (char*) out.wr() + nr, remain);
#else
ssize_t nr2 = read(fdin, (char*) out.wr() + nr, remain);
#endif
ssize_t nr2 = read(fdin, (char *)out.wr() + nr, remain);
if (nr2 <= 0)
{
fatal("leansdr::file_reader::run: partial read");
return;
}
fatal("partial read");
nr += nr2;
remain -= nr2;
}
@ -88,18 +93,27 @@ struct file_reader: runnable
out.written(nr / sizeof(T));
}
bool loop;
private:
void set_realtime(T &_filler)
{
int flags = fcntl(fdin, F_GETFL);
if (fcntl(fdin, F_SETFL, flags | O_NONBLOCK))
fatal("fcntl");
filler = new T(_filler);
}
private:
T *filler;
int fdin;
pipewriter<T> out;
};
// [file_writer] writes raw data from a [pipebuf] to a file descriptor.
template<typename T>
struct file_writer: runnable
template <typename T>
struct file_writer : runnable
{
file_writer(scheduler *sch, pipebuf<T> &_in, int _fdout) :
runnable(sch, _in.name), in(_in), fdout(_fdout)
file_writer(scheduler *sch, pipebuf<T> &_in, int _fdout) : runnable(sch, _in.name),
in(_in), fdout(_fdout)
{
}
void run()
@ -107,29 +121,17 @@ struct file_writer: runnable
int size = in.readable() * sizeof(T);
if (!size)
return;
#ifdef _MSC_VER
int nw = _write(fdout, in.rd(), size);
#else
int nw = write(fdout, in.rd(), size);
#endif
if (!nw)
{
fatal("leansdr::file_writer::run: pipe");
return;
}
fatal("pipe");
if (nw < 0)
{
fatal("leansdr::file_writer::run: write");
return;
}
fatal("write");
if (nw % sizeof(T))
{
fatal("leansdr::file_writer::run:partial write");
return;
}
fatal("partial write");
in.read(nw / sizeof(T));
}
private:
private:
pipereader<T> in;
int fdout;
};
@ -137,11 +139,14 @@ private:
// [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
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)
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()
@ -156,27 +161,18 @@ struct file_printer: runnable
char buf[256];
int len = snprintf(buf, sizeof(buf), format, (*pin) * scale);
if (len < 0)
{
fatal("leansdr::file_printer::run: obsolete glibc");
return;
}
#ifdef _MSC_VER
int nw = _write(fdout, buf, len);
#else
fatal("obsolete glibc");
int nw = write(fdout, buf, len);
#endif
if (nw != len)
{
fatal("leansdr::file_printer::run: partial write");
return;
}
fatal("partial write");
}
}
in.read(n);
}
T scale;
int decimation;
private:
private:
pipereader<T> in;
const char *format;
int fdout;
@ -187,11 +183,18 @@ private:
// to a file descriptor on a single line.
// Special case for complex.
template<typename T>
struct file_carrayprinter: runnable
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"))
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()
@ -218,32 +221,36 @@ struct file_carrayprinter: runnable
}
}
T scale;
int fixed_size; // Number of elements per batch, or 0.
private:
pipereader<complex<T> > in;
int fixed_size; // Number of elements per batch, or 0.
private:
pipereader<complex<T>> in;
const char *head, *format, *sep, *tail;
FILE *fout;
};
template<typename T, int N>
struct file_vectorprinter: runnable
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)
file_vectorprinter(scheduler *sch,
const char *_head,
const char *_format,
const char *_sep,
const char *_tail,
pipebuf<T[N]> &_in, int _fdout, int _n = N) : runnable(sch, _in.name), scale(1), in(_in),
head(_head), format(_format), sep(_sep), tail(_tail), n(_n)
{
fout = fdopen(_fdout, "w");
if (!fout)
{
fatal("leansdr::file_vectorprinter::file_vectorprinter: fdopen");
}
fatal("fdopen");
}
void run()
{
while (in.readable() >= 1)
{
fprintf(fout, head, N);
T (*pin)[N] = in.rd();
for (int i = 0; i < N; ++i)
fprintf(fout, head, n);
T(*pin)
[N] = in.rd();
for (int i = 0; i < n; ++i)
{
if (i)
fprintf(fout, "%s", sep);
@ -255,20 +262,23 @@ struct file_vectorprinter: runnable
fflush(fout);
}
T scale;
private:
private:
pipereader<T[N]> in;
const char *head, *format, *sep, *tail;
FILE *fout;
int n;
};
// [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
template <typename Tin, typename Tout>
struct itemcounter : runnable
{
itemcounter(scheduler *sch, pipebuf<Tin> &_in, pipebuf<Tout> &_out) :
runnable(sch, "itemcounter"), in(_in), out(_out)
itemcounter(scheduler *sch, pipebuf<Tin> &_in, pipebuf<Tout> &_out)
: runnable(sch, "itemcounter"),
in(_in), out(_out)
{
}
void run()
@ -281,20 +291,23 @@ struct itemcounter: runnable
out.write(count);
in.read(count);
}
private:
private:
pipereader<Tin> in;
pipewriter<Tout> out;
};
// [decimator] forwards 1 in N sample.
template<typename T>
struct decimator: runnable
template <typename T>
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)
decimator(scheduler *sch, int _d, pipebuf<T> &_in, pipebuf<T> &_out)
: runnable(sch, "decimator"),
d(_d),
in(_in), out(_out)
{
}
void run()
@ -306,7 +319,8 @@ struct decimator: runnable
in.read(count * d);
out.written(count);
}
private:
private:
pipereader<T> in;
pipewriter<T> out;
};
@ -314,13 +328,18 @@ private:
// [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)
{
}
@ -339,12 +358,12 @@ struct rate_estimator: runnable
den.read(count);
if (acc_den >= sample_size)
{
rate.write((float) acc_num / acc_den);
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;
@ -352,18 +371,16 @@ private:
// SERIALIZER
template<typename Tin, typename Tout>
struct serializer: runnable
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)
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)
{
(void) sch;
if (nin * sizeof(Tin) != nout * sizeof(Tout))
{
fail("serializer::serializer", "incompatible sizes");
return;
}
fail("serializer: incompatible sizes");
}
void run()
{
@ -374,61 +391,63 @@ struct serializer: runnable
out.written(nout);
}
}
private:
private:
int nin, nout;
pipereader<Tin> in;
pipewriter<Tout> out;
};
// serializer
}; // serializer
// [buffer_reader] reads from a user-supplied buffer.
template<typename T>
struct buffer_reader: runnable
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)
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));
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.
template<typename T>
struct buffer_writer: runnable
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)
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));
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 leansdr
#endif // LEANSDR_GENERIC_H
#endif // LEANSDR_GENERIC_H

File diff suppressed because it is too large Load Diff

View File

@ -1,3 +1,19 @@
// This file is part of LeanSDR Copyright (C) 2016-2018 <pabr@pabr.org>.
// See the toplevel README for more information.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
#ifndef LEANSDR_HDLC_H
#define LEANSDR_HDLC_H
@ -11,19 +27,14 @@ namespace leansdr
struct hdlc_dec
{
hdlc_dec(int _minframesize, // Including CRC, excluding HDLC flags.
int _maxframesize, bool _invert) :
minframesize(_minframesize),
maxframesize(_maxframesize),
invertmask(_invert ? 0xff : 0),
framebuf(new u8[maxframesize]),
debug(false)
hdlc_dec(int _minframesize, // Including CRC, excluding HDLC flags.
int _maxframesize,
bool _invert) : minframesize(_minframesize),
maxframesize(_maxframesize),
invertmask(_invert ? 0xff : 0),
framebuf(new u8[maxframesize]),
debug(false)
{
byte_out = 0;
nbits_out = 0;
framesize = 0;
crc16 = 0;
reset();
}
@ -45,24 +56,26 @@ struct hdlc_dec
// Return number of checksum errors in *fcs_errors.
// *ppin will have increased by at least 1 (unless count==0).
u8 *decode(u8 **ppin, int count, int *pdatasize, int *hdlc_errors,
int *fcs_errors)
u8 *decode(u8 **ppin, int count, int *pdatasize, int *hdlc_errors, int *fcs_errors)
{
*hdlc_errors = 0;
*fcs_errors = 0;
*pdatasize = -1;
u8 *pin = *ppin, *pend = pin + count;
for (; pin < pend; ++pin)
{
u8 byte_in = (*pin) ^ invertmask;
for (int bits = 8; bits--; byte_in <<= 1)
{
u8 bit_in = byte_in & 128;
shiftreg = (shiftreg >> 1) | bit_in;
if (!inframe)
{
if (shiftreg == 0x7e)
{ // HDLC flag 01111110
{ // HDLC flag 01111110
inframe = true;
nbits_out = 0;
begin_frame();
@ -71,11 +84,11 @@ struct hdlc_dec
else
{
if ((shiftreg & 0xfe) == 0x7c)
{ // 0111110x HDLC stuffing
{ // 0111110x HDLC stuffing
// Unstuff this 0
}
else if (shiftreg == 0x7e)
{ // 01111110 HDLC flag
{ // 01111110 HDLC flag
if (nbits_out != 7)
{
// Not at byte boundary
@ -87,8 +100,7 @@ struct hdlc_dec
{
// Checksum
crc16 ^= 0xffff;
if (framesize < 2 || framesize < minframesize
|| crc16 != crc16_check)
if (framesize < 2 || framesize < minframesize || crc16 != crc16_check)
{
if (debug)
fprintf(stderr, "!");
@ -111,7 +123,7 @@ struct hdlc_dec
// Special cases 0111111 and 1111111 cannot affect *pdatasize.
}
else if (shiftreg == 0xfe)
{ // 11111110 HDLC invalid
{ // 11111110 HDLC invalid
if (framesize)
{
if (debug)
@ -121,7 +133,7 @@ struct hdlc_dec
inframe = false;
}
else
{ // Data bit
{ // Data bit
byte_out = (byte_out >> 1) | bit_in; // HDLC is LSB first
++nbits_out;
if (nbits_out == 8)
@ -134,8 +146,8 @@ struct hdlc_dec
nbits_out = 0;
}
}
} // inframe
} // bits
} // inframe
} // bits
if (*pdatasize != -1)
{
// Found a complete frame
@ -143,88 +155,98 @@ struct hdlc_dec
return framebuf;
}
}
*ppin = pin;
return NULL;
}
private:
private:
// Config
int minframesize, maxframesize;
u8 invertmask;
u8 *framebuf; // [maxframesize]
u8 *framebuf; // [maxframesize]
// State
u8 shiftreg; // Input bitstream
bool inframe; // Currently receiving a frame ?
u8 byte_out; // Accumulator for data bits
int nbits_out; // Number of data bits in byte_out
int framesize; // Number of bytes in framebuf, if inframe
u16 crc16; // CRC of framebuf[framesize]
u8 shiftreg; // Input bitstream
bool inframe; // Currently receiving a frame ?
u8 byte_out; // Accumulator for data bits
int nbits_out; // Number of data bits in byte_out
int framesize; // Number of bytes in framebuf, if inframe
u16 crc16; // CRC of framebuf[framesize]
// CRC
static const u16 crc16_init = 0xffff;
static const u16 crc16_poly = 0x8408; // 0x1021 MSB-first
static const u16 crc16_poly = 0x8408; // 0x1021 MSB-first
static const u16 crc16_check = 0x0f47;
void crc16_byte(u8 data)
{
crc16 ^= data;
for (int bit = 8; bit--;)
{
crc16 = (crc16 & 1) ? (crc16 >> 1) ^ crc16_poly : (crc16 >> 1);
}
}
public:
public:
bool debug;
};
// hdlc_dec
// HDLC synchronizer with polarity detection
struct hdlc_sync: runnable
struct hdlc_sync : runnable
{
hdlc_sync(scheduler *sch,
pipebuf<u8> &_in, // Packed bits
pipebuf<u8> &_out, // Bytes
int _minframesize, // Including CRC, excluding HDLC flags.
int _maxframesize,
// Status
pipebuf<int> *_lock_out = NULL,
pipebuf<int> *_framecount_out = NULL,
pipebuf<int> *_fcserrcount_out = NULL,
pipebuf<int> *_hdlcbytecount_out = NULL,
pipebuf<int> *_databytecount_out = NULL) :
runnable(sch, "hdlc_sync"), minframesize(_minframesize), maxframesize(
_maxframesize), chunk_size(maxframesize + 2), in(_in), out(
_out, _maxframesize + chunk_size), lock_out(
opt_writer(_lock_out)), framecount_out(
opt_writer(_framecount_out)), fcserrcount_out(
opt_writer(_fcserrcount_out)), hdlcbytecount_out(
opt_writer(_hdlcbytecount_out)), databytecount_out(
opt_writer(_databytecount_out)), cur_sync(0), resync_phase(
0), lock_state(false), resync_period(32), header16(false)
hdlc_sync(scheduler *sch, pipebuf<u8> &_in, // Packed bits
pipebuf<u8> &_out, // Bytes
int _minframesize, // Including CRC, excluding HDLC flags.
int _maxframesize,
// Status
pipebuf<int> *_lock_out = NULL,
pipebuf<int> *_framecount_out = NULL,
pipebuf<int> *_fcserrcount_out = NULL,
pipebuf<int> *_hdlcbytecount_out = NULL,
pipebuf<int> *_databytecount_out = NULL) : runnable(sch, "hdlc_sync"),
minframesize(_minframesize),
maxframesize(_maxframesize),
chunk_size(maxframesize + 2),
in(_in),
out(_out, _maxframesize + chunk_size),
lock_out(opt_writer(_lock_out)),
framecount_out(opt_writer(_framecount_out)),
fcserrcount_out(opt_writer(_fcserrcount_out)),
hdlcbytecount_out(opt_writer(_hdlcbytecount_out)),
databytecount_out(opt_writer(_databytecount_out)),
cur_sync(0),
resync_phase(0),
lock_state(false),
resync_period(32),
header16(false)
{
for (int s = 0; s < NSYNCS; ++s)
{
syncs[s].dec = new hdlc_dec(minframesize, maxframesize, s != 0);
for (int h = 0; h < NERRHIST; ++h)
syncs[s].errhist[h] = 0;
}
syncs[cur_sync].dec->debug = sch->debug;
errslot = 0;
}
void run()
{
if (!opt_writable(lock_out) || !opt_writable(framecount_out)
|| !opt_writable(fcserrcount_out)
|| !opt_writable(hdlcbytecount_out)
|| !opt_writable(databytecount_out))
if (!opt_writable(lock_out) || !opt_writable(framecount_out) || !opt_writable(fcserrcount_out) || !opt_writable(hdlcbytecount_out) || !opt_writable(databytecount_out))
{
return;
}
bool previous_lock_state = lock_state;
int fcserrcount = 0, framecount = 0;
int hdlcbytecount = 0, databytecount = 0;
// Note: hdlc_dec may already hold one frame ready for output.
while ((long) in.readable() >= chunk_size
&& (long) out.writable() >= maxframesize + chunk_size)
while ((long)in.readable() >= chunk_size && (long)out.writable() >= maxframesize + chunk_size)
{
if (!resync_phase)
{
@ -233,14 +255,15 @@ struct hdlc_sync: runnable
{
if (s != cur_sync)
syncs[s].dec->reset();
syncs[s].errhist[errslot] = 0;
for (u8 *pin = in.rd(), *pend = pin + chunk_size;
pin < pend;)
for (u8 *pin = in.rd(), *pend = pin + chunk_size; pin < pend;)
{
int datasize, hdlc_errors, fcs_errors;
u8 *f = syncs[s].dec->decode(&pin, pend - pin,
&datasize, &hdlc_errors, &fcs_errors);
u8 *f = syncs[s].dec->decode(&pin, pend - pin, &datasize, &hdlc_errors, &fcs_errors);
syncs[s].errhist[errslot] += hdlc_errors;
if (s == cur_sync)
{
if (f)
@ -250,32 +273,39 @@ struct hdlc_sync: runnable
databytecount += datasize;
++framecount;
}
fcserrcount += fcs_errors;
framecount += fcs_errors;
}
}
}
errslot = (errslot + 1) % NERRHIST;
// Switch to another sync option ?
// Compare total error counts over about NERRHIST frames.
int total_errors[NSYNCS];
for (int s = 0; s < NSYNCS; ++s)
{
total_errors[s] = 0;
for (int h = 0; h < NERRHIST; ++h)
total_errors[s] += syncs[s].errhist[h];
}
int best = cur_sync;
for (int s = 0; s < NSYNCS; ++s)
if (total_errors[s] < total_errors[best])
best = s;
if (best != cur_sync)
{
lock_state = false;
if (sch->debug)
fprintf(stderr, "[%d:%d->%d:%d]", cur_sync,
total_errors[cur_sync], best,
total_errors[best]);
fprintf(stderr, "[%d:%d->%d:%d]", cur_sync, total_errors[cur_sync], best, total_errors[best]);
// No verbose messages on candidate syncs
syncs[cur_sync].dec->debug = false;
cur_sync = best;
@ -288,8 +318,8 @@ struct hdlc_sync: runnable
for (u8 *pin = in.rd(), *pend = pin + chunk_size; pin < pend;)
{
int datasize, hdlc_errors, fcs_errors;
u8 *f = syncs[cur_sync].dec->decode(&pin, pend - pin,
&datasize, &hdlc_errors, &fcs_errors);
u8 *f = syncs[cur_sync].dec->decode(&pin, pend - pin, &datasize, &hdlc_errors, &fcs_errors);
if (f)
{
lock_state = true;
@ -297,25 +327,29 @@ struct hdlc_sync: runnable
databytecount += datasize;
++framecount;
}
fcserrcount += fcs_errors;
framecount += fcs_errors;
}
} // resync_phase
} // resync_phase
in.read(chunk_size);
hdlcbytecount += chunk_size;
if (++resync_phase >= resync_period)
resync_phase = 0;
} // Work to do
} // Work to do
if (lock_state != previous_lock_state)
opt_write(lock_out, lock_state ? 1 : 0);
opt_write(framecount_out, framecount);
opt_write(fcserrcount_out, fcserrcount);
opt_write(hdlcbytecount_out, hdlcbytecount);
opt_write(databytecount_out, databytecount);
}
private:
private:
void output_frame(u8 *f, int size)
{
if (header16)
@ -324,6 +358,7 @@ private:
out.write(size >> 8);
out.write(size & 255);
}
memcpy(out.wr(), f, size);
out.written(size);
opt_write(framecount_out, 1);
@ -336,23 +371,26 @@ private:
pipewriter<int> *lock_out;
pipewriter<int> *framecount_out, *fcserrcount_out;
pipewriter<int> *hdlcbytecount_out, *databytecount_out;
static const int NSYNCS = 2; // Two possible polarities
static const int NERRHIST = 2; // Compare error counts over two frames
static const int NSYNCS = 2; // Two possible polarities
static const int NERRHIST = 2; // Compare error counts over two frames
struct
{
hdlc_dec *dec;
int errhist[NERRHIST];
} syncs[NSYNCS];
int errslot;
int cur_sync;
int resync_phase;
bool lock_state;
public:
public:
int resync_period;
bool header16; // Output length prefix
bool header16; // Output length prefix
};
// hdlc_sync
}// namespace
} // namespace leansdr
#endif // LEANSDR_HDLC_H
#endif // LEANSDR_HDLC_H

View File

@ -1,58 +1,79 @@
// This file is part of LeanSDR Copyright (C) 2016-2018 <pabr@pabr.org>.
// See the toplevel README for more information.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
#ifndef LEANSDR_IESS_H
#define LEANSDR_IESS_H
#include "leansdr/framework.h"
namespace leansdr {
namespace leansdr
{
// SELF-SYNCHRONIZING DESCRAMBLER
// Per ETSI TR 192 figure 8 (except Q20/ not connected to CLOCK).
// This implementation operates on packed bits, MSB first.
// SELF-SYNCHRONIZING DESCRAMBLER
// Per ETSI TR 192 figure 8 (except Q20/ not connected to CLOCK).
// This implementation operates on packed bits, MSB first.
struct etr192_descrambler : runnable {
struct etr192_descrambler : runnable
{
etr192_descrambler(scheduler *sch,
pipebuf<u8> &_in, // Packed scrambled bits
pipebuf<u8> &_out) // Packed bits
: runnable(sch, "etr192_dec"),
in(_in), out(_out),
shiftreg(0), counter(0)
pipebuf<u8> &_in, // Packed scrambled bits
pipebuf<u8> &_out) // Packed bits
: runnable(sch, "etr192_dec"),
in(_in), out(_out),
shiftreg(0), counter(0)
{
}
void run() {
int count = min(in.readable(), out.writable());
for ( u8 *pin=in.rd(), *pend=pin+count, *pout=out.wr();
pin<pend; ++pin,++pout) {
u8 byte_in=*pin, byte_out=0;
for ( int b=8; b--; byte_in<<=1 ) {
// Levels before clock transition
int bit_in = (byte_in&128) ? 1 : 0;
int reset_counter = (shiftreg ^ (shiftreg>>8)) & 1;
int counter_overflow = (counter==31) ? 1 : 0;
int taps = (shiftreg>>2) ^ (shiftreg>>19);
int bit_out = (taps ^ counter_overflow ^ bit_in ^ 1) & 1;
// Execute clock transition
#if 1 // Descramble
shiftreg = (shiftreg<<1) | bit_in;
#else // Scramble
shiftreg = (shiftreg<<1) | bit_out;
void run()
{
int count = min(in.readable(), out.writable());
for (u8 *pin = in.rd(), *pend = pin + count, *pout = out.wr();
pin < pend; ++pin, ++pout)
{
u8 byte_in = *pin, byte_out = 0;
for (int b = 8; b--; byte_in <<= 1)
{
// Levels before clock transition
int bit_in = (byte_in & 128) ? 1 : 0;
int reset_counter = (shiftreg ^ (shiftreg >> 8)) & 1;
int counter_overflow = (counter == 31) ? 1 : 0;
int taps = (shiftreg >> 2) ^ (shiftreg >> 19);
int bit_out = (taps ^ counter_overflow ^ bit_in ^ 1) & 1;
// Execute clock transition
#if 1 // Descramble
shiftreg = (shiftreg << 1) | bit_in;
#else // Scramble
shiftreg = (shiftreg << 1) | bit_out;
#endif
counter = reset_counter ? 0 : (counter+1)&31;
byte_out = (byte_out<<1) | bit_out;
}
*pout = byte_out;
}
in.read(count);
out.written(count);
counter = reset_counter ? 0 : (counter + 1) & 31;
byte_out = (byte_out << 1) | bit_out;
}
*pout = byte_out;
}
in.read(count);
out.written(count);
}
private:
pipereader<u8> in;
pipewriter<u8> out;
u32 shiftreg; // 20 bits
u8 counter; // 5 bits
}; // etr192_descrambler
u32 shiftreg; // 20 bits
u8 counter; // 5 bits
}; // etr192_descrambler
} // namespace
} // namespace leansdr
#endif // LEANSDR_IESS_H
#endif // LEANSDR_IESS_H

View File

@ -1,64 +0,0 @@
///////////////////////////////////////////////////////////////////////////////////
// Copyright (C) 2016 Edouard Griffiths, F4EXB //
// //
// This program is free software; you can redistribute it and/or modify //
// it under the terms of the GNU General Public License as published by //
// the Free Software Foundation as version 3 of the License, or //
// //
// This program is distributed in the hope that it will be useful, //
// but WITHOUT ANY WARRANTY; without even the implied warranty of //
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the //
// GNU General Public License V3 for more details. //
// //
// You should have received a copy of the GNU General Public License //
// along with this program. If not, see <http://www.gnu.org/licenses/>. //
///////////////////////////////////////////////////////////////////////////////////
#ifndef SDRBASE_UTIL_INCREMENTALARRAY_H_
#define SDRBASE_UTIL_INCREMENTALARRAY_H_
#include <stdint.h>
namespace leansdr
{
template<typename T>
class IncrementalArray
{
public:
T *m_array;
IncrementalArray();
~IncrementalArray();
void allocate(uint32_t size);
private:
uint32_t m_size;
};
template<typename T>
IncrementalArray<T>::IncrementalArray() :
m_array(0),
m_size(0)
{
}
template<typename T>
IncrementalArray<T>::~IncrementalArray()
{
if (m_array) { delete[] m_array; }
}
template<typename T>
void IncrementalArray<T>::allocate(uint32_t size)
{
if (size <= m_size) { return; }
if (m_array) { delete[] m_array; }
m_array = new T[size];
m_size = size;
}
} // namespace
#endif /* SDRBASE_UTIL_INCREMENTALARRAY_H_ */

View File

@ -0,0 +1,520 @@
// This file is part of LeanSDR Copyright (C) 2016-2018 <pabr@pabr.org>.
// See the toplevel README for more information.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
#ifndef LEANSDR_LDPC_H
#define LEANSDR_LDPC_H
#define lfprintf(...) \
{ \
}
namespace leansdr
{
// LDPC sparse matrix specified like in the DVB-S2 standard
// Taddr must be wide enough to index message bits and address bits.
template <typename Taddr>
struct ldpc_table
{
// TBD Save space
static const int MAX_ROWS = 162; // 64800 * (9/10) / 360
static const int MAX_COLS = 13;
int q;
int nrows;
struct row
{
int ncols;
Taddr cols[MAX_COLS];
} rows[MAX_ROWS];
};
// LDPC ENGINE
// SOFTBITs can be hard (e.g. bool) or soft (e.g. llr_t).
// They are stored as SOFTWORDs containing SWSIZE SOFTBITs.
// See interface in softword.h.
template <typename SOFTBIT, typename SOFTWORD, int SWSIZE, typename Taddr>
struct ldpc_engine
{
ldpc_engine()
: vnodes(NULL), cnodes(NULL)
{
}
// vnodes: Value/variable nodes (message bits)
// cnodes: Check nodes (parity bits)
int k; // Message size in bits
int n; // Codeword size in bits
struct node
{
Taddr *edges;
int nedges;
static const int CHUNK = 4; // Grow edges[] in steps of CHUNK.
void append(Taddr a)
{
if (nedges % CHUNK == 0)
{ // Full ?
edges = (Taddr *)realloc(edges, (nedges + CHUNK) * sizeof(Taddr));
if (!edges)
fatal("realloc");
}
edges[nedges++] = a;
}
};
node *vnodes; // [k]
node *cnodes; // [n-k]
// Initialize from a S2-style table.
ldpc_engine(const ldpc_table<Taddr> *table, int _k, int _n)
: k(_k), n(_n)
{
// Sanity checks
if (360 % SWSIZE)
fatal("Bad LDPC word size");
if (k % SWSIZE)
fatal("Bad LDPC k");
if (n % SWSIZE)
fatal("Bad LDPC n");
if (k != table->nrows * 360)
fatal("Bad table");
int n_k = n - k;
if (table->q * 360 != n_k)
fatal("Bad q");
vnodes = new node[k];
memset(vnodes, 0, sizeof(node) * k);
cnodes = new node[n_k];
memset(cnodes, 0, sizeof(node) * n_k);
// Expand the graph.
int m = 0;
// Iterate over rows
for (const typename ldpc_table<Taddr>::row *prow = table->rows;
prow < table->rows + table->nrows;
++prow)
{
// Process 360 bits per row.
int q = table->q;
int qoffs = 0;
for (int mw = 360; mw--; ++m, qoffs += q)
{
const Taddr *pa = prow->cols;
for (int nc = prow->ncols; nc--; ++pa)
{
int a = (int)*pa + qoffs;
if (a >= n_k)
a -= n_k; // Modulo n-k. Note qoffs<360*q.
if (a >= n_k)
fail("Invalid LDPC table");
vnodes[m].append(a);
cnodes[a].append(m);
}
}
}
}
void print_node_stats()
{
int nedges = count_edges(vnodes, k);
fprintf(stderr, "LDPC(%5d,%5d)(%.2f)"
" %5.2f edges/vnode, %5.2f edges/cnode\n",
k, n - k, (float)k / n, (float)nedges / k, (float)nedges / (n - k));
}
int count_edges(node *nodes, int nnodes)
{
int c = 0;
for (int i = 0; i < nnodes; ++i)
c += nodes[i].nedges;
return c;
}
// k: Message size in bits
// n: Codeword size in bits
// integrate: Optional S2-style post-processing
#if 0
void encode_hard(const ldpc_table<Taddr> *table, const uint8_t *msg,
int k, int n, uint8_t *parity, bool integrate=true) {
// Sanity checks
if ( 360 % SWSIZE ) fatal("Bad LDPC word size");
if ( k % SWSIZE ) fatal("Bad LDPC k");
if ( n % SWSIZE ) fatal("Bad LDPC n");
if ( k != table->nrows*360 ) fatal("Bad table");
int n_k = n - k;
if ( table->q*360 != n_k ) fatal("Bad q");
for ( int i=0; i<n_k/SWSIZE; ++i ) softword_zero(&parity[i]);
// Iterate over rows
for ( const typename ldpc_table<Taddr>::row *prow = table->rows; // quirk
prow < table->rows+table->nrows;
++prow ) {
// Process 360 bits per row, in words of SWSIZE bits
int q = table->q;
int qoffs = 0;
for ( int mw=360/SWSIZE; mw--; ++msg ) {
SOFTWORD msgword = *msg;
for ( int wbit=0; wbit<SWSIZE; ++wbit,qoffs+=q ) {
SOFTBIT msgbit = softword_get(msgword, wbit);
if ( ! msgbit ) continue; // TBD Generic soft version
const Taddr *pa = prow->cols;
for ( int nc=prow->ncols; nc--; ++pa ) {
// Don't wrap modulo range of Taddr
int a = (int)*pa + qoffs;
// Note: qoffs < 360*q=n-k
if ( a >= n_k ) a -= n_k; // TBD not predictable
softwords_flip(parity, a);
}
}
}
}
if ( integrate )
integrate_bits(parity, parity, n_k/SWSIZE);
}
#endif
void encode(const ldpc_table<Taddr> *table, const SOFTWORD *msg,
int k, int n, SOFTWORD *parity, int integrate = true)
{
// Sanity checks
if (360 % SWSIZE)
fatal("Bad LDPC word size");
if (k % SWSIZE)
fatal("Bad LDPC k");
if (n % SWSIZE)
fatal("Bad LDPC n");
if (k != table->nrows * 360)
fatal("Bad table");
int n_k = n - k;
if (table->q * 360 != n_k)
fatal("Bad q");
for (int i = 0; i < n_k / SWSIZE; ++i)
softword_zero(&parity[i]);
// Iterate over rows
for (const typename ldpc_table<Taddr>::row *prow = table->rows; // quirk
prow < table->rows + table->nrows;
++prow)
{
// Process 360 bits per row, in words of SWSIZE bits
int q = table->q;
int qoffs = 0;
for (int mw = 360 / SWSIZE; mw--; ++msg)
{
SOFTWORD msgword = *msg;
for (int wbit = 0; wbit < SWSIZE; ++wbit, qoffs += q)
{
SOFTBIT msgbit = softword_get(msgword, wbit);
if (!softbit_harden(msgbit))
continue;
const Taddr *pa = prow->cols;
for (int nc = prow->ncols; nc--; ++pa)
{
int a = (int)*pa + qoffs;
// Note: qoffs < 360*q=n-k
if (a >= n_k)
a -= n_k; // TBD not predictable
softwords_flip(parity, a);
}
}
}
}
if (integrate)
integrate_bits(parity, parity, n_k / SWSIZE);
}
// Flip bits connected to parity errors, one at a time,
// as long as things improve and max_bitflips is not exceeded.
// cw: codeword (k value bits followed by n-k check bits)
static const int PPCM = 39;
typedef int64_t score_t;
score_t compute_scores(SOFTWORD *m, SOFTWORD *p, SOFTWORD *q, int nc,
score_t *score, int k)
{
int total = 0;
memset(score, 0, k * sizeof(*score));
for (int c = 0; c < nc; ++c)
{
SOFTBIT err = softwords_xor(p, q, c);
if (softbit_harden(err))
{
Taddr *pe = cnodes[c].edges;
for (int e = cnodes[c].nedges; e--; ++pe)
{
int v = *pe;
int s = err * softwords_weight<SOFTBIT, SOFTWORD>(m, v) * PPCM / vnodes[v].nedges;
//fprintf(stderr, "c[%d] bad => v[%d] += %d (%d*%d)\n",
///c, v, s, err, softwords_weight<SOFTBIT,SOFTWORD>(m,*pe));
score[v] += s;
total += s;
}
}
}
return total;
}
int decode_bitflip(const ldpc_table<Taddr> *table, SOFTWORD *cw,
int k, int n,
int max_bitflips)
{
if (!vnodes)
fail("LDPC graph not initialized");
int n_k = n - k;
// Compute the expected check bits (without the final mixing)
SOFTWORD expected[n_k / SWSIZE];
encode(table, cw, k, n, expected, false);
// Reverse the integrator mixing from the received check bits
SOFTWORD received[n_k / SWSIZE];
diff_bits(cw + k / SWSIZE, received, n_k / SWSIZE);
// Compute initial scores
score_t score[k];
score_t tots = compute_scores(cw, expected, received, n_k, score, k);
lfprintf(stderr, "Initial score %d\n", (int)tots);
int nflipped = 0;
score_t score_threshold;
{
SOFTBIT one;
softbit_set(&one, true);
score_threshold = (int)one * 2;
}
bool progress = true;
while (progress && nflipped < max_bitflips)
{
progress = false;
// Try to flip parity bits.
// Due to differential decoding, they appear as consecutive errors.
SOFTBIT prev_err = softwords_xor(expected, received, 0);
for (int b = 0; b < n - k - 1; ++b)
{
prev_err = softwords_xor(expected, received, b); //TBD
SOFTBIT err = softwords_xor(expected, received, b + 1);
if (softbit_harden(prev_err) && softbit_harden(err))
{
lfprintf(stderr, "flip parity %d\n", b);
softwords_flip(received, b);
softwords_flip(received, b + 1);
++nflipped; // Counts as one flip before differential decoding.
progress = true;
int dtot = 0;
// Depenalize adjacent message bits.
{
Taddr *pe = cnodes[b].edges;
for (int e = cnodes[b].nedges; e--; ++pe)
{
int d = prev_err * softwords_weight<SOFTBIT, SOFTWORD>(cw, *pe) * PPCM / vnodes[*pe].nedges;
score[*pe] -= d;
dtot -= d;
}
}
{
Taddr *pe = cnodes[b + 1].edges;
for (int e = cnodes[b + 1].nedges; e--; ++pe)
{
int d = err * softwords_weight<SOFTBIT, SOFTWORD>(cw, *pe) * PPCM / vnodes[*pe].nedges;
score[*pe] -= d;
dtot -= d;
}
}
tots += dtot;
#if 1
// Also update the codeword in-place.
// TBD Useful for debugging only.
softwords_flip(cw, k + b);
#endif
// TBD version soft. err = ! err;
}
prev_err = err;
} // c nodes
score_t maxs = -(1 << 30);
for (int v = 0; v < k; ++v)
if (score[v] > maxs)
maxs = score[v];
if (!maxs)
break;
lfprintf(stderr, "maxs %d\n", (int)maxs);
// Try to flip each message bits with maximal score
for (int v = 0; v < k; ++v)
{
if (score[v] < score_threshold)
continue;
// if ( score[v] < maxs*9/10 ) continue;
if (score[v] < maxs - 4)
continue;
lfprintf(stderr, " flip %d score=%d\n", (int)v, (int)score[v]);
// Update expected parities and scores that depend on them.
score_t dtot = 0;
for (int commit = 0; commit <= 1; ++commit)
{
Taddr *pe = vnodes[v].edges;
for (int e = vnodes[v].nedges; e--; ++pe)
{
Taddr c = *pe;
SOFTBIT was_bad = softwords_xor(expected, received, c);
if (softbit_harden(was_bad))
{
Taddr *pe = cnodes[c].edges;
for (int e = cnodes[c].nedges; e--; ++pe)
{
int d = was_bad * softwords_weight<SOFTBIT, SOFTWORD>(cw, *pe) * PPCM / vnodes[*pe].nedges;
if (commit)
score[*pe] -= d;
else
dtot -= d;
}
}
softwords_flip(expected, c);
SOFTBIT is_bad = softwords_xor(expected, received, c);
if (softbit_harden(is_bad))
{
Taddr *pe = cnodes[c].edges;
for (int e = cnodes[c].nedges; e--; ++pe)
{
int d = is_bad * softwords_weight<SOFTBIT, SOFTWORD>(cw, *pe) * PPCM / vnodes[*pe].nedges;
if (commit)
score[*pe] += d;
else
dtot += d;
}
}
if (!commit)
softwords_flip(expected, c);
}
if (!commit)
{
if (dtot >= 0)
{
lfprintf(stderr, " abort %d\n", v);
break; // Next v
}
}
else
{
softwords_flip(cw, v);
++nflipped;
tots += dtot;
progress = true;
v = k - 1; // Force exit to update maxs ?
}
} // commit
} // v
lfprintf(stderr, "progress %d\n", progress);
#if 0
fprintf(stderr, "CHECKING TOTS INCREMENT (slow) %d\n", tots);
score_t tots2 = compute_scores(cw, expected, received, n_k, score, k);
if ( tots2 != tots ) fail("bad tots update");
#endif
}
return nflipped;
}
// EN 302 307-1 5.3.2.1 post-processing of parity bits.
// In-place allowed.
#if 1
static void integrate_bits(const SOFTWORD *in, SOFTWORD *out, int nwords)
{
SOFTBIT sum;
softbit_clear(&sum);
for (int i = 0; i < nwords; ++i)
{
SOFTWORD w = in[i];
for (int b = 0; b < SWSIZE; ++b)
{
sum = softbit_xor(sum, softword_get(w, b));
softword_write(w, b, sum);
}
out[i] = w;
}
}
#else
// Optimized for hard_sb
static void integrate_bits(const uint8_t *in, uint8_t *out, int nwords)
{
// TBD Optimize
uint8_t prev = 0;
for (int i = 0; i < nwords; ++i)
{
uint8_t c = in[i];
for (int j = SWSIZE; j--;)
{
c ^= prev << j;
prev = (c >> j) & 1;
}
out[i] = c;
}
}
#endif
// Undo EN 302 307-1 5.3.2.1, post-processing of parity bits.
// In-place allowed.
#if 1
static void diff_bits(const SOFTWORD *in, SOFTWORD *out, int nwords)
{
SOFTBIT prev;
softbit_clear(&prev);
for (int i = 0; i < nwords; ++i, ++in, ++out)
{
SOFTWORD w = *in;
for (int b = 0; b < SWSIZE; ++b)
{
SOFTBIT n = softword_get(w, b);
softword_write(w, b, softbit_xor(prev, n));
prev = n;
}
*out = w;
}
}
#else
// Optimized for hard_sb
static void diff_bits(const uint8_t *in, uint8_t *out, int nwords)
{
uint8_t prev = 0;
for (int i = 0; i < nwords; ++i)
{
uint8_t c = in[i];
out[i] = c ^ (prev | (c >> 1));
prev = (c & 1) << (SWSIZE - 1);
}
}
#endif
}; // ldpc_engine
} // namespace leansdr
#endif // LEANSDR_LDPC_H

View File

@ -0,0 +1,56 @@
#include "math.h"
namespace leansdr
{
int hamming_weight(uint8_t x)
{
static const int lut[16] = {0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4};
return lut[x & 15] + lut[x >> 4];
}
int hamming_weight(uint16_t x)
{
return hamming_weight((uint8_t)x) + hamming_weight((uint8_t)(x >> 8));
}
int hamming_weight(uint32_t x)
{
return hamming_weight((uint16_t)x) + hamming_weight((uint16_t)(x >> 16));
}
int hamming_weight(uint64_t x)
{
return hamming_weight((uint32_t)x) + hamming_weight((uint32_t)(x >> 32));
}
unsigned char parity(uint8_t x)
{
x ^= x >> 4;
return (0x6996 >> (x & 15)) & 1; // 16-entry look-up table
}
unsigned char parity(uint16_t x)
{
return parity((uint8_t)(x ^ (x >> 8)));
}
unsigned char parity(uint32_t x)
{
return parity((uint16_t)(x ^ (x >> 16)));
}
unsigned char parity(uint64_t x)
{
return parity((uint32_t)(x ^ (x >> 32)));
}
int log2i(uint64_t x)
{
int n = -1;
for (; x; ++n, x >>= 1)
;
return n;
}
} // leansdr

View File

@ -1,99 +1,155 @@
// This file is part of LeanSDR Copyright (C) 2016-2018 <pabr@pabr.org>.
// See the toplevel README for more information.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
#ifndef LEANSDR_MATH_H
#define LEANSDR_MATH_H
#define _USE_MATH_DEFINES
#include <math.h>
#include <stdint.h>
namespace leansdr {
namespace leansdr
{
template<typename T>
struct complex {
template <typename T>
struct complex
{
T re, im;
complex() : re(0), im(0) { }
complex(T x) : re(x), im(0) { }
complex(T x, T y) : re(x), im(y) { }
inline void operator +=(const complex<T> &x) { re+=x.re; im+=x.im; }
};
template<typename T>
complex<T> operator +(const complex<T> &a, const complex<T> &b) {
return complex<T>(a.re+b.re, a.im+b.im);
}
template<typename T>
complex<T> operator *(const complex<T> &a, const complex<T> &b) {
return complex<T>(a.re*b.re-a.im*b.im, a.re*b.im+a.im*b.re);
}
template<typename T>
complex<T> operator *(const complex<T> &a, const T &k) {
return complex<T>(a.re*k, a.im*k);
}
template<typename T>
complex<T> operator *(const T &k, const complex<T> &a) {
return complex<T>(k*a.re, k*a.im);
}
// TBD Optimize with dedicated instructions
inline int hamming_weight(uint8_t x) {
static const int lut[16] = { 0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4 };
return lut[x&15] + lut[x>>4];
}
inline int hamming_weight(uint16_t x) {
return hamming_weight((uint8_t)x)
+ hamming_weight((uint8_t)(x>>8));
}
inline int hamming_weight(uint32_t x) {
return hamming_weight((uint16_t)x)
+ hamming_weight((uint16_t)(x>>16));
}
inline int hamming_weight(uint64_t x) {
return hamming_weight((uint32_t)x)
+ hamming_weight((uint32_t)(x>>32));
}
inline unsigned char parity(uint8_t x) {
x ^= x>>4;
return (0x6996 >> (x&15)) & 1; // 16-entry look-up table
}
inline unsigned char parity(uint16_t x) {
return parity((uint8_t)(x^(x>>8)));
}
inline unsigned char parity(uint32_t x) {
return parity((uint16_t)(x^(x>>16)));
}
inline unsigned char parity(uint64_t x) {
return parity((uint32_t)(x^(x>>32)));
}
inline int log2i(uint64_t x) {
int n = -1;
for ( ; x; ++n,x>>=1 ) ;
return n;
}
// Pre-computed sin/cos for 16-bit angles
struct trig16 {
complex<float> lut[65536]; // TBD static and shared
trig16() {
for ( int a=0; a<65536; ++a ) {
float af = a * 2*M_PI / 65536;
lut[a].re = cosf(af);
lut[a].im = sinf(af);
}
complex() {}
complex(T x) : re(x), im(0) {}
complex(T x, T y) : re(x), im(y) {}
inline void operator+=(const complex<T> &x)
{
re += x.re;
im += x.im;
}
inline const complex<float> &expi(uint16_t a) const {
return lut[a];
inline void operator*=(const complex<T> &c)
{
T tre = re * c.re - im * c.im;
im = re * c.im + im * c.re;
re = tre;
}
inline void operator*=(const T &k)
{
re *= k;
im *= k;
}
};
template <typename T>
complex<T> operator+(const complex<T> &a, const complex<T> &b)
{
return complex<T>(a.re + b.re, a.im + b.im);
}
template <typename T>
complex<T> operator*(const complex<T> &a, const complex<T> &b)
{
return complex<T>(a.re * b.re - a.im * b.im, a.re * b.im + a.im * b.re);
}
template <typename T>
complex<T> operator*(const complex<T> &a, const T &k)
{
return complex<T>(a.re * k, a.im * k);
}
template <typename T>
complex<T> operator*(const T &k, const complex<T> &a)
{
return complex<T>(k * a.re, k * a.im);
}
template <typename T>
T dotprod(const T *u, const T *v, int n)
{
T acc = 0;
while (n--)
acc += (*u++) * (*v++);
return acc;
}
template <typename T>
inline T cnorm2(const complex<T> &u)
{
return u.re * u.re + u.im * u.im;
}
template <typename T>
T cnorm2(const complex<T> *p, int n)
{
T res = 0;
for (; n--; ++p)
res += cnorm2(*p);
return res;
}
// Return conj(u)*v
template <typename T>
inline complex<T> conjprod(const complex<T> &u, const complex<T> &v)
{
return complex<T>(u.re * v.re + u.im * v.im,
u.re * v.im - u.im * v.re);
}
// Return sum(conj(u[i])*v[i])
template <typename T>
complex<T> conjprod(const complex<T> *u, const complex<T> *v, int n)
{
complex<T> acc = 0;
while (n--)
acc += conjprod(*u++, *v++);
return acc;
}
// TBD Optimize with dedicated instructions
int hamming_weight(uint8_t x);
int hamming_weight(uint16_t x);
int hamming_weight(uint32_t x);
int hamming_weight(uint64_t x);
unsigned char parity(uint8_t x);
unsigned char parity(uint16_t x);
unsigned char parity(uint32_t x);
unsigned char parity(uint64_t x);
int log2i(uint64_t x);
// Pre-computed sin/cos for 16-bit angles
struct trig16
{
complex<float> lut[65536]; // TBD static and shared
trig16()
{
for (int a = 0; a < 65536; ++a)
{
float af = a * 2 * M_PI / 65536;
lut[a].re = cosf(af);
lut[a].im = sinf(af);
}
}
inline const complex<float> &expi(uint16_t a) const
{
return lut[a];
}
// a must fit in a int32_t, otherwise behaviour is undefined
inline const complex<float> &expi(float a) const {
return expi((uint16_t)(int16_t)(int32_t)a);
inline const complex<float> &expi(float a) const
{
return expi((uint16_t)(int16_t)(int32_t)a);
}
};
};
} // namespace
} // namespace leansdr
#endif // LEANSDR_MATH_H
#endif // LEANSDR_MATH_H

View File

@ -1,3 +1,19 @@
// This file is part of LeanSDR Copyright (C) 2016-2018 <pabr@pabr.org>.
// See the toplevel README for more information.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
#ifndef LEANSDR_RS_H
#define LEANSDR_RS_H
@ -28,16 +44,13 @@ namespace leansdr
// 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>
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;
}
fail("alpha!=2 not implemented");
// Precompute log and exp tables.
Tp alpha_i = 1;
for (Tp i = 0; i < (1 << N); ++i)
@ -45,20 +58,14 @@ struct gf2x_p
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
alpha_i <<= 1; // Multiply by alpha=[X] i.e. increase degrees
if (alpha_i & (1 << N))
alpha_i ^= P; // Modulo P iteratively
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 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)
@ -77,16 +84,11 @@ struct gf2x_p
// 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
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];
};
@ -98,14 +100,14 @@ struct rs_engine
// 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 }
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
G[i] = (i == 16) ? 1 : 0; // Init G=1
for (int d = 0; d < 16; ++d)
{
// Multiply by (X-alpha^d)
@ -115,7 +117,8 @@ struct rs_engine
}
#if DEBUG_RS
fprintf(stderr, "RS generator:");
for ( int i=0; i<=16; ++i ) fprintf(stderr, " %02x", G[i]);
for (int i = 0; i <= 16; ++i)
fprintf(stderr, " %02x", G[i]);
fprintf(stderr, "\n");
#endif
}
@ -162,12 +165,13 @@ struct rs_engine
{
// TBD Avoid copying
u8 p[204];
memcpy(p, msg, 204); // was 188 but causing underflow (PVS https://www.viva64.com/en/w/v512/)
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]);
for (int i = 0; i < 204; ++i)
fprintf(stderr, " %d", p[i]);
fprintf(stderr, "\n");
#endif
// Compute remainder modulo G
@ -183,7 +187,8 @@ struct rs_engine
}
#if DEBUG_RS
fprintf(stderr, "coded:");
for ( int i=0; i<204; ++i ) fprintf(stderr, " %d", p[i]);
for (int i = 0; i < 204; ++i)
fprintf(stderr, " %d", p[i]);
fprintf(stderr, "\n");
#endif
memcpy(msg + 188, p + 188, 16);
@ -193,14 +198,13 @@ struct rs_engine
// 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)
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 };
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;
@ -234,11 +238,13 @@ struct rs_engine
// 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, "[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]);
for (int i = 0; i < 16; ++i)
fprintf(stderr, " %d", synd[i]);
fprintf(stderr, "]\n");
#endif
@ -255,7 +261,8 @@ struct rs_engine
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]);
for (int i = 0; i < 16; ++i)
fprintf(stderr, " %d", omega[i]);
fprintf(stderr, "\n");
#endif
@ -265,7 +272,8 @@ struct rs_engine
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]);
for (int i = 0; i < 15; ++i)
fprintf(stderr, " %d", Cprime[i]);
fprintf(stderr, "\n");
#endif
@ -274,15 +282,15 @@ struct rs_engine
int roots_found = 0;
for (int i = 0; i < 255; ++i)
{
u8 r = gf.exp(i); // Candidate root alpha^0..alpha^254
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)
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)
{
@ -308,9 +316,8 @@ struct rs_engine
else
return false;
}
};
} // namespace
} // namespace leansdr
#endif // LEANSDR_RS_H
#endif // LEANSDR_RS_H

View File

@ -0,0 +1,98 @@
#include "sdr.h"
namespace leansdr
{
const char *cstln_base::names[] =
{
[BPSK] = "BPSK",
[QPSK] = "QPSK",
[PSK8] = "8PSK",
[APSK16] = "16APSK",
[APSK32] = "32APSK",
[APSK64E] = "64APSKe",
[QAM16] = "16QAM",
[QAM64] = "64QAM",
[QAM256] = "256QAM"
};
void softsymb_harden(llr_ss *ss)
{
for (int b = 0; b < 8; ++b)
ss->bits[b] = (ss->bits[b] < 0) ? -127 : 127;
}
void softsymb_harden(hard_ss *ss)
{
(void)ss; // NOP
}
void softsymb_harden(eucl_ss *ss)
{
for (int s = 0; s < ss->MAX_SYMBOLS; ++s)
ss->dists2[s] = (s == ss->nearest) ? 0 : 1;
}
uint8_t softsymb_to_dump(const llr_ss &ss, int bit)
{
return 128 - ss.bits[bit];
}
uint8_t softsymb_to_dump(const hard_ss &ss, int bit)
{
return ((ss >> bit) & 1) ? 255 : 0;
}
uint8_t softsymb_to_dump(const eucl_ss &ss, int bit)
{
(void)bit;
return ss.dists2[ss.nearest] >> 8;
}
void to_softsymb(const full_ss *fss, hard_ss *ss)
{
*ss = fss->nearest;
}
void to_softsymb(const full_ss *fss, eucl_ss *ss)
{
for (int s = 0; s < ss->MAX_SYMBOLS; ++s)
ss->dists2[s] = fss->dists2[s];
uint16_t best = 65535, best2 = 65535;
for (int s = 0; s < ss->MAX_SYMBOLS; ++s)
{
if (fss->dists2[s] < best)
{
best2 = best;
best = fss->dists2[s];
}
else if (fss->dists2[s] < best2)
{
best2 = fss->dists2[s];
}
}
ss->discr2 = best2 - best;
ss->nearest = fss->nearest;
}
void to_softsymb(const full_ss *fss, llr_ss *ss)
{
for (int b = 0; b < 8; ++b)
{
float v = (1.0f - fss->p[b]) / (fss->p[b] + 1e-6);
int r = logf(v) * 5; // TBD Optimal scaling vs saturation ?
if (r < -127)
r = -127;
if (r > 127)
r = 127;
ss->bits[b] = r;
}
}
} // leansdr

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,188 @@
#ifndef LEANSDR_SOFTWORD_H
#define LEANSDR_SOFTWORD_H
namespace leansdr
{
// This file implements options for SOFTBYTE in s2_frame_receiver
// (representation of bytes after deinterleaving).
// Also used as SOFTWORD in ldpc_engine
// (representation of packed SOFTBITs).
// Interface for ldpc_engine:
// SOFTBIT softword_get(const SOFTWORD &p, int b)
// SOFTBIT softwords_xor(const SOFTWORD p1[], const SOFTWORD p2[], int b)
// void softword_zero(SOFTWORD *p, int b)
// void softwords_set(SOFTWORD p[], int b)
// void softwords_flip(SOFTWORD p[], int b)
// HARD SOFTWORD / SOFTBYTE
typedef uint8_t hard_sb; // Bit 7 is transmitted first.
inline bool softword_get(const hard_sb &p, int b)
{
return p & (1 << (7 - b));
}
inline void softword_set(hard_sb *p, int b, bool v)
{
hard_sb mask = 1 << (7 - b);
*p = ((*p) & ~mask) | (v << (7 - b));
}
inline void softword_clear(hard_sb *p) { *p = 0; }
inline bool softword_weight(const bool &l)
{
return true;
}
inline void softbit_set(bool *p, bool v) { *p = v; }
inline bool softbit_harden(bool b) { return b; }
inline uint8_t softbyte_harden(const hard_sb &b) { return b; }
inline bool softbit_xor(bool x, bool y) { return x ^ y; }
inline void softbit_clear(bool *p) { *p = false; }
inline bool softwords_xor(const hard_sb p1[], const hard_sb p2[], int b)
{
return (p1[b / 8] ^ p2[b / 8]) & (1 << (7 - (b & 7)));
}
inline void softword_zero(hard_sb *p)
{
*p = 0;
}
inline void softwords_set(hard_sb p[], int b)
{
p[b / 8] |= 1 << (7 - (b & 7));
}
inline void softword_write(hard_sb &p, int b, bool v)
{
hard_sb mask = 1 << (7 - b);
p = (p & ~mask) | (hard_sb)v << (7 - b);
}
inline void softwords_write(hard_sb p[], int b, bool v)
{
softword_write(p[b / 8], b & 7, v);
}
inline void softwords_flip(hard_sb p[], int b)
{
p[b / 8] ^= 1 << (7 - (b & 7));
}
uint8_t *softbytes_harden(hard_sb p[], int nbytes, uint8_t storage[])
{
return p;
}
// LLR SOFTWORD
struct llr_sb
{
llr_t bits[8]; // bits[0] is transmitted first.
};
float prob(llr_t l)
{
return (127.0 + l) / 254;
}
llr_t llr(float p)
{
int r = -127 + 254 * p;
if (r < -127)
r = -127;
if (r > 127)
r = 127;
return r;
}
inline llr_t llr_xor(llr_t lx, llr_t ly)
{
float px = prob(lx), py = prob(ly);
return llr(px * (1 - py) + (1 - px) * py);
}
inline llr_t softword_get(const llr_sb &p, int b)
{
return p.bits[b];
}
// Special case to avoid computing b/8*8+b%7. Assumes llr_sb[] packed.
inline llr_t softwords_get(const llr_sb p[], int b)
{
return p[0].bits[b]; // Beyond bounds on purpose
}
inline void softword_set(llr_sb *p, int b, llr_t v)
{
p->bits[b] = v;
}
inline void softword_clear(llr_sb *p)
{
memset(p->bits, 0, sizeof(p->bits));
}
// inline llr_t softwords_get(const llr_sb p[], int b) {
// return softword_get(p[b/8], b&7);
// }
inline llr_t softword_weight(llr_t l) { return abs(l); }
inline void softbit_set(llr_t *p, bool v) { *p = v ? -127 : 127; }
inline bool softbit_harden(llr_t l) { return (l < 0); }
inline uint8_t softbyte_harden(const llr_sb &b)
{
// Without conditional jumps
uint8_t r = (((b.bits[0] & 128) >> 0) |
((b.bits[1] & 128) >> 1) |
((b.bits[2] & 128) >> 2) |
((b.bits[3] & 128) >> 3) |
((b.bits[4] & 128) >> 4) |
((b.bits[5] & 128) >> 5) |
((b.bits[6] & 128) >> 6) |
((b.bits[7] & 128) >> 7));
return r;
}
inline llr_t softbit_xor(llr_t x, llr_t y) { return llr_xor(x, y); }
inline void softbit_clear(llr_t *p) { *p = -127; }
inline llr_t softwords_xor(const llr_sb p1[], const llr_sb p2[], int b)
{
return llr_xor(p1[b / 8].bits[b & 7], p2[b / 8].bits[b & 7]);
}
inline void softword_zero(llr_sb *p)
{
memset(p, -127, sizeof(*p));
}
inline void softwords_set(llr_sb p[], int b)
{
p[b / 8].bits[b & 7] = 127;
}
inline void softword_write(llr_sb &p, int b, llr_t v)
{
p.bits[b] = v;
}
inline void softwords_write(llr_sb p[], int b, llr_t v)
{
softword_write(p[b / 8], b & 7, v);
}
inline void softwords_flip(llr_sb p[], int b)
{
llr_t *l = &p[b / 8].bits[b & 7];
*l = -*l;
}
uint8_t *softbytes_harden(llr_sb p[], int nbytes, uint8_t storage[])
{
for (uint8_t *q = storage; nbytes--; ++p, ++q)
*q = softbyte_harden(*p);
return storage;
}
// Generic wrappers
template <typename SOFTBIT, typename SOFTBYTE>
inline SOFTBIT softwords_get(const SOFTBYTE p[], int b)
{
return softword_get(p[b / 8], b & 7);
}
template <typename SOFTBIT, typename SOFTBYTE>
inline void softwords_set(SOFTBYTE p[], int b, SOFTBIT v)
{
softword_set(&p[b / 8], b & 7, v);
}
template <typename SOFTBIT, typename SOFTBYTE>
inline SOFTBIT softwords_weight(const SOFTBYTE p[], int b)
{
return softword_weight(softwords_get<SOFTBIT, SOFTBYTE>(p, b));
}
} // namespace leansdr
#endif // LEANSDR_SOFTWORD_H

View File

@ -1,3 +1,19 @@
// This file is part of LeanSDR Copyright (C) 2016-2018 <pabr@pabr.org>.
// See the toplevel README for more information.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
#ifndef LEANSDR_VITERBI_H
#define LEANSDR_VITERBI_H
@ -24,7 +40,7 @@ namespace leansdr
// 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>
template <typename TS, int NSTATES, typename TUS, int NUS, int NCS>
struct trellis
{
static const int NOSTATE = NSTATES + 1;
@ -33,9 +49,9 @@ struct trellis
{
struct branch
{
TS pred; // Predecessor state or NOSTATE
TUS us; // Uncoded symbol
} branches[NCS]; // Incoming branches indexed by coded symbol
TS pred; // Predecessor state or NOSTATE
TUS us; // Uncoded symbol
} branches[NCS]; // Incoming branches indexed by coded symbol
} states[NSTATES];
trellis()
@ -50,8 +66,7 @@ struct trellis
{
if (NCS & (NCS - 1))
{
fprintf(stderr, "leansdr::trellis::init_convolutional: NCS must be a power of 2\n");
return;
fprintf(stderr, "NCS must be a power of 2\n");
}
// Derive number of polynomials from NCS.
int nG = log2i(NCS);
@ -61,29 +76,27 @@ struct trellis
for (TUS us = 0; us < NUS; ++us)
{
// Run the convolutional encoder from state s with input us
uint64_t shiftreg = s; // TBD type
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
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
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;
fprintf(stderr, "Invalid convolutional code\n");
}
b->pred = s;
b->us = us;
}
}
}
void dump()
@ -102,42 +115,46 @@ struct trellis
fprintf(stderr, "\n");
}
}
};
// Interface that hides the templated internals.
template<typename TUS, typename TCS, typename TBM, typename TPM>
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;
virtual TUS update(TBM *costs, TPM *quality = NULL) = 0;
virtual TUS update(TCS s, TBM cost, TPM *quality = NULL) = 0;
};
template<typename TS, int NSTATES, typename TUS, int NUS, typename TCS, int NCS, typename TBM, typename TPM, typename TP>
struct viterbi_dec: viterbi_dec_interface<TUS, TCS, TBM, TPM>
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
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
statebank *states, *newstates; // Alternate between banks
viterbi_dec(trellis<TS, NSTATES, TUS, NUS, NCS> *_trellis) :
trell(_trellis)
viterbi_dec(trellis<TS, NSTATES, TUS, NUS, NCS> *_trellis) : trell(_trellis)
{
states = &statebanks[0];
newstates = &statebanks[1];
for (TS s = 0; s < NSTATES; ++s)
(*states)[s].cost = 0;
// Determine max value that can fit in TPM
max_tpm = (TPM) 0 - 1;
max_tpm = (TPM)0 - 1;
if (max_tpm < 0)
{
// TPM is signed
@ -160,12 +177,13 @@ struct viterbi_dec: viterbi_dec_interface<TUS, TCS, TBM, TPM>
// 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];
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
{ // <= guarantees one match
best_m = m;
best_b = b;
}
@ -193,10 +211,10 @@ struct viterbi_dec: viterbi_dec_interface<TUS, TCS, TBM, TPM>
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)
@ -221,12 +239,13 @@ struct viterbi_dec: viterbi_dec_interface<TUS, TCS, TBM, 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]];
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
{ // <= guarantees one match
best_m = m;
best_b = b;
}
@ -238,7 +257,8 @@ struct viterbi_dec: viterbi_dec_interface<TUS, TCS, TBM, TPM>
// 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];
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;
@ -272,10 +292,10 @@ struct viterbi_dec: viterbi_dec_interface<TUS, TCS, TBM, TPM>
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)
@ -301,7 +321,7 @@ struct viterbi_dec: viterbi_dec_interface<TUS, TCS, TBM, TPM>
fprintf(stderr, "\n");
}
private:
private:
TPM max_tpm;
};
@ -310,24 +330,15 @@ private:
// 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>
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
} // namespace leansdr
#endif // LEANSDR_VITERBI_H
#endif // LEANSDR_VITERBI_H