diff --git a/auto_rx/autorx/decode.py b/auto_rx/autorx/decode.py index a3c75f7..f66e6b7 100644 --- a/auto_rx/autorx/decode.py +++ b/auto_rx/autorx/decode.py @@ -244,26 +244,14 @@ class SondeDecoder(object): elif self.sonde_type == "DFM": # DFM06/DFM09 Sondes. - - # We need to handle inversion of DFM sondes in a bit of an odd way. - # Using our current receive chain (rtl_fm), rs_detect will detect: - # DFM06's as non-inverted ('DFM') - # DFM09's as inverted ('-DFM') - # HOWEVER, dfm09dm_ecc makes the assumption that the incoming signal is a DFM09, and - # inverts by default. - # So, to be able to support DFM06s, we need to flip the invert flag. - self.inverted = not self.inverted - - if self.inverted: - _invert_flag = "-i" - else: - _invert_flag = "" + # As of 2019-02-10, dfm09ecc auto-detects if the signal is inverted, + # so we don't need to specify an invert flag. # Note: Have removed a 'highpass 20' filter from the sox line, will need to re-evaluate if adding that is useful in the future. decode_cmd = "%s %s-p %d -d %s %s-M fm -F9 -s 15k -f %d 2>/dev/null |" % (self.sdr_fm, bias_option, int(self.ppm), str(self.device_idx), gain_param, self.sonde_freq) decode_cmd += "sox -t raw -r 15k -e s -b 16 -c 1 - -r 48000 -b 8 -t wav - highpass 20 lowpass 2000 2>/dev/null |" # DFM decoder - decode_cmd += "./dfm09ecc -vv --ecc --json %s 2>/dev/null" % _invert_flag + decode_cmd += "./dfm09ecc -vv --ecc --json --auto 2>/dev/null" elif self.sonde_type == "M10": # M10 Sondes diff --git a/demod/dfm09dm_dft.c b/demod/dfm09dm_dft.c index b8fb24a..fbbe064 100644 --- a/demod/dfm09dm_dft.c +++ b/demod/dfm09dm_dft.c @@ -40,22 +40,33 @@ typedef struct { double dir; double horiV; double vertV; float meas24[5]; float status[2]; + float _frmcnt; + char sonde_id[16]; } gpx_t; gpx_t gpx; +typedef struct { + int ec; + float ts; +} pcksts_t; + +pcksts_t pck[9]; + char dat_str[9][13+1]; -// Buffer to store sonde ID -char sonde_id[] = "DFMxx-xxxxxxxx"; +// JSON Buffer to store sonde ID +char json_sonde_id[] = "DFMxx-xxxxxxxxyy"; int option_verbose = 0, // ausfuehrliche Anzeige option_raw = 0, // rohe Frames option_inv = 0, // invertiert Signal + option_auto = 0, + option_dist = 0, // continuous pcks 0..8 option_ecc = 0, option_ptu = 0, option_ths = 0, - option_json = 0, // JSON blob output (for auto_rx) + option_json = 0, // JSON blob output (for auto_rx) wavloaded = 0; int wav_channel = 0; // audio channel: left @@ -154,9 +165,13 @@ int check(ui8_t code[8]) { int hamming(ui8_t *ham, int L, ui8_t *sym) { int i, j; - int ret = 0; // L = 7, 13 + int ecc = 0, ret = 0; // L = 7, 13 for (i = 0; i < L; i++) { // L * 2 nibble (data+parity) - if (option_ecc) ret |= check(ham+B*i); + if (option_ecc) { + ecc = check(ham+B*i); + if (ecc > 0) ret |= (1<= 0 && fr_id <= 8) { + if (fr_id >= 0 && fr_id <= 8) { for (i = 0; i < 13; i++) { nib = bits2val(dat_bits+4*i, 4); dat_str[fr_id][i] = nib2chr(nib); } dat_str[fr_id][13] = '\0'; + + pck[fr_id].ts = gpx._frmcnt; // time_stamp,frame_count,... + if (option_ecc) { + pck[fr_id].ec = ec; // option_ecc laesst -1 garnicht durch + if (ec > 0) { + ui8_t ecn = 0; + for (i = 0; i < 15; i++) { + if ( (ec>>i)&1 ) ecn++; + } + pck[fr_id].ec = ecn; + if ((option_dist || option_json) && ecn > 4) pck[fr_id].ec = -2; // threshold: #errors > 4 + } + } } if (fr_id == 0) { - start = 1; + start = 0x1000; frnr = bits2val(dat_bits+24, 8); gpx.frnr = frnr; } @@ -230,10 +257,10 @@ int dat_out(ui8_t *dat_bits) { if (fr_id == 5) { } - if (fr_id == 6) { + if (fr_id == 6) { // sat data } - if (fr_id == 7) { + if (fr_id == 7) { // sat data } if (fr_id == 8) { @@ -333,7 +360,7 @@ float get_Temp2(float *meas) { // meas[0..4] R = (f-f1)/g; // meas[0,3,4] > 0 ? if (R > 0) T = 1/(1/T0 + 1/B0 * log(R/R0)); - if (option_ptu && ptu_out && option_verbose == 2) { + if (option_ptu && ptu_out && option_verbose == 3) { printf(" (Rso: %.1f , Rb: %.1f)", Rs_o/1e3, Rb/1e3); } @@ -380,9 +407,8 @@ float get_Temp4(float *meas) { // meas[0..4] } -#define RSNbit 0x0100 // radiosonde DFM-06,DFM-09 -#define PSNbit 0x0200 // pilotsonde PS-15 -int conf_out(ui8_t *conf_bits) { +#define SNbit 0x0100 +int conf_out(ui8_t *conf_bits, int ec) { int conf_id; int ret = 0; int val, hl; @@ -400,18 +426,20 @@ int conf_out(ui8_t *conf_bits) { // // SN/ID immer im letzten Kanal? if ((gpx.sonde_typ & 0xF) < 7 && conf_id == 6) { - SN6 = bits2val(conf_bits+4, 4*6); // DFM-06: Kanal 6 - if ( SN6 == gpx.SN6 && SN6 != 0) { // nur Nibble-Werte 0..9 - gpx.sonde_typ = RSNbit | 6; + SN6 = bits2val(conf_bits+4, 4*6); // DFM-06: Kanal 6 + if (SN6 == gpx.SN6 && SN6 != 0) { // nur Nibble-Werte 0..9 + gpx.sonde_typ = SNbit | 6; ptu_out = 6; ret = 6; + sprintf(gpx.sonde_id, "ID06:%6X", gpx.SN6); + sprintf(json_sonde_id, "DFM06-%6X", gpx.SN6); } else { gpx.sonde_typ = 0; } gpx.SN6 = SN6; } - if (conf_id == 0xA) { // 0xACxxxxy + if (conf_id == 0xA) { // 0xACxxxxy , DFM-09 val = bits2val(conf_bits+8, 4*5); hl = (val & 1); // val&0xF 0,1? chA[hl] = (val >> 4) & 0xFFFF; @@ -419,10 +447,12 @@ int conf_out(ui8_t *conf_bits) { if (chAbit == 3) { // DFM-09: Kanal A SN = (chA[0] << 16) | chA[1]; if ( SN == SN_A ) { - gpx.sonde_typ = RSNbit | 0xA; + gpx.sonde_typ = SNbit | 0xA; gpx.SN = SN; ptu_out = 9; ret = 9; + sprintf(gpx.sonde_id, "ID09:%6u", gpx.SN); + sprintf(json_sonde_id, "DFM09-%6u", gpx.SN); } else { gpx.sonde_typ = 0; @@ -431,7 +461,7 @@ int conf_out(ui8_t *conf_bits) { chAbit = 0; } } - if (conf_id == 0xC) { // 0xCCxxxxy + if (conf_id == 0xC) { // 0xCCxxxxy , DFM-17? val = bits2val(conf_bits+8, 4*5); hl = (val & 1); chC[hl] = (val >> 4) & 0xFFFF; @@ -439,10 +469,12 @@ int conf_out(ui8_t *conf_bits) { if (chCbit == 3) { // DFM-17? Kanal C SN = (chC[0] << 16) | chC[1]; if ( SN == SN_C ) { - gpx.sonde_typ = RSNbit | 0xC; + gpx.sonde_typ = SNbit | 0xC; gpx.SN = SN; ptu_out = 9; ret = 17; + sprintf(gpx.sonde_id, "ID-%1X:%6u", gpx.sonde_typ & 0xF, gpx.SN); + sprintf(json_sonde_id, "DFM17-%6u", gpx.SN); } else { gpx.sonde_typ = 0; @@ -451,7 +483,7 @@ int conf_out(ui8_t *conf_bits) { chCbit = 0; } } - if (conf_id == 0xD) { // 0xDCxxxxy + if (conf_id == 0xD) { // 0xDCxxxxy , DFM-17? val = bits2val(conf_bits+8, 4*5); hl = (val & 1); chD[hl] = (val >> 4) & 0xFFFF; @@ -459,10 +491,12 @@ int conf_out(ui8_t *conf_bits) { if (chDbit == 3) { // DFM-17? Kanal D SN = (chD[0] << 16) | chD[1]; if ( SN == SN_D ) { - gpx.sonde_typ = RSNbit | 0xD; + gpx.sonde_typ = SNbit | 0xD; gpx.SN = SN; ptu_out = 9; ret = 18; + sprintf(gpx.sonde_id, "ID-%1X:%6u", gpx.sonde_typ & 0xF, gpx.SN); + sprintf(json_sonde_id, "DFM17-%6u", gpx.SN); } else { gpx.sonde_typ = 0; @@ -471,7 +505,7 @@ int conf_out(ui8_t *conf_bits) { chDbit = 0; } } - if (conf_id == 0x7) { // 0x70xxxxy + if (conf_id == 0x7) { // 0x70xxxxy , pilotsonde PS-15? val = bits2val(conf_bits+8, 4*5); hl = (val & 1); ch7[hl] = (val >> 4) & 0xFFFF; @@ -479,10 +513,12 @@ int conf_out(ui8_t *conf_bits) { if (ch7bit == 3) { // PS-15: Kanal 7 SN = (ch7[0] << 16) | ch7[1]; if ( SN == SN_7 ) { - gpx.sonde_typ = PSNbit | 0x7; + gpx.sonde_typ = SNbit | 0x7; gpx.SN = SN; ptu_out = 0; ret = 15; + sprintf(gpx.sonde_id, "ID15:%6u", gpx.SN); + sprintf(json_sonde_id, "DFM15-%6u", gpx.SN); } else { gpx.sonde_typ = 0; @@ -492,6 +528,7 @@ int conf_out(ui8_t *conf_bits) { } } + if (conf_id >= 0 && conf_id <= 4) { val = bits2val(conf_bits+4, 4*6); gpx.meas24[conf_id] = fl24(val); @@ -519,95 +556,110 @@ int conf_out(ui8_t *conf_bits) { } void print_gpx() { - int i, j; + int i, j; + int contgps = 0; + int output = 0; + int jsonout = 0; - if (start) { + output |= start; - if (option_raw == 2) { - for (i = 0; i < 9; i++) { - printf(" %s", dat_str[i]); - } - for (i = 0; i < 9; i++) { - for (j = 0; j < 13; j++) dat_str[i][j] = ' '; - } - } - else { - //if (option_auto && option_verbose) printf("[%c] ", option_inv?'-':'+'); - printf("[%3d] ", gpx.frnr); - printf("%4d-%02d-%02d ", gpx.jahr, gpx.monat, gpx.tag); - printf("%02d:%02d:%04.1f ", gpx.std, gpx.min, gpx.sek); - printf(" "); - printf("lat: %.6f ", gpx.lat); - printf("lon: %.6f ", gpx.lon); - printf("alt: %.1f ", gpx.alt); - printf(" vH: %5.2f ", gpx.horiV); - printf(" D: %5.1f ", gpx.dir); - printf(" vV: %5.2f ", gpx.vertV); - if (option_ptu && ptu_out) { - float t = get_Temp(gpx.meas24); - if (t > -270.0) printf(" T=%.1fC ", t); - if (option_verbose == 2) { - float t2 = get_Temp2(gpx.meas24); - float t4 = get_Temp4(gpx.meas24); - if (t2 > -270.0) printf(" T2=%.1fC ", t2); - if (t4 > -270.0) printf(" T4=%.1fC ", t4); - printf(" f0: %.2f ", gpx.meas24[0]); - printf(" f3: %.2f ", gpx.meas24[3]); - printf(" f4: %.2f ", gpx.meas24[4]); - } - } - if (option_verbose == 2 && (gpx.sonde_typ & 0xF) == 0xA) { - printf(" U: %.2fV ", gpx.status[0]); - printf(" Ti: %.1fK ", gpx.status[1]); - } - if ( option_verbose ) - { - if (gpx.sonde_typ & RSNbit) - { - if ((gpx.sonde_typ & 0xF) == 6) { // DFM-06 - printf(" (ID6:%06X) ", gpx.SN6); - sprintf(sonde_id, "DFM06-%06X", gpx.SN6); - } - if ((gpx.sonde_typ & 0xF) == 0xA) { // DFM-09 - printf(" (ID9:%06u) ", gpx.SN); - sprintf(sonde_id, "DFM09-%06u", gpx.SN); - } - if ((gpx.sonde_typ & 0xF) == 0xC || // DFM-17? - (gpx.sonde_typ & 0xF) == 0xD ) { - printf(" (ID-%1X:%06u) ", gpx.sonde_typ & 0xF, gpx.SN); - sprintf(sonde_id, "DFM17-%06u", gpx.SN); - } - gpx.sonde_typ ^= RSNbit; - } - if (gpx.sonde_typ & PSNbit) { - if ((gpx.sonde_typ & 0xF) == 0x7) { // PS-15? - printf(" (ID15:%06u) ", gpx.SN); - sprintf(sonde_id, "DFM15-%06u", gpx.SN); - } - gpx.sonde_typ ^= PSNbit; - } - } - } - printf("\n"); + for (i = 0; i < 9/*8*/; i++) { // trigger: pck8 + if ( !( (option_dist || option_json) && pck[i].ec < 0) ) + { + if (pck[8].ts - pck[i].ts < 6.0) { output |= (1<= 2) printf("[%c] ", option_inv?'-':'+'); + printf("[%3d] ", gpx.frnr); + printf("%4d-%02d-%02d ", gpx.jahr, gpx.monat, gpx.tag); + printf("%02d:%02d:%04.1f ", gpx.std, gpx.min, gpx.sek); + if (option_verbose >= 2 && option_ecc) printf("[%1X,%1X,%1X] ", pck[0].ec&0xF, pck[8].ec&0xF, pck[1].ec&0xF); + printf(" "); + printf(" lat: %.5f ", gpx.lat); if (option_verbose >= 2 && option_ecc) printf("[%1X] ", pck[2].ec&0xF); + printf(" lon: %.5f ", gpx.lon); if (option_verbose >= 2 && option_ecc) printf("[%1X] ", pck[3].ec&0xF); + printf(" alt: %.1f ", gpx.alt); if (option_verbose >= 2 && option_ecc) printf("[%1X] ", pck[4].ec&0xF); + printf(" vH: %5.2f ", gpx.horiV); + printf(" D: %5.1f ", gpx.dir); + printf(" vV: %5.2f ", gpx.vertV); + if (option_ptu && ptu_out) { + float t = get_Temp(gpx.meas24); + if (t > -270.0) printf(" T=%.1fC ", t); + if (option_verbose == 3) { + float t2 = get_Temp2(gpx.meas24); + float t4 = get_Temp4(gpx.meas24); + if (t2 > -270.0) printf(" T2=%.1fC ", t2); + if (t4 > -270.0) printf(" T4=%.1fC ", t4); + printf(" f0: %.2f ", gpx.meas24[0]); + printf(" f3: %.2f ", gpx.meas24[3]); + printf(" f4: %.2f ", gpx.meas24[4]); + } + } + if (option_verbose == 3 && (gpx.sonde_typ & 0xF) == 0xA) { + printf(" U: %.2fV ", gpx.status[0]); + printf(" Ti: %.1fK ", gpx.status[1]); + } + if (option_verbose) + { + if (gpx.sonde_typ & SNbit) { + printf(" (%s) ", gpx.sonde_id); + gpx.sonde_typ ^= SNbit; + } + } + } + printf("\n"); + + if (option_json && jsonout) + { + // Print JSON blob // valid sonde_ID? + printf("{ \"frame\": %d, \"id\": \"%s\", \"datetime\": \"%04d-%02d-%02dT%02d:%02d:%06.3fZ\", \"lat\": %.5f, \"lon\": %.5f, \"alt\": %.5f, \"vel_h\": %.5f, \"heading\": %.5f, \"vel_v\": %.5f", + gpx.frnr, json_sonde_id, gpx.jahr, gpx.monat, gpx.tag, gpx.std, gpx.min, gpx.sek, gpx.lat, gpx.lon, gpx.alt, gpx.horiV, gpx.dir, gpx.vertV); + if (ptu_out) { // get temperature + float t = get_Temp(gpx.meas24); // ecc-valid temperature? + if (t > -270.0) printf(", \"temp\": %.1f", t); + } + printf(" }\n"); + printf("\n"); + } + + } + + for (i = 0; i < 9; i++) pck[i].ec = -1; } -int print_frame() { +int print_frame(float frmcnt) { int i; int nib = 0; int frid = -1; int ret0, ret1, ret2; int ret = 0; + gpx._frmcnt = frmcnt; deinterleave(frame_bits+CONF, 7, hamming_conf); deinterleave(frame_bits+DAT1, 13, hamming_dat1); @@ -655,24 +707,24 @@ int print_frame() { else if (option_ecc) { if (ret0 == 0 || ret0 > 0) { - conf_out(block_conf); + conf_out(block_conf, ret0); } if (ret1 == 0 || ret1 > 0) { - frid = dat_out(block_dat1); + frid = dat_out(block_dat1, ret1); if (frid == 8) print_gpx(); } if (ret2 == 0 || ret2 > 0) { - frid = dat_out(block_dat2); + frid = dat_out(block_dat2, ret2); if (frid == 8) print_gpx(); } } else { - conf_out(block_conf); - frid = dat_out(block_dat1); + conf_out(block_conf, ret0); + frid = dat_out(block_dat1, ret1); if (frid == 8) print_gpx(); - frid = dat_out(block_dat2); + frid = dat_out(block_dat2, ret2); if (frid == 8) print_gpx(); } @@ -699,12 +751,14 @@ int main(int argc, char **argv) { int headerlen = 0; int frm = 0, nfrms = 8; // nfrms=1,2,4,8 - int k,K; + int k, K; float mv; unsigned int mv_pos, mv0_pos; int mp = 0; - float thres = 0.6; + float frm_cnt = 0.0; + + float thres = 0.65; int bitofs = 0; int symlen = 2; @@ -732,7 +786,8 @@ int main(int argc, char **argv) { else if ( (strcmp(*argv, "-v") == 0) || (strcmp(*argv, "--verbose") == 0) ) { option_verbose = 1; } - else if ( (strcmp(*argv, "-vv") == 0) ) { option_verbose = 2; } + else if ( (strcmp(*argv, "-vv" ) == 0) ) { option_verbose = 2; } + else if ( (strcmp(*argv, "-vvv") == 0) ) { option_verbose = 3; } else if ( (strcmp(*argv, "-r") == 0) || (strcmp(*argv, "--raw") == 0) ) { option_raw = 1; } @@ -743,8 +798,10 @@ int main(int argc, char **argv) { option_inv = 0x1; } else if ( (strcmp(*argv, "--ecc") == 0) ) { option_ecc = 1; } - else if ( (strcmp(*argv, "--json") == 0) ) { option_json = 1; } - else if ( (strcmp(*argv, "--ptu") == 0) ) { option_ptu = 1; ptu_out = 1; } + else if ( (strcmp(*argv, "--ptu") == 0) ) { option_ptu = 1; ptu_out = 1; } + else if ( (strcmp(*argv, "--auto") == 0) ) { option_auto = 1; } + else if ( (strcmp(*argv, "--dist") == 0) ) { option_dist = 1; option_ecc = 1; } + else if ( (strcmp(*argv, "--json") == 0) ) { option_json = 1; option_ecc = 1; } else if ( (strcmp(*argv, "--ch2") == 0) ) { wav_channel = 1; } // right channel (default: 0=left) else if ( (strcmp(*argv, "--ths") == 0) ) { ++argv; @@ -777,6 +834,9 @@ int main(int argc, char **argv) { } + for (k = 0; k < 9; k++) pck[k].ec = -1; // init ecc-status + + symlen = 2; headerlen = strlen(rawheader); bitofs = 2; // +1 .. +2 @@ -786,16 +846,15 @@ int main(int argc, char **argv) { return -1; }; - k = 0; - mv = -1; mv_pos = 0; + mv = 0; mv_pos = 0; while ( f32buf_sample(fp, option_inv, 1) != EOF ) { k += 1; if (k >= K-4) { mv0_pos = mv_pos; - mp = getCorrDFT(0, K, 0, &mv, &mv_pos); + mp = getCorrDFT(option_auto, K, 0, &mv, &mv_pos); k = 0; } else { @@ -803,7 +862,8 @@ int main(int argc, char **argv) { continue; } - if (mv > thres && mp > 0) { + if ( mp > 0 && (mv > thres || mv < -thres)) { + if (mv_pos > mv0_pos) { header_found = 0; @@ -816,7 +876,16 @@ int main(int argc, char **argv) { herr1 = 1; } } - if (herrs <= 1) header_found = 1; // herrs <= 1 bitfehler in header + if (herrs <= 1) { + header_found = 1; // herrs <= 1 bitfehler in header + if (mv < 0) header_found = -header_found; + } + + if (header_found < 0) { + // read_sbit(option_inv) buffer reset? + if (option_auto) option_inv ^= 0x1; + else header_found = 0; + } if (header_found) { @@ -824,8 +893,11 @@ int main(int argc, char **argv) { pos = headerlen; pos /= 2; + //if (fabs(mv) > 0.85) nfrms = 8; else nfrms = 4; // test OK/KO/NO count + frm = 0; while ( frm < nfrms ) { // nfrms=1,2,4,8 + frm_cnt = mv_pos/(spb*2.0*BITFRAME_LEN) + frm; while ( pos < BITFRAME_LEN ) { header_found = !(frm==nfrms-1 && pos>=BITFRAME_LEN-10); bitQ = read_sbit(fp, symlen, &bit, option_inv, bitofs, bitpos==0, !header_found); // symlen=2, return: zeroX/bit @@ -835,7 +907,7 @@ int main(int argc, char **argv) { bitpos += 1; } frame_bits[pos] = '\0'; - ret = print_frame(); + ret = print_frame(frm_cnt); if (pos < BITFRAME_LEN) break; pos = 0; frm += 1; diff --git a/demod/rs41dm_dft.c b/demod/rs41dm_dft.c index 7f42e5f..8982b2b 100644 --- a/demod/rs41dm_dft.c +++ b/demod/rs41dm_dft.c @@ -647,9 +647,8 @@ int get_GPS1() { err = check_CRC(pos_GPS1, pck_GPS1); if (err) gpx.crc |= crc_GPS1; - //err = 0; - err |= get_GPSweek(); - err |= get_GPStime(); + err |= get_GPSweek(); // no plausibility-check + err |= get_GPStime(); // no plausibility-check return err; } @@ -726,7 +725,7 @@ int get_GPSkoord() { gpx.lat = lat; gpx.lon = lon; gpx.alt = alt; - if ((alt < -1000) || (alt > 80000)) return -3; + if ((alt < -1000) || (alt > 80000)) return -3; // plausibility-check: altitude, if ecef=(0,0,0) // ECEF-Velocities @@ -765,7 +764,7 @@ int get_GPS3() { err = check_CRC(pos_GPS3, pck_GPS3); if (err) gpx.crc |= crc_GPS3; - err |= get_GPSkoord(); + err |= get_GPSkoord(); // plausibility-check: altitude, if ecef=(0,0,0) return err; } @@ -1038,64 +1037,62 @@ int print_position(int ec) { //fprintf(stdout, " (%.1f %.1f %.1f) ", gpx.vN, gpx.vE, gpx.vU); fprintf(stdout," vH: %4.1f D: %5.1f° vV: %3.1f ", gpx.vH, gpx.vD, gpx.vU); } - - if (option_json){ - // Print JSON output required by auto_rx. - if (!err1 && !err2 && !err3){ - if (option_ptu && !err0 && gpx.T > -273.0) { - printf("\n{ \"frame\": %d, \"id\": \"%s\", \"datetime\": \"%04d-%02d-%02dT%02d:%02d:%06.3fZ\", \"lat\": %.5f, \"lon\": %.5f, \"alt\": %.5f, \"vel_h\": %.5f, \"heading\": %.5f, \"vel_v\": %.5f, \"temp\":%.1f }\n", gpx.frnr, gpx.id, gpx.jahr, gpx.monat, gpx.tag, gpx.std, gpx.min, gpx.sek, gpx.lat, gpx.lon, gpx.alt, gpx.vH, gpx.vD, gpx.vU, gpx.T ); - } else { - printf("\n{ \"frame\": %d, \"id\": \"%s\", \"datetime\": \"%04d-%02d-%02dT%02d:%02d:%06.3fZ\", \"lat\": %.5f, \"lon\": %.5f, \"alt\": %.5f, \"vel_h\": %.5f, \"heading\": %.5f, \"vel_v\": %.5f }\n", gpx.frnr, gpx.id, gpx.jahr, gpx.monat, gpx.tag, gpx.std, gpx.min, gpx.sek, gpx.lat, gpx.lon, gpx.alt, gpx.vH, gpx.vD, gpx.vU ); - } - } - } } if (option_ptu && !err0) { if (gpx.T > -273.0) printf(" T=%.1fC ", gpx.T); } - //if (output) - { - if (option_crc) { - fprintf(stdout, " # "); - if (option_ecc && ec >= 0 && (gpx.crc & 0x1F) != 0) { - int pos, blk, len, crc; // unexpected blocks - int flen = NDATA_LEN; - if (frametype() < 0) flen += XDATA_LEN; - pos = pos_FRAME; - while (pos < flen-1) { - blk = frame[pos]; // 0x80XX: encrypted block - len = frame[pos+1]; // 0x76XX: 00-padding block - crc = check_CRC(pos, blk<<8); - fprintf(stdout, " %02X%02X", frame[pos], frame[pos+1]); - fprintf(stdout, "[%d]", crc&1); - pos = pos+2+len+2; - } + if (option_crc) { // show CRC-checks (and RS-check) + fprintf(stdout, " # "); + if (option_ecc && ec >= 0 && (gpx.crc & 0x1F) != 0) { + int pos, blk, len, crc; // unexpected blocks + int flen = NDATA_LEN; + if (frametype() < 0) flen += XDATA_LEN; + pos = pos_FRAME; + while (pos < flen-1) { + blk = frame[pos]; // 0x80XX: encrypted block + len = frame[pos+1]; // 0x76XX: 00-padding block + crc = check_CRC(pos, blk<<8); + fprintf(stdout, " %02X%02X", frame[pos], frame[pos+1]); + fprintf(stdout, "[%d]", crc&1); + pos = pos+2+len+2; } - else { - fprintf(stdout, "["); - for (i=0; i<5; i++) fprintf(stdout, "%d", (gpx.crc>>i)&1); - fprintf(stdout, "]"); - } - if (option_ecc == 2) { - if (ec > 0) fprintf(stdout, " (%d)", ec); - if (ec < 0) { - if (ec == -1) fprintf(stdout, " (-+)"); - else if (ec == -2) fprintf(stdout, " (+-)"); - else /*ec == -3*/ fprintf(stdout, " (--)"); - } + } + else { + fprintf(stdout, "["); + for (i=0; i<5; i++) fprintf(stdout, "%d", (gpx.crc>>i)&1); + fprintf(stdout, "]"); + } + if (option_ecc == 2) { + if (ec > 0) fprintf(stdout, " (%d)", ec); + if (ec < 0) { + if (ec == -1) fprintf(stdout, " (-+)"); + else if (ec == -2) fprintf(stdout, " (+-)"); + else /*ec == -3*/ fprintf(stdout, " (--)"); } } } get_Calconf(output); - //if (output) - { - if (option_verbose > 1) get_Aux(); - fprintf(stdout, "\n"); // fflush(stdout); + if (option_verbose > 1) get_Aux(); + + fprintf(stdout, "\n"); // fflush(stdout); + + + if (option_json) { + // Print JSON output required by auto_rx. + if (!err && !err1 && !err3) { // frame-nb/id && gps-time && gps-position (crc-)ok; 3 CRCs, RS not needed + if (option_ptu && !err0 && gpx.T > -273.0) { + printf("{ \"frame\": %d, \"id\": \"%s\", \"datetime\": \"%04d-%02d-%02dT%02d:%02d:%06.3fZ\", \"lat\": %.5f, \"lon\": %.5f, \"alt\": %.5f, \"vel_h\": %.5f, \"heading\": %.5f, \"vel_v\": %.5f, \"temp\":%.1f }\n", gpx.frnr, gpx.id, gpx.jahr, gpx.monat, gpx.tag, gpx.std, gpx.min, gpx.sek, gpx.lat, gpx.lon, gpx.alt, gpx.vH, gpx.vD, gpx.vU, gpx.T ); + } else { + printf("{ \"frame\": %d, \"id\": \"%s\", \"datetime\": \"%04d-%02d-%02dT%02d:%02d:%06.3fZ\", \"lat\": %.5f, \"lon\": %.5f, \"alt\": %.5f, \"vel_h\": %.5f, \"heading\": %.5f, \"vel_v\": %.5f }\n", gpx.frnr, gpx.id, gpx.jahr, gpx.monat, gpx.tag, gpx.std, gpx.min, gpx.sek, gpx.lat, gpx.lon, gpx.alt, gpx.vH, gpx.vD, gpx.vU ); + } + printf("\n"); + } } + } err |= err1 | err3; @@ -1227,7 +1224,7 @@ int main(int argc, char *argv[]) { else if (strcmp(*argv, "--std2") == 0) { frmlen = 518; } // NDATA_LEN+XDATA_LEN else if (strcmp(*argv, "--sat") == 0) { option_sat = 1; } else if (strcmp(*argv, "--ptu") == 0) { option_ptu = 1; } - else if (strcmp(*argv, "--json") == 0) { option_json = 1; } + else if (strcmp(*argv, "--json") == 0) { option_json = 1; option_ecc = 2; option_crc = 1; } else if (strcmp(*argv, "--ch2") == 0) { wav_channel = 1; } // right channel (default: 0=left) else if (strcmp(*argv, "--ths") == 0) { ++argv; diff --git a/ecc/bch_ecc.c b/ecc/bch_ecc.c index 2a96788..d892296 100644 --- a/ecc/bch_ecc.c +++ b/ecc/bch_ecc.c @@ -104,6 +104,9 @@ static RS_t RS256 = { 255, 12, 24, 231, 0, 1, 1, {0}}; static RS_t RS256ccsds = { 255, 16, 32, 223, 112, 11, 116, {0}}; static RS_t BCH64 = { 63, 2, 12, 51, 1, 1, 1, {0}}; +// static RS_t RS16_0 = { 15, 3, 6, 9, 0, 1, 1, {0}}; +static RS_t RS16ccsds = { 15, 2, 4, 11, 6, 1, 1, {0}}; + static GF_t GF; static RS_t RS; @@ -810,15 +813,40 @@ int rs_init_BCH64() { return check_gen; } +int rs_init_RS15ccsds() { + int i, check_gen; + ui8_t Xalp[MAX_DEG+1]; + + GF = GF16RS; + check_gen = GF_genTab( GF, exp_a, log_a); + + //RS = RS16_0; // N=15, t=3, b=0, p=1 + RS = RS16ccsds; // N=15, t=2, b=6, p=1 + for (i = 0; i <= MAX_DEG; i++) RS.g[i] = 0; + for (i = 0; i <= MAX_DEG; i++) Xalp[i] = 0; + + // g(X)=(X-alpha^b)...(X-alpha^(b+2t-1)) + RS.g[0] = 0x01; + Xalp[1] = 0x01; // X + for (i = 0; i < 2*RS.t; i++) { + Xalp[0] = exp_a[(RS.b+i) % (GF.ord-1)]; // Xalp[0..1]: X - alpha^(b+i) + poly_mul(RS.g, Xalp, RS.g); + } + + return check_gen; +} + int rs_encode(ui8_t cw[]) { int j; - ui8_t parity[MAX_DEG+1], + ui8_t __cw[MAX_DEG+1], + parity[MAX_DEG+1], d[MAX_DEG+1]; - for (j = 0; j < RS.R; j++) cw[j] = 0; for (j = 0; j <= MAX_DEG; j++) parity[j] = 0; - poly_divmod(cw, RS.g, d, parity); + for (j = 0; j <= MAX_DEG; j++) __cw[j] = 0; + for (j = RS.R; j < RS.N; j++) __cw[j] = cw[j]; + poly_divmod(__cw, RS.g, d, parity); //if (poly_deg(parity) >= RS.R) return -1; - for (j = 0; j <= poly_deg(parity); j++) cw[j] = parity[j]; + for (j = 0; j < RS.R; j++) cw[j] = parity[j]; return 0; } diff --git a/ecc/ecc-rs_gf16.c b/ecc/ecc-rs_gf16.c new file mode 100644 index 0000000..5134a5e --- /dev/null +++ b/ecc/ecc-rs_gf16.c @@ -0,0 +1,288 @@ + +/* + GF(2^4) - RS(15,11) + (no bit-packing, i.e. 1 nibble <-> 1 byte) + +a) + ka9q-fec: + gcc -c init_rs_char.c + gcc -c encode_rs_char.c + gcc -c decode_rs_char.c + + gcc -DKA9Q ecc-rs_gf16.c init_rs_char.o encode_rs_char.o decode_rs_char.o -o ecc-rs15ccsds + +b) + gcc ecc-rs_gf16.c -o ecc-rs15ccsds + +*/ + + +#include +#include + +#ifdef KA9Q + #include "fec.h" +#endif + + +typedef unsigned char ui8_t; +typedef unsigned int ui32_t; + +#include "bch_ecc.c" + + +#define BFSIZE 512 +#define rs_N 15 +#define rs_R 4 +#define rs_K (rs_N-rs_R) + +ui8_t data[BFSIZE]; + +ui8_t cw[rs_N+1]; // cw[MAX_DEG+1]; // fixed in encode() ... +ui8_t par[rs_N], msg[rs_N]; +int errs; +ui8_t err_pos[rs_R], err_val[rs_R]; + +#ifdef KA9Q + void *rs; + ui8_t codeword[rs_N]; + int errors, errpos[rs_R]; +#endif + + +int main(int argc, char *argv[]) { + int i, l; + + int len1, len2; + char *str1 = NULL, + *str2 = NULL; + + + rs_init_RS15ccsds(); // 0x13: X^4 + X + 1, t=2, b=6 + +#ifdef KA9Q + //rs = init_rs_char( 4, 0x13, 0, 1, 6, 0); // RS16_0 + rs = init_rs_char( 4, 0x13, 6, 1, 4, 0); // RS16ccsds +#endif + + if (argv[0] == NULL) return -1; + + for (i = 0; i < BFSIZE; i++) data[i] = 0; + for (i = 0; i < rs_N; i++) cw[i] = 0; +#ifdef KA9Q + for (i = 0; i < rs_N; i++) codeword[i] = 0; +#endif + + + str1 = argv[1]; + + len1 = strlen(str1); + if (len1 > BFSIZE) return -1; + + for (i = 0; i < len1; i++) { + l = sscanf(str1+i, "%1hhx", data+i); if (l < 1) {/*len1 = i;*/} + } + // 1byte/nibble + + + for (i = 0; i < rs_K; i++) cw[rs_R+i] = data[i]; // codeword[rs_N-1-i] = cw[i]; + + for (i = 0; i < rs_N; i++) msg[i] = 0; + for (i = 0; i < rs_N; i++) par[i] = 0; + for (i = 0; i < rs_K; i++) msg[i] = data[rs_K-1-i]; + + +// +// GF(16) %02X -> %1X (1 nibble / 1 byte) +// printf("%1X", cw[i]); // dbg: printf("%02X", cw[i]); +// printf("%1X", cw[i]&0xF); // dbg: printf("%02X", cw[i]); +// + + + printf("\n"); + + if (argv[2]) { + + str2 = argv[2]; + len2 = strlen(str2); + if (len2 > rs_R) len2 = rs_R; + + for (i = 0; i < len2; i++) { + l = sscanf(str2+i, "%1hhx", par+i); if (l < 1) {/*len2 = i;*/} + } + while (i < rs_R) par[i++] = 0; + + for (i = 0; i < rs_N; i++) cw[i] = 0; + for (i = 0; i < rs_R; i++) cw[i] = par[i]; + for (i = 0; i < len1; i++) cw[rs_R+i] = msg[rs_K-1-i]; + + +#ifdef KA9Q + for (i = 0; i < rs_N; i++) codeword[rs_N-1-i] = cw[i]; +#endif + + printf("(received)\n"); + printf("msg: "); + for (i = rs_R; i < rs_N; i++) printf("%1X", cw[i]); // dbg: printf("%02X", cw[i]); + printf("\n"); + printf("par: "); + for (i = 0; i < rs_R; i++) printf("%1X", cw[i]); + printf("\n"); + printf("msg+par:\n"); + for (i = 0; i < rs_N; i++) printf("%1X", cw[i]); + printf("\n"); + + printf("\n"); + + + printf("cw\n"); + errs = rs_decode(cw, err_pos, err_val); + printf("errs: %d\n", errs); + if (errs > 0) { + printf("pos: "); + for (i = 0; i < errs; i++) printf(" %d", err_pos[i]); + printf("\n"); + } + for (i = 0; i < rs_N; i++) printf("%1X", cw[i]); printf("\n"); + + +#ifdef KA9Q + printf("\n"); + + printf("codeword\n"); + errors = decode_rs_char(rs, codeword, errpos, 0); + printf("errors: %d\n", errors); + if (errors > 0) { + printf("pos: "); + for (i = 0; i < errors; i++) printf(" %d", errpos[i]); + printf("\n"); + } + for (i = 0; i < rs_N; i++) printf("%1x", codeword[i]); printf("\n"); +#endif + + } else { + + printf("msg: "); + for (i = rs_R; i < rs_N; i++) printf("%1X", cw[i]); // dbg: printf("%02X", cw[i]); + printf("\n"); + printf("\n"); + + + printf("cw\n"); + rs_encode(cw); + printf("par: "); + for (i = 0; i < rs_R; i++) printf("%1X", cw[i]); printf("\n"); + printf("cw-enc:\n"); + for (i = 0; i < rs_N; i++) { + //if (i == rs_R) printf(" "); + printf("%1X", cw[i]); + } + printf("\n"); + + // check + errs = rs_decode(cw, err_pos, err_val); + if (errs) { + printf("errs: %d\n", errs); + printf("cw-dec:\n"); + for (i = 0; i < rs_N; i++) { + //if (i == rs_R) printf(" "); + printf("%1X", cw[i]); + } + printf("\n"); + } + + +#ifdef KA9Q + printf("\n"); + + printf("codeword\n"); + printf("message: "); + for (i = 0; i < rs_K; i++) printf("%1x", msg[i]); printf("\n"); + encode_rs_char(rs, msg, par); + printf("parity : "); + for (i = 0; i < rs_R; i++) printf("%1x", par[i]); printf("\n"); + + for (i = 0; i < rs_K; i++) codeword[i] = msg[i]; + for (i = 0; i < rs_R; i++) codeword[rs_K+i] = par[i]; + printf("codeword:\n"); + for (i = 0; i < rs_N; i++) printf("%1x", codeword[i]); printf("\n"); +#endif + + } + + + printf("\n"); + + return 0; +} + + +/* + +RS(15,11): + codeword length 15 + message length 11 + parity length 4 + +./ecc-rs15 msg [par] + msg: 11 nibbles + par: 4 nibbles + +ecc-rs15 input/output: nibbles +cw[]: 1 byte / 1 nibble + + +examples: + +$ ./ecc-rs15ccsds 00000000001 + +msg: 00000000001 + +cw +par: 8281 +cw-enc: +828100000000001 + +codeword +message: 10000000000 +parity : 1828 +codeword: +100000000001828 + +$ ./ecc-rs15ccsds 00000000001 8281 + +(received) +msg: 00000000001 +par: 8281 +msg+par: +828100000000001 + +cw +errs: 0 +828100000000001 + +codeword +errors: 0 +100000000001828 + +$ ./ecc-rs15ccsds 00000000001 8283 + +(received) +msg: 00000000001 +par: 8283 +msg+par: +828300000000001 + +cw +errs: 1 +pos: 3 +828100000000001 + +codeword +errors: 1 +pos: 11 +100000000001828 + +*/ + + diff --git a/scan/dft_detect.c b/scan/dft_detect.c new file mode 100644 index 0000000..5c681ba --- /dev/null +++ b/scan/dft_detect.c @@ -0,0 +1,881 @@ + +#include +#include +#include +#include +#include + + +typedef unsigned char ui8_t; +typedef unsigned short ui16_t; +typedef unsigned int ui32_t; +typedef short i16_t; +typedef int i32_t; + + +static int option_verbose = 0, // ausfuehrliche Anzeige + option_inv = 0, // invertiert Signal + option_dc = 0, + option_silent = 0, + wavloaded = 0; +static int wav_channel = 0; // audio channel: left + + +//int dfm_bps = 2500; +static char dfm_header[] = "01100101011001101010010110101010"; + +//int vai_bps = 4800; +static char rs41_header[] = "00001000011011010101001110001000" + "01000100011010010100100000011111"; +static char rs92_header[] = //"10100110011001101001" + //"10100110011001101001" + "10100110011001101001" + "10100110011001101001" + "1010011001100110100110101010100110101001"; + +//int lms_bps = 4800; +static char lms6_header[] = "0101011000001000""0001110010010111" + "0001101010100111""0011110100111110"; + +//int m10_bps = 9600; +static char m10_header[] = "10011001100110010100110010011001"; +// frame byte[0..1]: byte[0]=framelen-1, byte[1]=type(8F=M2K2,9F=M10,AF=M10+) +// M2K2 : 64 8F : 0110010010001111 +// M10 : 64 9F : 0110010010011111 (framelen 0x64+1) +// M10-aux: 76 9F : 0111011010011111 (framelen 0x76+1) +// M10+ : 64 AF : 0110010010101111 (w/ gtop-GPS) + +// imet_9600 / 1200 Hz; +static char imet_preamble[] = "11110000111100001111000011110000" + "11110000111100001111000011110000" + "11110000111100001111000011110000" + "11110000111100001111000011110000"; // 1200 Hz preamble + + +//int imet1ab_bps = 9600; // 1200 bits/sec +static char imet1ab_header[] = "11110000111100001111000011110000" + // "11110000""10101100110010101100101010101100" + "11110000""10101100110010101100101010101100"; + + +// 11110000:1 , 001100110:0 // 11/4=2.1818.. +static char imet1rs_header[] = + "0000""1111""0000""1111""0000""1111" // preamble + "0000""1111"; + +// imet1rs/imet4 1200Hz preamble , lead_out , 8N1 byte: lead-in 8bits lead-out , ... +// 1:1200Hz/0:2200Hz tones, bit-duration 1/1200 sec, phase ... +// bits: 1111111111111111111 10 10000000 10 ..; + + +// C34/C50: 2400 baud, 1:2900Hz/0:4800Hz +static char c34_preheader[] = +"01010101010101010101010101010101"; // 2900 Hz tone +// dft, dB-max(1000Hz..5000Hz) = 2900Hz ? + +typedef struct { + int bps; // header: here bps means baudrate ... + int hLen; + int N; + char *header; + float BT; + float spb; + float thres; + float complex *Fm; + char *type; + unsigned char tn; +} rsheader_t; + +#define Nrs 9 +#define idxAB 7 +#define idxRS 8 +static rsheader_t rs_hdr[Nrs] = { + { 2500, 0, 0, dfm_header, 1.0, 0.0, 0.65, NULL, "DFM", 2}, + { 4800, 0, 0, rs41_header, 0.5, 0.0, 0.70, NULL, "RS41", 3}, + { 4800, 0, 0, rs92_header, 0.5, 0.0, 0.70, NULL, "RS92", 4}, + { 4800, 0, 0, lms6_header, 1.0, 0.0, 0.70, NULL, "LMS6", 8}, + { 9600, 0, 0, m10_header, 1.0, 0.0, 0.76, NULL, "M10", 5}, + { 5800, 0, 0, c34_preheader, 1.5, 0.0, 0.80, NULL, "C34C50", 9}, // C34/C50 2900 Hz tone + { 9600, 0, 0, imet_preamble, 0.5, 0.0, 0.80, NULL, "IMET", 6}, // IMET1AB=7, IMET1RS=8 + { 9600, 0, 0, imet1ab_header, 1.0, 0.0, 0.80, NULL, "IMET1AB", 6}, + { 9600, 0, 0, imet1rs_header, 0.5, 0.0, 0.80, NULL, "IMET1RS", 7} // IMET4 +}; + + +/* +// m10-false-positive: +// m10-preamble similar to rs41-preamble, parts of rs92/imet1ab; diffs: +// - iq: - modulation-index rs41 < rs92 < m10, +// - power level / frame < 1s, noise +// - fm: - frame duration <-> noise (variance/standard deviation) +// - pulse-shaping +// m10: 00110011 at 9600 bps +// rs41: 0 1 0 1 at 4800 bps +// - m10 top-carrier, fm-mean/average +// - m10-header ..110(1)0110011()011.. bit shuffle +// - m10 frame byte[1]=type(M2K2,M10,M10+) +*/ + +/* +// rs92 +// imet1ab-false-positive +// ... +*/ + + +static int sample_rate = 0, bits_sample = 0, channels = 0; +static int wav_ch = 0; // 0: links bzw. mono; 1: rechts + +static unsigned int sample_in, sample_out, delay; + +static int M; // N + +static float *bufs = NULL; + +static char *rawbits = NULL; + +static int Nvar = 0; // < M +static double xsum = 0; +static float *xs = NULL; +/* +static double xsum=0, qsum=0; +static float *xs = NULL, + *qs = NULL; +*/ + +static float dc_ofs = 0.0; +static float dc = 0.0; + +/* ------------------------------------------------------------------------------------ */ + + +static int LOG2N, N_DFT; + +static float complex *ew; + +static float complex *X, *Z, *cx; +static float *xn; +static float *db; + +static void dft_raw(float complex *Z) { + int s, l, l2, i, j, k; + float complex w1, w2, T; + + j = 1; + for (i = 1; i < N_DFT; i++) { + if (i < j) { + T = Z[j-1]; + Z[j-1] = Z[i-1]; + Z[i-1] = T; + } + k = N_DFT/2; + while (k < j) { + j = j - k; + k = k/2; + } + j = j + k; + } + + for (s = 0; s < LOG2N; s++) { + l2 = 1 << s; + l = l2 << 1; + w1 = (float complex)1.0; + w2 = ew[s]; // cexp(-I*M_PI/(float)l2) + for (j = 1; j <= l2; j++) { + for (i = j; i <= N_DFT; i += l) { + k = i + l2; + T = Z[k-1] * w1; + Z[k-1] = Z[i-1] - T; + Z[i-1] = Z[i-1] + T; + } + w1 = w1 * w2; + } + } +} + +static void dft(float *x, float complex *Z) { + int i; + for (i = 0; i < N_DFT; i++) Z[i] = (float complex)x[i]; + dft_raw(Z); +} + +static void Nidft(float complex *Z, float complex *z) { + int i; + for (i = 0; i < N_DFT; i++) z[i] = conj(Z[i]); + dft_raw(z); + // idft(): + // for (i = 0; i < N_DFT; i++) z[i] = conj(z[i])/(float)N_DFT; // hier: z reell +} + +static float freq2bin(int f) { + return f * N_DFT / (float)sample_rate; +} + +static float bin2freq(int k) { + return sample_rate * k / (float)N_DFT; +} + +/* ------------------------------------------------------------------------------------ */ +/* +static float get_bufvar(int ofs) { + float mu = xs[(sample_out+M + ofs) % M]/Nvar; + float var = qs[(sample_out+M + ofs) % M]/Nvar - mu*mu; + return var; +} +*/ +static float get_bufmu(int ofs) { + float mu = xs[(sample_out+M + ofs) % M]/Nvar; + return mu; +} + + +static int getCorrDFT(int abs, int K, unsigned int pos, float *maxv, unsigned int *maxvpos, rsheader_t rshd) { + int i; + int mp = -1; + float mx = 0.0; + double xnorm = 1; + unsigned int mpos = 0; + + dc = 0.0; + + if (rshd.N + K > N_DFT/2 - 2) return -1; + if (sample_in < delay+rshd.N+K) return -2; + + if (pos == 0) pos = sample_out; + + + for (i = 0; i < rshd.N+K; i++) xn[i] = bufs[(pos+M -(rshd.N+K-1) + i) % M]; + while (i < N_DFT) xn[i++] = 0.0; + + dft(xn, X); + + dc = get_bufmu(pos-sample_out); //oder: dc = creal(X[0])/N_DFT; + + for (i = 0; i < N_DFT; i++) Z[i] = X[i]*rshd.Fm[i]; + + Nidft(Z, cx); + + + if (abs) { + for (i = rshd.N; i < rshd.N+K; i++) { + if (fabs(creal(cx[i])) > fabs(mx)) { // imag(cx)=0 + mx = creal(cx[i]); + mp = i; + } + } + } + else { + for (i = rshd.N; i < rshd.N+K; i++) { + if (creal(cx[i]) > mx) { // imag(cx)=0 + mx = creal(cx[i]); + mp = i; + } + } + } + if (mp == rshd.N || mp == rshd.N+K-1) return -4; // Randwert + + mpos = pos - ( rshd.N+K-1 - mp ); + + //xnorm = sqrt(qs[(mpos + 2*M) % M]); + xnorm = 0.0; + for (i = 0; i < rshd.N; i++) xnorm += bufs[(mpos-i + M) % M]*bufs[(mpos-i + M) % M]; + xnorm = sqrt(xnorm); + + mx /= xnorm*N_DFT; + + *maxv = mx; + *maxvpos = mpos; + + + return mp; +} + +/* ------------------------------------------------------------------------------------ */ + +static int findstr(char *buff, char *str, int pos) { + int i; + for (i = 0; i < 4; i++) { + if (buff[(pos+i)%4] != str[i]) break; + } + return i; +} + +static int read_wav_header(FILE *fp, int wav_channel) { + char txt[4+1] = "\0\0\0\0"; + unsigned char dat[4]; + int byte, p=0; + + if (fread(txt, 1, 4, fp) < 4) return -1; + if (strncmp(txt, "RIFF", 4)) return -1; + if (fread(txt, 1, 4, fp) < 4) return -1; + // pos_WAVE = 8L + if (fread(txt, 1, 4, fp) < 4) return -1; + if (strncmp(txt, "WAVE", 4)) return -1; + // pos_fmt = 12L + for ( ; ; ) { + if ( (byte=fgetc(fp)) == EOF ) return -1; + txt[p % 4] = byte; + p++; if (p==4) p=0; + if (findstr(txt, "fmt ", p) == 4) break; + } + if (fread(dat, 1, 4, fp) < 4) return -1; + if (fread(dat, 1, 2, fp) < 2) return -1; + + if (fread(dat, 1, 2, fp) < 2) return -1; + channels = dat[0] + (dat[1] << 8); + + if (fread(dat, 1, 4, fp) < 4) return -1; + memcpy(&sample_rate, dat, 4); //sample_rate = dat[0]|(dat[1]<<8)|(dat[2]<<16)|(dat[3]<<24); + + if (fread(dat, 1, 4, fp) < 4) return -1; + if (fread(dat, 1, 2, fp) < 2) return -1; + //byte = dat[0] + (dat[1] << 8); + + if (fread(dat, 1, 2, fp) < 2) return -1; + bits_sample = dat[0] + (dat[1] << 8); + + // pos_dat = 36L + info + for ( ; ; ) { + if ( (byte=fgetc(fp)) == EOF ) return -1; + txt[p % 4] = byte; + p++; if (p==4) p=0; + if (findstr(txt, "data", p) == 4) break; + } + if (fread(dat, 1, 4, fp) < 4) return -1; + + + fprintf(stderr, "sample_rate: %d\n", sample_rate); + fprintf(stderr, "bits : %d\n", bits_sample); + fprintf(stderr, "channels : %d\n", channels); + + if (wav_channel >= 0 && wav_channel < channels) wav_ch = wav_channel; + else wav_ch = 0; + fprintf(stderr, "channel-In : %d\n", wav_ch+1); + + if ((bits_sample != 8) && (bits_sample != 16)) return -1; + + return 0; +} + +static int f32read_sample(FILE *fp, float *s) { + int i; + short b = 0; + + for (i = 0; i < channels; i++) { + + if (fread( &b, bits_sample/8, 1, fp) != 1) return EOF; + + if (i == wav_ch) { // i = 0: links bzw. mono + //if (bits_sample == 8) sint = b-128; // 8bit: 00..FF, centerpoint 0x80=128 + //if (bits_sample == 16) sint = (short)b; + + if (bits_sample == 8) { b -= 128; } + *s = b/128.0; + if (bits_sample == 16) { *s /= 256.0; } + } + } + + return 0; +} + + +static int f32buf_sample(FILE *fp, int inv, int cm) { + float s = 0.0; + float xneu, xalt; + + + if (f32read_sample(fp, &s) == EOF) return EOF; + + if (inv) s = -s; + bufs[sample_in % M] = s - dc_ofs; + + xneu = bufs[(sample_in ) % M]; + xalt = bufs[(sample_in+M - Nvar) % M]; + xsum += xneu - xalt; // + xneu - xalt + xs[sample_in % M] = xsum; +/* + qsum += (xneu - xalt)*(xneu + xalt); // + xneu*xneu - xalt*xalt + qs[sample_in % M] = qsum; +*/ + + if (0 && cm) { + // direct correlation + } + + + sample_out = sample_in - delay; + + sample_in += 1; + + return 0; +} + +static int read_bufbit(int symlen, char *bits, unsigned int mvp, int reset, float spb) { +// symlen==2: manchester2 0->10,1->01->1: 2.bit + + static unsigned int rcount; + static float rbitgrenze; + + double sum = 0.0; + + if (reset) { + rcount = 0; + rbitgrenze = 0; + } + + + rbitgrenze += spb; + do { + sum += bufs[(rcount + mvp + M) % M]; + rcount++; + } while (rcount < rbitgrenze); // n < spb + + if (symlen == 2) { + rbitgrenze += spb; + do { + sum -= bufs[(rcount + mvp + M) % M]; + rcount++; + } while (rcount < rbitgrenze); // n < spb + } + + + if (symlen != 2) { + if (sum >= 0) *bits = '1'; + else *bits = '0'; + } + else { + if (sum >= 0) strncpy(bits, "10", 2); + else strncpy(bits, "01", 2); + } + + return 0; +} + +static int headcmp(int symlen, char *hdr, int len, unsigned int mvp, int inv, int option_dc, float spb) { + int errs = 0; + int pos; + int step = 1; + char sign = 0; + + if (symlen != 1) step = 2; + if (inv) sign=1; + + for (pos = 0; pos < len; pos += step) { + read_bufbit(symlen, rawbits+pos, mvp+1-(int)(len*spb), pos==0, spb); + } + rawbits[pos] = '\0'; + + while (len > 0) { + if ((rawbits[len-1]^sign) != hdr[len-1]) errs += 1; + len--; + } + + if (option_dc && errs < 3) { + dc_ofs += dc; + } + + return errs; +} + +/* -------------------------------------------------------------------------- */ + + +#define SQRT2 1.4142135624 // sqrt(2) +// sigma = sqrt(log(2)) / (2*PI*BT): +//#define SIGMA 0.2650103635 // BT=0.5: 0.2650103635 , BT=0.3: 0.4416839392 + +// Gaussian FM-pulse +static double Q(double x) { + return 0.5 - 0.5*erf(x/SQRT2); +} +static double pulse(double t, double sigma) { + return Q((t-0.5)/sigma) - Q((t+0.5)/sigma); +} + + +static double norm2_match(float *match, int n) { + int i; + double x, y = 0.0; + for (i = 0; i < n; i++) { + x = match[i]; + y += x*x; + } + return y; +} + +static int init_buffers() { + + int i, j, pos; + double t; + double b0, b1, b2, b; + float normMatch; + + int K, NN; + int n, k; + float *match = NULL; + float *m = NULL; + + double BT = 0.5; + double sigma = sqrt(log(2)) / (2*M_PI*BT); + + char *bits = NULL; + float spb = 0.0; + + int hLen = 0; + + for (j = 0; j < Nrs; j++) { + rs_hdr[j].spb = sample_rate/(float)rs_hdr[j].bps; + rs_hdr[j].hLen = strlen(rs_hdr[j].header); + rs_hdr[j].N = rs_hdr[j].hLen * rs_hdr[j].spb + 0.5; + if (rs_hdr[j].hLen > hLen) hLen = rs_hdr[j].hLen; + } + + NN = hLen * sample_rate/2500.0 + 0.5; // max(hLen*spb) + + M = 3*NN; + //if (samples_per_bit < 6) M = 6*N; + + delay = NN/16; + sample_in = 0; + + K = M-NN - delay; // N+K < M + + LOG2N = 2 + (int)(log(NN+K)/log(2)); + N_DFT = 1 << LOG2N; + + while (NN + K > N_DFT/2 - 2) { + LOG2N += 1; + N_DFT <<= 1; + } + + Nvar = NN; // wenn Nvar fuer xnorm, dann Nvar=rshd.N + + rawbits = (char *)calloc( hLen+1, sizeof(char)); if (rawbits == NULL) return -100; + bufs = (float *)calloc( M+1, sizeof(float)); if (bufs == NULL) return -100; + xs = (float *)calloc( M+1, sizeof(float)); if (xs == NULL) return -100; +/* + qs = (float *)calloc( M+1, sizeof(float)); if (qs == NULL) return -100; +*/ + + xn = calloc(N_DFT+1, sizeof(float)); if (xn == NULL) return -1; + db = calloc(N_DFT+1, sizeof(float)); if (db == NULL) return -1; + + ew = calloc(LOG2N+1, sizeof(complex float)); if (ew == NULL) return -1; + X = calloc(N_DFT+1, sizeof(complex float)); if (X == NULL) return -1; + Z = calloc(N_DFT+1, sizeof(complex float)); if (Z == NULL) return -1; + cx = calloc(N_DFT+1, sizeof(complex float)); if (cx == NULL) return -1; + + for (n = 0; n < LOG2N; n++) { + k = 1 << n; + ew[n] = cexp(-I*M_PI/(float)k); + } + + match = (float *)calloc( NN+1, sizeof(float)); if (match == NULL) return -1; + m = (float *)calloc(N_DFT+1, sizeof(float)); if (m == NULL) return -1; + + + for (j = 0; j < Nrs-1; j++) + { + rs_hdr[j].Fm = (float complex *)calloc(N_DFT+1, sizeof(complex float)); if (rs_hdr[j].Fm == NULL) return -1; + bits = rs_hdr[j].header; + spb = rs_hdr[j].spb; + sigma = sqrt(log(2)) / (2*M_PI*rs_hdr[j].BT); + + for (i = 0; i < rs_hdr[j].N; i++) { + + pos = i/spb; + t = (i - pos*spb)/spb - 0.5; + + b1 = ((bits[pos] & 0x1) - 0.5)*2.0; + b = b1*pulse(t, sigma); + + if (pos > 0) { + b0 = ((bits[pos-1] & 0x1) - 0.5)*2.0; + b += b0*pulse(t+1, sigma); + } + + if (pos < hLen) { + b2 = ((bits[pos+1] & 0x1) - 0.5)*2.0; + b += b2*pulse(t-1, sigma); + } + + match[i] = b; + } + + normMatch = sqrt(norm2_match(match, rs_hdr[j].N)); + for (i = 0; i < rs_hdr[j].N; i++) { + match[i] /= normMatch; + } + + for (i = 0; i < rs_hdr[j].N; i++) m[rs_hdr[j].N-1 - i] = match[i]; + while (i < N_DFT) m[i++] = 0.0; + dft(m, rs_hdr[j].Fm); + + } + + + free(match); match = NULL; + free(m); m = NULL; + + return K; +} + +static int free_buffers() { + int j; + + if (bufs) { free(bufs); bufs = NULL; } + if (xs) { free(xs); xs = NULL; } +/* + if (qs) { free(qs); qs = NULL; } +*/ + if (rawbits) { free(rawbits); rawbits = NULL; } + + if (xn) { free(xn); xn = NULL; } + if (db) { free(xn); xn = NULL; } + if (ew) { free(ew); ew = NULL; } + if (X) { free(X); X = NULL; } + if (Z) { free(Z); Z = NULL; } + if (cx) { free(cx); cx = NULL; } + + for (j = 0; j < Nrs-1; j++) { + if (rs_hdr[j].Fm) { free(rs_hdr[j].Fm); rs_hdr[j].Fm = NULL; } + } + + return 0; +} + +/* ------------------------------------------------------------------------------------ */ + + +int main(int argc, char **argv) { + + FILE *fp = NULL; + char *fpname = NULL; + + int j; + int k, K; + float mv[Nrs]; + unsigned int mv_pos[Nrs], mv0_pos[Nrs]; + int mp[Nrs]; + + int header_found = 0; + int herrs; + float thres = 0.76; + + int j_max; + float mv_max; + +#ifdef CYGWIN + _setmode(fileno(stdin), _O_BINARY); // _setmode(_fileno(stdin), _O_BINARY); +#endif + setbuf(stdout, NULL); + + fpname = argv[0]; + ++argv; + while ((*argv) && (!wavloaded)) { + if ( (strcmp(*argv, "-h") == 0) || (strcmp(*argv, "--help") == 0) ) { + fprintf(stderr, "%s [options] audio.wav\n", fpname); + fprintf(stderr, " options:\n"); + //fprintf(stderr, " -v, --verbose\n"); + return 0; + } + else if ( (strcmp(*argv, "-v") == 0) || (strcmp(*argv, "--verbose") == 0) ) { + option_verbose = 1; + } + else if ( (strcmp(*argv, "--dc") == 0) ) { + option_dc = 1; + } + else if ( (strcmp(*argv, "-s") == 0) || (strcmp(*argv, "--silent") == 0) ) { + option_silent = 1; + } + else if ( (strcmp(*argv, "--ch2") == 0) ) { wav_channel = 1; } // right channel (default: 0=left) + else if ( (strcmp(*argv, "--ths") == 0) ) { + ++argv; + if (*argv) { + thres = atof(*argv); + for (j = 0; j < Nrs; j++) rs_hdr[j].thres = thres; + } + else return -1; + } + else { + fp = fopen(*argv, "rb"); + if (fp == NULL) { + fprintf(stderr, "%s konnte nicht geoeffnet werden\n", *argv); + return -1; + } + wavloaded = 1; + } + ++argv; + } + if (!wavloaded) fp = stdin; + + + j = read_wav_header(fp, wav_channel); + if ( j < 0 ) { + fclose(fp); + fprintf(stderr, "error: wav header\n"); + return -1; + } + + + K = init_buffers(); + if ( K < 0 ) { + fprintf(stderr, "error: init buffers\n"); + return -1; + }; + + for (j = 0; j < Nrs; j++) { + mv[j] = 0; + mv_pos[j] = 0; + mp[j] = 0; + } + + k = 0; + + while ( f32buf_sample(fp, option_inv, 1) != EOF ) { + + k += 1; + + if (k >= K-4) { + for (j = 0; j < Nrs-2; j++) { + mv0_pos[j] = mv_pos[j]; + mp[j] = getCorrDFT(-1, K, 0, mv+j, mv_pos+j, rs_hdr[j]); + } + k = 0; + } + else { + //for (j = 0; j < Nrs; j++) mv[j] = 0.0; + continue; + } + + header_found = 0; + for (j = 0; j < Nrs-2; j++) + { + if (mp[j] > 0 && (mv[j] > rs_hdr[j].thres || mv[j] < -rs_hdr[j].thres)) { + if (mv_pos[j] > mv0_pos[j]) { + + herrs = headcmp(1, rs_hdr[j].header, rs_hdr[j].hLen, mv_pos[j], mv[j]<0, option_dc, rs_hdr[j].spb); + if (herrs < 2) { // max 1 bitfehler in header + + if ( strncmp(rs_hdr[j].type, "IMET", 4) == 0 ) + { + int n, m; + int D = N_DFT/2 - 3; + float df; + float pow2200, pow2400; + int bin2200, bin2400; + + for (n = 0; n < N_DFT; n++) { + xn[n] = 0.0; + db[n] = 0.0; + } + + n = 0; + while (n < sample_rate) { // 1 sec + + if (f32buf_sample(fp, option_inv, 1) == EOF) break;//goto ende; + + xn[n % D] = bufs[sample_out % M]; + n++; + + if (n % D == 0) { + dft(xn, X); + for (m = 0; m < N_DFT; m++) db[m] += cabs(X[m]); + } + } + + df = bin2freq(1); + m = 50.0/df; + if (m < 1) m = 1; + if (freq2bin(2500) > N_DFT/2) goto ende; + + bin2200 = freq2bin(2200); + pow2200 = 0.0; + for (n = 0; n < m; n++) pow2200 += db[ bin2200 - m/4 + n ]; + + bin2400 = freq2bin(2400); + pow2400 = 0.0; + for (n = 0; n < m; n++) pow2400 += db[ bin2400 - m/4 + n ]; + + + mv[j] = fabs(mv[j]); + + if (pow2200 > pow2400) { // IMET1RS + mv[idxRS] = mv[j]; + mv[j] = 0; // IMET1 -> IMET1RS + mv_pos[idxRS] = mv_pos[j]; + j = idxRS; + header_found = 1; + } + else { // IMET1AB + mv[j] = 0; + j = idxAB; + mv_pos[j] = sample_out; + n = 0; + + // detect header/polarity + k = 0; + while ( n < 4*sample_rate && f32buf_sample(fp, option_inv, 1) != EOF ) { + + n += 1; + k += 1; + + if (k >= K-4) { + mv0_pos[j] = mv_pos[j]; + mp[j] = getCorrDFT(-1, K, 0, mv+j, mv_pos+j, rs_hdr[j]); + k = 0; + } + else { + //mv[j] = 0.0; + continue; + } + + if (mp[j] > 0 && (mv[j] > rs_hdr[j].thres || mv[j] < -rs_hdr[j].thres)) { + header_found = 1; + if (mv[j] < 0) header_found = -1; + break; + } + mv[j] = 0.0; + } + } + } + else header_found = 1; + + if (header_found) { + if (!option_silent) { + fprintf(stdout, "sample: %d\n", mv_pos[j]); + fprintf(stdout, "%s: %.4f\n", rs_hdr[j].type, mv[j]); + } + if ((j < 3) && mv[j] < 0) header_found = -1; + } + } + } + } + } + + if (header_found) break; + header_found = 0; + for (j = 0; j < Nrs; j++) mv[j] = 0.0; + } + +ende: + free_buffers(); + fclose(fp); + + // return only best result + // latest: j + k = j; + j_max = 0; + mv_max = 0.0; + for (j = 0; j < Nrs; j++) { + if ( fabs(mv_max) < fabs(mv[j]) ) { + mv_max = mv[j]; + j_max = j; + } + } + + + // rs_hdr[k].tn + return (header_found * rs_hdr[j_max].tn); +} +