// // ebnaut.c // // Copyright (c) 2015 Paul Nicholson // All rights reserved. // // Redistribution and use in source and binary forms, with or without // modification, are permitted provided that the following conditions // are met: // 1. Redistributions of source code must retain the above copyright // notice, this list of conditions and the following disclaimer. // 2. Redistributions in binary form must reproduce the above copyright // notice, this list of conditions and the following disclaimer in the // documentation and/or other materials provided with the distribution. // // THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR // IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES // OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. // IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, // INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT // NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF // THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. // // http://abelian.org/ebnaut/ #define VERSION "1.1" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #define FALSE 0 #define TRUE (!FALSE) // // Encoding or decoding. // #define ENCODE 1 #define DECODE 2 static int mode = 0; // // Output options. // #define OP_CARR 2 // Carrier output #define OP_OFF 4 // Carrier off, noise only static uint32_t opmode = 0; static int TFLAG = FALSE; // -t option: output encoded bits in text format static int SFLAG = 0; // -s option: do source coding only // // General options. // static double sample_rate = 0; // I/O sample rate, from -r static double sym_period = 0; // Signal bit period, -S option static int Nchars = 0; // -N option, message length, number of characters static int anorm = 0; // Set by -a option static char *mtext = NULL; // Argument of -M option static double start_offset = 0; // From -T, seconds into the input file of // first symbol static double freq_offset = 0; // Frequency offset, -F option, Hertz static int ncpu = 1; // -c option: number of CPUs to use #define MAXCPU 64 // // Noise settings. // static int NFLAG = 0; // Non-zero if noise is to be added static double request_ebn0; // Carrier to noise ratio in info bit bandwidth // // Phase options. // static int PSTEP = 0; // -PS option: do phase search static int PRAND = FALSE; // -PR: send with random phase static int PU = FALSE; // -PU option: uniform phase only static int PSTEP_INIT = 0; static int PFLAG = FALSE; // TRUE if -P supplies a specific phase to use static double refP = 0.0; // Phase specified by -P option // // Diagnostic flags, various special functions. // static int FFLAG2 = FALSE; // TRUE with -f2 static int FFLAG3 = FALSE; // TRUE with -f3 static int FFLAG4 = FALSE; // TRUE with -f4 static int FFLAG6 = FALSE; // TRUE with -f6 static int FFLAG7 = FALSE; // TRUE with -f7 static int FFLAG8 = FALSE; // TRUE with -f8 static int FFLAG9 = FALSE; // TRUE with -f9 static int FFLAG10 = FALSE; // TRUE with -f10 static int FFLAG13 = FALSE; // TRUE with -f13 static int FFLAG14 = FALSE; // TRUE with -f14 static int FFLAG15 = FALSE; // TRUE with -f15 static int FFLAG16 = FALSE; // TRUE with -f16 static int FFLAG17 = FALSE; // TRUE with -f17 static int FFLAG18 = FALSE; // TRUE with -f18 static int FFLAG19 = FALSE; // TRUE with -f19 static int FFLAG20 = FALSE; // TRUE with -f20 static int XFLAG = FALSE; // -X option: sometimes used to turn on experimental // code // // // static int QD = 0; // Number of output bits per input bit static int LB = 0; // Viterbi block size, Nchars * 6 + CRCLEN static int Nsymbols = 0; // Number of symbols per block static time_t Tstart = 0; // Decode start time, for elapsed time calculation static int IW = 0; // Interleave stride static double *sigbits = NULL; // Signal bit NRZ amplitudes static complex double *symamps = NULL; // Complex amplitude of each symbol // after frequency offset applied static complex double *rotamps = NULL; // Complex amplitude of each symbol // after frequency offset, phase // pattern, and de-interleave have // been applied static complex double *refph = NULL; // Phase adjustment for each symbol static double *normf = NULL; // Normalisation factor for each symbol #define MAX_CHARS 2000 // Maximum number of characters static char source_text[MAX_CHARS+1]; // Message text #define MAX_BITS 10000 static uint8_t source_bits[MAX_BITS+1]; // Source-encoded message bits #define MAX_RANK 100000000 /////////////////////////////////////////////////////////////////////////////// // Predefined Polynomial Sets // /////////////////////////////////////////////////////////////////////////////// static struct POLYLIST { char *alias; char *poly; } polylist[] = { { "1K0", "1" }, // Uncoded messages { "2K3A", "7,5" }, { "2K4A", "0xd,0xf" }, { "2K5A", "0x13,0x1d" }, { "2K6A", "0x2b,0x3d" }, { "2K7A", "0x5b,0x79" }, { "2K8A", "0xa7,0xf9" }, { "2K9A", "0x131,0x1eb" }, { "2K10A", "0x277,0x365" }, { "2K11A", "0x4dd,0x7b1" }, { "2K12A", "04335,05723" }, { "2K13A", "010533,017661" }, { "2K13B", "011651,014677" }, { "2K14A", "021675,027123" }, { "2K14B", "026651,036477" }, { "2K15A", "047244,065231" }, { "2K15B", "046253,077361" }, { "2K16A", "0514576,0715022" }, { "2K16B", "0114727,0176121" }, { "2K17A", "0207225,0330747" }, { "2K17B", "0213631,0353323" }, { "2K18A", "0507517,0654315" }, { "2K21A", "05016515,06770547" }, { "2K23A", "051202215,066575563" }, { "3K3A", "5,7,7" }, { "3K4A", "013,015,017" }, { "3K5A", "025,033,037" }, { "3K6A", "047,053,075" }, { "3K7A", "0x5b,0x65,0x7d" }, { "3K8A", "0225,0331,0367" }, { "3K9A", "0557,0663,0711" }, { "3K10A", "01117,01365,01633" }, { "3K11A", "02353,02671,03175" }, { "3K12A", "04767,05723,06265" }, { "3K13A", "010533,010675,017661" }, { "3K14A", "021645,035661,037133" }, { "4K13A", "011145,012477,015537,016727" }, { "4K14A", "021113,023175,035527,035537" }, { "4K15A", "046321,051271,063667,070535" }, { "4K15B", "050575,051447,056533,066371" }, { "4K16A", "0123175,0135233,0156627,0176151" }, { "4K17A", "0235433,0247631,0264335,0311727" }, { "4K19A", "01122645,01373343,01531175,01634157" }, { "4K21A", "04542365,05342433,06347677,07232671" }, { "4K23A", "022346335,024275433,035520777,036271251" }, { "4K25A", "0106042635,0125445117,0152646773,0167561761" }, { "8K17A", "0222537,0226345,0255077,0267667,0306347,0315117,0326721," "0372751" }, { "8K19A", "01154451,01214561,01221517,01226537,01276775,01523533," "01543563,01632737" }, { "8K21A", "04470745,04714453,05175161,05366615,06375351,06752427," "06775363,07310533" }, { "8K23A", "023372665,024534765,025561775,027325743,031716223,032210543," "035451321,036744713" }, { "8K25A", "0105341521,0111747353,0114371121,0122355637,0134552773," "0150627331,0164507577,0172554315" }, { "16K19A", "01047113,01072363,01131677,01153145,01214573,01306665," "01334255,01346335,01456535,01555113,01676751,01723045," "01737471,01743251,01761527,01765237" }, { "16K21A", "04232547,04505631,05114365,05234555,05275577,05333523," "05634353,05663763,05735521,06110751,06273753,06523067," "06723455,07454677,07647511,07710617" }, { "16K23A", "021570463,022334751,023471145,023564713,024225255," "025336745,026522675,030653713,030664223,031367237," "032361377,033355313,036377267,037212147,037352421," "037553071" }, { "16K25A", "0104722251,0107672671,0122541663,0122545535,0123470447," "0125364663,0131553555,0131615665,0145114761,0145522105," "0146174323,0157705733,0162370765,0166576517,0167572771," "0176274217" }, { NULL, NULL } // Sentinel }; /////////////////////////////////////////////////////////////////////////////// // Utilities // /////////////////////////////////////////////////////////////////////////////// // // Terminate with a fatal error message. // static void bailout( char *format, ...) { va_list ap; char temp[ 200]; va_start( ap, format); vsprintf( temp, format, ap); va_end( ap); fprintf( stderr, "bailout: %s\n", temp); exit( 1); } #if 0 // Only used for testing static double rusage_time( void) { struct rusage r; if( getrusage( RUSAGE_SELF, &r)) bailout( "rusage, %s", strerror( errno)); return r.ru_utime.tv_sec + 1e-6 * r.ru_utime.tv_usec + r.ru_stime.tv_sec + 1e-6 * r.ru_stime.tv_usec; } #endif // // Report something to standard error stream, subject to log level. // static int loglevel = 0; // Incremented by -v options static void report( int n, char *format, ...) { va_list ap; char temp[ 200]; if( n > loglevel) return; va_start( ap, format); vsprintf( temp, format, ap); va_end( ap); fprintf( stderr, "%s\n", temp); fflush( stderr); } static void * __attribute__ ((malloc)) safe_malloc( size_t size) { void *p = malloc( size); if( !p) bailout( "not enough memory [%zd]", size); return p; } /////////////////////////////////////////////////////////////////////////////// // Source Coding // /////////////////////////////////////////////////////////////////////////////// // // * (dec 42) to ; (dec 59) coded as 0 to 17 // ! (dec 33) coded to 18 // = (dec 61) coded to 19 // & (dec 38) coded to 20 // ? (dec 63) to Z (dec 90) coded as 21 to 48 // '\n' coded as 49 but decoded to a space // space and tab both coded to 50 // underscore codes as 51 // static void encode_char( char c, int *nt) { int u = -1; c = toupper( c); if( c >= 42 && c <= 59) u = c - 42; else if( c >= 63 && c <= 90) u = c - 42; else switch( c) { case '!': u = 18; break; case '=': u = 19; break; case '&': u = 20; break; case '\n': u = 49; break; case ' ': case '\t': u = 50; break; case '_': u = 51; break; } if( u >= 0) // A valid character to encode? { for( int i = 0; i < 6; i++) source_bits[*nt + i] = u & (1 << i) ? 1 : 0; *nt += 6; } } static char decode_char( uint8_t *bits) { int u = 0; for( int i = 0; i < 6; i++) if( bits[i]) u |= 1 << i; if( u >= 0 && u <= 17) return u + 42; if( u >= 21 && u <= 48) return u + 42; switch( u) { case 18: return '!'; case 19: return '='; case 20: return '&'; case 49: return ' '; case 50: return ' '; case 51: return '_'; } return '#'; // Indicates an invalid character } /////////////////////////////////////////////////////////////////////////////// // CRC // /////////////////////////////////////////////////////////////////////////////// static int CRCLEN = 16; // http://users.ece.cmu.edu/~koopman/crc/ // http://users.ece.cmu.edu/~koopman/crc/hw_data.html #define CRCBASE 0 // The CRC in the first entry of crclist[] #define CRCMAX 32 // Max CRC in crclist[] static uint64_t crclist[] = { 0x1, // 0 0x3, // 1 0x7, // 2 0xb, // 3 CRC-3K 0x13, // 4 CCITT-4 0x2b, // 5 CRC-5-ITU 0x47, // 6 CRC-6-CDMA2000-B 0xe5, // 7 CRC-7-MVB 0x137, // 8 CRC-8F/6 0x279, // 9 CRC-9F/6.2 0x537, // 10 CRC-10F/7 0x9eb, // 11 CRC-11F/8 0x149f, // 12 CRC-12F/8.2 0x216f, // 13 CRC-10F/8.2 0x46e3, // 14 CRC-14F/8.2 0xb48f, // 15 CRC-15K/9 0x11021, // 16 CCITT-16 0x2ed4f, // 17 CRC-17K/10 0x4d47b, // 18 CRC-18K/11 0xa3af3, // 19 CRC-19K/12 0x13ab0f, // 20 CRC-20K/12 0x2caea3, // 21 CRC-21K/12 0x5b9d23, // 22 CRC-22K/13 0x96f3a3, // 23 CRC-23K/14 0x16e7d23, // 24 CRC-24K/14 0x252399f, // 25 CRC-25K/14 0x4a3ecd7, // 26 CRC-26K/16 0xcb7aa27, // 27 0x19e2372b, // 28 0x2bc2cb4d, // 29 0x453be359, // 30 0x8dcad4f9, // 31 0x1741b8cd7 // 32 CRC-32K/6.1 }; // // 'data' and 'crc' are arrays containing one bit per byte. // static void make_crc( uint8_t *data, int len, uint8_t *crc) { if( !CRCLEN) return; uint8_t poly[33]; memset( poly, 0, 33); uint64_t c = crclist[CRCLEN - CRCBASE]; for( int i = 0; i <= CRCLEN; i++) poly[i] = c & (1 << (CRCLEN - i)) ? 1 : 0; uint8_t *buf = safe_malloc( len + CRCLEN); memcpy( buf, data, len); memset( buf + len, 0, CRCLEN); uint8_t *p = buf; for( int i = 0; i < len; i++, p++) if( *p) for( int j = 0; j <= CRCLEN; j++) p[j] ^= poly[j]; memcpy( crc, buf+len, CRCLEN); free( buf); } // // Check that a candidate decode has a valid CRC. // static int check_crc( void) { if( !CRCLEN) return TRUE; if( FFLAG19) return TRUE; uint8_t crc[CRCLEN]; make_crc( source_bits, LB - CRCLEN, crc); for( int i = 0; i < CRCLEN; i++) if( source_bits[LB - CRCLEN + i] != crc[i]) return FALSE; return TRUE; } // // Check that each character of a candidate decode has a valid source // encoding. // static int valid_message( void) { if( FFLAG19) return TRUE; for( int i = 0; i < LB - CRCLEN; i += 6) if( decode_char( source_bits + i) == '#') return FALSE; return TRUE; } /////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////// // // Load the source text from stdin. // static void load_source( void) { int n = 0, nt = 0; if( mtext) // Message supplied by -M option? { strcpy( source_text, mtext); n = strlen( source_text); } else n = fread( source_text, 1, MAX_CHARS, stdin); if( n <= 0) bailout( "no source"); source_text[n] = 0; for( int i = 0; i < n; i++) if( i != n-1 || source_text[i] != '\n') encode_char( source_text[i], &nt); if( nt / 6 != Nchars) bailout( "incorrect source length %d/%d", nt/6, Nchars); make_crc( source_bits, nt, source_bits + nt); } // // Generate a random source text. // static void gen_random_text( void) { for( int i = 0; i < Nchars; i++) { int c = 51 * (int64_t) rand() / RAND_MAX; for( int j = 0; j < 6; j++) source_bits[i*6 + j] = c & (1 << j) ? 1 : 0; } make_crc( source_bits, LB - CRCLEN, source_bits + LB - CRCLEN); } // // Generate the (N+1)th source text. // static void gen_seq_text( uint64_t N) { for( int i = 0; i < Nchars; i++) { int c = N % 52; N /= 52; for( int j = 0; j < 6; j++) source_bits[i*6 + j] = c & (1 << j) ? 1 : 0; } make_crc( source_bits, LB - CRCLEN, source_bits + LB - CRCLEN); } // // Generate a random source. // static void gen_random_source( void) { uint32_t r = rand(); int k = 0; for( int i = 0; i < LB; i++, k++) { if( k == 30) { k = 0; r = rand(); } source_bits[i] = r & (1 << k) ? 1 : 0; } } // // Generate random binary source. // // // Generate the (N+1)th binary source text // static void gen_seq_source( uint64_t N) { for( int i = 0; i < LB; i++) source_bits[i] = N & ((uint64_t)1 << i) ? 1 : 0; } /////////////////////////////////////////////////////////////////////////////// // Interleaving // /////////////////////////////////////////////////////////////////////////////// // // For offset 'in' in the transmission order, return the offset into the // non-interleaved data. // static int interleave( int in) { int N = Nsymbols; if( in >= N) bailout( "interleave input out of range"); int H = N / IW; // Number of rows int L = N % IW; if( L) H++; // Extra incomplete row of length L int x = in % IW; int y = in / IW; if( !L || x < L) return x * H + y; return L * H + (x - L) * (H-1) + y; } /////////////////////////////////////////////////////////////////////////////// // Modulation Layer // /////////////////////////////////////////////////////////////////////////////// // // Noise generator. Produce a pair of Gaussian noise samples by the // polar method. // static double noise_sigma = 1.0; static void noise_gen( double *n1, double *n2) { double v1, v2, s; do { double u1 = drand48(); double u2 = drand48(); v1 = 2 * u1 - 1; v2 = 2 * u2 - 1; s = v1*v1 + v2*v2; } while( s >= 1); double a = sqrt( -2 * log(s) / s); *n1 = a * v1 * noise_sigma; *n2 = a * v2 * noise_sigma; } static void send_symbol( int n) { double bit = sigbits[interleave( n)]; if( TFLAG) // -t option: output the message bits in text format { printf( "%d", bit > 0 ? 1 : 0); return; } if( opmode & OP_CARR) bit = 1; // Base band signal of unit amplitude and unit power double vc = bit * cos( refP); double vs = bit * sin( refP); if( opmode & OP_OFF) vc = vs = 0; for( int i = 0; i < (int)(sample_rate * sym_period); i++) { double noise1 = 0, noise2 = 0; if( NFLAG) noise_gen( &noise1, &noise2); printf( "%.7e %.7e\n", vc + noise1, vs + noise2); } } /////////////////////////////////////////////////////////////////////////////// // Demodulation Layer // /////////////////////////////////////////////////////////////////////////////// static double input_sum_squared = 0; // // Read a pair of I/Q samples. Input is ASCII two or three column. // If three column, the first column is assumed to be a timestamp and ignored. // static int recv_analogue( complex double *v) { char temp[100]; if( !fgets( temp, 99, stdin)) return FALSE; double v1, v2, v3; int n = sscanf( temp, "%lf %lf %lf\n", &v1, &v2, &v3); if( n < 2 || n > 3) bailout( "input fields %d [%s]", n, temp); if( n == 2) *v = v1 + I* v2; else *v = v2 + I * v3; return TRUE; } static int offset_samples = 0; static int nsbase = 0; // The current input sample, offset into input file data static int nspad1 = 0; // Number of zero samples added for padding at start static int nspad2 = 0; // Number of zero samples added for padding at end static int nsskip = 0; // Number of input samples skipped at beginning of file // // Receive symbol number n, 0 <= n < Nsymbols. // static void receive_symbol( int n) { // Rotation to implement the frequency offset double ph_inc = freq_offset/sample_rate * 2 * M_PI; // Radians per sample complex double a = 0; int ns2 = sample_rate * sym_period * (n + 1); for( ; nsbase < ns2; nsbase++) { complex double v; if( offset_samples < 0) { offset_samples++; nspad1++; } else if( recv_analogue( &v)) { double ph = ph_inc * nsbase; a += v * (cos(ph) + I * sin(ph)); input_sum_squared += cabs(v) * cabs(v); } else nspad2++; } symamps[n] = a; } // // Construct a pattern of reference phases. Returns FALSE if there are no // more phase patterns to try. // static char p_phase_txt[100]; static int prepare_phase( int n) { #define PMOD (360 / PSTEP) static int klist[12] = { 0, 1, -1, 2, -2, 3, -3, 4, -4, 5, -5, 6 }; int ns = n; int k1 = n % PMOD; n /= PMOD; int kp = k1 & 1; k1 >>= 1; // Polarity 0 or 180 int ph_a = klist[k1] + (kp ? PMOD/2 : 0); if( PU && n) return FALSE; int k2 = n % 3; n /= 3; int ph_b[4]; for( int i = 0; i < 4; i++) ph_b[i] = 0; if( k2 == 1) ph_b[0] = ph_b[1] = 1; if( k2 == 2) ph_b[2] = ph_b[3] = 1; int k3 = n % 3; n /= 3; int ph_c[4]; for( int i = 0; i < 4; i++) ph_c[i] = 0; if( k3 == 1) ph_c[0] = ph_c[2] = 1; if( k3 == 2) ph_c[1] = ph_c[3] = 1; if( ns && !kp && !k1 && !k2 && !k3) return FALSE; int ph[4]; for( int i = 0; i < 4; i++) { ph[i] = ph_a + ph_b[i] + ph_c[i]; ph[i] *= PSTEP; ph[i] %= 360; if( ph[i] > 180) ph[i] -= 360; } sprintf( p_phase_txt, "%d %d,%d,%d,%d", ns, ph[0], ph[1], ph[2], ph[3]); for( int i = 0; i < Nsymbols; i++) { int j = (4 * i)/(double) Nsymbols; double pht = ph[j] * M_PI/180; complex double a = cos(pht) + I * sin(pht); sigbits[interleave( i)] = creal( symamps[i] * a); rotamps[interleave( i)] = symamps[i] * a; } report( 1, "phase %3d %4d %4d %4d %4d", ns, ph[0], ph[1], ph[2], ph[3]); return TRUE; } // // Initialise the demodulator and receive Nsymbols from the input file. // static void receive_block( void) { complex double dummy; // // Initialise and read in all the samples. // offset_samples = round( start_offset * sample_rate); nsbase = nsskip = nspad1 = nspad2 = 0; if( offset_samples > 0) { for( int i = 0; i < offset_samples; i++) recv_analogue( &dummy); nsskip = offset_samples; } input_sum_squared = 0; for( int i = 0; i < Nsymbols; i++) receive_symbol( i); if( nsskip) report( 0, "skipped %.6f seconds to start", nsskip/sample_rate); if( nspad1) report( 0, "padded %.3f seconds at start", nspad1/sample_rate); if( nspad2) report( 0, "padded %.3f seconds at end", nspad2/sample_rate); // // If an explicit reference phase has been given with a -P option, apply // that now by rotating all the symbol amplitudes. // if( PFLAG) { double ph = refP; complex double a = cos(ph) - I * sin(ph); for( int i = 0; i < Nsymbols; i++) symamps[i] /= a; } // // Do amplitude normalisation: -a option // if( anorm) { int inorm = anorm / 2; for( int i = 0; i < Nsymbols; i++) { double sum = 0; int n = 0; for( int j = i - inorm; j <= i + inorm; j++) if( j >= 0 && j < Nsymbols) { sum += cabs( symamps[j]); n++; } normf[i] = sum / n; } for( int i = 0; i < Nsymbols; i++) if( normf[i]) symamps[i] /= normf[i]; } // // Output a map of the input symbols. // FILE *fd = fopen( "rawsyms", "w"); for( int i = 0; i < Nsymbols; i++) { complex double a = symamps[i]; fprintf( fd, "%d %.6e %.6e %.2f\n", i, creal(a), cimag(a), carg(a) * 180/M_PI); } fclose( fd); // // Try to determine a reference phase by summing the squared complex // symbol amplitudes. // complex double sumsq = 0; for( int i = 0; i < Nsymbols; i++) sumsq += symamps[i] * symamps[i]; complex double ref = csqrt(sumsq); // // Only report the initial reference phase if we're going to use it. // With -f4 we report the phase and we're finished. // if( !PFLAG && !FFLAG15) report( 1, "initial reference phase %.1f amplitude %.3e", -180 * carg(ref)/M_PI, sqrt( cabs(sumsq))); if( FFLAG4) { printf( "%.2f %.3e\n", -180 * carg(ref)/M_PI, sqrt( cabs(sumsq))); exit( 0); } // // Rotate input signal to zero the reference phase. Only if a reference // phase has not been specified with -P. Only if no -f15 option. // if( !PFLAG && !FFLAG15 && cabs( ref)) { ref /= cabs( ref); for( int i = 0; i < Nsymbols; i++) symamps[i] /= ref; } for( int i = 0; i < Nsymbols; i++) { rotamps[interleave( i)] = symamps[i]; sigbits[interleave( i)] = creal( symamps[i]); } } // // Receive in text mode when -t option is given. Each symbol is an ASCII // '0' or '1' character. Line terminators are ignored, any other character // is an error. // static void receive_block_text( void) { for( int i = 0; i < Nsymbols; i++) { int c = getchar(); if( c == EOF) bailout( "short input, read %d, expecting %d", i, Nsymbols); if( c == '1') sigbits[interleave( i)] = +1.0; else if( c == '0') sigbits[interleave( i)] = -1.0; else if( c == '\n' || c == '\r') continue; else bailout( "invalid input char [%d]", c); } } /////////////////////////////////////////////////////////////////////////////// // Convolutional Encoder // /////////////////////////////////////////////////////////////////////////////// static int QK = 0; // Encoder length: 1 + number of memory bits static int QS; // Number of encoder states static int QP[16]; // A set of polynomials, as set by -p option static uint32_t encstate; // Current encoder state static uint32_t outword = 0; static uint32_t weight = 0; static inline void reset_encoder( uint32_t val) { encstate = val; weight = 0; } static void encode( int b) { outword = 0; encstate |= b << (QK-1); for( int i = 0; i < QD; i++) outword = (outword << 1) | __builtin_parity( QP[i] & encstate); encstate >>= 1; // Discarding the oldest memory bit } // Like encode(), but just counts the output code weight static inline void encode_count( int b) { uint32_t u = encstate | b; for( int i = 0; i < QD; i++) if( __builtin_popcount( QP[i] & u) & 1) weight++; encstate = u << 1; } /////////////////////////////////////////////////////////////////////////////// // Viterbi Decoder // /////////////////////////////////////////////////////////////////////////////// typedef float VTYPE; // Data type for Viterbi typedef VTYPE TUPLE[16]; typedef uint16_t ETYPE; // // Table to list the two previous states that branch into each state, and their // associated info bit values. // static struct BMAP { uint32_t prev[2]; int bit[2]; } *BMAP; // // Table to list the output pattern and next state for each encoder state // transition. // struct STATE { ETYPE emit[2]; int next[2]; } *STATE; static int T = 0; // Number of code tuples in the input buffer static int trellis_length; // Length of trellis struct TRELLIS { VTYPE M1, // Accumulated branch metric of the most likely path to this state M2; // Accumulated branch metric of the other path to this state int v1, v2; // Previous states corresponding to paths with M1 and M2 } **trellis; // // Return the squared Euclidean distance between the received vector 'v', // and the expected vector 'e'. We drop the two terms which are constant // within each column of the trellis. What remains is just an inner product. // static inline VTYPE distance( ETYPE e, TUPLE v) { VTYPE ret = 0; for( int i = 0; i < QD; i++) ret -= e & (1 << i) ? v[i] : -v[i]; return ret; } static TUPLE *decoder_input = NULL; static int vt = 0; // The current column in the trellis during the forward pass static void forward_pass( int cpu) { // // Determine which trellis rows (states) this thread is going to process. // int j1 = cpu * QS/ncpu; int j2 = (cpu+1) * QS/ncpu; if( cpu == ncpu - 1) j2 = QS; // // Process the states in column 'vt' of the trellis. // for( int j = j1; j < j2; j++) // For each current state { trellis[j][vt].v1 = -1; trellis[j][vt].v2 = -1; // // Find the transition to state j which has minimum accumulated branch // metric. The previous state gets stored in v1 and its ABM in M1. // // The other transition (if there is one) gets put into v2 with ABM M2. // for( int n = 0; n < 2; n++) // For the two prev states that branch into j { int k = BMAP[j].prev[n]; if( trellis[k][vt-1].v1 < 0) // Previous state k never reached? continue; int bit = BMAP[j].bit[n]; VTYPE abm = trellis[k][vt-1].M1 + distance( STATE[k].emit[bit], decoder_input[vt-1]); if( trellis[j][vt].v1 < 0) { trellis[j][vt].M1 = abm; trellis[j][vt].v1 = k; } else if( abm >= trellis[j][vt].M1) { trellis[j][vt].M2 = abm; trellis[j][vt].v2 = k; } else { trellis[j][vt].M2 = trellis[j][vt].M1; trellis[j][vt].M1 = abm; trellis[j][vt].v2 = trellis[j][vt].v1; trellis[j][vt].v1 = k; } } } } // Semaphores for synchronising forward pass threads static sem_t vt_sem[MAXCPU]; // To synchronise launch of background work static sem_t vt_cnt; // To wait for background work complete static pthread_t pids[MAXCPU]; static void *decode3_start( void *arg) { int cpu = (int) (long) arg; while( 1) { sem_wait( &vt_sem[cpu]); if( vt > T) break; forward_pass( cpu); sem_post( &vt_cnt); } return 0; } static void fwd_viterbi( void) { vt = 1; sem_init( &vt_cnt, 0, 0); for( int i = 1; i < ncpu; i++) sem_init( &vt_sem[i], 0, 0); for( int i = 1; i < ncpu; i++) pthread_create( pids + i, NULL, decode3_start, (void *)(long)i); for( vt = 1; vt <= T; vt++) // For each trellis column { for( int i = 1; i < ncpu; i++) sem_post( &vt_sem[i]); forward_pass( 0); for( int i = 1; i < ncpu; i++) sem_wait( &vt_cnt); } for( int i = 1; i < ncpu; i++) sem_post( &vt_sem[i]); for( int i = 1; i < ncpu; i++) pthread_join( pids[i], NULL); } //////////////////////////////////////////////////////////////////////////////// // List Decoder Stack // //////////////////////////////////////////////////////////////////////////////// // // The stack for the backward passes. // // This is just a simple binary tree. Maybe it should be a more balanced // structure such as a red-black tree. // // Smaller branch metrics take the left side of the tree, larger metrics to // the right. // static int max_rank = 20000; // Set by -L option struct VSTACK { VTYPE m; // Path metric int q; // Backward pass number at insert int t; // Trellis column VTYPE u; struct VSTACK *left, *right; }; static struct VSTACK *vpool; static int nvstack = 0; static void vstack_init( void) { vpool = safe_malloc( max_rank * sizeof( struct VSTACK)); } static struct VSTACK *vroot = NULL, // The root of the tree *vfree; // Linked list of free nodes static int maxheight = 0; // // Prepare the pool of VSTACK structures for use. String them all into a // free list using the left pointer. // static void vstack_reset( void) { memset( vpool, 0, max_rank * sizeof( struct VSTACK)); for( int i = 0; i < max_rank - 1; i++) vpool[i].left = vpool + i + 1; vroot = NULL; vfree = vpool; nvstack = 0; } // // Return a node to the head of the free list. // static inline void vstack_free( struct VSTACK *vp) { vp->left = vfree; vfree = vp; } // // Remove and discard the node with largest branch metric from the stack. // static void vstack_pop_end( void) { struct VSTACK **pp = &vroot; while( (*pp)->right) pp = &(*pp)->right; struct VSTACK *vp = *pp; *pp = vp->left; nvstack--; vstack_free( vp); } // // Allocate a free VSTACK node from the head of the free list. // static struct VSTACK *vstack_alloc( void) { if( !vfree) vstack_pop_end(); struct VSTACK *vp = vfree; vfree = vp->left; return vp; } // // DSW rebalancing. // static void dsw_compress( struct VSTACK *s, int count) { for( int i = 0; i < count; i++) { struct VSTACK *c = s->right; s->right = c->right; s = s->right; c->right = s->left; s->left = c; } } static void dsw_rebalance( void) { struct VSTACK psroot = { 0, 0, 0, 0, NULL, vroot }; struct VSTACK *tail = &psroot; struct VSTACK *rest = tail->right; while( rest) if( !rest->left) { tail = rest; rest = rest->right; } else { struct VSTACK *temp = rest->left; rest->left = temp->right; temp->right = rest; rest = temp; tail->right = temp; } int size = nvstack; int leaves = size + 1 - (1 << (int) floor( log10(size + 1)/log10(2.0))); dsw_compress( &psroot, leaves); size = size - leaves; while( size > 1) { dsw_compress( &psroot, size/2); size /= 2; } vroot = psroot.right; } // // Insert a VSTACK node into the tree at the appropriate place for its // branch metric vp->m. // static void vstack_insert( struct VSTACK *vp) { struct VSTACK **pp = &vroot; int h = 0; // For measuring the tree height while( *pp) { pp = vp->m > (*pp)->m ? &(*pp)->right : &(*pp)->left; h++; } *pp = vp; vp->left = vp->right = NULL; nvstack++; if( h > maxheight) // Track the maximum height { maxheight = h; if( !FFLAG20 && // Rebalancing not disabled by -f20 ? maxheight > log10(nvstack)/log10(2.0) * 15.0) // Worth rebalancing? { dsw_rebalance(); report( 2, "rebalance %d", nvstack); maxheight = 0; } } } // // Pop the node with smallest branch metric from the stack and return it. // The node remains allocated and the caller must free it after use. // static struct VSTACK *vstack_pop( void) { struct VSTACK **pp = &vroot; while( (*pp)->left) pp = &(*pp)->left; struct VSTACK *vp = *pp; *pp = vp->right; nvstack--; return vp; } /////////////////////////////////////////////////////////////////////////////// // List Algorithm // /////////////////////////////////////////////////////////////////////////////// // Notes and steps refer to: // // "Fast Tree-Trellis List Viterbi Decoding" // Martin Roeder, Raouf Hamzaoui // // https://www.uni-konstanz.de/mmsp/pubsys/publishedFiles/RoHa05.pdf // Coded as per the above, but with rank (variable k) starting at zero // instead of one. static int **pk = NULL; static VTYPE backward_pass( int rank) { int t, i; VTYPE m; if( rank == 0) { // Step 1 t = T; i = 0; m = 0; pk[rank][t] = 0; } else { // Steps 2 and 3 struct VSTACK *vp = vstack_pop(); t = vp->t; i = pk[vp->q][t]; m = vp->u; for( int y = t; y <= T; y++) pk[rank][y] = pk[vp->q][y]; vstack_free( vp); // Step 4 int j = trellis[i][t].v2; // Step 5 m = m + trellis[i][t].M2 - trellis[j][t-1].M1; t--; i = pk[rank][t] = j; } do { // Step6 if( trellis[i][t].v2 >= 0) { struct VSTACK *vp = vstack_alloc(); vp->m = trellis[i][t].M2 + m; vp->q = rank; vp->t = t; vp->u = m; vstack_insert( vp); while( nvstack > max_rank - rank) vstack_pop_end(); } // Step 7 int j = trellis[i][t].v1; // Step 8 m = m + trellis[i][t].M1 - trellis[j][t-1].M1; t--; i = pk[rank][t] = j; } while( t > 0); return m; } // // Convert a discovered path (list of states) into the source bitstream. // Skip the final QK steps of the path - these are just the flushing bits. // static void backtrace( int rank) { if( XFLAG) // Which is quicker? { int *p = pk[rank]; int j = p[0]; for( int i = 0; i <= T - QK; i++) { int n = p[i+1]; source_bits[i] = STATE[j].next[0] == n ? 0 : 1; j = n; } } else for( int i = 0; i <= T - QK; i++) source_bits[i] = STATE[pk[rank][i]].next[0] == pk[rank][i+1] ? 0 : 1; } // // Once-only initialisation of the Viterbi decoder. // void init_viterbi( void) { // Trellis contains one entry for the initial state plus LB info bits // plus QK-1 flushing bits trellis_length = LB + QK; report( 2, "decoder length: %d, states %d, QK %d QD %d", trellis_length, QS, QK, QD); if( !FFLAG3) // Trellis not required for catastrophic code check { trellis = safe_malloc( sizeof( struct TRELLIS *) * QS); for( int i = 0; i < QS; i++) trellis[i] = safe_malloc( sizeof( struct TRELLIS) * trellis_length); decoder_input = safe_malloc( sizeof( TUPLE) * trellis_length); } // STATE[prev_state].emit[bit] gives the code tuple produced // by each encoder state transition. // STATE[from_state].next[bit] = to_state indicates which state transitions // occur with each info bit STATE = safe_malloc( sizeof( struct STATE) * QS); memset( STATE, 0, sizeof( struct STATE) * QS); // BMAP[to_state].prev[0..1] = from_state: Branch map lists the two // previous states that branch into to_state; BMAP = (struct BMAP *) safe_malloc( sizeof( struct BMAP) * QS); memset( BMAP, 0, sizeof( struct BMAP) * QS); uint8_t *tmp = safe_malloc( QS); memset( tmp, 0, QS); // Build the above tables simply by running the encoder through all of its // states, exercising each state with input bits 0 and 1 for( int i = 0; i < QS; i++) // For each encoder state { reset_encoder( i); encode( 0); STATE[i].emit[0] = outword; STATE[i].next[0] = encstate; BMAP[encstate].bit[tmp[encstate]] = 0; BMAP[encstate].prev[tmp[encstate]++] = i; reset_encoder( i); encode( 1); STATE[i].emit[1] = outword; STATE[i].next[1] = encstate; BMAP[encstate].bit[tmp[encstate]] = 1; BMAP[encstate].prev[tmp[encstate]++] = i; } for( int i = 0; i < QS; i++) // For each encoder state if( tmp[i] != 2) bailout( "prev branch count not 2"); free( tmp); if( !FFLAG3) { pk = safe_malloc( max_rank * sizeof( int *)); for( int i = 0; i < max_rank; i++) pk[i] = safe_malloc( sizeof( int) * trellis_length); } } static int ti = 0; void decode_append( VTYPE val) { if( T >= trellis_length) bailout( "max trellis length reached"); // Collect signal values until we have a complete tuple static TUPLE TP; TP[QD - 1 - ti] = val; // Bits arrive most-significant first if( ++ti < QD) return; // Not yet a full tuple // Insert tuple into the decoder input memcpy( &decoder_input[T], TP, sizeof( TUPLE)); T++; // Reset for the next tuple ti = 0; memset( TP, 0, sizeof( TUPLE)); } static void reset_decoder( void) { T = 0; ti = 0; // Set up the first column of the trellis. The -1 for v1 and v2 indicate // that those states are never reached. trellis[0][0].M1 = 0; // Encoder always begins in state 0 for( int i = 1; i < QS; i++) { trellis[i][0].v1 = -1; trellis[i][0].v2 = -1; } // Reset the list decoder stack vstack_reset(); } /////////////////////////////////////////////////////////////////////////////// // Unencoded detection // /////////////////////////////////////////////////////////////////////////////// // If the signal is unencoded (option -p1) we have QK = QD = 1 and Viterbi // decoding (which relies on state transitions) doesn't work because the // encoder has no state (or rather, it has a constant state). // // With -p1 this function is called to convert the signal to source bits. // // A 'path metric' is still computed and returned. This is used by the phase // search to select the most likely message from the phase search results. // // No list decoding is possible but a CRC can be used if desired. This can // be helpful during the phase search. // // You can use repetition codes, eg -p1,1 or -p1,1,1 and so on. // static double uncoded_detection( void) { double m = 0; double tsum = 0; // Temporary sum over the repetitions (if any, QD > 1) // of the source bits int tcnt = 0; for( int i = 0; i < Nsymbols; i++) { tsum += sigbits[i]; if( ++tcnt == QD) { int bit = tsum > 0 ? +1 : -1; m -= bit * tsum; source_bits[i/QD] = bit > 0 ? 1 : 0; tsum = 0; tcnt = 0; } } return m; } /////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////// static struct RESULT { int valid; char message[MAX_CHARS+1]; double ebn0; double esn0; int rank; int errs; int syms; int weight; double m; double ber; double effrate; double shcap; double carrier_phase; double carrier_amp; double carrier_sn; double carrier_uhz; double carrier_esn0; double carrier_ebn0; double info_rate; char phasetxt[100]; } result; static void write_errmap( char *errmap) { FILE *f = fopen( "errmap", "w"); for( int i = 0; i < Nsymbols; i++) { char e = errmap[interleave( i)]; fprintf( f, "%d %d\n", i, e); } fclose( f); } static void analyse_signal( void) { // // Undo any normalisation of the input signal // if( anorm) for( int i = 0; i < Nsymbols; i++) rotamps[interleave( i)] *= normf[i]; // // Re-encode the message and reverse the modulation of the received signal. // Sum the carrier amplitudes, counter symbol errors and build the error // map array. // int n = 0; char *errmap = safe_malloc( Nsymbols); complex double carrier_sum = 0; // For coherent sum of the complete carrier result.errs = result.weight = 0; reset_encoder( 0); for( int i = 0; i < LB; i++) { encode( source_bits[i]); for( int j = QD-1; j >= 0; j--) { int bit = outword & (1<= 0) { errmap[n] = 1; result.errs++; } if( bit) result.weight++; n++; } } for( int i = 0; i < QK-1; i++) { encode( 0); for( int j = QD-1; j >= 0; j--) { int bit = outword & (1<= 0) { errmap[n] = 1; result.errs++; } if( bit) result.weight++; n++; } } result.carrier_phase = -carg( carrier_sum); result.carrier_amp = cabs( carrier_sum)/Nsymbols/sample_rate/sym_period; double duration = Nsymbols * sym_period; double path_metric = creal(carrier_sum); report( 2, "path metric: %.7e", path_metric); double carrier_amplitude = cabs(carrier_sum)/(duration * sample_rate); double carrier_power = carrier_amplitude * carrier_amplitude; report( 2, "carrier amplitude: %.3e", carrier_amplitude); report( 2, "carrier power: %.3e", carrier_power); // Total power of the input signal double total_power = input_sum_squared / (duration * sample_rate); report( 2, "total power: %.3e", total_power); // Total noise power of the input signal double total_noise_power = total_power - carrier_power; report( 2, "total noise power: %.3e", total_noise_power); // Total input BW: two channels, each with bandwidth sample_rate/2 double input_bw = 2 * sample_rate/2; // One sided noise spectral density double noise_density = total_noise_power/input_bw; double noise_bw = 1/duration; double noise_in_carrier_bw = total_noise_power * noise_bw/input_bw; double info_bit_period = duration / (Nchars * 5.7); result.carrier_sn = 10 * log10( carrier_power/noise_in_carrier_bw); result.carrier_uhz = 1e6/duration; double ebn0_sn = carrier_power * info_bit_period/noise_density; double esn0_sn = carrier_power * sym_period/noise_density; result.carrier_ebn0 = 10 * log10( ebn0_sn); result.carrier_esn0 = 10 * log10( esn0_sn); // Estimate Es/N0 from the symbol error rate result.ber = result.errs/(double) Nsymbols; // Symbol error rate double esn0 = 0; // Es/N0 in dB for( esn0 = -30; esn0 < 60 ; esn0 += 0.1) { double xber = 0.5 * erfc( sqrt( pow( 10, esn0/10))); if( xber < result.ber) break; } // Overall code rate based on 5.7 bits of info per character result.effrate = Nchars * 5.7/(double) n; result.esn0 = esn0; result.ebn0 = result.esn0 + 10 * log10( 1/result.effrate); double snr = pow( 10, result.carrier_esn0/10); result.shcap = 3600 * 1/(double)sym_period * log2(1 + snr); result.info_rate = Nchars * 5.7/duration * 3600; result.valid = TRUE; if( !FFLAG9 && !FFLAG10 && !FFLAG18) write_errmap( errmap); free( errmap); } static void report_message( int rank, double m) { // // If the current decode is less likely than an earlier one, // then discard it. // if( result.valid && !FFLAG9 && !FFLAG10 && m > result.m) return; // // Save the result, overwriting any previous, and analyse the signal. // for( int i = 0; i < Nchars; i++) result.message[i] = decode_char( source_bits + 6*i); result.message[Nchars] = 0; result.rank = rank; result.m = m; if( !TFLAG) analyse_signal(); // // Output the result. With -f10 (report all list outputs) use a simplified // one line format. // if( !FFLAG10 && !TFLAG) { printf( "found rank %d ber %.4e Eb/N0 %.1f M %.9e", result.rank, result.ber, result.ebn0, result.m); if( PSTEP) printf( " ph %s", p_phase_txt); printf( " [%s]\n", result.message); report( 1, "carrier phase: %.1f deg", result.carrier_phase * 180/M_PI); report( 1, "carrier amplitude: %.3e", result.carrier_amp); report( 1, "carrier Eb/N0: %.1f dB", result.carrier_ebn0); report( 1, "carrier Es/N0: %.2f dB", result.carrier_esn0); report( 1, "carrier S/N: %.2f dB in %.1f uHz," " %.2f dB in 1Hz, %.2f dB in 2.5kHz", result.carrier_sn, result.carrier_uhz, result.carrier_sn - 10 * log10(1e6/result.carrier_uhz), result.carrier_sn - 10 * log10(2.5e9/result.carrier_uhz)); report( 2, "bits per hour: %.1f", result.info_rate); report( 2, "Shannon capacity: %.1f% bits/hour", result.shcap); report( 2, "Shannon efficiency: %.1f%%", 100.0 * result.info_rate/result.shcap); } if( !FFLAG10) { report( 1, "elapsed %d", (int)(time(NULL) - Tstart)); report( 2, "Code weight: %d", result.weight); } if( FFLAG10) printf( "%d %d %d %.5f %.8e\n", rank, result.weight, result.errs, result.ber, -m); if( fflush( stdout)) bailout( "output file closed"); } static void run_decode( void) { // // Make sure max_rank doesn't exceed the number of possible // inner code messages. // if( LB < 31 && max_rank > 1 << LB) { max_rank = 1 << LB; report( 0, "list length reduced to %d", max_rank); } Tstart = time(NULL); init_viterbi(); vstack_init(); if( TFLAG) receive_block_text(); else receive_block(); result.valid = FALSE; for( int ph_idx = PSTEP_INIT; !PSTEP || prepare_phase( ph_idx); ph_idx++) { if( QS == 1) // No encoding so run uncoded detection { double m = uncoded_detection(); if( check_crc() && valid_message()) if( !result.valid || m < result.m) report_message( 0, m); } else // Run Viterbi decoding { reset_decoder(); for( int i = 0; i < Nsymbols; i++) decode_append( sigbits[i]); fwd_viterbi(); for( int rank = 0; rank < max_rank; rank++) { double m = backward_pass( rank); if( result.valid && !FFLAG9 && !FFLAG10 && m > result.m) break; backtrace( rank); if( FFLAG10 || (check_crc() && valid_message())) { report_message( rank, m); if( !FFLAG9 && !FFLAG10) break; } } } if( !PSTEP) break; } report( 1, "vstack: max height %d", maxheight); } // // Check for a catastrophic polynomial set. // static int cat_check( void) { init_viterbi(); for( int k = 1; k < QS; k++) // For each encoder state { int b = 0; // Branch count int i = k; while( b++ < QS) { if( STATE[i].emit[0] == 0) i = STATE[i].next[0]; else if( STATE[i].emit[1] == 0) i = STATE[i].next[1]; else break; if( i == k) { printf( "catastrophic cycle from state %d b=%d\n", k, b - 1); return 1; } } } return 0; } static void run_encode( void) { // // If noise is requested, set the noise amplitude to produce the // requested Eb/N0. Signal power and RMS amplitude are always unity, // and the noise sigma is calculated to give the required noise density. // if( NFLAG) { double duration = Nsymbols * sym_period; double info_bits = LB - CRCLEN; // Allow for the unused 0.3 bits per character if using source encoding info_bits *= 5.7/6.0; double info_bit_period = duration / info_bits; double signal_power = 1; double energy_per_bit = signal_power * info_bit_period; double noise_density = energy_per_bit / request_ebn0; double total_bw = sample_rate/2.0; // Bandwidth of each channel, I and Q double noise_power = noise_density * total_bw; // Noise power per chan noise_sigma = sqrt( noise_power); double signal_sn = signal_power / noise_density; report( 1, "Carrier S/N %.2f dB in 1Hz", 10 * log10(signal_sn)); report( 1, "noise sigma %.4e", noise_sigma); } if( PRAND) { refP = drand48() * 2 * M_PI; report( 1, "random phase %.1f", refP * 180/M_PI); } load_source(); reset_encoder( 0); int nb = 0; for( int i = 0; i < LB; i++) { encode( source_bits[i]); for( int j = QD-1; j >= 0; j--) sigbits[nb++] = outword & (1<= 0; j--) sigbits[nb++] = outword & (1<alias && strcasecmp( list->alias, arg)) list++; if( list->alias) p = list->poly; QD = 0; QK = 0; while( 1) { uint32_t poly = (uint32_t) strtol( p, &q, 0); if( p == q) break; QP[QD++] = poly; int i = 0; do { poly >>= 1; i++; } while( poly != 0); if( i > QK) QK = i; p = q; if( *p == 0) break; if( *p != ',') { fprintf( stderr, "error parsing [%s]\n", arg); exit( 1); } p++; } } static void usage( void) { fprintf( stderr, "ebnaut version %s\n" "usage: ebnaut [options]\n" "\n" "options:\n" "\n" " -v Increase verbosity\n" " -c num Number of CPUs to use, default 1\n" "\n" " -e Encode\n" " -d Decode\n" " -N nchars Message length in characters (no default)\n" " -p poly Specify convolution polynomial set or alias\n" " -k crclen CRC length in bits, 8 to 32 (default 16)\n" "\n" " -t Output encoded bits in text format\n" " -s Do source coding only\n" " -f2 Output the minimum distance spectrum of the inner code\n" " -f3 Check for catastrophic code\n" " -f4 No decode, just measure reference phase\n" " -f6 Output interleave map\n" " -f7 Output code weight and length of cascaded code\n" " -f8 Distance spectrum, exhaustive, cascaded code\n" " -f9 Report all decodes (default the most likely)\n" " -f10 Report all list outputs\n" " -f13 Output code weight and length of inner code\n" " -f14 Distance spectrum, exhaustive, inner code\n" " -f15 Disable use of initial reference phase\n" " -f16 Analyse received signal\n" " -f17 Output polynomial(s)\n" " -f18 Do not write errmap file\n" " -f19 Disable use of redundancy\n" " -f20 Disable DSW rebalancing of list stack\n" "\n" " -M 'message text' Alternative to stdin\n" "\n" " -n dB Set Eb/N0\n" " -PS Do phase search\n" "\n" " -oc Output unmodulated carrier and noise\n" " -on Output noise only\n" "\n" " -r rate Sample rate\n" " -S secs Symbol period\n" " -F hertz Frequency offset, Hz\n" " -T secs Start offset, seconds\n" "\n" " -L max List decode up to max trials (default 20000)\n", VERSION ); exit( 1); } int main( int argc, char *argv[]) { while( 1) { int c = getopt( argc, argv, "vedstp:n:N:L:o:r:S:F:T:f:c:P:k:a:M:?X"); if( c == 'X') XFLAG = TRUE; else if( c == 'v') loglevel++; else if( c == 's') SFLAG = 1; else if( c == 'N') Nchars = atoi( optarg); else if( c == 'c') ncpu = atoi( optarg); else if( c == 'L') { max_rank = atoi( optarg); if( max_rank > MAX_RANK) bailout( "-L too large, max %d", MAX_RANK); if( max_rank <= 0) max_rank = 1; } else if( c == 'r') sample_rate = atof( optarg); else if( c == 'P') { if( !strcmp( optarg, "U") || !strcmp( optarg, "U30")) { PSTEP = 30; PU = TRUE; } else if( !strcmp( optarg, "U15")) { PSTEP = 15; PU = TRUE; } else if( !strcmp( optarg, "S") || !strcmp( optarg, "S30")) PSTEP = 30; else if( !strcmp( optarg, "S15")) PSTEP = 15; else if( !strcmp( optarg, "R")) PRAND = TRUE; else if( !strncmp( optarg, "N", 1)) PSTEP_INIT = atoi( optarg+1); else { refP = atof( optarg) * M_PI/180; PFLAG = 1; } } else if( c == 'f') { int flag = atoi( optarg); switch( flag) { case 2: FFLAG2 = TRUE; break; case 3: FFLAG3 = TRUE; break; case 4: FFLAG4 = TRUE; break; case 6: FFLAG6 = TRUE; break; case 7: FFLAG7 = TRUE; break; case 8: FFLAG8 = TRUE; break; case 9: FFLAG9 = TRUE; break; case 10: FFLAG10 = TRUE; break; case 13: FFLAG13 = TRUE; break; case 14: FFLAG14 = TRUE; break; case 15: FFLAG15 = TRUE; break; case 16: FFLAG16 = TRUE; break; case 17: FFLAG17 = TRUE; break; case 18: FFLAG18 = TRUE; break; case 19: FFLAG19 = TRUE; break; case 20: FFLAG20 = TRUE; break; default: bailout( "unrecognised -f option %d", flag); } } else if( c == 'o') { if( !strcmp( optarg, "c")) opmode |= OP_CARR; else if( !strcmp( optarg, "n")) opmode |= OP_OFF; else bailout( "unrecognised output option"); } else if( c == 'n') { char *p; if( (p = strchr( optarg, 's')) != NULL) { long int n = atol( p+1); if( !n) n = time(NULL) * getpid(); srand48( n); report( 1, "seeding drand48 %ld", n); } request_ebn0 = pow( 10, atof(optarg)/10); NFLAG = 1; } else if( c == 'F') freq_offset = atof( optarg); else if( c == 'T') start_offset = atof( optarg); else if( c == 'e') mode = ENCODE; else if( c == 'd') mode = DECODE; else if( c == 'p') parse_polys( optarg); else if( c == 'S') sym_period = atof( optarg); else if( c == 'k') CRCLEN = atoi( optarg); else if( c == 't') TFLAG = TRUE; else if( c == 'a') anorm = atoi( optarg); else if( c == 'M') mtext = strdup( optarg); else if( c == -1) break; else { report( 0, "unknown option %c", c); usage(); } } if( ncpu > MAXCPU) bailout( "too many CPU with -c, max %d", MAXCPU); if( !Nchars && !FFLAG2 && !FFLAG3 && !FFLAG17) bailout( "needs -N option to specify number of characters"); if( CRCLEN) if( CRCLEN < CRCBASE || CRCLEN > CRCMAX) bailout( "invalid CRC length"); if( mode != ENCODE && PRAND) bailout( "-PR only valid with -e"); if( SFLAG) // Do source encoding only? { load_source(); for( int i = 0; i < Nchars * 6; i++) printf( "%d", source_bits[i]); printf( "\n"); return 0; } if( FFLAG17) { if( QD) { for( int i = 0; i < QD; i++) printf( "%s0%o", i ? "," : "", QP[i]); printf( "\n"); } else for( int i = 0; polylist[i].alias; i++) printf( "%s %s\n", polylist[i].alias, polylist[i].poly); return 0; } if( !QK) bailout( "no polynomials given, needs -p"); QS = 1 << (QK-1); // Number of encoder states; LB = Nchars * 6 + CRCLEN; // Viterbi block size if( LB > MAX_BITS) bailout( "too many bits to encode"); Nsymbols = QD * (LB + QK - 1); // Output block size report( 2, "Nsymbols = %d", Nsymbols); report( 3, "QK=%d QS=%d QD=%d LB=%d Nsymbols=%d", QK, QS, QD, LB, Nsymbols); IW = ceil( sqrt( Nsymbols)); if( FFLAG6) // -f6: Output interleave map { for( int i = 0; i < Nsymbols; i++) printf( "%d %d\n", i, interleave( i)); return 0; } // // -f7: Encode random text and output the weight of the cascaded code. // if( FFLAG7) { while( 1) { gen_random_text(); reset_encoder( 0); for( int i = 0; i < LB; i++) encode_count( source_bits[i]); for( int i = 0; i < QK-1; i++) encode_count( 0); printf( "%d\n", weight); } return 0; // Never reached } // // -f8: Output the complete distance (weight) spectrum of the cascaded code // if( FFLAG8) { uint64_t nmax = round( pow( 52, Nchars)); report( 0, "nmax %Lu", nmax); uint64_t counts[5000]; memset( counts, 0, sizeof( counts)); int max_weight = 0, min_weight = 0; for( uint64_t n = 0; n < nmax; n++) { gen_seq_text( n); reset_encoder( 0); for( int i = 0; i < LB; i++) encode_count( source_bits[i]); for( int i = 0; i < QK-1; i++) encode_count( 0); counts[weight]++; if( weight > max_weight) max_weight = weight; if( weight < min_weight || !min_weight) min_weight = weight; } for( int i = min_weight; i <= max_weight; i++) printf( "%d %lu\n", i, (long) counts[i]); return 0; } // // -f13: Encode random source and output the weight of the inner code. // if( FFLAG13) { while( 1) { gen_random_source(); reset_encoder( 0); for( int i = 0; i < LB; i++) encode_count( source_bits[i]); for( int i = 0; i < QK-1; i++) encode_count( 0); printf( "%d\n", weight); } return 0; // Never reached } // // -f14: Output the complete distance (weight) spectrum of the inner code // if( FFLAG14) { uint64_t nmax = (uint64_t) 1 << LB; report( 0, "Nchars %d LB %d nmax %Lu", Nchars, LB, nmax); uint64_t counts[5000]; memset( counts, 0, sizeof( counts)); int max_weight = 0, min_weight = 0; for( uint64_t n = 0; n < nmax; n++) { gen_seq_source( n); reset_encoder( 0); for( int i = 0; i < LB; i++) encode_count( source_bits[i]); for( int i = 0; i < QK-1; i++) encode_count( 0); counts[weight]++; if( weight > max_weight) max_weight = weight; if( weight < min_weight || !min_weight) min_weight = weight; } for( int i = min_weight; i <= max_weight; i++) printf( "%d %lu\n", i, (long) counts[i]); return 0; } if( FFLAG3) return cat_check(); // -f3: check for catastrophic code if( FFLAG2) // -f2: Output weight spectrum { report( 1, "QK=%d", QK); int lim = 1 << QK; uint64_t counts[5000]; memset( counts, 0, sizeof( counts)); int max_weight = 0, min_weight = 0; for( int i = 1; i < lim; i++) { reset_encoder( 0); uint32_t u = i; for( int j = 0; j < QK; j++, u >>= 1) encode_count( u & 1); for( int j = 0; j < QK-1; j++) encode_count( 0); counts[weight]++; if( weight > max_weight) max_weight = weight; if( weight < min_weight || !min_weight) min_weight = weight; } for( int i = min_weight; i <= max_weight; i++) printf( "%d %lu\n", i, (long) counts[i]); return 0; } if( !mode && !FFLAG2 && !FFLAG3 && !FFLAG6 && !FFLAG7 && !FFLAG8 && !FFLAG13 && !FFLAG14 && !FFLAG16) bailout( "needs -d or -e"); if( !TFLAG) { if( !sym_period) bailout( "no symbol period given - needs -S"); if( !sample_rate) bailout( "no sample rate given - needs -r"); } sigbits = (double *) malloc( sizeof( double) * Nsymbols); if( mode == ENCODE) { run_encode(); return 0; } symamps = (complex double *) malloc( sizeof( complex double) * Nsymbols); rotamps = (complex double *) malloc( sizeof( complex double) * Nsymbols); refph = (complex double *) malloc( sizeof( complex double) * Nsymbols); normf = (double *) malloc( sizeof( double) * Nsymbols); if( FFLAG16) // -f16: Analyse signal for known message { // // -M supplies the known message. // Reverse the received modulation and measure the reconstructed // carrier. // if( TFLAG) report( 0, "-t ignored with -f16"); load_source(); receive_block(); analyse_signal(); report( 1, "carrier phase: %.1f", result.carrier_phase * 180/M_PI); report( 1, "carrier Eb/N0: %.1f dB", result.carrier_ebn0); report( 1, "carrier S/N: %.2f dB in %.1f uHz," " %.2f dB in 1Hz, %.2f dB in 2.5kHz", result.carrier_sn, result.carrier_uhz, result.carrier_sn - 10 * log10(1e6/result.carrier_uhz), result.carrier_sn - 10 * log10(2.5e9/result.carrier_uhz)); return 0; } if( mode == DECODE) run_decode(); return 0; }