// // ebnaut-rx.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 "0.9a" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include // Missing prototypes in the mingw headers int strcasecmp(const char *s1, const char *s2); // Missing from string.h char *strdup(const char *s); // Missing from string.h /////////////////////////////////////////////////////////////////////////////// // Globals // /////////////////////////////////////////////////////////////////////////////// static HWND hwnd, hwnd_period, hwnd_list, hwnd_chars, hwnd_ncpu, hwnd_sigfile, hwnd_result, hwnd_start, hwnd_freq, hwnd_pstep, hwnd_runstop, hwnd_status, hwnd_signal, hwnd_codesel, hwnd_browse, hwnd_crc; static int stop = FALSE; // Set TRUE to signal the decoder tasks to stop static double sym_period = 1; // Signal bit period, seconds static int Nchars = 5; // Message length, number of characters static char sigfile[500] = ""; // Input WAV filename static double sample_rate = 120.0; static double start_offset = 0; // Seconds into the input file of first symbol static double freq_offset = 0; static int autostart = FALSE; // Set TRUE by -as command line option static int autoexit = FALSE; // Set TRUE by -ax command line option static char *inf1_ut = NULL; // 'ut' field from WAV file info header static char *inf1_rf = NULL; // 'rf' field from WAV file info header static FILE *fin = NULL; // Handle for input file static char *poly_string = "8K19A"; // Selected polynomial set or alias static time_t Tstart = 0; // Decoder start time, for elapsed time reporting static int QD = 0; // Inner code, number of symbols per data bit, eg 8, 16, etc static int LB = 0; // Viterbi block size, Nchars * 6 + CRCLEN static int Nsymbols = 0; // Number of symbols per block static double *sigbits = NULL; // Receive signal bit NRZ amplitudes static complex double *symamps = NULL; // Complex amplitude of each symbol // after frequeny offset applied static complex double *rotamps = NULL; // Complex amplitude of each symbol // after frequency offset, phase // pattern, and de-interleave are // applied static char *logfile = NULL; // Filename of log file static int IW = 0; // Interleave stride #define MAX_CHARS 1000 // Maximum number of characters #define MAX_BITS 6050 static uint8_t source_bits[MAX_BITS+1]; // Source-encoded message bits static int ncpu = 1; // Number of CPUs to use #define MAXCPU 64 typedef float VTYPE; // Floating point type for the Viterbi trellis /////////////////////////////////////////////////////////////////////////////// // Predefined Polynomial Sets // /////////////////////////////////////////////////////////////////////////////// struct POLYLIST { char *alias; char *poly; } polylist[] = { { "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 // /////////////////////////////////////////////////////////////////////////////// // The 32 bit library doesn't contain asprintf or vasprintf. int vasprintf( char **s, const char * format, va_list va) { *s = NULL; int n = _vsnprintf( *s, 0, format, va); if( n < 0) return -1; *s = malloc( n + 100); if( !*s) return -1; return vsprintf( *s, format, va); } int asprintf( char **, const char *, ...) __attribute__ ((format (printf, 2, 3))); int asprintf( char **s, const char * format, ...) { va_list va; va_start( va, format); int r = vasprintf( s, format, va); va_end( va); return r; } // // 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); } // // Append a message to the log file. The file is opened and closed for // each message. // static int loglevel = 0; static void report( int, char *, ...) __attribute__ ((format (printf, 2, 3))); static void report( int n, char *format, ...) { if( n > loglevel) return; char *text = NULL; va_list ap; va_start( ap, format); vasprintf( &text, format, ap); va_end( ap); fprintf( stderr, "%s\n", text); if( logfile) { FILE *f = fopen( logfile, "at"); if( f) { fprintf( f, "%s\n", text); fclose( f); } } free( text); } // // Stop the decode and show an error message. // static void show_run( void); static void error_stop( char *format, ...) { if( !stop) { va_list ap; char *text = NULL; va_start( ap, format); vasprintf( &text, format, ap); va_end( ap); char *message = NULL; asprintf( &message, "Error: %s", text); SetWindowText( hwnd_status, message); report( 0, "ERROR: %s", text); free( message); free( text); show_run(); } ExitThread( 0); } // // Display one or two line decoder status message. // static void show_status( char *format, ...) { if( stop) return; va_list ap; char *text = NULL; va_start( ap, format); vasprintf( &text, format, ap); va_end( ap); SetWindowText( hwnd_status, text); free( text); } // // Enable or disable the setup controls. // static void enable_setup( int enable) { EnableWindow( hwnd_codesel, enable); EnableWindow( hwnd_crc, enable); EnableWindow( hwnd_period, enable); EnableWindow( hwnd_list, enable); EnableWindow( hwnd_sigfile, enable); EnableWindow( hwnd_chars, enable); EnableWindow( hwnd_start, enable); EnableWindow( hwnd_freq, enable); EnableWindow( hwnd_ncpu, enable); EnableWindow( hwnd_browse, enable); EnableWindow( hwnd_pstep, enable); } /////////////////////////////////////////////////////////////////////////////// // Source Coding // /////////////////////////////////////////////////////////////////////////////// // // 6 source coded payload bits per character. // // * (dec 42) to Z (dec 90) coded as 0 to 48 // '\n' coded as 49 but decoded to a space // space and tab both coded to 50 // underscore codes as 51 // // Codes 52 to 63 are un-assigned and contribute strength to the outer error // detection layer. // 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 } /////////////////////////////////////////////////////////////////////////////// // Interleaving // /////////////////////////////////////////////////////////////////////////////// // // A simple rectangular (approximately square) interleaving. // // For offset 'in' in the transmission order, return the offset into the // non-interleaved data. // static int interleave( int in) { if( in >= Nsymbols) bailout( "interleave input out of range"); int H = Nsymbols / IW; // Number of rows int L = Nsymbols % 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; } /////////////////////////////////////////////////////////////////////////////// // 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 = 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; 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) { for( int i = 0; i < LB - CRCLEN; i += 6) if( decode_char( source_bits + i) == '#') return FALSE; return TRUE; } /////////////////////////////////////////////////////////////////////////////// // WAV File // /////////////////////////////////////////////////////////////////////////////// static char * format_timestamp( long double T) { static char temp[50]; time_t sec = (time_t) T; struct tm *tm = gmtime( &sec); int d = 0.5 + 1e3 * (T - sec); if( d == 1000) { d = 0; sec++; } sprintf( temp, "%04d-%02d-%02d %02d:%02d:%02d.%03d", tm->tm_year + 1900, tm->tm_mon+1, tm->tm_mday, tm->tm_hour, tm->tm_min, tm->tm_sec, d); return temp; } static void open_wav( void) { fin = fopen( sigfile, "rb"); if( !fin) error_stop( "Cannot open %s, %s", sigfile, strerror( errno)); size_t n; char buff[500]; if( fread( buff, 1, 20, fin) != 20) error_stop( "input file read error %s", strerror(errno)); if( strncmp( buff, "RIFF", 4)) error_stop( "input file incorrect format (1)"); if( strncmp( buff+8, "WAVEfmt ", 8)) error_stop( "input file incorrect format (2)"); if( *(int32_t *)(buff+16) != 40) error_stop( "input file read error %s", strerror(errno)); if( fread( buff, 1, 40, fin) != 40) error_stop( "input file read error %s", strerror(errno)); // if( *(uint16_t *) buff != 1) // Compression code // bailout( "input file incorrect format (4)"); int chans = *(uint16_t *)(buff + 2); if( chans != 2) error_stop( "Input file must have two channels"); sample_rate = *(uint32_t *)(buff + 4); // Skip over 'fact' chunk if( fread( buff, 1, 8, fin) != 8) error_stop( "input file read error %s", strerror(errno)); if( strncmp( buff, "fact", 4)) error_stop( "input file incorrect format (5)"); n = *(uint32_t *)(buff + 4); if( fread( buff, 1, n, fin) != n) error_stop( "input file read error %s", strerror(errno)); if( fread( buff, 1, 8, fin) != 8) error_stop( "input file read error %s", strerror(errno)); // Read 'inf1' chunk if( strncmp( buff, "inf1", 4)) error_stop( "input file incorrect format (6)"); n = *(uint32_t *)(buff+4); if( fread( buff, 1, n, fin) != n) error_stop( "input file incorrect format (7)"); report( 1, "inf1 [%s]", buff); if( inf1_ut) free( inf1_ut); if( inf1_rf) free( inf1_rf); char *p, *q, *ut = NULL; for( p = buff; p && *p; ) { q = strchr( p, ' '); if( q) *q++ = 0; if( !strncmp( p, "ut=", 3)) ut = p + 3; else if( !strncmp( p, "rf=", 3)) inf1_rf = strdup( p + 3); else if( !strncmp( p, "sr=", 3)) sample_rate = atof( p + 3); p = q; } // Read data chunk header if( fread( buff, 1, 8, fin) != 8) error_stop( "input file read error %s", strerror(errno)); if( strncmp( buff, "data", 4)) error_stop( "input file incorrect format (8)"); char temp[500] = ""; sprintf( temp, "Sample rate: %.6f/sec ", sample_rate); if( inf1_rf) sprintf( strchr( temp, 0), "Rx freq: %s Hz\n", inf1_rf); if( ut) { long double T = strtold( ut, NULL); inf1_ut = strdup( format_timestamp( T)); sprintf( strchr( temp, 0), "File start: %s\n", inf1_ut); } if( !stop) SetWindowText( hwnd_signal, temp); } /////////////////////////////////////////////////////////////////////////////// // Demodulation Layer // /////////////////////////////////////////////////////////////////////////////// static double input_sum_squared = 0; // // Read a pair of I/Q samples. Trusts that the input file has 32 bit IEEE // format samples. // static int recv_analogue( complex double *v) { char buff[8]; if( fread( buff, 1, 8, fin) != 8) return 0; double vi, vq; vi = *(float *)(buff + 0); vq = *(float *)(buff + 4); *v = vi + I * vq; return 1; } 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, where 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; } // // Initialise the demodulator and receive Nsymbols from the input file. // static int 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); fclose( fin); fin = NULL; 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); // // Output a map of the raw received symbols. // FILE *fd = fopen( "rawsyms.txt", "wt"); 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. // complex double sumsq = 0; for( int i = 0; i < Nsymbols; i++) sumsq += symamps[i] * symamps[i]; complex double ref = csqrt(sumsq); report( 0, "initial reference phase %.1f", -180 * carg(ref)/M_PI); // // Rotate input signal to zero the reference phase. // if( cabs( ref)) { ref /= cabs( ref); for( int i = 0; i < Nsymbols; i++) symamps[i] /= ref; } return TRUE; } /////////////////////////////////////////////////////////////////////////////// // Reference Phase Search // /////////////////////////////////////////////////////////////////////////////// static char p_phase_txt[100]; // Text reporting the current phase pattern // // Set up reference phase pattern n, 0 <= n < 108. Returns FALSE if there // are no more patterns to try. // static int PSTEP = 30; // Or 15 for a more fine-grained search static int PU = FALSE; // TRUE with -PU option for uniform phase search 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; // All done 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", ph[0], ph[1], ph[2], ph[3]); // // Apply the reference phase pattern to the raw symbols symamps[] and // de-interleave to produce rotamps[]. Extract the real component into // sigbits[] for the decoder to use. // 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); rotamps[interleave( i)] = symamps[i] * a; sigbits[interleave( i)] = creal( symamps[i] * a); } double progress = ns * 100.0/(216 * PSTEP/15); show_status( "Running... %.0f%% complete\nReference phase pattern: %s", progress, p_phase_txt); return TRUE; } /////////////////////////////////////////////////////////////////////////////// // Convolutional Encoder // /////////////////////////////////////////////////////////////////////////////// // // The encoder is used to create tables for the decoder to use, and // to re-encode a decoded message to compare with the received signal. // static int QK = 0; // Encoder length: 1 + number of memory bits static int QS; // Number of encoder states static int QP[16]; // Up to 16 polynomials, set by -p option static uint32_t encstate; // Current encoder state static uint32_t outword = 0; // Decoder output word static inline void reset_encoder( uint32_t val) { encstate = val; } 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; } /////////////////////////////////////////////////////////////////////////////// // Viterbi Decoder // /////////////////////////////////////////////////////////////////////////////// // Semaphores for synchronising forward pass threads static HANDLE vt_cnt; // To synchronise launch of background work static HANDLE vt_sem[MAXCPU]; // To wait for background work complete static struct BMAP { uint32_t prev[2]; int bit[2]; } *BMAP; typedef VTYPE TUPLE[16]; typedef uint16_t ETYPE; // // Table to list the output pattern and next state for each encoder state // transition. // struct STATE { ETYPE emit[2]; int32_t next[2]; } *STATE; 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; static int T = 0; // Number of code tuples in the input buffer // // 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 && !stop; 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 { // Demote v1,M1 to v2,M2 and insert the new best into v1,M1 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; } } } } // // Work along the trellis calling the forward pass for each column. // // After each column is done, all the threads have to rendezvous before // advancing to the next column. // static DWORD WINAPI decode3_start( LPVOID arg) { int cpu = (intptr_t) arg; while( 1) { WaitForSingleObject( vt_sem[cpu], INFINITE); // Wait to begin next col if( vt > T) break; forward_pass( cpu); ReleaseSemaphore( vt_cnt, 1, NULL); // Signal that this thread is done if( stop) break; } return 0; } static void decode2( void) { vt = 1; HANDLE pids[MAXCPU]; // Create ncpu-1 extra threads to share the work of the forward pass for( int i = 1; i < ncpu; i++) pids[i] = CreateThread( NULL, 0, (LPTHREAD_START_ROUTINE) decode3_start, (LPVOID)(intptr_t) i, 0, 0); for( vt = 1; !stop && vt <= T; vt++) // For each trellis column { // Signal the other threads to start work on this column for( int i = 1; i < ncpu; i++) ReleaseSemaphore( vt_sem[i], 1, NULL); forward_pass( 0); // Wait for the other threads to complete this column for( int i = 1; i < ncpu; i++) WaitForSingleObject( vt_cnt, INFINITE); } // Release the other threads so that they can exit for( int i = 1; i < ncpu; i++) ReleaseSemaphore( vt_sem[i], 1, NULL); // Wait for them to exit for( int i = 1; i < ncpu; i++) WaitForSingleObject( pids[i], INFINITE); } static int ti = 0; // Count of bits added so far to the current input tuple void decode_append( VTYPE val) { // 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)); if( T >= trellis_length) bailout( "max trellis length reached"); } static void reset_viterbi( 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; } } /////////////////////////////////////////////////////////////////////////////// // 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. // #define MAX_RANK 20000000 // Maximum list size static int max_rank = 20000; // The default list size 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 = NULL; static int nvstack = 0; static void vstack_init( void) { vpool = realloc( vpool, max_rank * sizeof( struct VSTACK)); if( !vpool) error_stop( "Not enough memory (vpool)"); } 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( 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 (their variable k) starting at zero // instead of one. // // Array of previously discovered paths. There's a way to avoid storing // the whole path but I haven't bothered to implement that. // // It would use less space if it was dynamically allocated. // static int32_t **pk = NULL; static int pk_alloc = 0; 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 { // Step 6 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) { for( int i = 0; i <= T - QK; i++) { int bit = 0; if( STATE[pk[rank][i]].next[0] == pk[rank][i+1]) bit = 0; else if( STATE[pk[rank][i]].next[1] == pk[rank][i+1]) bit = 1; else // Only occurs if there's a software error bailout( "illegal state transition"); source_bits[i] = bit; } } void init_viterbi( void) { static int trellis_QS = 0; static int trellis_length_alloc = 0; // Trellis contains one entry for the initial state plus LB info bits // plus QK-1 flushing bits trellis_length = LB + QK; double mem = trellis_length * (double) QS * sizeof( struct TRELLIS); mem += (sizeof( struct STATE) + sizeof( struct BMAP)) * (double) QS; mem += sizeof( struct TRELLIS) * (trellis_length + 1); mem += 380e6; report( 1, "decoder mem: %.3f Gb", mem/1e9); // If the trellis size has changed, release the old one if( trellis && (trellis_QS != QS || trellis_length_alloc != trellis_length)) { for( int i = 0; i < trellis_QS; i++) if( trellis[i]) free( trellis[i]); free( trellis); trellis = NULL; } trellis_length_alloc = trellis_length; // Allocate a new trellis if necessary if( !trellis) { trellis = malloc( sizeof( struct TRELLIS *) * QS); if( !trellis) error_stop( "Not enough memory (trellis)"); memset( trellis, 0, sizeof( struct TRELLIS *) * QS); for( int i = 0; i < QS; i++) { trellis[i] = malloc( sizeof( struct TRELLIS) * trellis_length); if( !trellis[i]) { while( --i >= 0) free( trellis[i]); free( trellis); trellis = NULL; error_stop( "Not enough memory (trellis)"); } } trellis_QS = QS; decoder_input = malloc( sizeof( TUPLE) * trellis_length); if( !decoder_input) error_stop( "Not enough memory (decoder input)"); } // 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 = realloc( STATE, 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 *) realloc( BMAP, sizeof( struct BMAP) * QS); memset( BMAP, 0, sizeof( struct BMAP) * QS); uint8_t *tmp = 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); // Reallocate the pk matrix if( pk) { for( int i = 0; i < pk_alloc; i++) if( pk[i]) free( pk[i]); free( pk); pk_alloc = 0; pk = NULL; } pk = malloc( max_rank * sizeof( int *)); if( !pk) error_stop( "Not enough memory (pk)"); memset( pk, 0, max_rank * sizeof( int *)); for( int i = 0; i < max_rank; i++) { pk[i] = malloc( sizeof( int) * trellis_length); if( !pk[i]) { while( --i >= 0) free( pk[i]); free( pk); pk = NULL; error_stop( "Not enough memory (pk)"); } } pk_alloc = max_rank; } // // A structure holding the current best decode. // static struct RESULT { int valid; char message[MAX_CHARS+1]; double ebn0; double esn0; int rank; int errs; int weight; double m; double ber; double effrate; double shcap; double carrier_phase; double carrier_sn; double carrier_uhz; double carrier_esn0; double carrier_ebn0; double info_rate; char phasetxt[100]; } result; static void show_result( void) { if( stop) return; if( !result.valid) { SetWindowText( hwnd_result, ""); return; } char *display_text = NULL; asprintf( &display_text, "Message: %s\n" "Rank: %d Es/N0: %.1fdB Eb/N0: %.1fdB\n" "Symbol errors: %d/%d BER: %.1f %%\n" "Reference phase: %s\n" "Carrier S/N %.2f dB in %.1f uHz, carrier Eb/N0: %.1f dB\n" "Info rate: %.2f bits/hour, %.1f %% of Shannon capacity\n", result.message, result.rank, result.esn0, result.ebn0, result.errs, Nsymbols, 100.0 * result.errs/(double)Nsymbols, result.phasetxt, result.carrier_sn, result.carrier_uhz, result.carrier_ebn0, result.info_rate, result.info_rate/result.shcap * 100.0); SetWindowText( hwnd_result, display_text); free( display_text); } /////////////////////////////////////////////////////////////////////////////// // Report Decode // /////////////////////////////////////////////////////////////////////////////// // // This analyses a decoded message by re-encoding it and comparing it with // the received signal. // static int reencode( char *msg, int rank, double m) { // // If the current decode is less likely than an earlier one, // then discard it. // if( result.valid && m > result.m) return FALSE; int n = 0; char *errmap = malloc( Nsymbols); char *bitmap = malloc( Nsymbols); complex double carrier_sum = 0; // For coherent sum of the complete carrier result.errs = result.weight = 0; // // Re-encode and compare with the signal in rotamps[] to determine the // total number of symbol errors. // // As we go through, sum the complex carrier amplitude after the modulation // is undone. // reset_encoder( 0); for( int i = 0; i < LB; i++) // For each payload bit... { encode( source_bits[i]); // Encode to produce QD bits in outword for( int j = QD-1; j >= 0; j--) { int bit = outword & (1<= 0) { errmap[n] = 1; result.errs++; } n++; } } for( int i = 0; i < QK-1; i++) // For each flushing bit... { encode( 0); for( int j = QD-1; j >= 0; j--) { int bit = outword & (1<= 0) { errmap[n] = 1; result.errs++; } n++; } } // // Calculate the performance of the decode and update the results // structure. // strcpy( result.message, msg); result.rank = rank; result.m = m; report( 0, "---------------------------------------------"); report( 0, "found %s", result.message); report( 0, "list rank %d", rank); strcpy( result.phasetxt, p_phase_txt); report( 0, "reference phase %s", result.phasetxt); result.carrier_phase = -carg( carrier_sum); double duration = Nsymbols * sym_period; 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); report( 0, "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( 0, "carrier Es/N0 %.2f dB", result.carrier_esn0); report( 0, "carrier Eb/N0 %.2f dB", result.carrier_ebn0); report( 0, "info bit period %.2f seconds", info_bit_period); result.ber = result.errs/(double) Nsymbols; // Symbol error rate // // Estimate Es/N0 from symbol error rate. Grossly inaccurate for low bit // errors and fails completely if there are no errors. But for typical // signals it does quite well. // 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; } result.esn0 = esn0; // Overall code rate based on 5.7 bits of info per character result.effrate = Nchars * 5.7/(double) Nsymbols; result.ebn0 = result.esn0 + 10 * log10( 1/result.effrate); result.valid = TRUE; 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; show_result(); report( 0, "symbol error rate %d/%d = %.3f %%", result.errs, Nsymbols, 100.0 * result.errs/(double) Nsymbols); if( result.errs) { report( 0, "Es/N0 from symbol errors %.1f dB", result.esn0); report( 0, "Eb/N0 from symbol errors %.1f dB", result.ebn0); } report( 0, "Shannon capacity %.1f bits/hour", result.shcap); report( 0, "Shannon efficiency %.1f %%", 100.0 * result.info_rate/result.shcap); // // Output the symbol file. // char *name_prefix = strdup( sigfile); char *p = strrchr( name_prefix, '.'); if( p) *p = 0; char *symbols_filename = malloc( strlen(name_prefix) + 12 + 1); sprintf( symbols_filename, "%s-symbols.csv", name_prefix); report( 0, "symbols file %s", symbols_filename); FILE *f = fopen( symbols_filename, "w"); for( int i = 0; i < Nsymbols; i++) { fprintf( f, "%d,", i); // Symbol number int j = interleave( i); complex double a = symamps[i]; fprintf( f, "%.6e,%.6e,%.1f,", creal(a), cimag(a), carg(a) * 180/M_PI); fprintf( f, "%d,", errmap[j]); fprintf( f, "%d\n", bitmap[j]); } fclose( f); report( 0, "elapsed %d seconds", (int)(time(NULL) - Tstart)); free( errmap); free( bitmap); free( symbols_filename); return TRUE; } // // Handle a successful decode: Reverse the source coding and call reencode() // to evaluate the performance and store the result. // static void report_message( int rank, double m) { char message[MAX_CHARS + 1]; int i, n; for( i = n = 0; i < LB - CRCLEN; i += 6, n++) message[n] = decode_char( source_bits + i); message[n] = 0; if( n == 0) return; reencode( message, rank, m); } // // Main sequence of the decoder. // static void run_decoder( void) { report( 0, "polynomial %s", poly_string); report( 0, "crc size %d", CRCLEN); report( 0, "number of chars %d", Nchars); report( 0, "block size %d", LB); report( 0, "symbol period %.3f", sym_period); report( 0, "number of symbols %d", Nsymbols); report( 0, "start offset %.4f", start_offset); report( 0, "freq offset %.6f", freq_offset); report( 0, "cores %d", ncpu); // // 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); } else report( 0, "list length %d", max_rank); // // Initialise everything. // reset_encoder( 0); init_viterbi(); vstack_init(); if( stop) return; // // Read in the audio from the input file and demodulate. // show_status( "Reading signal..."); if( !receive_block()) { error_stop( "No input or input too short"); return; } if( stop) return; // // Work through all the reference phase trial patterns. // show_status( "Running..."); for( int ph_idx = 0; prepare_phase( ph_idx); ph_idx++) { // Reset for the next pass reset_viterbi(); vstack_reset(); // Forward pass for( int i = 0; i < Nsymbols; i++) decode_append( sigbits[i]); decode2(); if( stop) break; // Backward passes, up to the list length for( int rank = 0; rank < max_rank && !stop; rank++) { double m = backward_pass( rank); // Give up if there are no more in the list better than the best // decode so far if( result.valid && m > result.m) break; backtrace( rank); if( check_crc() && valid_message()) report_message( rank, m); } if( stop) break; } if( logfile) { free( logfile); logfile = NULL; } } // // Parse the -p argument or the dialog box entry. It can be a comma- // separated list of polynomials, or an alias for one of the built-in // codes. // static void parse_polys( void) { char *p = poly_string, *q; // See if we have an alias for a built-in polynomial. If so, substitute. struct POLYLIST *list = polylist; while( list->alias && strcasecmp( list->alias, poly_string)) list++; if( list->alias) p = list->poly; // Now parse the comma-separate polynomial 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 != ',') bailout( "error parsing polynomial set\n"); p++; } } // // The decoder task. A separate thread so that the main thread continues // to service the user interface. // static DWORD WINAPI decoder_task( __attribute__ ((unused)) LPVOID arg) { show_status( "Initialising..."); int len; char temp[100]; len = GetWindowTextLength( hwnd_ncpu); if( len) { GetWindowText( hwnd_ncpu, temp, len + 1); ncpu = atoi( temp); } len = GetWindowTextLength( hwnd_start); if( len) { GetWindowText( hwnd_start, temp, len + 1); start_offset = atof( temp); } len = GetWindowTextLength( hwnd_freq); if( len) { GetWindowText( hwnd_freq, temp, len + 1); freq_offset = atof( temp); } len = GetWindowTextLength( hwnd_period); if( len) { GetWindowText( hwnd_period, temp, len + 1); sym_period = atof( temp); } len = GetWindowTextLength( hwnd_list); if( len) { GetWindowText( hwnd_list, temp, len + 1); max_rank = atoi( temp); } len = GetWindowTextLength( hwnd_chars); if( len) { GetWindowText( hwnd_chars, temp, len + 1); Nchars = atoi( temp); } len = GetWindowTextLength( hwnd_sigfile); if( len) GetWindowText( hwnd_sigfile, sigfile, len + 1); else sigfile[0] = 0; result.valid = FALSE; show_result(); if( ncpu < 1 || ncpu > MAXCPU) error_stop( "Invalid number of CPUs"); if( max_rank < 1 || max_rank > MAX_RANK) error_stop( "Invalid list length"); if( Nchars < 1 || Nchars > MAX_CHARS) error_stop( "Invalid message length"); if( !sigfile[0]) error_stop( "Signal filename not given"); if( !poly_string) error_stop( "Coding not set"); if( sym_period < 0.01) error_stop( "Invalid symbol period"); vt_cnt = CreateSemaphore( NULL, 0, MAXCPU, NULL); for( int i = 1; i < ncpu; i++) vt_sem[i] = CreateSemaphore( NULL, 0, MAXCPU, NULL); parse_polys(); if( QD < 2 || QD > 16) bailout( "invalid polynomials"); LB = Nchars * 6 + CRCLEN; if( LB > MAX_BITS) error_stop( "Too many bits (%d) to decode", LB); QS = 1 << (QK-1); // Number of encoder states; Nsymbols = QD * (LB + QK - 1); symamps = (complex double *) realloc( symamps, sizeof( complex double) * Nsymbols); sigbits = (double *) realloc( sigbits, sizeof( double) * Nsymbols); rotamps = (complex double *) realloc( rotamps, sizeof( complex double) * Nsymbols); if( !symamps || !sigbits || !rotamps) error_stop( "Not enough memory"); IW = ceil( sqrt( Nsymbols)); // Set up the interleave stride if( stop) return 0; Tstart = time( NULL); // Reference for elapsed time reporting // // Begin the log file. // if( logfile) free( logfile); char *name_prefix = strdup( sigfile); char *p = strrchr( name_prefix, '.'); if( p) *p = 0; logfile = malloc( strlen(name_prefix) + 8 + 1); if( !logfile) error_stop( "Not enough memory (logfile)"); sprintf( logfile, "%s-log.txt", name_prefix); FILE *f = fopen( logfile, "w"); if( !f) { free( logfile); logfile = NULL; } else fclose( f); report( 0, "ebnaut-rx V%s", VERSION); // // Open the input file and read its header. open_wav() will not return // if it finds a problem. // report( 0, "input file %s", sigfile); open_wav(); report( 0, "sample rate %.6f per second", sample_rate); if( inf1_rf) report( 0, "rx frequency %s", inf1_rf); if( inf1_ut) report( 0, "file start time %s", inf1_ut); // // Run the decoder, until it has finished or been signalled to stop. // run_decoder(); int elapsed = (int)(time(NULL) - Tstart); if( stop) return 0; show_status( "Finished: elapsed %d seconds", elapsed); report( 0, "finished elapsed %d seconds", elapsed); show_run(); if( autoexit) { report( 0, "auto exit"); exit( 0); } return 0; } static int runstop = FALSE; static void show_run( void) { if( stop) return; runstop = FALSE; SetWindowText( hwnd_runstop, "Run"); enable_setup( TRUE); } static void show_stop( void) { if( stop) return; runstop = TRUE; SetWindowText( hwnd_runstop, "Stop"); enable_setup( FALSE); } static HANDLE decoder_thread = NULL; static void decoder_start( void) { show_stop(); decoder_thread = CreateThread( NULL, 0, (LPTHREAD_START_ROUTINE) decoder_task, 0, 0, 0); } static void decoder_stop( void) { stop = TRUE; WaitForSingleObject( decoder_thread, INFINITE); stop = FALSE; show_status( "Stopped"); show_run(); } static void select_wav( void) { OPENFILENAME ofn; TCHAR szFile[MAX_PATH]; ZeroMemory(&ofn, sizeof(ofn)); ofn.lStructSize = sizeof(ofn); ofn.lpstrFile = szFile; ofn.lpstrFile[0] = '\0'; ofn.hwndOwner = hwnd; ofn.nMaxFile = sizeof(szFile); ofn.lpstrFilter = TEXT( "All files(*.*)\0*.*\0"); ofn.nFilterIndex = 1; ofn.lpstrInitialDir = NULL; ofn.lpstrFileTitle = NULL; ofn.Flags = OFN_PATHMUSTEXIST | OFN_FILEMUSTEXIST; if( GetOpenFileName( &ofn)) { strcpy( sigfile, szFile); SetWindowText( hwnd_sigfile, sigfile); } } /////////////////////////////////////////////////////////////////////////////// // Window/Display // /////////////////////////////////////////////////////////////////////////////// // // IDs for the various display widgets. // #define ID_PERIOD 8 #define ID_LIST 10 #define ID_RUN 11 #define ID_SIGFILE 20 #define ID_NCPU 21 #define ID_START 22 #define ID_CHARS 23 #define ID_FREQ 24 #define ID_PSTEP 25 #define ID_BROWSE 26 #define ID_CODESEL 27 #define ID_CRC 28 // // One-time initialisation of the main window and all the controls. // static void setup_display( void) { char temp[200]; int Y = 10; // // Code selection. // CreateWindow( "button", "Select Coding", WS_CHILD | WS_VISIBLE | BS_GROUPBOX, 10, Y, 425, 65, hwnd, (HMENU) 1, NULL, NULL); hwnd_codesel = CreateWindow( "combobox", NULL, WS_CHILD | WS_VISIBLE | WS_VSCROLL | CBS_DROPDOWNLIST, 20, Y+20, 130, 80, hwnd, (HMENU) ID_CODESEL, NULL, NULL); for( struct POLYLIST *p = polylist; p->alias; p++) SendMessage( hwnd_codesel, CB_ADDSTRING, 0, (LPARAM) p->alias); for( struct POLYLIST *p = polylist; p->alias; p++) if( !strcmp( p->alias, poly_string)) SendMessage( hwnd_codesel, CB_SETCURSEL, (WPARAM) (p - polylist), (LPARAM) 0); // // CRC length. // hwnd_crc = CreateWindow( "combobox", NULL, WS_CHILD | WS_VISIBLE | WS_VSCROLL | CBS_DROPDOWNLIST, 160, Y+20, 80, 80, hwnd, (HMENU) ID_CRC, NULL, NULL); for( int i = CRCBASE; i <= CRCMAX; i++) { sprintf( temp, "CRC %d", i); SendMessage( hwnd_crc, CB_ADDSTRING, 0, (LPARAM) temp); } SendMessage( hwnd_crc, CB_SETCURSEL, (WPARAM) (CRCLEN - CRCBASE), (LPARAM) 0); // // Symbol period. // CreateWindow( "STATIC", "Symbol period:", WS_CHILD | WS_VISIBLE | SS_LEFT, 260, Y + 23, 100, 20, hwnd, (HMENU) 1, NULL, NULL); hwnd_period = CreateWindow( "Edit", NULL, WS_CHILD | WS_VISIBLE | WS_BORDER, 370, Y + 20, 40, 20, hwnd, (HMENU) ID_PERIOD, NULL, NULL); sprintf( temp, "%.1f", sym_period); SetWindowText( hwnd_period, temp); // // Decoder Settings. // Y += 80; CreateWindow( "button", "Decoder Settings", WS_CHILD | WS_VISIBLE | BS_GROUPBOX, 10, Y, 425, 150, hwnd, (HMENU) 1, NULL, NULL); CreateWindow( "STATIC", "File:", WS_CHILD | WS_VISIBLE | SS_LEFT, 20, Y+23, 40, 20, hwnd, (HMENU) 1, NULL, NULL); hwnd_sigfile = CreateWindow( "Edit", NULL, WS_CHILD | WS_VISIBLE | WS_BORDER, 54, Y+20, 300, 20, hwnd, (HMENU) ID_SIGFILE, NULL, NULL); SetWindowText( hwnd_sigfile, sigfile); hwnd_browse = CreateWindow( "button", "Browse", WS_VISIBLE | WS_CHILD , 364, Y+17, 60, 23, hwnd, (HMENU) ID_BROWSE, NULL, NULL); CreateWindow( "STATIC", "Message length:", WS_CHILD | WS_VISIBLE | SS_LEFT, 20, Y+53, 150, 20, hwnd, (HMENU) 1, NULL, NULL); hwnd_chars = CreateWindow( "Edit", NULL, WS_CHILD | WS_VISIBLE | WS_BORDER, 140, Y+50, 45, 20, hwnd, (HMENU) ID_CHARS, NULL, NULL); sprintf( temp, "%d", Nchars); SetWindowText( hwnd_chars, temp); CreateWindow( "STATIC", "Start offset:", WS_CHILD | WS_VISIBLE | SS_LEFT, 195, Y+53, 100, 20, hwnd, (HMENU) 1, NULL, NULL); hwnd_start = CreateWindow( "Edit", NULL, WS_CHILD | WS_VISIBLE | WS_BORDER, 275, Y+50, 100, 20, hwnd, (HMENU) ID_START, NULL, NULL); sprintf( temp, "%.4f", start_offset); SetWindowText( hwnd_start, temp); CreateWindow( "STATIC", "Freq offset:", WS_CHILD | WS_VISIBLE | SS_LEFT, 195, Y+83, 100, 20, hwnd, (HMENU) 1, NULL, NULL); hwnd_freq = CreateWindow( "Edit", NULL, WS_CHILD | WS_VISIBLE | WS_BORDER, 275, Y+80, 100, 20, hwnd, (HMENU) ID_FREQ, NULL, NULL); sprintf( temp, "%.6f", freq_offset); SetWindowText( hwnd_freq, temp); CreateWindow( "STATIC", "List length:", WS_CHILD | WS_VISIBLE | SS_LEFT, 20, Y+83, 100, 20, hwnd, (HMENU) 1, NULL, NULL); hwnd_list = CreateWindow( "Edit", NULL, WS_CHILD | WS_VISIBLE | WS_BORDER, 100, Y+80, 70, 20, hwnd, (HMENU) ID_LIST, NULL, NULL); sprintf( temp, "%d", max_rank); SetWindowText( hwnd_list, temp); CreateWindow( "STATIC", "CPUs:", WS_CHILD | WS_VISIBLE | SS_LEFT, 20, Y+113, 60, 20, hwnd, (HMENU) 1, NULL, NULL); hwnd_ncpu = CreateWindow( "Edit", NULL, WS_CHILD | WS_VISIBLE | WS_BORDER, 70, Y+110, 25, 20, hwnd, (HMENU) ID_NCPU, NULL, NULL); sprintf( temp, "%d", ncpu); SetWindowText( hwnd_ncpu, temp); CreateWindow( "STATIC", "Phase step:", WS_CHILD | WS_VISIBLE | SS_LEFT, 195, Y+113, 100, 20, hwnd, (HMENU) 1, NULL, NULL); hwnd_pstep = CreateWindow( "combobox", NULL, WS_CHILD | WS_VISIBLE | WS_VSCROLL | CBS_DROPDOWNLIST, 275, Y+110, 130, 80, hwnd, (HMENU) ID_PSTEP, NULL, NULL); SendMessage( hwnd_pstep, CB_ADDSTRING, 0, (LPARAM) "30 degrees"); SendMessage( hwnd_pstep, CB_ADDSTRING, 0, (LPARAM) "15 degrees"); if( PSTEP == 30) SendMessage( hwnd_pstep, CB_SETCURSEL, (WPARAM) 0, (LPARAM) 0); else if( PSTEP == 15) SendMessage( hwnd_pstep, CB_SETCURSEL, (WPARAM) 1, (LPARAM) 0); // // Signal file information box. // Y += 165; CreateWindow( "button", "Signal File", WS_CHILD | WS_VISIBLE | BS_GROUPBOX, 10, Y, 425, 70, hwnd, (HMENU) 1, NULL, NULL); hwnd_signal = CreateWindow( "STATIC", "", WS_CHILD | WS_VISIBLE | SS_LEFT, 20, Y+20, 390, 40, hwnd, (HMENU) 1, NULL, NULL); // // Decoder status box. // Y += 80; CreateWindow( "button", "Decoder Status", WS_CHILD | WS_VISIBLE | BS_GROUPBOX, 10, Y, 425, 70, hwnd, (HMENU) 1, NULL, NULL); hwnd_status = CreateWindow( "STATIC", "", WS_CHILD | WS_VISIBLE | SS_LEFT, 20, Y+20, 390, 40, hwnd, (HMENU) 1, NULL, NULL); // // Decoder output box. // Y += 80; CreateWindow( "button", "Decoder Output", WS_CHILD | WS_VISIBLE | BS_GROUPBOX, 10, Y, 425, 150, hwnd, (HMENU) 1, NULL, NULL); hwnd_result = CreateWindow( "STATIC", "", WS_CHILD | WS_VISIBLE | SS_LEFT | SS_NOPREFIX, 20, Y+20, 390, 120, hwnd, (HMENU) 1, NULL, NULL); Y += 160; hwnd_runstop = CreateWindow( "button", "", WS_VISIBLE | WS_CHILD , 20, Y, 80, 25, hwnd, (HMENU) ID_RUN, NULL, NULL); show_run(); } // // Main window event handler. // LRESULT CALLBACK main_event(HWND h, UINT msg, WPARAM wParam, LPARAM lParam) { switch( msg) { case WM_CREATE: hwnd = h; setup_display(); if( autostart) decoder_start(); return 0; case WM_COMMAND: switch( LOWORD(wParam)) { case ID_CRC: if( HIWORD(wParam) == CBN_SELCHANGE) { int i = SendMessage( hwnd_crc, CB_GETCURSEL, 0, 0); CRCLEN = CRCBASE + i; } return 0; case ID_CODESEL: if( HIWORD(wParam) == CBN_SELCHANGE) { int i = SendMessage( hwnd_codesel, CB_GETCURSEL, 0, 0); poly_string = polylist[i].alias; } return 0; case ID_PSTEP: if( HIWORD(wParam) == CBN_SELCHANGE) { int i = SendMessage( hwnd_pstep, CB_GETCURSEL, 0, 0); PSTEP = i == 1 ? 15 : 30; } return 0; } if( HIWORD(wParam) == BN_CLICKED) switch( LOWORD(wParam)) { case ID_RUN: runstop ? decoder_stop() : decoder_start(); return 0; case ID_BROWSE: select_wav(); return 0; } break; case WM_DESTROY: PostQuitMessage(0); return 0; } return DefWindowProc( h, msg, wParam, lParam); } /////////////////////////////////////////////////////////////////////////////// // Command Line // /////////////////////////////////////////////////////////////////////////////// // // Split the command line into space separated arguments. // #define MAX_ARGS 100 static void split_cmdline( char *cmdline, int *argc, char *argv[]) { argv[0] = "ebnaut-rx.exe"; *argc = 1; char *p = cmdline; while( *argc < MAX_ARGS && *p) { if( isspace( *p)) { p++; continue; } if( *p == '"') { p++; argv[(*argc)++] = p; while( *p && *p != '"') p++; } else { argv[(*argc)++] = p; while( *p && !isspace( *p)) p++; } if( *p) *p++ = 0; } } static void usage( void) { fprintf( stderr, "options:\n" "\n" " -N nchars Message length in characters (no default)\n" " -f filename Input file name\n" " -c ncpu Number of CPU cores to use (default 1)\n" " -L size Viterbi list size (default 20000)\n" " -p poly Polynomial set, or alias\n" " -k crclen CRC length in bits, 8 to 32 (default 16)\n" " -F freq Frequency offset, Hz\n" " -T secs Start time offset, seconds\n" " -S secs Symbol period (default 1.0)\n" "\n" " -PS15 Phase search step 15 degrees\n" " -PS30 Phase search step 30 degrees (default)\n" " -PU Search uniform phase only (default full search)\n" "\n" " -as Start automatically (default start button)\n" " -ax Exit automatically when run complete\n" ); exit( 0); } int WINAPI WinMain( HINSTANCE hInstance, __attribute__ ((unused)) HINSTANCE hPrevInst, LPTSTR cmdline, int nShowCmd) { char *argv[MAX_ARGS]; int argc; split_cmdline( cmdline, &argc, argv); while( 1) { int c = getopt( argc, argv, "N:L:f:c:p:F:T:S:P:k:a:v?"); if( c == -1) break; if( c == 'N') Nchars = atoi( optarg); else if( c == 'L') max_rank = atoi( optarg); else if( c == 'f') strcpy( sigfile, optarg); else if( c == 'c') ncpu = atoi( optarg); else if( c == 'p') poly_string = strdup( optarg); else if( c == 'F') freq_offset = atof( optarg); else if( c == 'T') start_offset = atof( optarg); else if( c == 'a') { if( !strcmp( optarg, "s")) autostart = TRUE; else if( !strcmp( optarg, "x")) autoexit = TRUE; else if( !strcmp( optarg, "sx") || !strcmp( optarg, "xs")) { autostart = autoexit = TRUE; } else bailout( "unrecognised -a argument [%s]", optarg); } else if( c == 'S') sym_period = 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 bailout( "invalid -P argument [%s], optarg"); } else if( c == 'v') loglevel++; else if( c == 'k') CRCLEN = atoi( optarg); else if( c == '?') usage(); else bailout( "unrecognised option %c", c); } if( optind != argc) usage(); if( CRCLEN) if( CRCLEN < CRCBASE || CRCLEN > CRCMAX) bailout( "invalid CRC length"); WNDCLASS wc; memset( &wc, 0, sizeof( wc)); wc.lpszClassName = "Static Control"; wc.hInstance = hInstance; wc.hbrBackground = GetSysColorBrush( COLOR_3DFACE); wc.lpfnWndProc = main_event; wc.hCursor = LoadCursor( 0, IDC_ARROW); RegisterClass( &wc); char title[100]; sprintf( title, "EbNaut Decoder V%s", VERSION); HWND h = CreateWindow( wc.lpszClassName, title, WS_OVERLAPPEDWINDOW | WS_VISIBLE, CW_USEDEFAULT, CW_USEDEFAULT, 455, 640, 0, 0, hInstance, 0); ShowWindow( h, nShowCmd); UpdateWindow( h); MSG msg; while( GetMessage( &msg, NULL, 0, 0)) { TranslateMessage( &msg); DispatchMessage( &msg); } return (int) msg.wParam; }