From 1d04aeaa5ae7db0f4113c4bf27cecb3b4487ed8d Mon Sep 17 00:00:00 2001 From: Mark Jessop Date: Sat, 2 Feb 2019 15:23:26 +1030 Subject: [PATCH] Update to latest upstream decoders. Use --json options on decoders. Respect inverted DFM flag. --- auto_rx/auto_rx.py | 23 +- auto_rx/autorx/decode.py | 34 +- auto_rx/autorx/scan.py | 8 +- demod/README.md | 27 + demod/demod_dft.c | 95 ++- demod/demod_dft.h | 4 +- demod/demod_iq.c | 883 ++++++++++++++++++++++++ demod/demod_iq.h | 19 + demod/dfm09dm_dft.c | 163 ++++- demod/lms6dm_dft.c | 1065 +++++++++++++++++++++++++++++ demod/m10dm_dft.c | 80 ++- demod/rs41dm_dft.c | 64 +- demod/rs41dm_iq.c | 1390 ++++++++++++++++++++++++++++++++++++++ demod/rs92dm_dft.c | 23 +- ecc/bch_ecc.c | 84 ++- 15 files changed, 3854 insertions(+), 108 deletions(-) create mode 100644 demod/README.md create mode 100644 demod/demod_iq.c create mode 100644 demod/demod_iq.h create mode 100644 demod/lms6dm_dft.c create mode 100644 demod/rs41dm_iq.c diff --git a/auto_rx/auto_rx.py b/auto_rx/auto_rx.py index 15f27af..fc8a32f 100644 --- a/auto_rx/auto_rx.py +++ b/auto_rx/auto_rx.py @@ -231,11 +231,21 @@ def handle_scan_results(): # Already decoding this sonde, continue. continue else: - logging.info("Detected new %s sonde on %.3f MHz!" % (_type, _freq/1e6)) + + # Handle an inverted sonde detection. + if _type.startswith('-'): + _inverted = " (Inverted)" + _check_type = _type[1:] + else: + _check_type = _type + _inverted = "" + + # Note: We don't indicate if it's been detected as inverted here. + logging.info("Detected new %s sonde on %.3f MHz!" % (_check_type, _freq/1e6)) # Break if we don't support this sonde type. - if (_type not in VALID_SONDE_TYPES): - logging.error("Unsupported sonde type: %s" % _type) + if (_check_type not in VALID_SONDE_TYPES): + logging.error("Unsupported sonde type: %s" % _check_type) continue if allocate_sdr(check_only=True) is not None : @@ -355,7 +365,12 @@ def telemetry_filter(telemetry): if vaisala_callsign_valid or dfm_callsign_valid: return True else: - logging.warning("Payload ID %s does not match regex. Discarding." % telemetry['id']) + _id_msg = "Payload ID %s is invalid." % telemetry['id'] + # Add in a note about DFM sondes and their oddness... + if 'DFM' in telemetry['id']: + _id_msg += " Note: DFM sondes may take a while to get an ID." + + logging.warning(_id_msg) return False diff --git a/auto_rx/autorx/decode.py b/auto_rx/autorx/decode.py index 108683a..cdef477 100644 --- a/auto_rx/autorx/decode.py +++ b/auto_rx/autorx/decode.py @@ -120,6 +120,14 @@ class SondeDecoder(object): # This will become our decoder thread. self.decoder = None + # Detect if we have an 'inverted' sonde. + if self.sonde_type.startswith('-'): + self.inverted = True + # Strip off the leading '-' character' + self.sonde_type = self.sonde_type[1:] + else: + self.inverted = False + # Check if the sonde type is valid. if self.sonde_type not in self.VALID_SONDE_TYPES: self.log_error("Unsupported sonde type: %s" % self.sonde_type) @@ -201,7 +209,7 @@ class SondeDecoder(object): # 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 - lowpass 2600 2>/dev/null |" - decode_cmd += "./rs41ecc --crc --ecc --ptu 2>/dev/null" + decode_cmd += "./rs41ecc --crc --ecc --ptu --json 2>/dev/null" elif self.sonde_type == "RS92": # Decoding a RS92 requires either an ephemeris or an almanac file. @@ -232,16 +240,30 @@ class SondeDecoder(object): # rtl_fm -p 0 -g 26.0 -M fm -F9 -s 12k -f 400500000 | sox -t raw -r 12k -e s -b 16 -c 1 - -r 48000 -b 8 -t wav - highpass 20 lowpass 2500 2>/dev/null | ./rs92ecc -vx -v --crc --ecc --vel -e ephemeris.dat decode_cmd = "%s %s-p %d -d %s %s-M fm -F9 -s 12k -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 12k -e s -b 16 -c 1 - -r 48000 -b 8 -t wav - lowpass 2500 highpass 20 2>/dev/null |" - decode_cmd += "./rs92ecc -vx -v --crc --ecc --vel %s 2>/dev/null" % _rs92_gps_data + decode_cmd += "./rs92ecc -vx -v --crc --ecc --vel --json %s 2>/dev/null" % _rs92_gps_data elif self.sonde_type == "DFM": - # DFM06/DFM09 Sondes + # 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 = "" # 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 2>/dev/null" + decode_cmd += "./dfm09ecc -vv --ecc --json %s 2>/dev/null" % _invert_flag else: @@ -257,6 +279,8 @@ class SondeDecoder(object): # Timeout Counter. _last_packet = time.time() + self.log_debug("Decoder Command: %s" % self.decoder_command ) + # Start the thread. self.decode_process = subprocess.Popen(self.decoder_command, shell=True, stdin=None, stdout=subprocess.PIPE, preexec_fn=os.setsid) self.async_reader = AsynchronousFileReader(self.decode_process.stdout, autostart=True) @@ -362,7 +386,7 @@ class SondeDecoder(object): try: _telemetry['datetime_dt'] = parse(_telemetry['datetime']) except Exception as e: - self.log_error("Invalid date/time in telemetry dict - %s" % str(e)) + self.log_error("Invalid date/time in telemetry dict - %s (Sonde may not have GPS lock" % str(e)) return False # Add in the sonde frequency and type fields. diff --git a/auto_rx/autorx/scan.py b/auto_rx/autorx/scan.py index 7b60c39..4f0082c 100644 --- a/auto_rx/autorx/scan.py +++ b/auto_rx/autorx/scan.py @@ -191,8 +191,8 @@ def detect_sonde(frequency, rs_path="./", dwell_time=10, sdr_fm='rtl_fm', device else: gain_param = '' - rx_test_command = "timeout %ds %s %s-p %d -d %s %s-M fm -F9 -s 15k -f %d 2>/dev/null |" % (dwell_time, sdr_fm, bias_option, int(ppm), str(device_idx), gain_param, frequency) - rx_test_command += "sox -t raw -r 15k -e s -b 16 -c 1 - -r 48000 -t wav - highpass 20 2>/dev/null |" + rx_test_command = "timeout %ds %s %s-p %d -d %s %s-M fm -F9 -s 22k -f %d 2>/dev/null |" % (dwell_time, sdr_fm, bias_option, int(ppm), str(device_idx), gain_param, frequency) + rx_test_command += "sox -t raw -r 22k -e s -b 16 -c 1 - -r 48000 -t wav - highpass 20 2>/dev/null |" rx_test_command += os.path.join(rs_path,"rs_detect") + " -z -t 8 2>/dev/null >/dev/null" logging.debug("Scanner #%s - Attempting sonde detection on %.3f MHz" % (str(device_idx), frequency/1e6)) @@ -218,8 +218,8 @@ def detect_sonde(frequency, rs_path="./", dwell_time=10, sdr_fm='rtl_fm', device if (ret_code & 0x80) > 0: # If the inverted bit is set, we have to do some munging of the return code to get the sonde type. ret_code = abs(-1 * (0x100 - ret_code)) - # Currently ignoring the inverted flag, as rs_detect appears to detect some sondes as inverted incorrectly. - #inv = "-" + + inv = "-" else: ret_code = abs(ret_code) diff --git a/demod/README.md b/demod/README.md new file mode 100644 index 0000000..a638107 --- /dev/null +++ b/demod/README.md @@ -0,0 +1,27 @@ + +## Radiosonde decoders + +alternative decoders using cross-correlation for better header-synchronization + +#### Files + + * `demod_dft.c`, `demod_dft.h`,
+ `rs41dm_dft.c`, `rs92dm_dft.c`, `dfm09dm_dft.c`, `m10dm_dft.c`, `lms6dm_dft.c`,
+ `RS/ecc/bch_ecc.c` + +#### Compile + (copy `bch_ecc.c`)
+ `gcc -c demod_dft.c`
+ `gcc rs41dm_dft.c demod_dft.o -lm -o rs41dm_dft`
+ `gcc dfm09dm_dft.c demod_dft.o -lm -o dfm09dm_dft`
+ `gcc m10dm_dft.c demod_dft.o -lm -o m10dm_dft`
+ `gcc lms6dm_dft.c demod_dft.o -lm -o lms6dm_dft`
+ `gcc rs92dm_dft.c demod_dft.o -lm -o rs92dm_dft` (needs `RS/rs92/nav_gps_vel.c`) + +#### Usage/Examples + `./rs41dm_dft --ecc2 --crc -vx --ptu `
+ `./dfm09dm_dft --ecc -v --ptu ` (add `-i` for dfm06)
+ `./m10dm_dft --dc -vv --ptu -c `
+ `./lms6dm_dft --vit --ecc -v `
+ + diff --git a/demod/demod_dft.c b/demod/demod_dft.c index cc872cc..e90288b 100644 --- a/demod/demod_dft.c +++ b/demod/demod_dft.c @@ -167,6 +167,7 @@ int getCorrDFT(int abs, int K, unsigned int pos, float *maxv, unsigned int *maxv static int sample_rate = 0, bits_sample = 0, channels = 0; static float samples_per_bit = 0; +static int wav_ch = 0; // 0: links bzw. mono; 1: rechts static int findstr(char *buff, char *str, int pos) { int i; @@ -176,7 +177,7 @@ static int findstr(char *buff, char *str, int pos) { return i; } -float read_wav_header(FILE *fp, float baudrate) { +float read_wav_header(FILE *fp, float baudrate, int wav_channel) { char txt[4+1] = "\0\0\0\0"; unsigned char dat[4]; int byte, p=0; @@ -224,6 +225,10 @@ float read_wav_header(FILE *fp, float baudrate) { 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; samples_per_bit = sample_rate/baudrate; @@ -241,7 +246,7 @@ static int f32read_sample(FILE *fp, float *s) { if (fread( &b, bits_sample/8, 1, fp) != 1) return EOF; - if (i == 0) { // i = 0: links bzw. mono + 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; @@ -411,6 +416,92 @@ int read_sbit(FILE *fp, int symlen, int *bit, int inv, int ofs, int reset, int c /* -------------------------------------------------------------------------- */ +int read_softbit(FILE *fp, int symlen, int *bit, float *sb, float level, int inv, int ofs, int reset, int cm) { +// symlen==2: manchester2 10->0,01->1: 2.bit + + static double bitgrenze; + static unsigned long scount; + + float sample; + + double sum = 0.0; + int n = 0; + + if (reset) { + scount = 0; + bitgrenze = 0; + } + + if (symlen == 2) { + bitgrenze += samples_per_bit; + do { + if (buffered > 0) buffered -= 1; + else if (f32buf_sample(fp, inv, cm) == EOF) return EOF; + + sample = bufs[(sample_out-buffered + ofs + M) % M]; + if (scount > bitgrenze-samples_per_bit && scount < bitgrenze-2) + { + sum -= sample; + n++; + } + scount++; + } while (scount < bitgrenze); // n < samples_per_bit + } + + bitgrenze += samples_per_bit; + do { + if (buffered > 0) buffered -= 1; + else if (f32buf_sample(fp, inv, cm) == EOF) return EOF; + + sample = bufs[(sample_out-buffered + ofs + M) % M]; + if (scount > bitgrenze-samples_per_bit && scount < bitgrenze-2) + { + sum += sample; + n++; + } + scount++; + } while (scount < bitgrenze); // n < samples_per_bit + + if (sum >= 0) *bit = 1; + else *bit = 0; + + *sb = sum / n; + + if (*sb > +2.5*level) *sb = +0.8*level; + if (*sb > +level) *sb = +level; + + if (*sb < -2.5*level) *sb = -0.8*level; + if (*sb < -level) *sb = -level; + + *sb /= level; + + return 0; +} + +float header_level(char hdr[], int hLen, unsigned int pos, int inv) { + int n, bitn; + int sgn = 0; + double s = 0.0; + double sum = 0.0; + + n = 0; + bitn = 0; + while ( bitn < hLen && (n < N) ) { + sgn = (hdr[bitn]&1)*2-1; // {'0','1'} -> {-1,1} + s = bufs[(pos-N + n + M) % M]; + if (inv) s = -s; + sum += s * sgn; + n++; + bitn = n / samples_per_bit; + } + sum /= n; + + return sum; +} + +/* -------------------------------------------------------------------------- */ + + static double norm2_match() { int i; double x, y = 0.0; diff --git a/demod/demod_dft.h b/demod/demod_dft.h index f29d463..b54221c 100644 --- a/demod/demod_dft.h +++ b/demod/demod_dft.h @@ -1,7 +1,9 @@ -float read_wav_header(FILE*, float); +float read_wav_header(FILE*, float, int); int f32buf_sample(FILE*, int, int); int read_sbit(FILE*, int, int*, int, int, int, int); +int read_softbit(FILE *fp, int symlen, int *bit, float *sb, float level, int inv, int ofs, int reset, int cm); +float header_level(char hdr[], int hLen, unsigned int pos, int inv); int getCorrDFT(int, int, unsigned int, float *, unsigned int *); int headcmp(int, char*, int, unsigned int, int, int); diff --git a/demod/demod_iq.c b/demod/demod_iq.c new file mode 100644 index 0000000..1d2c12d --- /dev/null +++ b/demod/demod_iq.c @@ -0,0 +1,883 @@ + +/* + * sync header: correlation/matched filter + * compile: + * gcc -c demod_iq.c + * + * author: zilog80 + */ + +/* ------------------------------------------------------------------------------------ */ + +#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; + +#include "demod_iq.h" + + +static int sample_rate = 0, bits_sample = 0, channels = 0; +static float samples_per_bit = 0; + +static unsigned int sample_in, sample_out, delay; +static int buffered = 0; + +static int N, M; + +static float *match = NULL, + *bufs = NULL; + +static char *rawbits = NULL; + +static int Nvar = 0; // < M +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 option_iq = 0; + +/* ------------------------------------------------------------------------------------ */ + +#include + +static int LOG2N, N_DFT; +static int M_DFT; + +static float complex *ew; + +static float complex *Fm, *X, *Z, *cx; +static float *xn; + +static float complex *Hann; + +static int N_IQBUF; +static float complex *raw_iqbuf = NULL; +static float complex *rot_iqbuf = NULL; + +static double df = 0.0; +static int len_sq = 0; + +static unsigned int sample_posframe = 0; +static unsigned int sample_posnoise = 0; + +static double V_noise = 0.0; +static double V_signal = 0.0; +static double SNRdB = 0.0; + + +static void cdft(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 rdft(float *x, float complex *Z) { + int i; + for (i = 0; i < N_DFT; i++) Z[i] = (float complex)x[i]; + cdft(Z); +} + +static void Nidft(float complex *Z, float complex *z) { + int i; + for (i = 0; i < N_DFT; i++) z[i] = conj(Z[i]); + cdft(z); + // idft(): + // for (i = 0; i < N_DFT; i++) z[i] = conj(z[i])/(float)N_DFT; // hier: z reell +} + +static float bin2freq(int k) { + float fq = sample_rate * k / N_DFT; + if (fq > sample_rate/2.0) fq -= sample_rate; + return fq; +} + +static int max_bin(float complex *Z) { + int k, kmax; + double max; + + max = 0; kmax = 0; + for (k = 0; k < N_DFT; k++) { + if (cabs(Z[k]) > max) { + max = cabs(Z[k]); + kmax = k; + } + } + + return kmax; +} + +/* ------------------------------------------------------------------------------------ */ + +int getCorrDFT(int abs, int K, unsigned int pos, float *maxv, unsigned int *maxvpos) { + int i; + int mp = -1; + float mx = 0.0; + float xnorm = 1; + unsigned int mpos = 0; + + dc = 0.0; + + if (N + K > N_DFT/2 - 2) return -1; + if (sample_in < delay+N+K) return -2; + + if (pos == 0) pos = sample_out; + + + for (i = 0; i < N+K; i++) xn[i] = bufs[(pos+M -(N+K-1) + i) % M]; + while (i < N_DFT) xn[i++] = 0.0; + + rdft(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]*Fm[i]; + + Nidft(Z, cx); + + + if (abs) { + for (i = N; i < N+K; i++) { + if (fabs(creal(cx[i])) > fabs(mx)) { // imag(cx)=0 + mx = creal(cx[i]); + mp = i; + } + } + } + else { + for (i = N; i < N+K; i++) { + if (creal(cx[i]) > mx) { // imag(cx)=0 + mx = creal(cx[i]); + mp = i; + } + } + } + if (mp == N || mp == N+K-1) return -4; // Randwert + + mpos = pos - ( N+K-1 - mp ); + xnorm = sqrt(qs[(mpos + 2*M) % M]); + mx /= xnorm*N_DFT; + + *maxv = mx; + *maxvpos = mpos; + + if (pos == sample_out) buffered = sample_out-mpos; + + return mp; +} + +/* ------------------------------------------------------------------------------------ */ + +static int wav_ch = 0; // 0: links bzw. mono; 1: rechts + +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; +} + +float read_wav_header(FILE *fp, float baudrate, 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; + + samples_per_bit = sample_rate/baudrate; + + fprintf(stderr, "samples/bit: %.2f\n", samples_per_bit); + + return samples_per_bit; +} + +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 f32read_csample(FILE *fp, float complex *z) { + short x = 0, y = 0; + + if (fread( &x, bits_sample/8, 1, fp) != 1) return EOF; + if (fread( &y, bits_sample/8, 1, fp) != 1) return EOF; + + *z = x + I*y; + + if (bits_sample == 8) { *z -= 128 + I*128; } + *z /= 128.0; + if (bits_sample == 16) { *z /= 256.0; } + + return 0; +} + +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; +} + +float get_bufmu(int ofs) { + float mu = xs[(sample_out+M + ofs) % M]/Nvar; + return mu; +} + +int f32buf_sample(FILE *fp, int inv, int cm) { + float s = 0.0; + float xneu, xalt; + + float complex z, w; + static float complex z0; //= 1.0; + double gain = 1.0; + + double t = sample_in / (double)sample_rate; + + + if (option_iq) { + + if ( f32read_csample(fp, &z) == EOF ) return EOF; + raw_iqbuf[sample_in % N_IQBUF] = z; + + z *= cexp(-t*2*M_PI*df*I); + w = z * conj(z0); + s = gain * carg(w)/M_PI; + z0 = z; + rot_iqbuf[sample_in % N_IQBUF] = z; + + if (sample_posnoise > 0) + { + if (sample_out >= sample_posframe && sample_out < sample_posframe+len_sq) { + if (sample_out == sample_posframe) V_signal = 0.0; + V_signal += cabs(rot_iqbuf[sample_out % N_IQBUF]); + } + if (sample_out == sample_posframe+len_sq) V_signal /= (double)len_sq; + + if (sample_out >= sample_posnoise && sample_out < sample_posnoise+len_sq) { + if (sample_out == sample_posnoise) V_noise = 0.0; + V_noise += cabs(rot_iqbuf[sample_out % N_IQBUF]); + } + if (sample_out == sample_posnoise+len_sq) { + V_noise /= (double)len_sq; + if (V_signal > 0 && V_noise > 0) { + // iq-samples/V [-1..1] + // dBw = 2*dBv, P=c*U*U + // dBw = 2*10*log10(V/V0) + SNRdB = 20.0 * log10(V_signal/V_noise+1e-20); + } + } + } + + + if (option_iq >= 2) + { + double xbit = 0.0; + double h = 1.0; // modulation index, GFSK; h(rs41)=0.8? // rs-depend... + //float complex xi = cexp(+I*M_PI*h/samples_per_bit); + double f1 = -h*sample_rate/(2*samples_per_bit); + double f2 = -f1; + + float complex X1 = 0; + float complex X2 = 0; + + int n = samples_per_bit; + while (n > 0) { + n--; + t = -n / (double)sample_rate; + z = rot_iqbuf[(sample_in - n + N_IQBUF) % N_IQBUF]; // +1 + X1 += z*cexp(-t*2*M_PI*f1*I); + X2 += z*cexp(-t*2*M_PI*f2*I); + } + + xbit = cabs(X2) - cabs(X1); + + s = xbit / samples_per_bit; + } + } + else { + if (f32read_sample(fp, &s) == EOF) return EOF; + } + + if (inv) s = -s; // swap IQ? + bufs[sample_in % M] = s - dc_ofs; + + xneu = bufs[(sample_in ) % M]; + xalt = bufs[(sample_in+M - Nvar) % M]; + xsum += xneu - xalt; // + xneu - xalt + qsum += (xneu - xalt)*(xneu + xalt); // + xneu*xneu - xalt*xalt + xs[sample_in % M] = xsum; + 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) { +// 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 += samples_per_bit; + do { + sum += bufs[(rcount + mvp + M) % M]; + rcount++; + } while (rcount < rbitgrenze); // n < samples_per_bit + + if (symlen == 2) { + rbitgrenze += samples_per_bit; + do { + sum -= bufs[(rcount + mvp + M) % M]; + rcount++; + } while (rcount < rbitgrenze); // n < samples_per_bit + } + + + 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; +} + +int headcmp(int symlen, char *hdr, int len, unsigned int mvp, int inv, int option_dc) { + 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*samples_per_bit), pos==0); + } + 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; +} + +int get_fqofs(int hdrlen, unsigned int mvp, float *freq, float *snr) { + int j; + int buf_start; + int presamples = 256*samples_per_bit; + + if (presamples > M_DFT) presamples = M_DFT; + + buf_start = mvp - hdrlen*samples_per_bit - presamples; + + while (buf_start < 0) buf_start += N_IQBUF; + + for (j = 0; j < M_DFT; j++) { + Z[j] = Hann[j]*raw_iqbuf[(buf_start+j) % N_IQBUF]; + } + while (j < N_DFT) Z[j++] = 0; + + cdft(Z); + df = bin2freq(max_bin(Z)); + + // if |df| 1000.0) df = 0.0; + + + sample_posframe = sample_in; //mvp - hdrlen*samples_per_bit; + sample_posnoise = mvp + sample_rate*7/8.0; + + + *freq = df; + *snr = SNRdB; + + return 0; +} + +/* -------------------------------------------------------------------------- */ + +int read_sbit(FILE *fp, int symlen, int *bit, int inv, int ofs, int reset, int cm) { +// symlen==2: manchester2 10->0,01->1: 2.bit + + static double bitgrenze; + static unsigned long scount; + + float sample; + + double sum = 0.0; + + if (reset) { + scount = 0; + bitgrenze = 0; + } + + if (symlen == 2) { + bitgrenze += samples_per_bit; + do { + if (buffered > 0) buffered -= 1; + else if (f32buf_sample(fp, inv, cm) == EOF) return EOF; + + sample = bufs[(sample_out-buffered + ofs + M) % M]; + sum -= sample; + + scount++; + } while (scount < bitgrenze); // n < samples_per_bit + } + + bitgrenze += samples_per_bit; + do { + if (buffered > 0) buffered -= 1; + else if (f32buf_sample(fp, inv, cm) == EOF) return EOF; + + sample = bufs[(sample_out-buffered + ofs + M) % M]; + sum += sample; + + scount++; + } while (scount < bitgrenze); // n < samples_per_bit + + if (sum >= 0) *bit = 1; + else *bit = 0; + + return 0; +} + +int read_IDsbit(FILE *fp, int symlen, int *bit, int inv, int ofs, int reset, int cm) { +// symlen==2: manchester2 10->0,01->1: 2.bit + + static double bitgrenze; + static unsigned long scount; + + float sample; + + double sum = 0.0; + double mid; + double l = 1.0; + + if (reset) { + scount = 0; + bitgrenze = 0; + } + + if (symlen == 2) { + mid = bitgrenze + (samples_per_bit-1)/2.0; + bitgrenze += samples_per_bit; + do { + if (buffered > 0) buffered -= 1; + else if (f32buf_sample(fp, inv, cm) == EOF) return EOF; + + sample = bufs[(sample_out-buffered + ofs + M) % M]; + if (mid-l < scount && scount < mid+l) sum -= sample; + + scount++; + } while (scount < bitgrenze); // n < samples_per_bit + } + + mid = bitgrenze + (samples_per_bit-1)/2.0; + bitgrenze += samples_per_bit; + do { + if (buffered > 0) buffered -= 1; + else if (f32buf_sample(fp, inv, cm) == EOF) return EOF; + + sample = bufs[(sample_out-buffered + ofs + M) % M]; + if (mid-l < scount && scount < mid+l) sum += sample; + + scount++; + } while (scount < bitgrenze); // n < samples_per_bit + + if (sum >= 0) *bit = 1; + else *bit = 0; + + return 0; +} + +/* -------------------------------------------------------------------------- */ + +int read_softbit(FILE *fp, int symlen, int *bit, float *sb, float level, int inv, int ofs, int reset, int cm) { +// symlen==2: manchester2 10->0,01->1: 2.bit + + static double bitgrenze; + static unsigned long scount; + + float sample; + + double sum = 0.0; + int n = 0; + + if (reset) { + scount = 0; + bitgrenze = 0; + } + + if (symlen == 2) { + bitgrenze += samples_per_bit; + do { + if (buffered > 0) buffered -= 1; + else if (f32buf_sample(fp, inv, cm) == EOF) return EOF; + + sample = bufs[(sample_out-buffered + ofs + M) % M]; + if (scount > bitgrenze-samples_per_bit && scount < bitgrenze-2) + { + sum -= sample; + n++; + } + scount++; + } while (scount < bitgrenze); // n < samples_per_bit + } + + bitgrenze += samples_per_bit; + do { + if (buffered > 0) buffered -= 1; + else if (f32buf_sample(fp, inv, cm) == EOF) return EOF; + + sample = bufs[(sample_out-buffered + ofs + M) % M]; + if (scount > bitgrenze-samples_per_bit && scount < bitgrenze-2) + { + sum += sample; + n++; + } + scount++; + } while (scount < bitgrenze); // n < samples_per_bit + + if (sum >= 0) *bit = 1; + else *bit = 0; + + *sb = sum / n; + + if (*sb > +2.5*level) *sb = +0.8*level; + if (*sb > +level) *sb = +level; + + if (*sb < -2.5*level) *sb = -0.8*level; + if (*sb < -level) *sb = -level; + + *sb /= level; + + return 0; +} + +float header_level(char hdr[], int hLen, unsigned int pos, int inv) { + int n, bitn; + int sgn = 0; + double s = 0.0; + double sum = 0.0; + + n = 0; + bitn = 0; + while ( bitn < hLen && (n < N) ) { + sgn = (hdr[bitn]&1)*2-1; // {'0','1'} -> {-1,1} + s = bufs[(pos-N + n + M) % M]; + if (inv) s = -s; + sum += s * sgn; + n++; + bitn = n / samples_per_bit; + } + sum /= n; + + return sum; +} + +/* -------------------------------------------------------------------------- */ + + +#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() { + int i; + double x, y = 0.0; + for (i = 0; i < N; i++) { + x = match[i]; + y += x*x; + } + return y; +} + +int init_buffers(char hdr[], int hLen, float BT, int opt_iq) { + //hLen = strlen(header) = HEADLEN; + + int i, pos; + float b0, b1, b2, b, t; + float normMatch; + double sigma = sqrt(log(2)) / (2*M_PI*BT); + + int K; + int n, k; + float *m = NULL; + + option_iq = opt_iq; + + N = hLen * samples_per_bit + 0.5; + M = 3*N; + if (samples_per_bit < 6) M = 6*N; + Nvar = N; //N/2; // = N/k + + bufs = (float *)calloc( M+1, sizeof(float)); if (bufs == NULL) return -100; + match = (float *)calloc( N+1, sizeof(float)); if (match == 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; + + + rawbits = (char *)calloc( N+1, sizeof(char)); if (rawbits == NULL) return -100; + + for (i = 0; i < M; i++) bufs[i] = 0.0; + + + for (i = 0; i < N; i++) { + pos = i/samples_per_bit; + t = (i - pos*samples_per_bit)/samples_per_bit - 0.5; + + b1 = ((hdr[pos] & 0x1) - 0.5)*2.0; + b = b1*pulse(t, sigma); + + if (pos > 0) { + b0 = ((hdr[pos-1] & 0x1) - 0.5)*2.0; + b += b0*pulse(t+1, sigma); + } + + if (pos < hLen) { + b2 = ((hdr[pos+1] & 0x1) - 0.5)*2.0; + b += b2*pulse(t-1, sigma); + } + + match[i] = b; + } + + normMatch = sqrt(norm2_match()); + for (i = 0; i < N; i++) { + match[i] /= normMatch; + } + + + delay = N/16; + sample_in = 0; + + K = M-N - delay; //N/2 - delay; // N+K < M + + LOG2N = 2 + (int)(log(N+K)/log(2)); + N_DFT = 1 << LOG2N; + + while (N + K > N_DFT/2 - 2) { + LOG2N += 1; + N_DFT <<= 1; + } + + + xn = calloc(N_DFT+1, sizeof(float)); if (xn == NULL) return -1; + + ew = calloc(LOG2N+1, sizeof(float complex)); if (ew == NULL) return -1; + Fm = calloc(N_DFT+1, sizeof(float complex)); if (Fm == NULL) return -1; + X = calloc(N_DFT+1, sizeof(float complex)); if (X == NULL) return -1; + Z = calloc(N_DFT+1, sizeof(float complex)); if (Z == NULL) return -1; + cx = calloc(N_DFT+1, sizeof(float complex)); if (cx == NULL) return -1; + + M_DFT = M; + Hann = calloc(N_DFT+1, sizeof(float complex)); if (Hann == NULL) return -1; + for (i = 0; i < M_DFT; i++) Hann[i] = 0.5 * (1 - cos( 2 * M_PI * i / (double)(M_DFT-1) ) ); + + for (n = 0; n < LOG2N; n++) { + k = 1 << n; + ew[n] = cexp(-I*M_PI/(float)k); + } + + m = calloc(N_DFT+1, sizeof(float)); if (m == NULL) return -1; + for (i = 0; i < N; i++) m[N-1 - i] = match[i]; + while (i < N_DFT) m[i++] = 0.0; + rdft(m, Fm); + + free(m); m = NULL; + + + if (option_iq) + { + if (channels < 2) return -1; +/* + M_DFT = samples_per_bit*256+0.5; + while ( (1 << LOG2N) < M_DFT ) LOG2N++; + LOG2N++; + N_DFT = (1 << LOG2N); + N_IQBUF = M_DFT + samples_per_bit*(64+16); +*/ + N_IQBUF = N_DFT; + raw_iqbuf = calloc(N_IQBUF+1, sizeof(float complex)); if (raw_iqbuf == NULL) return -1; + rot_iqbuf = calloc(N_IQBUF+1, sizeof(float complex)); if (rot_iqbuf == NULL) return -1; + + len_sq = samples_per_bit*8; + } + + + return K; +} + +int free_buffers() { + + if (match) { free(match); match = NULL; } + 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 (ew) { free(ew); ew = NULL; } + if (Fm) { free(Fm); Fm = NULL; } + if (X) { free(X); X = NULL; } + if (Z) { free(Z); Z = NULL; } + if (cx) { free(cx); cx = NULL; } + + if (Hann) { free(Hann); Hann = NULL; } + + if (option_iq) + { + if (raw_iqbuf) { free(raw_iqbuf); raw_iqbuf = NULL; } + if (rot_iqbuf) { free(rot_iqbuf); rot_iqbuf = NULL; } + } + + return 0; +} + +/* ------------------------------------------------------------------------------------ */ + +unsigned int get_sample() { + return sample_out; +} + diff --git a/demod/demod_iq.h b/demod/demod_iq.h new file mode 100644 index 0000000..6c5a2ae --- /dev/null +++ b/demod/demod_iq.h @@ -0,0 +1,19 @@ + +float read_wav_header(FILE*, float, int); +int f32buf_sample(FILE*, int, int); +int read_sbit(FILE*, int, int*, int, int, int, int); +int read_IDsbit(FILE*, int, int*, int, int, int, int); +int read_softbit(FILE *fp, int symlen, int *bit, float *sb, float level, int inv, int ofs, int reset, int cm); +float header_level(char hdr[], int hLen, unsigned int pos, int inv); + +int getCorrDFT(int, int, unsigned int, float *, unsigned int *); +int headcmp(int, char*, int, unsigned int, int, int); +int get_fqofs(int, unsigned int, float *, float *); +float get_bufvar(int); +float get_bufmu(int); + +int init_buffers(char*, int, float, int); +int free_buffers(void); + +unsigned int get_sample(void); + diff --git a/demod/dfm09dm_dft.c b/demod/dfm09dm_dft.c index 9f88936..9ac246b 100644 --- a/demod/dfm09dm_dft.c +++ b/demod/dfm09dm_dft.c @@ -32,7 +32,7 @@ typedef struct { int frnr; int sonde_typ; ui32_t SN6; - ui32_t SN9; + ui32_t SN; int week; int gpssec; int jahr; int monat; int tag; int std; int min; float sek; @@ -55,7 +55,11 @@ int option_verbose = 0, // ausfuehrliche Anzeige option_ecc = 0, option_ptu = 0, option_ths = 0, + option_json = 0, // JSON blob output (for auto_rx) wavloaded = 0; +int wav_channel = 0; // audio channel: left + +int ptu_out = 0; int start = 0; @@ -198,7 +202,7 @@ int dat_out(ui8_t *dat_bits) { if (fr_id == 1) { // 00..31: ? GPS-Sats in Sicht? - msek = bits2val(dat_bits+32, 16); + msek = bits2val(dat_bits+32, 16); // UTC (= GPS - 18sec ab 1.1.2017) gpx.sek = msek/1000.0; } @@ -329,7 +333,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 && option_verbose == 2) { + if (option_ptu && ptu_out && option_verbose == 2) { printf(" (Rso: %.1f , Rb: %.1f)", Rs_o/1e3, Rb/1e3); } @@ -376,22 +380,30 @@ float get_Temp4(float *meas) { // meas[0..4] } -#define SNbit 0x0100 +#define RSNbit 0x0100 // radiosonde DFM-06,DFM-09 +#define PSNbit 0x0200 // pilotsonde PS-15 int conf_out(ui8_t *conf_bits) { int conf_id; int ret = 0; int val, hl; + ui32_t SN6, SN; static int chAbit, chA[2]; - ui32_t SN6, SN9; + static int chCbit, chC[2]; + static int chDbit, chD[2]; + static int ch7bit, ch7[2]; + static ui32_t SN_A, SN_C, SN_D, SN_7; conf_id = bits2val(conf_bits, 4); - //if (conf_id > 6) gpx.SN6 = 0; //// gpx.sonde_typ & 0xF = 9; // SNbit? - - if ((gpx.sonde_typ & 0xFF) < 9 && conf_id == 6) { - SN6 = bits2val(conf_bits+4, 4*6); // DFM-06: Kanal 6 - if ( SN6 == gpx.SN6 ) { // nur Nibble-Werte 0..9 - gpx.sonde_typ = SNbit | 6; + // gibt es Kanaele > 6 (2-teilige ID)? + // if (conf_id > 6) gpx.SN6 = 0; // -> DFM-09,PS-15 // SNbit? + // + // 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; + ptu_out = 6; ret = 6; } else { @@ -401,22 +413,84 @@ int conf_out(ui8_t *conf_bits) { } if (conf_id == 0xA) { // 0xACxxxxy val = bits2val(conf_bits+8, 4*5); - hl = (val & 1) == 0; + hl = (val & 1); // val&0xF 0,1? chA[hl] = (val >> 4) & 0xFFFF; chAbit |= 1 << hl; if (chAbit == 3) { // DFM-09: Kanal A - SN9 = (chA[1] << 16) | chA[0]; - if ( SN9 == gpx.SN9 ) { - gpx.sonde_typ = SNbit | 9; + SN = (chA[0] << 16) | chA[1]; + if ( SN == SN_A ) { + gpx.sonde_typ = RSNbit | 0xA; + gpx.SN = SN; + ptu_out = 9; ret = 9; } else { gpx.sonde_typ = 0; } - gpx.SN9 = SN9; + SN_A = SN; chAbit = 0; } } + if (conf_id == 0xC) { // 0xCCxxxxy + val = bits2val(conf_bits+8, 4*5); + hl = (val & 1); + chC[hl] = (val >> 4) & 0xFFFF; + chCbit |= 1 << hl; + if (chCbit == 3) { // DFM-17? Kanal C + SN = (chC[0] << 16) | chC[1]; + if ( SN == SN_C ) { + gpx.sonde_typ = RSNbit | 0xC; + gpx.SN = SN; + ptu_out = 9; + ret = 17; + } + else { + gpx.sonde_typ = 0; + } + SN_C = SN; + chCbit = 0; + } + } + if (conf_id == 0xD) { // 0xDCxxxxy + val = bits2val(conf_bits+8, 4*5); + hl = (val & 1); + chD[hl] = (val >> 4) & 0xFFFF; + chDbit |= 1 << hl; + if (chDbit == 3) { // DFM-17? Kanal D + SN = (chD[0] << 16) | chD[1]; + if ( SN == SN_D ) { + gpx.sonde_typ = RSNbit | 0xD; + gpx.SN = SN; + ptu_out = 9; + ret = 18; + } + else { + gpx.sonde_typ = 0; + } + SN_D = SN; + chDbit = 0; + } + } + if (conf_id == 0x7) { // 0x70xxxxy + val = bits2val(conf_bits+8, 4*5); + hl = (val & 1); + ch7[hl] = (val >> 4) & 0xFFFF; + ch7bit |= 1 << hl; + if (ch7bit == 3) { // PS-15: Kanal 7 + SN = (ch7[0] << 16) | ch7[1]; + if ( SN == SN_7 ) { + gpx.sonde_typ = PSNbit | 0x7; + gpx.SN = SN; + ptu_out = 0; + ret = 15; + } + else { + gpx.sonde_typ = 0; + } + SN_7 = SN; + ch7bit = 0; + } + } if (conf_id >= 0 && conf_id <= 4) { val = bits2val(conf_bits+4, 4*6); @@ -430,7 +504,7 @@ int conf_out(ui8_t *conf_bits) { } // STM32-status: Bat, MCU-Temp - if ((gpx.sonde_typ & 0xFF) == 9) { // DFM-09 (STM32) + if ((gpx.sonde_typ & 0xF) == 0xA) { // DFM-09 (STM32) if (conf_id == 0x5) { // voltage val = bits2val(conf_bits+8, 4*4); gpx.status[0] = val/1000.0; @@ -469,7 +543,7 @@ void print_gpx() { printf(" vH: %5.2f ", gpx.horiV); printf(" D: %5.1f ", gpx.dir); printf(" vV: %5.2f ", gpx.vertV); - if (option_ptu) { + if (option_ptu && ptu_out) { float t = get_Temp(gpx.meas24); if (t > -270.0) printf(" T=%.1fC ", t); if (option_verbose == 2) { @@ -482,31 +556,47 @@ void print_gpx() { printf(" f4: %.2f ", gpx.meas24[4]); } } - if (option_verbose == 2 && (gpx.sonde_typ & 0xFF) == 9) { + if (option_verbose == 2 && (gpx.sonde_typ & 0xF) == 0xA) { printf(" U: %.2fV ", gpx.status[0]); printf(" Ti: %.1fK ", gpx.status[1]); } - if (option_verbose && (gpx.sonde_typ & SNbit)) + if ( option_verbose ) { - if ((gpx.sonde_typ & 0xFF) == 6) { - printf(" (ID%1d:%06X) ", gpx.sonde_typ & 0xF, gpx.SN6); - sprintf(sonde_id, "DFM06-%06X", gpx.SN6); + 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 & 0xFF) == 9) { - printf(" (ID%1d:%06u) ", gpx.sonde_typ & 0xF, gpx.SN9); - sprintf(sonde_id, "DFM09-%06u", gpx.SN9); + 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; } - gpx.sonde_typ ^= SNbit; } } printf("\n"); - // Print JSON blob - // Get temperature - float t = get_Temp(gpx.meas24); - 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, 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, t ); - - + if (option_json) + { + // Print JSON blob + // Get temperature + float t = get_Temp(gpx.meas24); + 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, 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, t ); + } } } @@ -636,6 +726,7 @@ int main(int argc, char **argv) { fprintf(stderr, " -i, --invert\n"); fprintf(stderr, " --ecc (Hamming ECC)\n"); fprintf(stderr, " --ths (peak threshold; default=%.1f)\n", thres); + fprintf(stderr, " --json (JSON output)\n"); return 0; } else if ( (strcmp(*argv, "-v") == 0) || (strcmp(*argv, "--verbose") == 0) ) { @@ -652,7 +743,9 @@ int main(int argc, char **argv) { option_inv = 0x1; } else if ( (strcmp(*argv, "--ecc") == 0) ) { option_ecc = 1; } - else if ( (strcmp(*argv, "--ptu") == 0) ) { option_ptu = 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, "--ch2") == 0) ) { wav_channel = 1; } // right channel (default: 0=left) else if ( (strcmp(*argv, "--ths") == 0) ) { ++argv; if (*argv) { @@ -673,7 +766,7 @@ int main(int argc, char **argv) { if (!wavloaded) fp = stdin; - spb = read_wav_header(fp, (float)BAUD_RATE); + spb = read_wav_header(fp, (float)BAUD_RATE, wav_channel); if ( spb < 0 ) { fclose(fp); fprintf(stderr, "error: wav header\n"); diff --git a/demod/lms6dm_dft.c b/demod/lms6dm_dft.c new file mode 100644 index 0000000..bc280bd --- /dev/null +++ b/demod/lms6dm_dft.c @@ -0,0 +1,1065 @@ + +/* + * LMS6 + * (403 MHz) + * + * sync header: correlation/matched filter + * files: lms6dm_dft.c demod_dft.h demod_dft.c bch_ecc.c + * compile: + * gcc -c demod_dft.c + * gcc lms6dm_dft.c demod_dft.o -lm -o lms6dm_dft + * usage: + * ./lms6dm_dft -v --vit --ecc + * + * author: zilog80 + */ + +#include +#include +#include +#include + +#ifdef CYGWIN + #include // cygwin: _setmode() + #include +#endif + + +typedef unsigned char ui8_t; +typedef unsigned short ui16_t; +typedef unsigned int ui32_t; + +#include "demod_dft.h" + +#include "bch_ecc.c" // RS/ecc/ + + +int option_verbose = 0, // ausfuehrliche Anzeige + option_raw = 0, // rohe Frames + option_ecc = 0, + option_vit = 0, + option_inv = 0, // invertiert Signal + option_dc = 0, + wavloaded = 0; +int wav_channel = 0; // audio channel: left + + +/* -------------------------------------------------------------------------- */ + +#define BAUD_RATE 4800 + +#define BITS 8 +#define HEADOFS 0 //16 +#define HEADLEN ((4*16)-HEADOFS) + +char rawheader[] = "0101011000001000""0001110010010111""0001101010100111""0011110100111110"; // (c0,inv(c1)) +// (00) 58 f3 3f b8 +//char header[] = "0000001101011101""0100100111000010""0100111111110010""0110100001101011"; // (c0,c1) +ui8_t rs_sync[] = { 0x00, 0x58, 0xf3, 0x3f, 0xb8}; +// 0x58f33fb8 little-endian <-> 0x1ACFFC1D big-endian bytes + +#define SYNC_LEN 5 +#define FRM_LEN (223) +#define PAR_LEN (32) +#define FRMBUF_LEN (3*FRM_LEN) +#define BLOCKSTART (SYNC_LEN*BITS*2) +#define BLOCK_LEN (FRM_LEN+PAR_LEN+SYNC_LEN) // 255+5 = 260 +#define RAWBITBLOCK_LEN ((BLOCK_LEN+1)*BITS*2) // (+1 tail) + +// (00) 58 f3 3f b8 +char blk_rawbits[RAWBITBLOCK_LEN+SYNC_LEN*BITS*2 +8] = "0000000000000000""0000001101011101""0100100111000010""0100111111110010""0110100001101011"; +//char *block_rawbits = blk_rawbits+SYNC_LEN*BITS*2; + +float soft_rawbits[RAWBITBLOCK_LEN+SYNC_LEN*BITS*2 +8] = + { -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, + -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0, 1.0, -1.0, 1.0, + -1.0, 1.0, -1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.0, -1.0, + -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0, 1.0, -1.0, + -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0 }; + +ui8_t block_bytes[BLOCK_LEN+8]; + + +ui8_t frm_sync[] = { 0x24, 0x54, 0x00, 0x00}; +ui8_t frame[FRM_LEN] = { 0x24, 0x54, 0x00, 0x00}; // dataheader + +ui8_t *p_frame = frame; + + +#define FRAME_LEN (300) // 4800baud, 16bits/byte +#define BITFRAME_LEN (FRAME_LEN*BITS) +#define RAWBITFRAME_LEN (BITFRAME_LEN*2) +#define OVERLAP 64 +#define OFS 4 + + +char frame_bits[BITFRAME_LEN+OVERLAP*BITS +8]; // init L-1 bits mit 0 + +#define L 7 // d_f=10 + char polyA[] = "1001111"; // 0x4f: x^6+x^3+x^2+x+1 + char polyB[] = "1101101"; // 0x6d: x^6+x^5+x^3+x^2+1 +/* +// d_f=6 +qA[] = "1110011"; // 0x73: x^6+x^5+x^4+x+1 +qB[] = "0011110"; // 0x1e: x^4+x^3+x^2+x +pA[] = "10010101"; // 0x95: x^7+x^4+x^2+1 = (x+1)(x^6+x^5+x^4+x+1) = (x+1)qA +pB[] = "00100010"; // 0x22: x^5+x = (x+1)(x^4+x^3+x^2+x)=x(x+1)^3 = (x+1)qB +polyA = qA + x*qB +polyB = qA + qB +*/ + +char vit_rawbits[RAWBITFRAME_LEN+OVERLAP*BITS*2 +8]; +char vits_rawbits[RAWBITFRAME_LEN+OVERLAP*BITS*2 +8]; +char vits_bits[BITFRAME_LEN+OVERLAP*BITS +8]; + +#define N (1 << L) +#define M (1 << (L-1)) + +typedef struct { + ui8_t bIn; + ui8_t codeIn; + int w; + int prevState; + float sw; +} states_t; + +states_t vit_state[RAWBITFRAME_LEN+OVERLAP +8][M]; + +states_t vit_d[N]; + +ui8_t vit_code[N]; + + +int vit_initCodes() { + int cA, cB; + int i, bits; + + for (bits = 0; bits < N; bits++) { + cA = 0; + cB = 0; + for (i = 0; i < L; i++) { + cA ^= (polyA[L-1-i]&1) & ((bits >> i)&1); + cB ^= (polyB[L-1-i]&1) & ((bits >> i)&1); + } + vit_code[bits] = (cA<<1) | cB; + } + + return 0; +} + +int vit_dist(int c, char *rc) { + return (((c>>1)^rc[0])&1) + ((c^rc[1])&1); +} + +int vit_start(char *rc) { + int t, m, j, c, d; + + t = L-1; + m = M; + while ( t > 0 ) { // t=0..L-2: nextState 0) { + c = vit_state[t][j].codeIn; + vit_rawbits[2*t -2] = 0x30 + ((c>>1) & 1); + vit_rawbits[2*t -1] = 0x30 + (c & 1); + j = vit_state[t][j].prevState; + t--; + } + + return 0; +} + +int viterbi(char *rc) { + int t, tmax; + int j, j_min, w_min; + + vit_start(rc); + + tmax = strlen(rc)/2; + + for (t = L-1; t < tmax; t++) + { + vit_next(t, rc+2*t); + } + + w_min = -1; + for (j = 0; j < M; j++) { + if (w_min < 0) { + w_min = vit_state[tmax][j].w; + j_min = j; + } + if (vit_state[tmax][j].w < w_min) { + w_min = vit_state[tmax][j].w; + j_min = j; + } + } + vit_path(j_min, tmax); + + return 0; +} + + +float vits_dist(int c, float *rc) { + int bit0 = ((c>>1)&1) * 2 - 1; + int bit1 = (c&1) * 2 - 1; + return sqrt( (bit0-rc[0])*(bit0-rc[0]) + (bit1-rc[1])*(bit1-rc[1]) ); +} + +int vits_start(float *rc) { + int t, m, j, c; + float d; + + t = L-1; + m = M; + while ( t > 0 ) { // t=0..L-2: nextState 0) { + dec = vit_state[t][j].bIn; + vits_bits[t-1] = 0x30 + dec; + c = vit_state[t][j].codeIn; + vits_rawbits[2*t -2] = 0x30 + ((c>>1) & 1); + vits_rawbits[2*t -1] = 0x30 + (c & 1); + j = vit_state[t][j].prevState; + t--; + } + + return 0; +} + +int viterbi_soft(float *rc, int len) { + int t, tmax; + int j, j_min; + float sw_min; + + vits_start(rc); + + tmax = len/2; + + for (t = L-1; t < tmax; t++) + { + vits_next(t, rc+2*t); + } + + sw_min = -1.0; + for (j = 0; j < M; j++) { + if (sw_min < 0.0) { + sw_min = vit_state[tmax][j].sw; + j_min = j; + } + if (vit_state[tmax][j].sw < sw_min) { + sw_min = vit_state[tmax][j].sw; + j_min = j; + } + } + vits_path(j_min, tmax); + + return 0; +} + +// ------------------------------------------------------------------------ + +int deconv(char* rawbits, char *bits) { + + int j, n, bitA, bitB; + char *p; + int len; + int errors = 0; + int m = L-1; + + len = strlen(rawbits); + for (j = 0; j < m; j++) bits[j] = '0'; + n = 0; + while ( 2*(m+n) < len ) { + p = rawbits+2*(m+n); + bitA = bitB = 0; + for (j = 0; j < m; j++) { + bitA ^= (bits[n+j]&1) & (polyA[j]&1); + bitB ^= (bits[n+j]&1) & (polyB[j]&1); + } + if ( (bitA^(p[0]&1))==(polyA[m]&1) && (bitB^(p[1]&1))==(polyB[m]&1) ) bits[n+m] = '1'; + else if ( (bitA^(p[0]&1))==0 && (bitB^(p[1]&1))==0 ) bits[n+m] = '0'; + else { + if ( (bitA^(p[0]&1))!=(polyA[m]&1) && (bitB^(p[1]&1))==(polyB[m]&1) ) bits[n+m] = 0x39; + else bits[n+m] = 0x38; + errors = n; + break; + } + n += 1; + } + bits[n+m] = '\0'; + + return errors; +} + +// ------------------------------------------------------------------------ + +int crc16_0(ui8_t frame[], int len) { + int crc16poly = 0x1021; + int rem = 0x0, i, j; + int byte; + + for (i = 0; i < len; i++) { + byte = frame[i]; + rem = rem ^ (byte << 8); + for (j = 0; j < 8; j++) { + if (rem & 0x8000) { + rem = (rem << 1) ^ crc16poly; + } + else { + rem = (rem << 1); + } + rem &= 0xFFFF; + } + } + return rem; +} + +int check_CRC(ui8_t frame[]) { + ui32_t crclen = 0, + crcdat = 0; + + crclen = 221; + crcdat = (frame[crclen]<<8) | frame[crclen+1]; + if ( crcdat != crc16_0(frame, crclen) ) { + return 1; // CRC NO + } + else return 0; // CRC OK +} + +// ------------------------------------------------------------------------ + +int bits2bytes(char *bitstr, ui8_t *bytes) { + int i, bit, d, byteval; + int len = strlen(bitstr)/8; + int bitpos, bytepos; + + bitpos = 0; + bytepos = 0; + + while (bytepos < len) { + + byteval = 0; + d = 1; + for (i = 0; i < BITS; i++) { + bit=*(bitstr+bitpos+i); /* little endian */ + //bit=*(bitstr+bitpos+7-i); /* big endian */ + if ((bit == '1') || (bit == '9')) byteval += d; + else /*if ((bit == '0') || (bit == '8'))*/ byteval += 0; + d <<= 1; + } + bitpos += BITS; + bytes[bytepos++] = byteval & 0xFF; + } + + //while (bytepos < FRAME_LEN+OVERLAP) bytes[bytepos++] = 0; + + return bytepos; +} + +/* -------------------------------------------------------------------------- */ + + +typedef struct { + int frnr; + int sn; + int week; int gpstow; + int jahr; int monat; int tag; + int wday; + int std; int min; float sek; + double lat; double lon; double h; + double vH; double vD; double vV; + double vE; double vN; double vU; + //int freq; +} gpx_t; + +gpx_t gpx; + +gpx_t gpx0 = { 0 }; + + +#define pos_SondeSN (OFS+0x00) // ?4 byte 00 7A.... +#define pos_FrameNb (OFS+0x04) // 2 byte +//GPS Position +#define pos_GPSTOW (OFS+0x06) // 4 byte +#define pos_GPSlat (OFS+0x0E) // 4 byte +#define pos_GPSlon (OFS+0x12) // 4 byte +#define pos_GPSalt (OFS+0x16) // 4 byte +//GPS Velocity East-North-Up (ENU) +#define pos_GPSvO (OFS+0x1A) // 3 byte +#define pos_GPSvN (OFS+0x1D) // 3 byte +#define pos_GPSvV (OFS+0x20) // 3 byte + + +int get_SondeSN() { + unsigned byte; + + byte = (p_frame[pos_SondeSN]<<24) | (p_frame[pos_SondeSN+1]<<16) + | (p_frame[pos_SondeSN+2]<<8) | p_frame[pos_SondeSN+3]; + gpx.sn = byte & 0xFFFFFF; + + return 0; +} + +int get_FrameNb() { + int i; + unsigned byte; + ui8_t frnr_bytes[2]; + int frnr; + + gpx = gpx0; + + for (i = 0; i < 2; i++) { + byte = p_frame[pos_FrameNb + i]; + frnr_bytes[i] = byte; + } + + frnr = (frnr_bytes[0] << 8) + frnr_bytes[1] ; + gpx.frnr = frnr; + + return 0; +} + + +char weekday[7][3] = { "So", "Mo", "Di", "Mi", "Do", "Fr", "Sa"}; +//char weekday[7][4] = { "Sun", "Mon", "Tue", "Wed", "Thu", "Fri", "Sat"}; + +int get_GPStime() { + int i; + unsigned byte; + ui8_t gpstime_bytes[4]; + int gpstime = 0, // 32bit + day; + float ms; + + for (i = 0; i < 4; i++) { + byte = p_frame[pos_GPSTOW + i]; + gpstime_bytes[i] = byte; + } + gpstime = 0; + for (i = 0; i < 4; i++) { + gpstime |= gpstime_bytes[i] << (8*(3-i)); + } + + gpx.gpstow = gpstime; + + ms = gpstime % 1000; + gpstime /= 1000; + + day = gpstime / (24 * 3600); + gpstime %= (24*3600); + + if ((day < 0) || (day > 6)) return -1; + + gpx.wday = day; + gpx.std = gpstime / 3600; + gpx.min = (gpstime % 3600) / 60; + gpx.sek = gpstime % 60 + ms/1000.0; + + return 0; +} + +double B60B60 = 0xB60B60; // 2^32/360 = 0xB60B60.xxx + +int get_GPSlat() { + int i; + unsigned byte; + ui8_t gpslat_bytes[4]; + int gpslat; + double lat; + + for (i = 0; i < 4; i++) { + byte = p_frame[pos_GPSlat + i]; + if (byte > 0xFF) return -1; + gpslat_bytes[i] = byte; + } + + gpslat = 0; + for (i = 0; i < 4; i++) { + gpslat |= gpslat_bytes[i] << (8*(3-i)); + } + lat = gpslat / B60B60; + gpx.lat = lat; + + return 0; +} + +int get_GPSlon() { + int i; + unsigned byte; + ui8_t gpslon_bytes[4]; + int gpslon; + double lon; + + for (i = 0; i < 4; i++) { + byte = p_frame[pos_GPSlon + i]; + if (byte > 0xFF) return -1; + gpslon_bytes[i] = byte; + } + + gpslon = 0; + for (i = 0; i < 4; i++) { + gpslon |= gpslon_bytes[i] << (8*(3-i)); + } + lon = gpslon / B60B60; + gpx.lon = lon; + + return 0; +} + +int get_GPSalt() { + int i; + unsigned byte; + ui8_t gpsheight_bytes[4]; + int gpsheight; + double height; + + for (i = 0; i < 4; i++) { + byte = p_frame[pos_GPSalt + i]; + if (byte > 0xFF) return -1; + gpsheight_bytes[i] = byte; + } + + gpsheight = 0; + for (i = 0; i < 4; i++) { + gpsheight |= gpsheight_bytes[i] << (8*(3-i)); + } + height = gpsheight / 1000.0; + gpx.h = height; + + if (height < -100 || height > 60000) return -1; + return 0; +} + +int get_GPSvel24() { + int i; + unsigned byte; + ui8_t gpsVel_bytes[3]; + int vel24; + double vx, vy, vz, dir; //, alpha; + + for (i = 0; i < 3; i++) { + byte = p_frame[pos_GPSvO + i]; + if (byte > 0xFF) return -1; + gpsVel_bytes[i] = byte; + } + vel24 = gpsVel_bytes[0] << 16 | gpsVel_bytes[1] << 8 | gpsVel_bytes[2]; + if (vel24 > (0x7FFFFF)) vel24 -= 0x1000000; + vx = vel24 / 1e3; // ost + + for (i = 0; i < 3; i++) { + byte = p_frame[pos_GPSvN + i]; + if (byte > 0xFF) return -1; + gpsVel_bytes[i] = byte; + } + vel24 = gpsVel_bytes[0] << 16 | gpsVel_bytes[1] << 8 | gpsVel_bytes[2]; + if (vel24 > (0x7FFFFF)) vel24 -= 0x1000000; + vy= vel24 / 1e3; // nord + + for (i = 0; i < 3; i++) { + byte = p_frame[pos_GPSvV + i]; + if (byte > 0xFF) return -1; + gpsVel_bytes[i] = byte; + } + vel24 = gpsVel_bytes[0] << 16 | gpsVel_bytes[1] << 8 | gpsVel_bytes[2]; + if (vel24 > (0x7FFFFF)) vel24 -= 0x1000000; + vz = vel24 / 1e3; // hoch + + gpx.vE = vx; + gpx.vN = vy; + gpx.vU = vz; + + + gpx.vH = sqrt(vx*vx+vy*vy); +/* + alpha = atan2(vy, vx)*180/M_PI; // ComplexPlane (von x-Achse nach links) - GeoMeteo (von y-Achse nach rechts) + dir = 90-alpha; // z=x+iy= -> i*conj(z)=y+ix=re(i(pi/2-t)), Achsen und Drehsinn vertauscht + if (dir < 0) dir += 360; // atan2(y,x)=atan(y/x)=pi/2-atan(x/y) , atan(1/t) = pi/2 - atan(t) + gpx.vD2 = dir; +*/ + dir = atan2(vx, vy) * 180 / M_PI; + if (dir < 0) dir += 360; + gpx.vD = dir; + + gpx.vV = vz; + + return 0; +} + + +// RS(255,223)-CCSDS +#define rs_N 255 +#define rs_K 223 +#define rs_R (rs_N-rs_K) // 32 +ui8_t rs_cw[rs_N]; + +int lms6_ecc(ui8_t *cw) { + int errors; + ui8_t err_pos[rs_R], + err_val[rs_R]; + + errors = rs_decode(cw, err_pos, err_val); + + return errors; +} + +void print_frame(int crc_err, int len) { + int err=0; + + if (p_frame[0] != 0) + { + //if ((p_frame[pos_SondeSN+1] & 0xF0) == 0x70) // ? beginnen alle SNs mit 0x7A.... bzw 80..... ? + if ( p_frame[pos_SondeSN+1] ) + { + get_FrameNb(); + get_GPStime(); + get_SondeSN(); + if (option_verbose) printf(" (%7d) ", gpx.sn); + printf(" [%5d] ", gpx.frnr); + printf("%s ", weekday[gpx.wday]); + printf("(%02d:%02d:%06.3f) ", gpx.std, gpx.min, gpx.sek); // falls Rundung auf 60s: Ueberlauf + + get_GPSlat(); + get_GPSlon(); + err = get_GPSalt(); + if (!err) { + printf(" lat: %.6f° ", gpx.lat); + printf(" lon: %.6f° ", gpx.lon); + printf(" alt: %.2fm ", gpx.h); + //if (option_verbose) + { + get_GPSvel24(); + //if (option_verbose == 2) printf(" (%.1f ,%.1f,%.1f) ", gpx.vE, gpx.vN, gpx.vU); + printf(" vH: %.1fm/s D: %.1f° vV: %.1fm/s ", gpx.vH, gpx.vD, gpx.vV); + } + } + if (crc_err==0) printf(" [OK]"); else printf(" [NO]"); + + printf("\n"); + } + } +} + +int blk_pos = SYNC_LEN; +int frm_pos = 0; +int sf = 0; + +void proc_frame(int len) { + + char *rawbits = NULL; + int i, j; + int err = 0; + int errs = 0; + int crc_err = 0; + int flen, blen; + + + if ((len % 8) > 4) { + while (len % 8) blk_rawbits[len++] = '0'; + } + //if (len > RAWBITFRAME_LEN+OVERLAP*BITS*2) len = RAWBITFRAME_LEN+OVERLAP*BITS*2; + //for (i = len; i < RAWBITFRAME_LEN+OVERLAP*BITS*2; i++) frame_rawbits[i] = 0; // oder: '0' + blk_rawbits[len] = '\0'; + + flen = len / (2*BITS); + + if (option_vit == 1) { + viterbi(blk_rawbits); + rawbits = vit_rawbits; + } + else if (option_vit == 2) { + viterbi_soft(soft_rawbits, len); + rawbits = vits_rawbits; + } + else rawbits = blk_rawbits; + + err = deconv(rawbits, frame_bits); + + if (err) { for (i=err; i < RAWBITBLOCK_LEN/2; i++) frame_bits[i] = 0; } + + + blen = bits2bytes(frame_bits, block_bytes); + for (j = blen; j < flen; j++) block_bytes[j] = 0; + + + if (option_ecc) { + for (j = 0; j < rs_N; j++) rs_cw[rs_N-1-j] = block_bytes[SYNC_LEN+j]; + errs = lms6_ecc(rs_cw); + for (j = 0; j < rs_N; j++) block_bytes[SYNC_LEN+j] = rs_cw[rs_N-1-j]; + } + + if (option_raw == 2) { + for (i = 0; i < flen; i++) printf("%02x ", block_bytes[i]); + if (option_ecc) printf("(%d)", errs); + printf("\n"); + } + else if (option_raw == 4 && option_ecc) { + for (i = 0; i < rs_N; i++) printf("%02x", block_bytes[SYNC_LEN+i]); + printf(" (%d)", errs); + printf("\n"); + } + else if (option_raw == 8) { + if (option_vit == 1) { + for (i = 0; i < len; i++) printf("%c", vit_rawbits[i]); printf("\n"); + } + else if (option_vit == 2) { + for (i = 0; i < len; i++) printf("%c", vits_rawbits[i]); printf("\n"); + } + else { + for (i = 0; i < len; i++) printf("%c", blk_rawbits[i]); printf("\n"); + } + } + + blk_pos = SYNC_LEN; + + while ( blk_pos-SYNC_LEN < FRM_LEN ) { + + if (sf == 0) { + while ( blk_pos-SYNC_LEN < FRM_LEN ) { + sf = 0; + for (j = 0; j < 4; j++) sf += (block_bytes[blk_pos+j] == frm_sync[j]); + if (sf == 4) { + frm_pos = 0; + break; + } + blk_pos++; + } + } + + if ( sf && frm_pos < FRM_LEN ) { + frame[frm_pos] = block_bytes[blk_pos]; + frm_pos++; + blk_pos++; + } + + if (frm_pos == FRM_LEN) { + + crc_err = check_CRC(p_frame); + + if (option_raw == 1) { + for (i = 0; i < FRM_LEN; i++) printf("%02x ", p_frame[i]); + if (crc_err==0) printf(" [OK]"); else printf(" [NO]"); + printf("\n"); + } + + if (option_raw == 0) print_frame(crc_err, len); + + frm_pos = 0; + sf = 0; + } + + } + +} + + +int main(int argc, char **argv) { + + FILE *fp = NULL; + char *fpname = NULL; + float spb = 0.0; + int header_found = 0; + + int bit, rbit; + int bitpos = 0; + int bitQ; + int pos; + int herrs, herr1; + int headerlen = 0; + + int k,K; + float mv; + unsigned int mv_pos, mv0_pos; + int mp = 0; + + float thres = 0.76; + + int bitofs = 0; + int symlen = 1; + unsigned int bc = 0; + + float sb = 0.0; + float sbit = 0.0; + float level = -1.0, ll = -1.0; + + +#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"); + fprintf(stderr, " -r, --raw\n"); + fprintf(stderr, " --vit,--vit2 (Viterbi/soft)\n"); + fprintf(stderr, " --ecc (Reed-Solomon)\n"); + return 0; + } + else if ( (strcmp(*argv, "-v") == 0) || (strcmp(*argv, "--verbose") == 0) ) { + option_verbose = 1; + } + else if ( (strcmp(*argv, "-r") == 0) || (strcmp(*argv, "--raw") == 0) ) { + option_raw = 1; // bytes - rs_ecc_codewords + } + else if ( (strcmp(*argv, "-r0") == 0) || (strcmp(*argv, "--raw0") == 0) ) { + option_raw = 2; // bytes: sync + codewords + } + else if ( (strcmp(*argv, "-rc") == 0) || (strcmp(*argv, "--rawecc") == 0) ) { + option_raw = 4; // rs_ecc_codewords + } + else if ( (strcmp(*argv, "-R") == 0) || (strcmp(*argv, "--RAW") == 0) ) { + option_raw = 8; // rawbits + } + else if (strcmp(*argv, "--ecc" ) == 0) { option_ecc = 1; } // RS-ECC + else if (strcmp(*argv, "--vit" ) == 0) { option_vit = 1; } // viterbi-hard + else if (strcmp(*argv, "--vit2") == 0) { option_vit = 2; } // viterbi-soft + else if ( (strcmp(*argv, "-i") == 0) || (strcmp(*argv, "--invert") == 0) ) { + option_inv = 1; + } + else if ( (strcmp(*argv, "--dc") == 0) ) { + option_dc = 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); + } + else return -1; + } + else if ( (strcmp(*argv, "--level") == 0) ) { + ++argv; + if (*argv) { + ll = atof(*argv); + } + 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; + + + spb = read_wav_header(fp, (float)BAUD_RATE, wav_channel); + if ( spb < 0 ) { + fclose(fp); + fprintf(stderr, "error: wav header\n"); + return -1; + } + if ( spb < 8 ) { + fprintf(stderr, "note: sample rate low\n"); + } + + + if (option_raw == 4) option_ecc = 1; + + if (option_vit) { + vit_initCodes(); + } + if (option_ecc) { + rs_init_RS255ccsds(); // bch_ecc.c + } + + + symlen = 1; + headerlen = strlen(rawheader); + bitofs = 1; // +1 .. +2 + K = init_buffers(rawheader, headerlen, 2); // shape=2 (alt. shape=1) + if ( K < 0 ) { + fprintf(stderr, "error: init buffers\n"); + return -1; + }; + + level = ll; + k = 0; + mv = -1; mv_pos = 0; + + while ( f32buf_sample(fp, option_inv, 1) != EOF ) { + + k += 1; + if (k >= K-4) { + mv0_pos = mv_pos; + mp = getCorrDFT(-1, K, 0, &mv, &mv_pos); + k = 0; + } + else { + mv = 0.0; + continue; + } + + if (mp > 0 && (mv > thres || mv < -thres)) { + if (mv_pos > mv0_pos) { + + header_found = 0; + herrs = headcmp(symlen, rawheader, headerlen, mv_pos, mv<0, option_dc); // (symlen=1) + herr1 = 0; + + if (herrs <= 3 && herrs > 0) { + herr1 = headcmp(symlen, rawheader, headerlen, mv_pos+1, mv<0, option_dc); + if (herr1 < herrs) { + herrs = herr1; + herr1 = 1; + } + } + if (herrs <= 3) header_found = 1; // herrs <= 3 bitfehler in header + + if (header_found) { + + if (ll <= 0) level = header_level(rawheader, headerlen, mv_pos, mv<0) * 0.6; + + bitpos = 0; + pos = BLOCKSTART; + + if (mv > 0) bc = 0; else bc = 1; + + while ( pos < RAWBITBLOCK_LEN ) { + header_found = !(pos>=RAWBITBLOCK_LEN-10); + //bitQ = read_sbit(fp, symlen, &rbit, option_inv, bitofs, bitpos==0, !header_found); // symlen=1, return: zeroX/bit + bitQ = read_softbit(fp, symlen, &rbit, &sb, level, option_inv, bitofs, bitpos==0, !header_found); // symlen=1, return: zeroX/bit + if (bitQ == EOF) { break; } + + bit = rbit ^ (bc%2); // (c0,inv(c1)) + blk_rawbits[pos] = 0x30 + bit; + + sbit = sb * (-(int)(bc%2)*2+1); + soft_rawbits[pos] = sbit; + + bc++; + pos++; + bitpos += 1; + } + + blk_rawbits[pos] = '\0'; + soft_rawbits[pos] = 0; + + proc_frame(pos); + + if (pos < RAWBITBLOCK_LEN) break; + + pos = BLOCKSTART; + header_found = 0; + } + } + } + + } + + + free_buffers(); + + fclose(fp); + + return 0; +} + diff --git a/demod/m10dm_dft.c b/demod/m10dm_dft.c index 70a842c..039449d 100644 --- a/demod/m10dm_dft.c +++ b/demod/m10dm_dft.c @@ -50,6 +50,8 @@ int option_verbose = 0, // ausfuehrliche Anzeige option_ptu = 0, option_dc = 0, wavloaded = 0; +int wav_channel = 0; // audio channel: left + /* -------------------------------------------------------------------------- */ /* @@ -97,17 +99,23 @@ dduudduudduudduu duduudduuduudduu ddududuudduduudd uduuddududududud uudduduuddu #define HEADLEN 32 // HEADLEN+HEADOFS=32 <= strlen(header) #define HEADOFS 0 // Sync-Header (raw) // Sonde-Header (bits) -//char head[] = "11001100110011001010011001001100"; //"011001001001111100100000"; // M10: 64 9F 20 , M2K2: 64 8F 20 - //"011101101001111100100000"; // M??: 76 9F 20 - //"011001000100100100001001"; // M10-dop: 64 49 09 +//char head[] = "11001100110011001010011001001100"; //"0110010010011111"; // M10: 64 9F , M2K2: 64 8F + //"0111011010011111"; // M10: 76 9F , w/ aux-data + //"0110010001001001"; // M10-dop: 64 49 09 + //"0110010010101111"; // M10+: 64 AF w/ gtop-GPS char rawheader[] = "10011001100110010100110010011001"; -#define FRAME_LEN 102 +#define FRAME_LEN (100+1) // 0x64+1 #define BITFRAME_LEN (FRAME_LEN*BITS) -ui8_t frame_bytes[FRAME_LEN+10]; +#define AUX_LEN 20 +#define BITAUX_LEN (AUX_LEN*BITS) -char frame_bits[BITFRAME_LEN+4]; +ui8_t frame_bytes[FRAME_LEN+AUX_LEN+4]; + +char frame_bits[BITFRAME_LEN+BITAUX_LEN+8]; + +int auxlen = 0; // 0 .. 0x76-0x64 int bits2bytes(char *bitstr, ui8_t *bytes) { @@ -117,7 +125,7 @@ int bits2bytes(char *bitstr, ui8_t *bytes) { bitpos = 0; bytepos = 0; - while (bytepos < FRAME_LEN) { + while (bytepos < FRAME_LEN+AUX_LEN) { byteval = 0; d = 1; @@ -134,24 +142,25 @@ int bits2bytes(char *bitstr, ui8_t *bytes) { } - //while (bytepos < FRAME_LEN) bytes[bytepos++] = 0; + //while (bytepos < FRAME_LEN+AUX_LEN) bytes[bytepos++] = 0; return 0; } /* -------------------------------------------------------------------------- */ +#define stdFLEN 0x64 // pos[0]=0x64 #define pos_GPSTOW 0x0A // 4 byte #define pos_GPSlat 0x0E // 4 byte #define pos_GPSlon 0x12 // 4 byte #define pos_GPSalt 0x16 // 4 byte #define pos_GPSweek 0x20 // 2 byte //Velocity East-North-Up (ENU) -#define pos_GPSvO 0x04 // 2 byte -#define pos_GPSvN 0x06 // 2 byte -#define pos_GPSvV 0x08 // 2 byte -#define pos_SN 0x5D // 2+3 byte -#define pos_Check 0x63 // 2 byte +#define pos_GPSvE 0x04 // 2 byte +#define pos_GPSvN 0x06 // 2 byte +#define pos_GPSvU 0x08 // 2 byte +#define pos_SN 0x5D // 2+3 byte +#define pos_Check (stdFLEN-1) // 2 byte #define ANSI_COLOR_RED "\x1b[31m" @@ -314,7 +323,7 @@ int get_GPSvel() { const double ms2kn100 = 2e2; // m/s -> knots: 1 m/s = 3.6/1.852 kn = 1.94 kn for (i = 0; i < 2; i++) { - byte = frame_bytes[pos_GPSvO + i]; + byte = frame_bytes[pos_GPSvE + i]; gpsVel_bytes[i] = byte; } vel16 = gpsVel_bytes[0] << 8 | gpsVel_bytes[1]; @@ -341,7 +350,7 @@ int get_GPSvel() { datum.vD = dir; for (i = 0; i < 2; i++) { - byte = frame_bytes[pos_GPSvV + i]; + byte = frame_bytes[pos_GPSvU + i]; gpsVel_bytes[i] = byte; } vel16 = gpsVel_bytes[0] << 8 | gpsVel_bytes[1]; @@ -637,8 +646,8 @@ int print_pos(int csOK) { if (option_verbose) { err |= get_GPSvel(); if (!err) { - //if (option_verbose == 2) fprintf(stdout, " "col_GPSvel"(%.1f , %.1f : %.1f°)"col_TXT" ", datum.vx, datum.vy, datum.vD2); - fprintf(stdout, " vH: "col_GPSvel"%.1f"col_TXT" D: "col_GPSvel"%.1f°"col_TXT" vV: "col_GPSvel"%.1f"col_TXT" ", datum.vH, datum.vD, datum.vV); + //if (option_verbose == 2) fprintf(stdout, " "col_GPSvel"(%.1f , %.1f : %.1f)"col_TXT" ", datum.vx, datum.vy, datum.vD2); + fprintf(stdout, " vH: "col_GPSvel"%.1f"col_TXT" D: "col_GPSvel"%.1f"col_TXT"° vV: "col_GPSvel"%.1f"col_TXT" ", datum.vH, datum.vD, datum.vV); } if (option_verbose >= 2) { get_SN(); @@ -705,26 +714,33 @@ int print_frame(int pos) { int i; ui8_t byte; int cs1, cs2; + int flen = stdFLEN; // stdFLEN=0x64, auxFLEN=0x76 bits2bytes(frame_bits, frame_bytes); + flen = frame_bytes[0]; + if (flen == stdFLEN) auxlen = 0; + else { + auxlen = flen - stdFLEN; + if (auxlen < 0 || auxlen > AUX_LEN) auxlen = 0; + } - cs1 = (frame_bytes[pos_Check] << 8) | frame_bytes[pos_Check+1]; - cs2 = checkM10(frame_bytes, pos_Check); + cs1 = (frame_bytes[pos_Check+auxlen] << 8) | frame_bytes[pos_Check+auxlen+1]; + cs2 = checkM10(frame_bytes, pos_Check+auxlen); if (option_raw) { if (option_color && frame_bytes[1] != 0x49) { fprintf(stdout, col_FRTXT); - for (i = 0; i < FRAME_LEN-1; i++) { + for (i = 0; i < FRAME_LEN+auxlen; i++) { byte = frame_bytes[i]; if ((i >= pos_GPSTOW) && (i < pos_GPSTOW+4)) fprintf(stdout, col_GPSTOW); if ((i >= pos_GPSlat) && (i < pos_GPSlat+4)) fprintf(stdout, col_GPSlat); if ((i >= pos_GPSlon) && (i < pos_GPSlon+4)) fprintf(stdout, col_GPSlon); if ((i >= pos_GPSalt) && (i < pos_GPSalt+4)) fprintf(stdout, col_GPSalt); if ((i >= pos_GPSweek) && (i < pos_GPSweek+2)) fprintf(stdout, col_GPSweek); - if ((i >= pos_GPSvO) && (i < pos_GPSvO+6)) fprintf(stdout, col_GPSvel); + if ((i >= pos_GPSvE) && (i < pos_GPSvE+6)) fprintf(stdout, col_GPSvel); if ((i >= pos_SN) && (i < pos_SN+5)) fprintf(stdout, col_SN); - if ((i >= pos_Check) && (i < pos_Check+2)) fprintf(stdout, col_Check); + if ((i >= pos_Check+auxlen) && (i < pos_Check+auxlen+2)) fprintf(stdout, col_Check); fprintf(stdout, "%02x", byte); fprintf(stdout, col_FRTXT); } @@ -736,7 +752,7 @@ int print_frame(int pos) { fprintf(stdout, ANSI_COLOR_RESET"\n"); } else { - for (i = 0; i < FRAME_LEN-1; i++) { + for (i = 0; i < FRAME_LEN+auxlen; i++) { byte = frame_bytes[i]; fprintf(stdout, "%02x", byte); } @@ -750,7 +766,7 @@ int print_frame(int pos) { } else if (frame_bytes[1] == 0x49) { if (option_verbose == 3) { - for (i = 0; i < FRAME_LEN-1; i++) { + for (i = 0; i < FRAME_LEN+auxlen; i++) { byte = frame_bytes[i]; fprintf(stdout, "%02x", byte); } @@ -777,7 +793,7 @@ int main(int argc, char **argv) { int herrs, herr1; int headerlen = 0; - int k,K; + int k, K; float mv; unsigned int mv_pos, mv0_pos; int mp = 0; @@ -825,7 +841,8 @@ int main(int argc, char **argv) { else if ( (strcmp(*argv, "--dc") == 0) ) { option_dc = 1; } - else if (strcmp(*argv, "--ths") == 0) { + 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); @@ -845,7 +862,7 @@ int main(int argc, char **argv) { if (!wavloaded) fp = stdin; - spb = read_wav_header(fp, (float)BAUD_RATE); + spb = read_wav_header(fp, (float)BAUD_RATE, wav_channel); if ( spb < 0 ) { fclose(fp); fprintf(stderr, "error: wav header\n"); @@ -858,7 +875,7 @@ int main(int argc, char **argv) { symlen = 2; headerlen = strlen(rawheader); - bitofs = 1; // +1 .. +2 + bitofs = 0; // 0 .. +2 K = init_buffers(rawheader, headerlen, 1); // shape=0 (alt. shape=1) if ( K < 0 ) { fprintf(stderr, "error: init buffers\n"); @@ -886,10 +903,11 @@ int main(int argc, char **argv) { if (mv_pos > mv0_pos) { header_found = 0; - herrs = headcmp(1, rawheader, headerlen, mv_pos, mv<0, option_dc); // (symlen=2) + herrs = headcmp(1, rawheader, headerlen, mv_pos, mv<0, option_dc); // header nicht manchester! herr1 = 0; if (herrs <= 3 && herrs > 0) { - herr1 = headcmp(symlen, rawheader, headerlen, mv_pos+1, mv<0, option_dc); + herr1 = headcmp(1, rawheader, headerlen, mv_pos+1, mv<0, option_dc); + //int herr2 = headcmp(1, rawheader, headerlen, mv_pos-1, mv<0, option_dc); if (herr1 < herrs) { herrs = herr1; herr1 = 1; @@ -904,7 +922,7 @@ int main(int argc, char **argv) { pos /= 2; bit0 = '0'; - while ( pos < BITFRAME_LEN ) { + while ( pos < BITFRAME_LEN+BITAUX_LEN ) { header_found = !(pos>=BITFRAME_LEN-10); bitQ = read_sbit(fp, symlen, &bit, option_inv, bitofs, bitpos==0, !header_found); // symlen=2, return: zeroX/bit if (bitQ == EOF) { break; } diff --git a/demod/rs41dm_dft.c b/demod/rs41dm_dft.c index 1fb24f5..7f42e5f 100644 --- a/demod/rs41dm_dft.c +++ b/demod/rs41dm_dft.c @@ -69,7 +69,9 @@ int option_verbose = 0, // ausfuehrliche Anzeige option_sat = 0, // GPS sat data option_ptu = 0, option_ths = 0, + option_json = 0, // JSON output (auto_rx) wavloaded = 0; +int wav_channel = 0; // audio channel: left #define BITS 8 @@ -454,7 +456,7 @@ int get_FrameConf() { if (crc == 0) { calfr = framebyte(pos_CalData); if (calfrchk[calfr] == 0) // const? - { + { // 0x32 not constant for (i = 0; i < 16; i++) { calibytes[calfr*16 + i] = framebyte(pos_CalData+1+i); } @@ -842,8 +844,12 @@ int get_Calconf(int out) { if (calfr == 0x02 && option_verbose /*== 2*/) { byte = framebyte(pos_Calburst); - burst = byte; // fw >= 0x4ef5, BK irrelevant? (killtimer in 0x31?) + burst = byte; // fw >= 0x4ef5, BK irrelevant? (burst-killtimer in 0x31?) fprintf(stdout, ": BK %02X ", burst); + if (option_verbose == 3) { // killtimer + int kt = frame[0x5A] + (frame[0x5B] << 8); // short? + if ( kt != 0xFFFF ) fprintf(stdout, ": kt 0x%04x = %dsec = %.1fmin ", kt, kt, kt/60.0); + } } if (calfr == 0x00 && option_verbose) { @@ -855,6 +861,12 @@ int get_Calconf(int out) { fprintf(stdout, ": fq %d ", freq); } + if (calfr == 0x31 && option_verbose == 3) { + int bt = frame[0x59] + (frame[0x5A] << 8); // short? + // fw >= 0x4ef5: default=[88 77]=0x7788sec=510min + if ( bt != 0x0000 ) fprintf(stdout, ": bt 0x%04x = %dsec = %.1fmin ", bt, bt, bt/60.0); + } + if (calfr == 0x21 && option_verbose /*== 2*/) { // eventuell noch zwei bytes in 0x22 for (i = 0; i < 9; i++) sondetyp[i] = 0; for (i = 0; i < 8; i++) { @@ -974,7 +986,12 @@ int rs41_ecc(int frmlen) { ret = errors1 + errors2; - if (errors1 < 0 || errors2 < 0) ret = -1; + if (errors1 < 0 || errors2 < 0) { + ret = 0; + if (errors1 < 0) ret |= 0x1; + if (errors2 < 0) ret |= 0x2; + ret = -ret; + } return ret; } @@ -1020,6 +1037,10 @@ 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 ); @@ -1037,10 +1058,34 @@ int print_position(int ec) { //if (output) { if (option_crc) { - fprintf(stdout, " # ["); - for (i=0; i<5; i++) fprintf(stdout, "%d", (gpx.crc>>i)&1); - fprintf(stdout, "]"); - if (option_ecc == 2 && ec > 0) fprintf(stdout, " (%d)", ec); + 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, " (--)"); + } + } } } @@ -1159,6 +1204,7 @@ int main(int argc, char *argv[]) { fprintf(stderr, " --ecc (Reed-Solomon)\n"); fprintf(stderr, " --std (std framelen)\n"); fprintf(stderr, " --ths (peak threshold; default=%.1f)\n", thres); + fprintf(stderr, " --json (JSON output)\n"); return 0; } else if ( (strcmp(*argv, "-v") == 0) || (strcmp(*argv, "--verbose") == 0) ) { @@ -1181,6 +1227,8 @@ 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, "--ch2") == 0) { wav_channel = 1; } // right channel (default: 0=left) else if (strcmp(*argv, "--ths") == 0) { ++argv; if (*argv) { @@ -1201,7 +1249,7 @@ int main(int argc, char *argv[]) { if (!wavloaded) fp = stdin; - spb = read_wav_header(fp, (float)BAUD_RATE); + spb = read_wav_header(fp, (float)BAUD_RATE, wav_channel); if ( spb < 0 ) { fclose(fp); fprintf(stderr, "error: wav header\n"); diff --git a/demod/rs41dm_iq.c b/demod/rs41dm_iq.c new file mode 100644 index 0000000..f232dfc --- /dev/null +++ b/demod/rs41dm_iq.c @@ -0,0 +1,1390 @@ + +/* + * rs41 + * sync header: correlation/matched filter + * files: rs41dm_iq.c bch_ecc.c demod_iq.h demod_iq.c + * compile: + * gcc -c demod_iq.c + * gcc rs41dm_iq.c demod_iq.o -lm -o rs41dm_iq + * + * author: zilog80 + */ + +#include +#include +#include +#include + +#ifdef CYGWIN + #include // cygwin: _setmode() + #include +#endif + +typedef unsigned char ui8_t; +typedef unsigned short ui16_t; +typedef unsigned int ui32_t; +typedef short i16_t; +typedef int i32_t; + +//#include "demod_dft.c" +#include "demod_iq.h" + +#include "bch_ecc.c" // RS/ecc/ + + +typedef struct { + int typ; + int msglen; + int msgpos; + int parpos; + int hdrlen; + int frmlen; +} rscfg_t; + +rscfg_t cfg_rs41 = { 41, (320-56)/2, 56, 8, 8, 320}; + + +typedef struct { + int frnr; + char id[9]; + int week; int gpssec; + int jahr; int monat; int tag; + int wday; + int std; int min; float sek; + double lat; double lon; double alt; + double vN; double vE; double vU; + double vH; double vD; double vD2; + float T; + ui32_t crc; +} gpx_t; + +gpx_t gpx; + +int option_verbose = 0, // ausfuehrliche Anzeige + option_raw = 0, // rohe Frames + option_inv = 0, // invertiert Signal + option_res = 0, // genauere Bitmessung + option_crc = 0, // check CRC + option_ecc = 0, // Reed-Solomon ECC + option_sat = 0, // GPS sat data + option_ptu = 0, + option_ths = 0, + option_iq = 0, + option_ofs = 0, + option_dbg = 0, + wavloaded = 0; +int wav_channel = 0; // audio channel: left + + +#define BITS 8 +#define HEADLEN 64 +#define FRAMESTART ((HEADLEN)/BITS) + +/* 10 B6 CA 11 22 96 12 F8 */ +char header[] = "0000100001101101010100111000100001000100011010010100100000011111"; + + +#define NDATA_LEN 320 // std framelen 320 +#define XDATA_LEN 198 +#define FRAME_LEN (NDATA_LEN+XDATA_LEN) // max framelen 518 +ui8_t //xframe[FRAME_LEN] = { 0x10, 0xB6, 0xCA, 0x11, 0x22, 0x96, 0x12, 0xF8}, = xorbyte( frame) + frame[FRAME_LEN] = { 0x86, 0x35, 0xf4, 0x40, 0x93, 0xdf, 0x1a, 0x60}; // = xorbyte(xframe) + +float byteQ[FRAME_LEN]; + +#define MASK_LEN 64 +ui8_t mask[MASK_LEN] = { 0x96, 0x83, 0x3E, 0x51, 0xB1, 0x49, 0x08, 0x98, + 0x32, 0x05, 0x59, 0x0E, 0xF9, 0x44, 0xC6, 0x26, + 0x21, 0x60, 0xC2, 0xEA, 0x79, 0x5D, 0x6D, 0xA1, + 0x54, 0x69, 0x47, 0x0C, 0xDC, 0xE8, 0x5C, 0xF1, + 0xF7, 0x76, 0x82, 0x7F, 0x07, 0x99, 0xA2, 0x2C, + 0x93, 0x7C, 0x30, 0x63, 0xF5, 0x10, 0x2E, 0x61, + 0xD0, 0xBC, 0xB4, 0xB6, 0x06, 0xAA, 0xF4, 0x23, + 0x78, 0x6E, 0x3B, 0xAE, 0xBF, 0x7B, 0x4C, 0xC1}; +/* LFSR: ab i=8 (mod 64): + * m[16+i] = m[i] ^ m[i+2] ^ m[i+4] ^ m[i+6] + * ________________3205590EF944C6262160C2EA795D6DA15469470CDCE85CF1 + * F776827F0799A22C937C3063F5102E61D0BCB4B606AAF423786E3BAEBF7B4CC196833E51B1490898 + */ + +/* ------------------------------------------------------------------------------------ */ + +#define BAUD_RATE 4800 + +/* ------------------------------------------------------------------------------------ */ + + +int bits2byte(char bits[]) { + int i, byteval=0, d=1; + for (i = 0; i < 8; i++) { // little endian + /* for (i = 7; i >= 0; i--) { // big endian */ + if (bits[i] == 1) byteval += d; + else if (bits[i] == 0) byteval += 0; + else return 0x100; + d <<= 1; + } + return byteval; +} + + +/* +ui8_t xorbyte(int pos) { + return xframe[pos] ^ mask[pos % MASK_LEN]; +} +*/ +ui8_t framebyte(int pos) { + return frame[pos]; +} + + +/* ------------------------------------------------------------------------------------ */ +/* + * Convert GPS Week and Seconds to Modified Julian Day. + * - Adapted from sci.astro FAQ. + * - Ignores UTC leap seconds. + */ +void Gps2Date(long GpsWeek, long GpsSeconds, int *Year, int *Month, int *Day) { + + long GpsDays, Mjd; + long J, C, Y, M; + + GpsDays = GpsWeek * 7 + (GpsSeconds / 86400); + Mjd = 44244 + GpsDays; + + J = Mjd + 2468570; + C = 4 * J / 146097; + J = J - (146097 * C + 3) / 4; + Y = 4000 * (J + 1) / 1461001; + J = J - 1461 * Y / 4 + 31; + M = 80 * J / 2447; + *Day = J - 2447 * M / 80; + J = M / 11; + *Month = M + 2 - (12 * J); + *Year = 100 * (C - 49) + Y + J; +} +/* ------------------------------------------------------------------------------------ */ + +ui32_t u4(ui8_t *bytes) { // 32bit unsigned int + ui32_t val = 0; + memcpy(&val, bytes, 4); + return val; +} + +ui32_t u3(ui8_t *bytes) { // 24bit unsigned int + int val24 = 0; + val24 = bytes[0] | (bytes[1]<<8) | (bytes[2]<<16); + // = memcpy(&val, bytes, 3), val &= 0x00FFFFFF; + return val24; +} + +int i3(ui8_t *bytes) { // 24bit signed int + int val = 0, + val24 = 0; + val = bytes[0] | (bytes[1]<<8) | (bytes[2]<<16); + val24 = val & 0xFFFFFF; if (val24 & 0x800000) val24 = val24 - 0x1000000; + return val24; +} + +ui32_t u2(ui8_t *bytes) { // 16bit unsigned int + return bytes[0] | (bytes[1]<<8); +} + +/* +double r8(ui8_t *bytes) { + double val = 0; + memcpy(&val, bytes, 8); + return val; +} + +float r4(ui8_t *bytes) { + float val = 0; + memcpy(&val, bytes, 4); + return val; +} +*/ + +/* +int crc16x(int start, int len) { + int crc16poly = 0x1021; + int rem = 0xFFFF, i, j; + int xbyte; + + if (start+len+2 > FRAME_LEN) return -1; + + for (i = 0; i < len; i++) { + xbyte = xorbyte(start+i); + rem = rem ^ (xbyte << 8); + for (j = 0; j < 8; j++) { + if (rem & 0x8000) { + rem = (rem << 1) ^ crc16poly; + } + else { + rem = (rem << 1); + } + rem &= 0xFFFF; + } + } + return rem; +} +*/ +int crc16(int start, int len) { + int crc16poly = 0x1021; + int rem = 0xFFFF, i, j; + int byte; + + if (start+len+2 > FRAME_LEN) return -1; + + for (i = 0; i < len; i++) { + byte = framebyte(start+i); + rem = rem ^ (byte << 8); + for (j = 0; j < 8; j++) { + if (rem & 0x8000) { + rem = (rem << 1) ^ crc16poly; + } + else { + rem = (rem << 1); + } + rem &= 0xFFFF; + } + } + return rem; +} + +int check_CRC(ui32_t pos, ui32_t pck) { + ui32_t crclen = 0, + crcdat = 0; + if (((pck>>8) & 0xFF) != frame[pos]) return -1; + crclen = frame[pos+1]; + if (pos + crclen + 4 > FRAME_LEN) return -1; + crcdat = u2(frame+pos+2+crclen); + if ( crcdat != crc16(pos+2, crclen) ) { + return 1; // CRC NO + } + else return 0; // CRC OK +} + + +/* + Pos: SubHeader, 1+1 byte (ID+LEN) +0x039: 7928 FrameNumber+SondeID + +(0x050: 0732 CalFrames 0x00..0x32) +0x065: 7A2A PTU +0x093: 7C1E GPS1: RXM-RAW (0x02 0x10) Week, TOW, Sats +0x0B5: 7D59 GPS2: RXM-RAW (0x02 0x10) pseudorange, doppler +0x112: 7B15 GPS3: NAV-SOL (0x01 0x06) ECEF-POS, ECEF-VEL +0x12B: 7611 00 +0x12B: 7Exx AUX-xdata +*/ + +#define crc_FRAME (1<<0) +#define xor_FRAME 0x1713 // ^0x6E3B=0x7928 +#define pck_FRAME 0x7928 +#define pos_FRAME 0x039 +#define pos_FrameNb 0x03B // 2 byte +#define pos_SondeID 0x03D // 8 byte +#define pos_CalData 0x052 // 1 byte, counter 0x00..0x32 +#define pos_Calfreq 0x055 // 2 byte, calfr 0x00 +#define pos_Calburst 0x05E // 1 byte, calfr 0x02 +// ? #define pos_Caltimer 0x05A // 2 byte, calfr 0x02 ? +#define pos_CalRSTyp 0x05B // 8 byte, calfr 0x21 (+2 byte in 0x22?) + // weitere chars in calfr 0x22/0x23; weitere ID + +#define crc_PTU (1<<1) +#define xor_PTU 0xE388 // ^0x99A2=0x0x7A2A +#define pck_PTU 0x7A2A // PTU +#define pos_PTU 0x065 + +#define crc_GPS1 (1<<2) +#define xor_GPS1 0x9667 // ^0xEA79=0x7C1E +#define pck_GPS1 0x7C1E // RXM-RAW (0x02 0x10) +#define pos_GPS1 0x093 +#define pos_GPSweek 0x095 // 2 byte +#define pos_GPSiTOW 0x097 // 4 byte +#define pos_satsN 0x09B // 12x2 byte (1: SV, 1: quality,strength) + +#define crc_GPS2 (1<<3) +#define xor_GPS2 0xD7AD // ^0xAAF4=0x7D59 +#define pck_GPS2 0x7D59 // RXM-RAW (0x02 0x10) +#define pos_GPS2 0x0B5 +#define pos_minPR 0x0B7 // 4 byte +#define pos_FF 0x0BB // 1 byte +#define pos_dataSats 0x0BC // 12x(4+3) byte (4: pseudorange, 3: doppler) + +#define crc_GPS3 (1<<4) +#define xor_GPS3 0xB9FF // ^0xC2EA=0x7B15 +#define pck_GPS3 0x7B15 // NAV-SOL (0x01 0x06) +#define pos_GPS3 0x112 +#define pos_GPSecefX 0x114 // 4 byte +#define pos_GPSecefY 0x118 // 4 byte +#define pos_GPSecefZ 0x11C // 4 byte +#define pos_GPSecefV 0x120 // 3*2 byte +#define pos_numSats 0x126 // 1 byte +#define pos_sAcc 0x127 // 1 byte +#define pos_pDOP 0x128 // 1 byte + +#define crc_AUX (1<<5) +#define pck_AUX 0x7E00 // LEN variable +#define pos_AUX 0x12B + +#define crc_ZERO (1<<6) // LEN variable +#define pck_ZERO 0x7600 + + +ui8_t calibytes[51*16]; +ui8_t calfrchk[51]; +float Rf1, // ref-resistor f1 (750 Ohm) + Rf2, // ref-resistor f2 (1100 Ohm) + co1[3], // { -243.911 , 0.187654 , 8.2e-06 } + calT1[3], // calibration T1 + co2[3], // { -243.911 , 0.187654 , 8.2e-06 } + calT2[3]; // calibration T2-Hum + + +double c = 299.792458e6; +double L1 = 1575.42e6; + +int get_SatData() { + int i, n; + int sv; + ui32_t minPR; + int Nfix; + double pDOP, sAcc; + + fprintf(stdout, "[%d]\n", u2(frame+pos_FrameNb)); + + fprintf(stdout, "iTOW: 0x%08X", u4(frame+pos_GPSiTOW)); + fprintf(stdout, " week: 0x%04X", u2(frame+pos_GPSweek)); + fprintf(stdout, "\n"); + minPR = u4(frame+pos_minPR); + fprintf(stdout, "minPR: %d", minPR); + fprintf(stdout, "\n"); + + for (i = 0; i < 12; i++) { + n = i*7; + sv = frame[pos_satsN+2*i]; + if (sv == 0xFF) break; + fprintf(stdout, " SV: %2d # ", sv); + fprintf(stdout, "prMes: %.1f", u4(frame+pos_dataSats+n)/100.0 + minPR); + fprintf(stdout, " "); + fprintf(stdout, "doMes: %.1f", -i3(frame+pos_dataSats+n+4)/100.0*L1/c); + fprintf(stdout, "\n"); + } + + fprintf(stdout, "ECEF-POS: (%d,%d,%d)\n", + (i32_t)u4(frame+pos_GPSecefX), + (i32_t)u4(frame+pos_GPSecefY), + (i32_t)u4(frame+pos_GPSecefZ)); + fprintf(stdout, "ECEF-VEL: (%d,%d,%d)\n", + (i16_t)u2(frame+pos_GPSecefV+0), + (i16_t)u2(frame+pos_GPSecefV+2), + (i16_t)u2(frame+pos_GPSecefV+4)); + + Nfix = frame[pos_numSats]; + sAcc = frame[pos_sAcc]/10.0; + pDOP = frame[pos_pDOP]/10.0; + fprintf(stdout, "numSatsFix: %2d sAcc: %.1f pDOP: %.1f\n", Nfix, sAcc, pDOP); + + + fprintf(stdout, "CRC: "); + fprintf(stdout, " %04X", pck_GPS1); + if (check_CRC(pos_GPS1, pck_GPS1)==0) fprintf(stdout, "[OK]"); else fprintf(stdout, "[NO]"); + //fprintf(stdout, "[%+d]", check_CRC(pos_GPS1, pck_GPS1)); + fprintf(stdout, " %04X", pck_GPS2); + if (check_CRC(pos_GPS2, pck_GPS2)==0) fprintf(stdout, "[OK]"); else fprintf(stdout, "[NO]"); + //fprintf(stdout, "[%+d]", check_CRC(pos_GPS2, pck_GPS2)); + fprintf(stdout, " %04X", pck_GPS3); + if (check_CRC(pos_GPS3, pck_GPS3)==0) fprintf(stdout, "[OK]"); else fprintf(stdout, "[NO]"); + //fprintf(stdout, "[%+d]", check_CRC(pos_GPS3, pck_GPS3)); + + fprintf(stdout, "\n"); + fprintf(stdout, "\n"); + + return 0; +} + + +int get_FrameNb() { + int i; + unsigned byte; + ui8_t frnr_bytes[2]; + int frnr; + + for (i = 0; i < 2; i++) { + byte = framebyte(pos_FrameNb + i); + frnr_bytes[i] = byte; + } + + frnr = frnr_bytes[0] + (frnr_bytes[1] << 8); + gpx.frnr = frnr; + + return 0; +} + +int get_SondeID(int crc) { + int i; + unsigned byte; + char sondeid_bytes[9]; + + if (crc == 0) { + for (i = 0; i < 8; i++) { + byte = framebyte(pos_SondeID + i); + //if ((byte < 0x20) || (byte > 0x7E)) return -1; + sondeid_bytes[i] = byte; + } + sondeid_bytes[8] = '\0'; + if ( strncmp(gpx.id, sondeid_bytes, 8) != 0 ) { + //for (i = 0; i < 51; i++) calfrchk[i] = 0; + memset(calfrchk, 0, 51); + memcpy(gpx.id, sondeid_bytes, 8); + gpx.id[8] = '\0'; + } + } + + return 0; +} + +int get_FrameConf() { + int crc, err; + ui8_t calfr; + int i; + + crc = check_CRC(pos_FRAME, pck_FRAME); + if (crc) gpx.crc |= crc_FRAME; + + err = crc; + err |= get_FrameNb(); + err |= get_SondeID(crc); + + if (crc == 0) { + calfr = framebyte(pos_CalData); + if (calfrchk[calfr] == 0) // const? + { // 0x32 not constant + for (i = 0; i < 16; i++) { + calibytes[calfr*16 + i] = framebyte(pos_CalData+1+i); + } + calfrchk[calfr] = 1; + } + } + + return err; +} + +int get_CalData() { + + memcpy(&Rf1, calibytes+61, 4); // 0x03*0x10+13 + memcpy(&Rf2, calibytes+65, 4); // 0x04*0x10+ 1 + + memcpy(co1+0, calibytes+77, 4); // 0x04*0x10+13 + memcpy(co1+1, calibytes+81, 4); // 0x05*0x10+ 1 + memcpy(co1+2, calibytes+85, 4); // 0x05*0x10+ 5 + + memcpy(calT1+0, calibytes+89, 4); // 0x05*0x10+ 9 + memcpy(calT1+1, calibytes+93, 4); // 0x05*0x10+13 + memcpy(calT1+2, calibytes+97, 4); // 0x06*0x10+ 1 + + memcpy(co2+0, calibytes+293, 4); // 0x12*0x10+ 5 + memcpy(co2+1, calibytes+297, 4); // 0x12*0x10+ 9 + memcpy(co2+2, calibytes+301, 4); // 0x12*0x10+13 + + memcpy(calT2+0, calibytes+305, 4); // 0x13*0x10+ 1 + memcpy(calT2+1, calibytes+309, 4); // 0x13*0x10+ 5 + memcpy(calT2+2, calibytes+313, 4); // 0x13*0x10+ 9 + + return 0; +} + +float get_Tc0(ui32_t f, ui32_t f1, ui32_t f2) { + // y = (f - f1) / (f2 - f1); + // y1 = (f - f1) / f2; // = (1 - f1/f2)*y + float a = 3.9083e-3, // Pt1000 platinum resistance + b = -5.775e-7, + c = -4.183e-12; // below 0C, else C=0 + float *cal = calT1; + float Rb = (f1*Rf2-f2*Rf1)/(f2-f1), // ofs + Ra = f * (Rf2-Rf1)/(f2-f1) - Rb, + raw = Ra/1000.0, + g_r = 0.8024*cal[0] + 0.0176, // empirisch + r_o = 0.0705*cal[1] + 0.0011, // empirisch + r = raw * g_r + r_o, + t = (-a + sqrt(a*a + 4*b*(r-1)))/(2*b); // t>0: c=0 + // R/R0 = 1 + at + bt^2 + c(t-100)t^3 , R0 = 1000 Ohm, t/Celsius + return t; +} +float get_Tc(ui32_t f, ui32_t f1, ui32_t f2) { + float *p = co1; + float *c = calT1; + float g = (float)(f2-f1)/(Rf2-Rf1), // gain + Rb = (f1*Rf2-f2*Rf1)/(float)(f2-f1), // ofs + Rc = f/g - Rb, + //R = (Rc + c[1]) * c[0], + //T = p[0] + p[1]*R + p[2]*R*R; + R = Rc * c[0], + T = (p[0] + p[1]*R + p[2]*R*R + c[1])*(1.0 + c[2]); + return T; +} + +int get_PTU() { + int err=0, i; + int bR, bc1, bT1, + bc2, bT2; + ui32_t meas[12]; + float Tc = -273.15; + float Tc0 = -273.15; + + get_CalData(); + + err = check_CRC(pos_PTU, pck_PTU); + if (err) gpx.crc |= crc_PTU; + + if (err == 0) + { + + for (i = 0; i < 12; i++) { + meas[i] = u3(frame+pos_PTU+2+3*i); + } + + bR = calfrchk[0x03] && calfrchk[0x04]; + bc1 = calfrchk[0x04] && calfrchk[0x05]; + bT1 = calfrchk[0x05] && calfrchk[0x06]; + bc2 = calfrchk[0x12] && calfrchk[0x13]; + bT2 = calfrchk[0x13]; + + if (bR && bc1 && bT1) { + Tc = get_Tc(meas[0], meas[1], meas[2]); + Tc0 = get_Tc0(meas[0], meas[1], meas[2]); + } + gpx.T = Tc; + + if (option_verbose == 4) + { + printf(" h: %8.2f # ", gpx.alt); // crc_GPS3 ? + + printf("1: %8d %8d %8d", meas[0], meas[1], meas[2]); + printf(" # "); + printf("2: %8d %8d %8d", meas[3], meas[4], meas[5]); + printf(" # "); + printf("3: %8d %8d %8d", meas[6], meas[7], meas[8]); + printf(" # "); + if (Tc > -273.0) { + printf(" T: %8.4f , T0: %8.4f ", Tc, Tc0); + } + printf("\n"); + + if (gpx.alt > -100.0) { + printf(" %9.2f ; %6.1f ; %6.1f ", gpx.alt, Rf1, Rf2); + printf("; %10.6f ; %10.6f ; %10.6f ;", calT1[0], calT1[1], calT1[2]); + printf(" %8d ; %8d ; %8d ", meas[0], meas[1], meas[2]); + printf("; %10.6f ; %10.6f ; %10.6f ;", calT2[0], calT2[1], calT2[2]); + printf(" %8d ; %8d ; %8d" , meas[6], meas[7], meas[8]); + printf("\n"); + } + } + + } + + return err; +} + +int get_GPSweek() { + int i; + unsigned byte; + ui8_t gpsweek_bytes[2]; + int gpsweek; + + for (i = 0; i < 2; i++) { + byte = framebyte(pos_GPSweek + i); + gpsweek_bytes[i] = byte; + } + + gpsweek = gpsweek_bytes[0] + (gpsweek_bytes[1] << 8); + //if (gpsweek < 0) { gpx.week = -1; return -1; } // (short int) + gpx.week = gpsweek; + + return 0; +} + +char weekday[7][3] = { "So", "Mo", "Di", "Mi", "Do", "Fr", "Sa"}; + +int get_GPStime() { + int i; + unsigned byte; + ui8_t gpstime_bytes[4]; + int gpstime = 0, // 32bit + day; + int ms; + + for (i = 0; i < 4; i++) { + byte = framebyte(pos_GPSiTOW + i); + gpstime_bytes[i] = byte; + } + + memcpy(&gpstime, gpstime_bytes, 4); + ms = gpstime % 1000; + gpstime /= 1000; + + gpx.gpssec = gpstime; + + day = (gpstime / (24 * 3600)) % 7; + //if ((day < 0) || (day > 6)) return -1; // besser CRC-check + + gpstime %= (24*3600); + + gpx.wday = day; + gpx.std = gpstime / 3600; + gpx.min = (gpstime % 3600) / 60; + gpx.sek = gpstime % 60 + ms/1000.0; + + return 0; +} + +int get_GPS1() { + int err=0; + + // ((framebyte(pos_GPS1)<<8) | framebyte(pos_GPS1+1)) != pck_GPS1 ? + if ( framebyte(pos_GPS1) != ((pck_GPS1>>8) & 0xFF) ) { + gpx.crc |= crc_GPS1; + return -1; + } + + err = check_CRC(pos_GPS1, pck_GPS1); + if (err) gpx.crc |= crc_GPS1; + + //err = 0; + err |= get_GPSweek(); + err |= get_GPStime(); + + return err; +} + +int get_GPS2() { + int err=0; + + err = check_CRC(pos_GPS2, pck_GPS2); + if (err) gpx.crc |= crc_GPS2; + + return err; +} + +#define EARTH_a 6378137.0 +#define EARTH_b 6356752.31424518 +#define EARTH_a2_b2 (EARTH_a*EARTH_a - EARTH_b*EARTH_b) + +double a = EARTH_a, + b = EARTH_b, + a_b = EARTH_a2_b2, + e2 = EARTH_a2_b2 / (EARTH_a*EARTH_a), + ee2 = EARTH_a2_b2 / (EARTH_b*EARTH_b); + +void ecef2elli(double X[], double *lat, double *lon, double *alt) { + double phi, lam, R, p, t; + + lam = atan2( X[1] , X[0] ); + + p = sqrt( X[0]*X[0] + X[1]*X[1] ); + t = atan2( X[2]*a , p*b ); + + phi = atan2( X[2] + ee2 * b * sin(t)*sin(t)*sin(t) , + p - e2 * a * cos(t)*cos(t)*cos(t) ); + + R = a / sqrt( 1 - e2*sin(phi)*sin(phi) ); + *alt = p / cos(phi) - R; + + *lat = phi*180/M_PI; + *lon = lam*180/M_PI; +} + +int get_GPSkoord() { + int i, k; + unsigned byte; + ui8_t XYZ_bytes[4]; + int XYZ; // 32bit + double X[3], lat, lon, alt; + ui8_t gpsVel_bytes[2]; + short vel16; // 16bit + double V[3], phi, lam, dir; + + + for (k = 0; k < 3; k++) { + + for (i = 0; i < 4; i++) { + byte = frame[pos_GPSecefX + 4*k + i]; + XYZ_bytes[i] = byte; + } + memcpy(&XYZ, XYZ_bytes, 4); + X[k] = XYZ / 100.0; + + for (i = 0; i < 2; i++) { + byte = frame[pos_GPSecefV + 2*k + i]; + gpsVel_bytes[i] = byte; + } + vel16 = gpsVel_bytes[0] | gpsVel_bytes[1] << 8; + V[k] = vel16 / 100.0; + + } + + + // ECEF-Position + ecef2elli(X, &lat, &lon, &alt); + gpx.lat = lat; + gpx.lon = lon; + gpx.alt = alt; + if ((alt < -1000) || (alt > 80000)) return -3; + + + // ECEF-Velocities + // ECEF-Vel -> NorthEastUp + phi = lat*M_PI/180.0; + lam = lon*M_PI/180.0; + gpx.vN = -V[0]*sin(phi)*cos(lam) - V[1]*sin(phi)*sin(lam) + V[2]*cos(phi); + gpx.vE = -V[0]*sin(lam) + V[1]*cos(lam); + gpx.vU = V[0]*cos(phi)*cos(lam) + V[1]*cos(phi)*sin(lam) + V[2]*sin(phi); + + // NEU -> HorDirVer + gpx.vH = sqrt(gpx.vN*gpx.vN+gpx.vE*gpx.vE); +/* + double alpha; + alpha = atan2(gpx.vN, gpx.vE)*180/M_PI; // ComplexPlane (von x-Achse nach links) - GeoMeteo (von y-Achse nach rechts) + dir = 90-alpha; // z=x+iy= -> i*conj(z)=y+ix=re(i(pi/2-t)), Achsen und Drehsinn vertauscht + if (dir < 0) dir += 360; // atan2(y,x)=atan(y/x)=pi/2-atan(x/y) , atan(1/t) = pi/2 - atan(t) + gpx.vD2 = dir; +*/ + dir = atan2(gpx.vE, gpx.vN) * 180 / M_PI; + if (dir < 0) dir += 360; + gpx.vD = dir; + + return 0; +} + +int get_GPS3() { + int err=0; + + // ((framebyte(pos_GPS3)<<8) | framebyte(pos_GPS3+1)) != pck_GPS3 ? + if ( framebyte(pos_GPS3) != ((pck_GPS3>>8) & 0xFF) ) { + gpx.crc |= crc_GPS3; + return -1; + } + + err = check_CRC(pos_GPS3, pck_GPS3); + if (err) gpx.crc |= crc_GPS3; + + err |= get_GPSkoord(); + + return err; +} + +int get_Aux() { +// +// "Ozone Sounding with Vaisala Radiosonde RS41" user's guide +// + int i, auxlen, auxcrc, count7E, pos7E; + + count7E = 0; + pos7E = pos_AUX; + + // 7Exx: xdata + while ( pos7E < FRAME_LEN && framebyte(pos7E) == 0x7E ) { + + auxlen = framebyte(pos7E+1); + auxcrc = framebyte(pos7E+2+auxlen) | (framebyte(pos7E+2+auxlen+1)<<8); + + if ( auxcrc == crc16(pos7E+2, auxlen) ) { + if (count7E == 0) fprintf(stdout, "\n # xdata = "); + else fprintf(stdout, " # "); + + //fprintf(stdout, " # %02x : ", framebyte(pos7E+2)); + for (i = 1; i < auxlen; i++) { + fprintf(stdout, "%c", framebyte(pos7E+2+i)); + } + count7E++; + pos7E += 2+auxlen+2; + } + else { + pos7E = FRAME_LEN; + gpx.crc |= crc_AUX; + } + } + + i = check_CRC(pos7E, 0x7600); // 0x76xx: 00-padding block + if (i) gpx.crc |= crc_ZERO; + + return count7E; +} + +int get_Calconf(int out) { + int i; + unsigned byte; + ui8_t calfr = 0; + ui8_t burst = 0; + ui16_t fw = 0; + int freq = 0, f0 = 0, f1 = 0; + char sondetyp[9]; + int err = 0; + + byte = framebyte(pos_CalData); + calfr = byte; + err = check_CRC(pos_FRAME, pck_FRAME); + + if (option_verbose == 3) { + fprintf(stdout, "\n"); // fflush(stdout); + fprintf(stdout, "[%5d] ", gpx.frnr); + fprintf(stdout, " 0x%02x: ", calfr); + for (i = 0; i < 16; i++) { + byte = framebyte(pos_CalData+1+i); + fprintf(stdout, "%02x ", byte); + } + if (err == 0) fprintf(stdout, "[OK]"); + else fprintf(stdout, "[NO]"); + fprintf(stdout, " "); + } + + if (out && err == 0) + { + if (calfr == 0x01 && option_verbose /*== 2*/) { + fw = framebyte(pos_CalData+6) | (framebyte(pos_CalData+7)<<8); + fprintf(stdout, ": fw 0x%04x ", fw); + } + + if (calfr == 0x02 && option_verbose /*== 2*/) { + byte = framebyte(pos_Calburst); + burst = byte; // fw >= 0x4ef5, BK irrelevant? (burst-killtimer in 0x31?) + fprintf(stdout, ": BK %02X ", burst); + if (option_verbose == 3) { // killtimer + int kt = frame[0x5A] + (frame[0x5B] << 8); // short? + if ( kt != 0xFFFF ) fprintf(stdout, ": kt 0x%04x = %dsec = %.1fmin ", kt, kt, kt/60.0); + } + } + + if (calfr == 0x00 && option_verbose) { + byte = framebyte(pos_Calfreq) & 0xC0; // erstmal nur oberste beiden bits + f0 = (byte * 10) / 64; // 0x80 -> 1/2, 0x40 -> 1/4 ; dann mal 40 + byte = framebyte(pos_Calfreq+1); + f1 = 40 * byte; + freq = 400000 + f1+f0; // kHz; + fprintf(stdout, ": fq %d ", freq); + } + + if (calfr == 0x31 && option_verbose == 3) { + int bt = frame[0x59] + (frame[0x5A] << 8); // short? + // fw >= 0x4ef5: default=[88 77]=0x7788sec=510min + if ( bt != 0x0000 ) fprintf(stdout, ": bt 0x%04x = %dsec = %.1fmin ", bt, bt, bt/60.0); + } + + if (calfr == 0x21 && option_verbose /*== 2*/) { // eventuell noch zwei bytes in 0x22 + for (i = 0; i < 9; i++) sondetyp[i] = 0; + for (i = 0; i < 8; i++) { + byte = framebyte(pos_CalRSTyp + i); + if ((byte >= 0x20) && (byte < 0x7F)) sondetyp[i] = byte; + else if (byte == 0x00) sondetyp[i] = '\0'; + } + fprintf(stdout, ": %s ", sondetyp); + } + } + + return 0; +} + +/* + frame[pos_FRAME-1] == 0x0F: len == NDATA_LEN(320) + frame[pos_FRAME-1] == 0xF0: len == FRAME_LEN(518) +*/ +int frametype() { // -4..+4: 0xF0 -> -4 , 0x0F -> +4 + int i; + ui8_t b = frame[pos_FRAME-1]; + int ft = 0; + for (i = 0; i < 4; i++) { + ft += ((b>>i)&1) - ((b>>(i+4))&1); + } + return ft; +} + +/* ------------------------------------------------------------------------------------ */ +/* + (uses fec-lib by KA9Q) + ka9q-fec: + gcc -c init_rs_char.c + gcc -c decode_rs_char.c + +#include "fec.h" // ka9q-fec + + +void *rs; +unsigned char codeword1[rs_N], codeword2[rs_N]; + + rs = init_rs_char( 8, 0x11d, 0, 1, rs_R, 0); + + // ka9q-fec301: p(x) = p[0]x^(N-1) + ... + p[N-2]x + p[N-1] + // -> cw[i] = codeword[RS.N-1-i] + +*/ + +#define rs_N 255 +#define rs_R 24 +#define rs_K (rs_N-rs_R) + +ui8_t cw1[rs_N], cw2[rs_N]; + +int rs41_ecc(int frmlen) { +// richtige framelen wichtig fuer 0-padding + + int i, leak, ret = 0; + int errors1, errors2; + ui8_t err_pos1[rs_R], err_pos2[rs_R], + err_val1[rs_R], err_val2[rs_R]; + + + if (frmlen > FRAME_LEN) frmlen = FRAME_LEN; + cfg_rs41.frmlen = frmlen; + cfg_rs41.msglen = (frmlen-56)/2; // msgpos=56; + leak = frmlen % 2; + + for (i = frmlen; i < FRAME_LEN; i++) frame[i] = 0; // FRAME_LEN-HDR = 510 = 2*255 + + + for (i = 0; i < rs_R; i++) cw1[i] = frame[cfg_rs41.parpos+i ]; + for (i = 0; i < rs_R; i++) cw2[i] = frame[cfg_rs41.parpos+i+rs_R]; + for (i = 0; i < rs_K; i++) cw1[rs_R+i] = frame[cfg_rs41.msgpos+2*i ]; + for (i = 0; i < rs_K; i++) cw2[rs_R+i] = frame[cfg_rs41.msgpos+2*i+1]; + + errors1 = rs_decode(cw1, err_pos1, err_val1); + errors2 = rs_decode(cw2, err_pos2, err_val2); + + + if (option_ecc == 2 && (errors1 < 0 || errors2 < 0)) { + frame[pos_FRAME] = (pck_FRAME>>8)&0xFF; frame[pos_FRAME+1] = pck_FRAME&0xFF; + frame[pos_PTU] = (pck_PTU >>8)&0xFF; frame[pos_PTU +1] = pck_PTU &0xFF; + frame[pos_GPS1] = (pck_GPS1 >>8)&0xFF; frame[pos_GPS1 +1] = pck_GPS1 &0xFF; + frame[pos_GPS2] = (pck_GPS2 >>8)&0xFF; frame[pos_GPS2 +1] = pck_GPS2 &0xFF; + frame[pos_GPS3] = (pck_GPS3 >>8)&0xFF; frame[pos_GPS3 +1] = pck_GPS3 &0xFF; + if (frametype() < -2) { + for (i = NDATA_LEN + 7; i < FRAME_LEN-2; i++) frame[i] = 0; + } + for (i = 0; i < rs_K; i++) cw1[rs_R+i] = frame[cfg_rs41.msgpos+2*i ]; + for (i = 0; i < rs_K; i++) cw2[rs_R+i] = frame[cfg_rs41.msgpos+2*i+1]; + errors1 = rs_decode(cw1, err_pos1, err_val1); + errors2 = rs_decode(cw2, err_pos2, err_val2); + } + + + // Wenn Fehler im 00-padding korrigiert wurden, + // war entweder der frame zu kurz, oder + // Fehler wurden falsch korrigiert; + // allerdings ist bei t=12 die Wahrscheinlichkeit, + // dass falsch korrigiert wurde mit 1/t! sehr gering. + + // check CRC32 + // CRC32 OK: + //for (i = 0; i < cfg_rs41.hdrlen; i++) frame[i] = data[i]; + for (i = 0; i < rs_R; i++) { + frame[cfg_rs41.parpos+ i] = cw1[i]; + frame[cfg_rs41.parpos+rs_R+i] = cw2[i]; + } + for (i = 0; i < rs_K; i++) { // cfg_rs41.msglen <= rs_K + frame[cfg_rs41.msgpos+ 2*i] = cw1[rs_R+i]; + frame[cfg_rs41.msgpos+1+2*i] = cw2[rs_R+i]; + } + if (leak) { + frame[cfg_rs41.msgpos+2*i] = cw1[rs_R+i]; + } + + + ret = errors1 + errors2; + if (errors1 < 0 || errors2 < 0) { + ret = 0; + if (errors1 < 0) ret |= 0x1; + if (errors2 < 0) ret |= 0x2; + ret = -ret; + } + + return ret; +} + +/* ------------------------------------------------------------------------------------ */ + + +int print_position(int ec) { + int i; + int err, err0, err1, err2, err3; + int output, out_mask; + + err = get_FrameConf(); + + err1 = get_GPS1(); + err2 = get_GPS2(); + err3 = get_GPS3(); + + err0 = get_PTU(); + + out_mask = crc_FRAME|crc_GPS1|crc_GPS3; + output = ((gpx.crc & out_mask) != out_mask); // (!err || !err1 || !err3); + + if (output) { + + if (!err) { + fprintf(stdout, "[%5d] ", gpx.frnr); + fprintf(stdout, "(%s) ", gpx.id); + } + if (!err1) { + Gps2Date(gpx.week, gpx.gpssec, &gpx.jahr, &gpx.monat, &gpx.tag); + fprintf(stdout, "%s ", weekday[gpx.wday]); + fprintf(stdout, "%04d-%02d-%02d %02d:%02d:%06.3f", + gpx.jahr, gpx.monat, gpx.tag, gpx.std, gpx.min, gpx.sek); + if (option_verbose == 3) fprintf(stdout, " (W %d)", gpx.week); + } + if (!err3) { + fprintf(stdout, " "); + fprintf(stdout, " lat: %.5f ", gpx.lat); + fprintf(stdout, " lon: %.5f ", gpx.lon); + fprintf(stdout, " alt: %.2f ", gpx.alt); + //if (option_verbose) + { + //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_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; + } + } + 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); + } + } + + err |= err1 | err3; + + return err; +} + +void print_frame(int len) { + int i, ec = 0, ft; + + gpx.crc = 0; + + //frame[pos_FRAME-1] == 0x0F: len == NDATA_LEN(320) + //frame[pos_FRAME-1] == 0xF0: len == FRAME_LEN(518) + ft = frametype(); + if (ft > 2) len = NDATA_LEN; + // STD-frames mit 00 auffuellen fuer Fehlerkorrektur + if (len > NDATA_LEN && len < NDATA_LEN+XDATA_LEN-10) { + if (ft < -2) { + len = NDATA_LEN + 7; // std-O3-AUX-frame + } + } + // AUX-frames mit vielen Fehlern besser mit 00 auffuellen + + for (i = len; i < FRAME_LEN-2; i++) { + frame[i] = 0; + } + if (ft > 2 || len == NDATA_LEN) { + frame[FRAME_LEN-2] = 0; + frame[FRAME_LEN-1] = 0; + } + if (len > NDATA_LEN) len = FRAME_LEN; + else len = NDATA_LEN; + + + if (option_ecc) { + ec = rs41_ecc(len); + } + + + if (option_raw) { + if (option_ecc == 2 && ec >= 0) { + if (len < FRAME_LEN && frame[FRAME_LEN-1] != 0) len = FRAME_LEN; + } + for (i = 0; i < len; i++) { + fprintf(stdout, "%02x", frame[i]); + } + if (option_ecc) { + if (ec >= 0) fprintf(stdout, " [OK]"); else fprintf(stdout, " [NO]"); + if (option_ecc == 2 && ec > 0) fprintf(stdout, " (%d)", ec); + } + fprintf(stdout, "\n"); + } + else if (option_sat) { + get_SatData(); + } + else { + print_position(ec); + } +} + + +int main(int argc, char *argv[]) { + + FILE *fp; + char *fpname = NULL; + float spb = 0.0; + char bitbuf[8]; + int bit_count = 0, + bitpos = 0, + byte_count = FRAMESTART, + ft_len = FRAME_LEN, + header_found = 0; + int bit, byte; + int frmlen = FRAME_LEN; + int headerlen; + int herrs, herr1; + int bitQ, Qerror_count; + + int k, K; + float mv; + unsigned int mv_pos, mv0_pos; + int mp = 0; + + float thres = 0.7; + + int bitofs = 0; + int symlen = 1; + int shift = 2; + + +#ifdef CYGWIN + _setmode(fileno(stdin), _O_BINARY); // _fileno(stdin) +#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, -vx, -vv (info, aux, info/conf)\n"); + fprintf(stderr, " -r, --raw\n"); + fprintf(stderr, " -i, --invert\n"); + fprintf(stderr, " --crc (check CRC)\n"); + fprintf(stderr, " --ecc (Reed-Solomon)\n"); + fprintf(stderr, " --std (std framelen)\n"); + fprintf(stderr, " --ths (peak threshold; default=%.1f)\n", thres); + fprintf(stderr, " --iq0,2,3 (IQ data)\n"); + return 0; + } + else if ( (strcmp(*argv, "-v") == 0) || (strcmp(*argv, "--verbose") == 0) ) { + option_verbose = 1; + } + else if (strcmp(*argv, "-vx") == 0) { option_verbose = 2; } + else if (strcmp(*argv, "-vv") == 0) { option_verbose = 3; } + else if (strcmp(*argv, "-vvv") == 0) { option_verbose = 4; } + else if (strcmp(*argv, "--crc") == 0) { option_crc = 1; } + else if (strcmp(*argv, "--res") == 0) { option_res = 1; } + else if ( (strcmp(*argv, "-r") == 0) || (strcmp(*argv, "--raw") == 0) ) { + option_raw = 1; + } + else if ( (strcmp(*argv, "-i") == 0) || (strcmp(*argv, "--invert") == 0) ) { + option_inv = 1; + } + else if (strcmp(*argv, "--ecc" ) == 0) { option_ecc = 1; } + else if (strcmp(*argv, "--ecc2") == 0) { option_ecc = 2; } + else if (strcmp(*argv, "--std" ) == 0) { frmlen = 320; } // NDATA_LEN + 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, "--ch2") == 0) { wav_channel = 1; } // right channel (default: 0=left) + else if (strcmp(*argv, "--ths") == 0) { + ++argv; + if (*argv) { + thres = atof(*argv); + } + else return -1; + } + else if ( (strcmp(*argv, "-d") == 0) ) { + ++argv; + if (*argv) { + shift = atoi(*argv); + if (shift > 6) shift = 6; + if (shift < -2) shift = -2; + } + else return -1; + } + else if (strcmp(*argv, "--iq0") == 0) { option_iq = 1; } // differential/FM-demod + else if (strcmp(*argv, "--iq2") == 0) { option_iq = 2; } + else if (strcmp(*argv, "--iq3") == 0) { option_iq = 3; } // iq2==iq3 + else if (strcmp(*argv, "--ofs") == 0) { option_ofs = 1; } + else if (strcmp(*argv, "--dbg") == 0) { option_dbg = 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; + + + if (option_iq) wav_channel = 0; + + + spb = read_wav_header(fp, (float)BAUD_RATE, wav_channel); + if ( spb < 0 ) { + fclose(fp); + fprintf(stderr, "error: wav header\n"); + return -1; + } + if ( spb < 8 ) { + fprintf(stderr, "note: sample rate low\n"); + } + + + if (option_ecc) { + rs_init_RS255(); + } + + // rs41: BT=0.5, h=0.8,1.0 ? + symlen = 1; + headerlen = strlen(header); + bitofs = shift; // +0 .. +3 // FM: +1 , IQ: +2 + K = init_buffers(header, headerlen, 0.5, option_iq); // BT=0.5 (IQ-Int: BT > 0.5 ?) + if ( K < 0 ) { + fprintf(stderr, "error: init buffers\n"); + return -1; + }; + + + k = 0; + mv = -1; 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); + k = 0; + } + else { + mv = 0.0; + continue; + } + + if (mv > thres && mp > 0) { + if (mv_pos > mv0_pos) { + + header_found = 0; + herrs = headcmp(symlen, header, headerlen, mv_pos, mv<0, 0); // symlen=1 + herr1 = 0; + if (herrs <= 3 && herrs > 0) { + herr1 = headcmp(symlen, header, headerlen, mv_pos+1, mv<0, 0); + if (herr1 < herrs) { + herrs = herr1; + herr1 = 1; + } + } + if (herrs <= 2) header_found = 1; // herrs <= 2 bitfehler in header + + if (header_found) { + + if (/*preamble &&*/ option_ofs) + { + float freq = 0.0; + float snr = 0.0; + get_fqofs(headerlen, mv_pos, &freq, &snr); + if (option_dbg) { + fprintf(stderr, "fq-ofs: %+.2f Hz snr: %.2f dB\n", freq, snr); + } + } + + byte_count = FRAMESTART; + bit_count = 0; // byte_count*8-HEADLEN + bitpos = 0; + + Qerror_count = 0; + ft_len = frmlen; + + while ( byte_count < frmlen ) { + if (option_iq) { + bitQ = read_IDsbit(fp, symlen, &bit, option_inv, bitofs, bit_count==0, 0); // symlen=1, return: zeroX/bit + } + else { + bitQ = read_sbit(fp, symlen, &bit, option_inv, bitofs, bit_count==0, 0); // symlen=1, return: zeroX/bit + } + if ( bitQ == EOF) break; + bit_count += 1; + bitbuf[bitpos] = bit; + bitpos++; + if (bitpos == BITS) { + bitpos = 0; + byte = bits2byte(bitbuf); + frame[byte_count] = byte ^ mask[byte_count % MASK_LEN]; + + byteQ[byte_count] = get_bufvar(0); + + if (byte_count > NDATA_LEN) { // Fehler erst ab minimaler framelen Zaehlen + if (byteQ[byte_count]*2 > byteQ[byte_count-300]*3) { // Var(frame)/Var(noise) ca. 1:2 + Qerror_count += 1; + } + } + + byte_count++; + } + if (Qerror_count == 4) { // framelen = 320 oder 518 + ft_len = byte_count; + Qerror_count += 1; + } + } + header_found = 0; + print_frame(ft_len); + + while ( bit_count < BITS*(FRAME_LEN-8+24) ) { + bitQ = read_sbit(fp, symlen, &bit, option_inv, bitofs, bit_count==0, 0); // symlen=1, return: zeroX/bit + if ( bitQ == EOF) break; + bit_count++; + } + + byte_count = FRAMESTART; + } + } + } + + } + + + free_buffers(); + + fclose(fp); + + return 0; +} + diff --git a/demod/rs92dm_dft.c b/demod/rs92dm_dft.c index 904bdf8..9af507e 100644 --- a/demod/rs92dm_dft.c +++ b/demod/rs92dm_dft.c @@ -79,7 +79,9 @@ int option_verbose = 0, // ausfuehrliche Anzeige option_aux = 0, // Aux/Ozon option_der = 0, // linErr option_ths = 0, + option_json = 0, // JSON output (auto_rx) rawin = 0; +int wav_channel = 0; // audio channel: left double dop_limit = 9.9; double d_err = 10000; @@ -1142,12 +1144,16 @@ int print_position() { // GPS-Hoehe ueber Ellipsoid } } } - // Print out telemetry data as JSON, even if we don't have a valid GPS lock. - if (!err1 && !err2){ - if (gpx.aux[0] != 0 || gpx.aux[1] != 0 || gpx.aux[2] != 0 || gpx.aux[3] != 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, \"aux\": \"%04x%04x%04x%04x\"}\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.aux[0], gpx.aux[1], gpx.aux[2], gpx.aux[3]); - } 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_json) + { + // Print out telemetry data as JSON, even if we don't have a valid GPS lock. + if (!err1 && !err2){ + if (gpx.aux[0] != 0 || gpx.aux[1] != 0 || gpx.aux[2] != 0 || gpx.aux[3] != 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, \"aux\": \"%04x%04x%04x%04x\"}\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.aux[0], gpx.aux[1], gpx.aux[2], gpx.aux[3]); + } 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 ); + } } } @@ -1262,6 +1268,7 @@ int main(int argc, char *argv[]) { fprintf(stderr, " --crc (CRC check GPS)\n"); fprintf(stderr, " --ecc (Reed-Solomon)\n"); fprintf(stderr, " --ths (peak threshold; default=%.1f)\n", thres); + fprintf(stderr, " --json (JSON output)\n"); return 0; } else if ( (strcmp(*argv, "--vel") == 0) ) { @@ -1335,6 +1342,8 @@ int main(int argc, char *argv[]) { else if (strcmp(*argv, "-g2") == 0) { option_vergps = 2; } // verbose2 GPS (bancroft) else if (strcmp(*argv, "-gg") == 0) { option_vergps = 8; } // vverbose GPS else if (strcmp(*argv, "--ecc") == 0) { option_ecc = 1; } + else if (strcmp(*argv, "--json") == 0) { option_json = 1; } // JSON output (for auto_rx) + else if (strcmp(*argv, "--ch2") == 0) { wav_channel = 1; } // right channel (default: 0=left) else if (strcmp(*argv, "--ths") == 0) { ++argv; if (*argv) { @@ -1379,7 +1388,7 @@ int main(int argc, char *argv[]) { } - spb = read_wav_header(fp, (float)BAUD_RATE); + spb = read_wav_header(fp, (float)BAUD_RATE, wav_channel); if ( spb < 0 ) { fclose(fp); fprintf(stderr, "error: wav header\n"); diff --git a/ecc/bch_ecc.c b/ecc/bch_ecc.c index 4815a42..2a96788 100644 --- a/ecc/bch_ecc.c +++ b/ecc/bch_ecc.c @@ -10,8 +10,14 @@ Vaisala RS92, RS41: RS(255, 231), t=12 + f=X^8+X^4+X^3+X^2+1, b=0 g(X) = (X-alpha^0)...(X-alpha^(2t-1)) + LMS6: + RS(255, 223), t=16 (CCSDS) + f=X^8+X^7+X^2+X+1, b=112 + g(X) = (X-(alpha^11)^112)...(X-(alpha^11)^(112+2t-1)) + Meisei: bin.BCH(63, 51), t=2 g(X) = (X^6+X+1)(X^6+X^4+X^2+X+1) @@ -66,6 +72,10 @@ static GF_t GF256RS = { 0x11D, // RS-GF(2^8): X^8 + X^4 + X^3 + X^2 + 1 : 0x11D 256, // 2^8 0x02 }; // generator: alpha = X +static GF_t GF256RSccsds = { 0x187, // RS-GF(2^8): X^8 + X^7 + X^2 + X + 1 : 0x187 + 256, // 2^8 + 0x02 }; // generator: alpha = X + static GF_t GF64BCH = { 0x43, // BCH-GF(2^6): X^6 + X + 1 : 0x43 64, // 2^6 0x02 }; // generator: alpha = X @@ -85,12 +95,14 @@ typedef struct { ui8_t R; // RS: R=2t, BCH: R<=mt ui8_t K; // K=N-R ui8_t b; + ui8_t p; ui8_t ip; // p*ip = 1 mod N ui8_t g[MAX_DEG+1]; // ohne g[] eventuell als init_return } RS_t; -static RS_t RS256 = { 255, 12, 24, 231, 0, {0}}; -static RS_t BCH64 = { 63, 2, 12, 51, 1, {0}}; +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 GF_t GF; @@ -627,7 +639,7 @@ int era_sigma(int n, ui8_t era_pos[], ui8_t *sigma) { sig[0] = 1; Xa[0] = 1; for (i = 0; i < n; i++) { // n <= 2*RS.t - a_i = exp_a[era_pos[i] % (GF.ord-1)]; + a_i = exp_a[(RS.p*era_pos[i]) % (GF.ord-1)]; Xa[1] = a_i; // Xalp[0..1]: (1-alpha^(j_)X) poly_mul(sig, Xa, sig); } @@ -642,9 +654,9 @@ int syndromes(ui8_t cw[], ui8_t *S) { int i, errors = 0; ui8_t a_i; - // syndromes: e_j=S(alpha^(b+i)) + // syndromes: e_j=S((alpha^p)^(b+i)) (wie in g(X)) for (i = 0; i < 2*RS.t; i++) { - a_i = exp_a[(RS.b+i) % (GF.ord-1)]; // alpha^(b+i) + a_i = exp_a[(RS.p*(RS.b+i)) % (GF.ord-1)]; // (alpha^p)^(b+i) S[i] = poly_eval(cw, a_i); if (S[i]) errors = 1; } @@ -718,7 +730,7 @@ int rs_init_RS255() { GF = GF256RS; check_gen = GF_genTab( GF, exp_a, log_a); - RS = RS256; // N=255, t=12, b=0 + RS = RS256; // N=255, t=12, b=0, p=1 for (i = 0; i <= MAX_DEG; i++) RS.g[i] = 0; for (i = 0; i <= MAX_DEG; i++) Xalp[i] = 0; @@ -733,6 +745,55 @@ int rs_init_RS255() { return check_gen; } +int rs_init_RS255ccsds() { + int i, check_gen; + ui8_t Xalp[MAX_DEG+1]; + + GF = GF256RSccsds; + check_gen = GF_genTab( GF, exp_a, log_a); + + RS = RS256ccsds; // N=255, t=16, b=112, p=11 + for (i = 0; i <= MAX_DEG; i++) RS.g[i] = 0; + for (i = 0; i <= MAX_DEG; i++) Xalp[i] = 0; + + // beta=alpha^p primitive root of g(X) + // beta^ip=alpha // N=255, p=11 -> ip=116 + for (i = 1; i < GF.ord-1; i++) { + if ( (RS.p * i) % (GF.ord-1) == 1 ) { + RS.ip = i; + break; + } + } + + // g(X)=(X-(alpha^p)^b)...(X-(alpha^p)^(b+2t-1)), b=112 + RS.g[0] = 0x01; + Xalp[1] = 0x01; // X + for (i = 0; i < 2*RS.t; i++) { + Xalp[0] = exp_a[(RS.p*(RS.b+i)) % (GF.ord-1)]; // Xalp[0..1]: X - (alpha^p)^(b+i) + poly_mul(RS.g, Xalp, RS.g); + } +/* + RS.g[ 0] = RS.g[32] = exp_a[0]; + RS.g[ 1] = RS.g[31] = exp_a[249]; + RS.g[ 2] = RS.g[30] = exp_a[59]; + RS.g[ 3] = RS.g[29] = exp_a[66]; + RS.g[ 4] = RS.g[28] = exp_a[4]; + RS.g[ 5] = RS.g[27] = exp_a[43]; + RS.g[ 6] = RS.g[26] = exp_a[126]; + RS.g[ 7] = RS.g[25] = exp_a[251]; + RS.g[ 8] = RS.g[24] = exp_a[97]; + RS.g[ 9] = RS.g[23] = exp_a[30]; + RS.g[10] = RS.g[22] = exp_a[3]; + RS.g[11] = RS.g[21] = exp_a[213]; + RS.g[12] = RS.g[20] = exp_a[50]; + RS.g[13] = RS.g[19] = exp_a[66]; + RS.g[14] = RS.g[18] = exp_a[170]; + RS.g[15] = RS.g[17] = exp_a[5]; + RS.g[16] = exp_a[24]; +*/ + return check_gen; +} + int rs_init_BCH64() { int i, check_gen; @@ -761,7 +822,7 @@ int rs_encode(ui8_t cw[]) { return 0; } -// 2*Errors + Erasure <= 2*RS.t +// 2*Errors + Erasure <= 2*t int rs_decode_ErrEra(ui8_t cw[], int nera, ui8_t era_pos[], ui8_t *err_pos, ui8_t *err_val) { ui8_t x, gamma; @@ -793,7 +854,7 @@ int rs_decode_ErrEra(ui8_t cw[], int nera, ui8_t era_pos[], if (nera > 0) { era_sigma(nera, era_pos, sigma); poly_mul(sigma, S, S); - for (i = 2*RS.t; i <= MAX_DEG; i++) S[i] = 0; // S = sig*S mod x^12 + for (i = 2*RS.t; i <= MAX_DEG; i++) S[i] = 0; // S = sig*S mod x^2t } if (errera) @@ -823,7 +884,8 @@ int rs_decode_ErrEra(ui8_t cw[], int nera, ui8_t era_pos[], x = (ui8_t)i; // roll-over if (poly_eval(sigLam, x) == 0) { // Lambda(x)=0 fuer x in erasures[] moeglich // error location index - err_pos[nerr] = log_a[GF_inv(x)]; + ui8_t x1 = GF_inv(x); + err_pos[nerr] = (log_a[x1]*RS.ip) % (GF.ord-1); // error value; bin-BCH: err_val=1 err_val[nerr] = forney(x, Omega, sigLam); //err_val[nerr] == 0, wenn era_val[pos]=0, d.h. cw[pos] schon korrekt @@ -832,7 +894,7 @@ int rs_decode_ErrEra(ui8_t cw[], int nera, ui8_t era_pos[], if (nerr >= deg_sigLam) break; } - // 2*Errors + Erasure <= 2*RS.t + // 2*Errors + Erasure <= 2*t if (nerr < deg_sigLam) errera = -1; // uncorrectable errors else { errera = nerr; @@ -843,7 +905,7 @@ int rs_decode_ErrEra(ui8_t cw[], int nera, ui8_t era_pos[], return errera; } -// Errors <= RS.t +// Errors <= t int rs_decode(ui8_t cw[], ui8_t *err_pos, ui8_t *err_val) { ui8_t tmp[1] = {0}; return rs_decode_ErrEra(cw, 0, tmp, err_pos, err_val);