From f228b43e02244d533a0b588ca49b0e3b87ffdda4 Mon Sep 17 00:00:00 2001 From: Mark Jessop Date: Sat, 29 Apr 2017 23:36:47 +0930 Subject: [PATCH] Modified rs92gps to be parsable easier, initial commit of auto sonde RX utilities. --- auto_rx/aprs_utils.py | 55 ++++++++ auto_rx/auto_rx.py | 310 ++++++++++++++++++++++++++++++++++-------- auto_rx/auto_rx.sh | 8 ++ rs92/rs92gps.c | 15 +- 4 files changed, 326 insertions(+), 62 deletions(-) create mode 100644 auto_rx/aprs_utils.py create mode 100755 auto_rx/auto_rx.sh diff --git a/auto_rx/aprs_utils.py b/auto_rx/aprs_utils.py new file mode 100644 index 0000000..947d17f --- /dev/null +++ b/auto_rx/aprs_utils.py @@ -0,0 +1,55 @@ +# APRS push utils for Sonde auto RX. + +from socket import * + + +# Push a Radiosonde data packet to APRS as an object. +def push_balloon_to_aprs(sonde_data, aprs_comment="BOM Balloon", aprsUser="N0CALL", aprsPass="00000", serverHost = 'rotate.aprs2.net', serverPort = 14580): + # Pad or limit the sonde ID to 9 characters. + object_name = sonde_data["id"].strip() + if len(object_name) > 9: + object_name = object_name[:9] + elif len(object_name) < 9: + object_name = object_name + " "*(9-len(object_name)) + + # Convert float latitude to APRS format (DDMM.MM) + lat = float(sonde_data["lat"]) + lat_degree = abs(int(lat)) + lat_minute = abs(lat - int(lat)) * 60.0 + lat_min_str = ("%02.2f" % lat_minute).zfill(5) + lat_dir = "S" + if lat>0.0: + lat_dir = "N" + lat_str = "%02d%s" % (lat_degree,lat_min_str) + lat_dir + + # Convert float longitude to APRS format (DDDMM.MM) + lon = float(sonde_data["lon"]) + lon_degree = abs(int(lon)) + lon_minute = abs(lon - int(lon)) * 60.0 + lon_min_str = ("%02.2f" % lon_minute).zfill(5) + lon_dir = "E" + if lon<0.0: + lon_dir = "W" + lon_str = "%03d%s" % (lon_degree,lon_min_str) + lon_dir + + # Convert Alt (in metres) to feet + alt = int(float(sonde_data["alt"])/0.3048) + + # Produce the APRS object string. + out_str = ";%s*111111z%s/%sO000/000/A=%06d %s" % (object_name,lat_str,lon_str,alt,aprs_comment) + + # Connect to an APRS-IS server, login, then push our object position in. + + # create socket & connect to server + sSock = socket(AF_INET, SOCK_STREAM) + sSock.connect((serverHost, serverPort)) + # logon + sSock.send('user %s pass %s vers VK5QI-Python 0.01\n' % (aprsUser, aprsPass) ) + # send packet + sSock.send('%s>APRS:%s\n' % (aprsUser, out_str) ) + + # close socket + sSock.shutdown(0) + sSock.close() + + return out_str \ No newline at end of file diff --git a/auto_rx/auto_rx.py b/auto_rx/auto_rx.py index de54451..8d619cd 100644 --- a/auto_rx/auto_rx.py +++ b/auto_rx/auto_rx.py @@ -4,89 +4,289 @@ # # 2017-04 Mark Jessop # +# The following binaries will need to be built and copied to this directory: +# rs92/rs92gps +# scan/rs_detect + import numpy as np +import sys +import logging +import datetime +import time +import Queue +import subprocess +import traceback +from aprs_utils import * +from threading import Thread from StringIO import StringIO from findpeaks import * from os import system -import sys + + +# Receiver Parameters +RX_PPM = 0 +RX_GAIN = 0 # 0 = Auto # Sonde Search Configuration Parameters -MIN_FREQ = 400.4e6 -MAX_FREQ = 403.5e6 -SEARCH_STEP = 500 -MIN_FREQ_DISTANCE = 10000 # Expect a minimum distance of 10 kHz between sondes. -MIN_SNR = 10 # Only takes peaks that are a minimum of 10dB above the noise floor. +MIN_FREQ = 400.4e6 # Search start frequency (Hz) +MAX_FREQ = 403.5e6 # Search stop frequency (Hz) +SEARCH_STEP = 800 # Search step (Hz) +FREQ_QUANTIZATION = 5000 # Quantize search results to 5 kHz steps. +MIN_FREQ_DISTANCE = 1000 # Expect a minimum distance of 10 kHz between sondes. +MIN_SNR = 10 # Only takes peaks that are a minimum of 10dB above the noise floor. +SEARCH_ATTEMPTS = 5 # Number of attempts to search before giving up +SEARCH_DELAY = 300 # Delay between search attempts (seconds) +# Other Receiver Parameters +MAX_RX_TIME = 3*60*60 + +# APRS Output +APRS_OUTPUT_ENABLED = True +APRS_UPDATE_RATE = 30 +APRS_USER = "N0CALL" +APRS_PASS = "000000" +aprs_queue = Queue.Queue(1) + + + +def run_rtl_power(start, stop, step, filename="log_power.csv", dwell = 20): + # Run rtl_power, with a timeout + # rtl_power -f 400400000:403500000:800 -i20 -1 log_power.csv + rtl_power_cmd = "timeout %d rtl_power -f %d:%d:%d -i %d -1 %s" % (dwell+10, start, stop, step, dwell, filename) + logging.info("Running frequency scan.") + ret_code = system(rtl_power_cmd) + if ret_code == 1: + logging.critical("rtl_power call failed!") + sys.exit(1) + else: + return True def read_rtl_power(filename): - ''' Read in frequency samples from a single-shot log file produced by rtl_power''' + ''' Read in frequency samples from a single-shot log file produced by rtl_power''' - # Output buffers. - freq = np.array([]) - power = np.array([]) + # Output buffers. + freq = np.array([]) + power = np.array([]) - freq_step = 0 + freq_step = 0 - # Open file. - f = open(filename,'r') + # Open file. + f = open(filename,'r') - # rtl_power log files are csv's, with the first 6 fields in each line describing the time and frequency scan parameters - # for the remaining fields, which contain the power samples. + # rtl_power log files are csv's, with the first 6 fields in each line describing the time and frequency scan parameters + # for the remaining fields, which contain the power samples. - for line in f: - # Split line into fields. - fields = line.split(',') + for line in f: + # Split line into fields. + fields = line.split(',') - if len(fields) < 6: - raise Exception("Invalid number of samples in input file - corrupt?") + if len(fields) < 6: + logging.error("Invalid number of samples in input file - corrupt?") + raise Exception("Invalid number of samples in input file - corrupt?") - start_date = fields[0] - start_time = fields[1] - start_freq = float(fields[2]) - stop_freq = float(fields[3]) - freq_step = float(fields[4]) - n_samples = int(fields[5]) + start_date = fields[0] + start_time = fields[1] + start_freq = float(fields[2]) + stop_freq = float(fields[3]) + freq_step = float(fields[4]) + n_samples = int(fields[5]) - freq_range = np.arange(start_freq,stop_freq,freq_step) - samples = np.loadtxt(StringIO(",".join(fields[6:])),delimiter=',') + freq_range = np.arange(start_freq,stop_freq,freq_step) + samples = np.loadtxt(StringIO(",".join(fields[6:])),delimiter=',') - # Add frequency range and samples to output buffers. - freq = np.append(freq, freq_range) - power = np.append(power, samples) + # Add frequency range and samples to output buffers. + freq = np.append(freq, freq_range) + power = np.append(power, samples) - f.close() - return (freq, power, freq_step) + f.close() + return (freq, power, freq_step) +def quantize_freq(freq_list, quantize=5000): + return np.round(freq_list/quantize)*quantize + +def detect_sonde(frequency): + """ Receive some FM and attempt to detect the presence of a radiosonde. """ + rx_test_command = "timeout 10s rtl_fm -p %d -M fm -s 15k -f %d 2>/dev/null |" % (RX_PPM, 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 += "./rs_detect -z -t 8 2>/dev/null" + + logging.info("Attempting sonde detection on %.3f MHz" % (frequency/1e6)) + ret_code = system(rx_test_command) + + ret_code = ret_code >> 8 + + if ret_code == 3: + logging.info("Detected a RS41!") + return "RS41" + elif ret_code == 4: + logging.info("Detected a RS92!") + return "RS92" + else: + return None + + +# $rs = ""; +# $WFM = ""; +# if (abs($ret) == 2) { $rs = "dfm"; $breite = "15k"; $dec = './dfm06 -vv --ecc'; $filter = "lowpass 2000 highpass 20"; } +# if (abs($ret) == 3) { $rs = "rs41"; $breite = "12k"; $dec = './rs41ecc --ecc -v'; $filter = "lowpass 2600"; } +# if (abs($ret) == 4) { $rs = "rs92"; $breite = "12k"; $dec = './rs92gps --vel2 -a almanac.txt'; $filter = "lowpass 2500 highpass 20"; } +# if (abs($ret) == 5) { $rs = "m10"; $breite = "24k"; $dec = './m10x -vv'; $filter = "highpass 20"; } +# if (abs($ret) == 6) { $rs = "imet"; $breite = "40k"; $dec = './imet1ab -v'; $filter = "highpass 20"; $WFM = "-o 4"; } +# if ($inv) { print "-";} print uc($rs)," ($utc)\n"; +# $utc = strftime('%Y%m%d_%H%M%S', gmtime); +# $wavfile = $rs."-".$utc."Z-".$freq."Hz.wav"; +# if ($rs) { +# system("timeout 30s rtl_fm -p $ppm -M fm $WFM -s $breite -f $freq 2>/dev/null |\ +# sox -t raw -r $breite -e s -b 16 -c 1 - -r 48000 -b 8 -t wav - $filter 2>/dev/null |\ +# tee $log_dir/$wavfile | $dec $inv 2>/dev/null"); + +def process_rs92_line(line): + try: + + params = line.split(',') + if len(params) < 9: + logging.error("Not enough parameters: %s" % line) + return + + # Attempt to extract parameters. + rs92_frame = {} + rs92_frame['frame'] = int(params[0]) + rs92_frame['id'] = str(params[1]) + rs92_frame['time'] = str(params[2]) + rs92_frame['lat'] = float(params[3]) + rs92_frame['lon'] = float(params[4]) + rs92_frame['alt'] = float(params[5]) + rs92_frame['vel_h'] = float(params[6]) + rs92_frame['heading'] = float(params[7]) + rs92_frame['vel_v'] = float(params[8]) + rs92_frame['ok'] = 'OK' + + logging.info("RS92: %s,%d,%s,%.5f,%.5f,%.1f" % (rs92_frame['id'], rs92_frame['frame'],rs92_frame['time'], rs92_frame['lat'], rs92_frame['lon'], rs92_frame['alt'])) + + return rs92_frame + + except: + logging.error("Could not parse string: %s" % line) + return None + +def decode_rs92(frequency, ppm=RX_PPM, rx_queue=None): + """ Decode a RS92 sonde """ + decode_cmd = "rtl_fm -p %d -M fm -s 12k -f %d 2>/dev/null |" % (ppm, frequency) + 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 += "./rs92gps --vel2 --crc -a almanac.txt" + + + rx_start_time = time.time() + + rx = subprocess.Popen(decode_cmd, shell=True, stdin=None, stdout=subprocess.PIPE) + + while True: + try: + line = rx.stdout.readline() + if (line != None) and (line != ""): + data = process_rs92_line(line) + + if data != None: + data['freq'] = "%.3f MHz" % (frequency/1e6) + + if rx_queue != None: + try: + rx_queue.put_nowait(data) + except: + pass + except: + traceback.print_exc() + logging.error("Could not read from rxer stdout?") + rx.kill() + return + + +def internet_push_thread(): + global aprs_queue, APRS_USER, APRS_PASS, APRS_UPDATE_RATE, APRS_OUTPUT_ENABLED + print("Started thread.") + while APRS_OUTPUT_ENABLED: + try: + data = aprs_queue.get_nowait() + except: + continue + + aprs_comment = "Sonde Auto-RX Test %s" % data['freq'] + if APRS_OUTPUT_ENABLED: + aprs_data = push_balloon_to_aprs(data,aprs_comment=aprs_comment,aprsUser=APRS_USER, aprsPass=APRS_PASS) + logging.debug("Data pushed to APRS-IS: %s" % aprs_data) + + time.sleep(APRS_UPDATE_RATE) + + print("Closing thread.") if __name__ == "__main__": - # Run rtl_power, with a timeout - # rtl_power -f 400400000:403500000:800 -i20 -1 log_power.csv - rtl_power_cmd = "timeout 30 rtl_power -f %d:%d:%d -i20 -1 log_power.csv" % (MIN_FREQ, MAX_FREQ, SEARCH_STEP) - print("Running: %s" % rtl_power_cmd) - ret_code = system(rtl_power_cmd) - if ret_code == 1: - print("rtl_power call failed!") - sys.exit(1) - # Read in result - (freq, power, step) = read_rtl_power('log_power.csv') + # Setup logging. + logging.basicConfig(format='%(asctime)s %(levelname)s:%(message)s', filename=datetime.datetime.utcnow().strftime("log/%Y%m%d-%H%M%S.log"), level=logging.DEBUG) + logging.getLogger().addHandler(logging.StreamHandler()) - # Rough approximation of the noise floor of the received power spectrum. - power_nf = np.mean(power) + sonde_freq = 0.0 + sonde_type = None - # Detect peaks. - peak_indices = detect_peaks(power, mph=(power_nf+MIN_SNR), mpd=(MIN_FREQ_DISTANCE/step)) + while SEARCH_ATTEMPTS>0: + # Scan Band + run_rtl_power(MIN_FREQ, MAX_FREQ, SEARCH_STEP) - if len(peak_indices) == 0: - print("No peaks found!") - sys.exit(1) - - peak_frequencies = freq[peak_indices] - - - print("Peaks found at: %s" % str(peak_frequencies/1e6)) + # Read in result + (freq, power, step) = read_rtl_power('log_power.csv') + # Rough approximation of the noise floor of the received power spectrum. + power_nf = np.mean(power) + + # Detect peaks. + peak_indices = detect_peaks(power, mph=(power_nf+MIN_SNR), mpd=(MIN_FREQ_DISTANCE/step), show = False) + + if len(peak_indices) == 0: + logging.info("No peaks found!") + continue + + # Sort peaks by power. + peak_powers = power[peak_indices] + peak_freqs = freq[peak_indices] + peak_frequencies = np.flip(peak_freqs[np.argsort(peak_powers)],0) + + # Quantize to nearest 5 kHz + peak_frequencies = quantize_freq(peak_frequencies, FREQ_QUANTIZATION) + + logging.info("Peaks found at (MHz): %s" % str(peak_frequencies/1e6)) + + for freq in peak_frequencies: + detected = detect_sonde(freq) + if detected != None: + sonde_freq = freq + sonde_type = detected + break + + if sonde_type != None: + break + else: + SEARCH_ATTEMPTS -= 1 + logging.warning("Search attempt failed, %d attempts remaining. Waiting %d seconds." % (SEARCH_ATTEMPTS, SEARCH_DELAY)) + time.sleep(SEARCH_DELAY) + + if SEARCH_ATTEMPTS == 0: + logging.error("No sondes detcted, exiting.") + sys.exit(0) + + logging.info("Starting decoding of %s on %.3f MHz" % (sonde_type, sonde_freq/1e6)) + + t = Thread(target=internet_push_thread) + t.start() + + if sonde_type == "RS92": + decode_rs92(sonde_freq, rx_queue=aprs_queue) + elif sonde_type == "RS41": + logging.error("Not implemented.") + else: + pass + + APRS_OUTPUT_ENABLED = False diff --git a/auto_rx/auto_rx.sh b/auto_rx/auto_rx.sh new file mode 100755 index 0000000..2ddf615 --- /dev/null +++ b/auto_rx/auto_rx.sh @@ -0,0 +1,8 @@ +#!/bin/bash +# Radiosonde Auto-RX Script + +# Start auto_rx process with a 3 hour timeout. +timeout 10800 python auto_rx.py + +# Clean up rtl_fm process. +killall rtl_fm \ No newline at end of file diff --git a/rs92/rs92gps.c b/rs92/rs92gps.c index 9e2e8ab..1626023 100644 --- a/rs92/rs92gps.c +++ b/rs92/rs92gps.c @@ -1287,8 +1287,8 @@ int print_position() { // GPS-Hoehe ueber Ellipsoid } if (!err1) { - fprintf(stdout, "[%5d] ", gpx.frnr); - fprintf(stdout, "(%s) ", gpx.id); + fprintf(stdout, "%5d,", gpx.frnr); + fprintf(stdout, "%s,", gpx.id); } if (!err2) { @@ -1297,20 +1297,21 @@ int print_position() { // GPS-Hoehe ueber Ellipsoid //fprintf(stdout, "(W %d) ", gpx.week); fprintf(stdout, "(%04d-%02d-%02d) ", gpx.jahr, gpx.monat, gpx.tag); } - fprintf(stdout, "%s ", weekday[gpx.wday]); // %04.1f: wenn sek >= 59.950, wird auf 60.0 gerundet - fprintf(stdout, "%02d:%02d:%06.3f", gpx.std, gpx.min, gpx.sek); + //fprintf(stdout, "%s ", weekday[gpx.wday]); // %04.1f: wenn sek >= 59.950, wird auf 60.0 gerundet + fprintf(stdout, "%02d:%02d:%06.3f,", gpx.std, gpx.min, gpx.sek); if (n > 0) { fprintf(stdout, " "); - if (almanac) fprintf(stdout, " lat: %.4f lon: %.4f alt: %.1f ", gpx.lat, gpx.lon, gpx.alt); - else fprintf(stdout, " lat: %.5f lon: %.5f alt: %.1f ", gpx.lat, gpx.lon, gpx.alt); + if (almanac) fprintf(stdout, "%.4f,%.4f,%.1f,", gpx.lat, gpx.lon, gpx.alt); + else fprintf(stdout, "%.5f,%.5f,%.1f,", gpx.lat, gpx.lon, gpx.alt); if (option_verbose && option_vergps != 8) { fprintf(stdout, " (d:%.1f)", gpx.diter); } if (option_vel /*&& option_vergps >= 2*/) { - fprintf(stdout," vH: %4.1f D: %5.1f° vV: %3.1f ", gpx.vH, gpx.vD, gpx.vU); + fprintf(stdout,"%4.1f,%5.1f,%3.1f,", gpx.vH, gpx.vD, gpx.vU); + //fprintf(stdout," VEL: vH: %4.1f D: %5.1f° vV: %3.1f ", gpx.vH, gpx.vD, gpx.vU); } if (option_verbose) { if (option_vergps != 2) {