kopia lustrzana https://github.com/projecthorus/radiosonde_auto_rx
Change dft_detect thresholds, commit working for new thresholds.
Mark Jessop
rodzic
a126f4a5dd
commit
cb984cf2b8
|
|
@ -23,7 +23,7 @@ GENERATED_DIR = "./generated"
|
|||
|
||||
# Range of Eb/N0 SNRs to produce.
|
||||
# 10-20 dB seems to be the range where the demodulators fall over.
|
||||
EBNO_RANGE = np.arange(10,20,0.5)
|
||||
EBNO_RANGE = np.arange(5,20.5,0.5)
|
||||
|
||||
# Normalise the samples to +/- 1.0!
|
||||
# If we don't do this, bad things can happen later down the track...
|
||||
|
|
@ -115,7 +115,7 @@ if __name__ == '__main__':
|
|||
for ebno in EBNO_RANGE:
|
||||
_data_noise = add_noise(_data, variance=_var, baud_rate=_baud_rate, ebno=ebno, fs=_fs)
|
||||
|
||||
_out_file = _source.split('.bin')[0] + "_%.1fdB"%ebno + ".bin"
|
||||
_out_file = _source.split('.bin')[0] + "_%04.1fdB"%ebno + ".bin"
|
||||
|
||||
save_sample(_data_noise, _out_file)
|
||||
print("Saved file: %s" % _out_file)
|
||||
|
|
|
|||
|
|
@ -0,0 +1,76 @@
|
|||
# dft_detect Threshold Optimzation
|
||||
|
||||
Mark Jessop - 2019-03-02
|
||||
|
||||
auto_rx uses the [dft_detect](https://github.com/rs1729/RS/blob/master/scan/dft_detect.c) utility, developed by rs1729. This utility accepts fm-demodulated samples (fs=48kHz, unsigned 8-bit) and correlates the incoming samples against a set of known radiosonde packet headers.
|
||||
|
||||
The incoming sample stream is considered to contain a radiosonde signal when the correlation score for a particular radiosonde type header exceeds a threshold, defined in the code [here](https://github.com/rs1729/RS/blob/master/scan/dft_detect.c#L92).
|
||||
These threshold values need to be set selected such that detection occurs at the lowest signal SNR as possible, though without producing false positive, or mis-detections.
|
||||
|
||||
To assist with selection of these thresholds, an investigation was performed where dft_detect was run over a set of radiosonde samples, with calibrated SNRs (Eb/N0) between 5dB and 20 dB (using generate_lowsnr.py). dft_detect was also run over noise samples.
|
||||
|
||||
The following analysis was performed using the rtlfm FM demodulator, using a 22 kHz output sample rate. The following command-line was used to process samples:
|
||||
```
|
||||
$ cat test_file.bin | csdr convert_f_s16 | ./tsrc - - 1.8333 | csdr convert_s16_f | csdr shift_addition_cc -0.25000 2>/dev/null | csdr convert_f_u8 | ./rtl_fm_stdin -M fm -f 401000000 -F9 -s 22000 2>/dev/null| sox -t raw -r 22000 -e s -b 16 -c 1 - -r 48000 -b 8 -t wav - highpass 20 2>/dev/null |../dft_detect --debug
|
||||
```
|
||||
|
||||
## dft_detect 'noise floor'
|
||||
Correlation of the known headers with noise (no radiosonde signals present) will result in a non-zero correlation score. These scores will set the minimum setting possible for the threshold, to avoid false detections.
|
||||
|
||||
dft_detect was run over a set of noise samples (noise1_96k_float.bin, noise2_96k_float.bin, noise3_96k_float.bin), which represent commonly observed background noise on my auto_rx station. The three samples have slightly different noise characteristics and contain different types of local interference.
|
||||
|
||||
The worst-case (highest) correlation scores out of the three samples are shown for each radiosonde type:
|
||||
|
||||
Sonde Type | Worse-Case Correlation Score
|
||||
-----------|-----------------------------
|
||||
DFM | 0.3050
|
||||
RS41 | 0.3261
|
||||
RS92 | 0.2961
|
||||
M10 | 0.6106
|
||||
LMS6 | 0.3157
|
||||
|
||||
We can see that the M10 radiosonde has a very high correlation score against noise. This is due to the M10 header being about half the size of the other radiosondes.
|
||||
|
||||
Note that as we don't have any LMS6 IQ samples at the moment, we will not be including LMS6 results in any further analysis.
|
||||
|
||||
## Threshold Experimentings
|
||||
Using these minimum threshold figures, the next aim is to determine thresholds which result in no false detections of other sonde types. We can do this by running dft_detect over all the calibrated-SNR samples and observing the resultant correlation scores.
|
||||
|
||||
We start out with thresholds set 0.1 higher than the 'correlation noise floor'. The following plots show the correlation score for each set of radiosonde samples (RS41, RS92, DFM, M10), vs Eb/N0. The dotted points indicate a point in which the correlation score exceeds the threshold, and hence we considered a sonde type to be detected.
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
We can see that while detection of the wanted radiosonde type occurs from a very low SNR, there are also a large number of mis-detections, particularly at higher SNRs. By looking at the maximum observed correlation scores for each dataset, we can set a new, higher threshold to avoid mis-detections.
|
||||
|
||||
Dataset / Sonde Type: | RS41 | RS92 | DFM | M10
|
||||
----------------------|------|------|-----|-----
|
||||
RS41 | 0.9162 | 0.4390 | 0.5299 | 0.7209
|
||||
RS92 | 0.4070 | 0.9050 | 0.4836 | 0.7151
|
||||
DFM | 0.4352 | 0.4421 | 0.9037 | 0.5781
|
||||
M10 | 0.3519 | 0.3585 | 0.3707 | 0.9810
|
||||
|
||||
Based on the above table, we can define new thresholds at just above the highest 'unwanted' detection score:
|
||||
* RS41: 0.53
|
||||
* RS92: 0.54
|
||||
* DFM: 0.62
|
||||
* M10: 0.75
|
||||
|
||||
We can now re-generate the above plots, with the new threshold values:
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
Much improvement! From the plots, we can now determine the minimum SNR at which detection is possible:
|
||||
* RS41: 8.5 dB
|
||||
* RS92: 10.5 dB
|
||||
* DFM: 10.5 dB
|
||||
* M10: 5.0 dB (There may be a SNR miscalibration issue with the M10 samples)
|
||||
|
||||
## Future Work
|
||||
* It would probably be advantagous to tweak these detection thresholds such that the minimum detection threshold matches the minimum SNR for valid telemetry decoding - but that's for another time!
|
||||
* Work out how the current peak-detection threshold SNR lines up with the above Eb/N0 figures for different radiosonde types. Can we even do reliable peak detection (i.e. not too many peaks) with signals with such a low SNR?
|
||||
|
|
@ -0,0 +1,61 @@
|
|||
#!/usr/bin/env python
|
||||
#
|
||||
# Investigate dft_detect correlation outputs - Plotting
|
||||
# Plots the CSV output generated by dft_detect_scores.py
|
||||
#
|
||||
# Copyright (C) 2018 Mark Jessop <vk5qi@rfhead.net>
|
||||
# Released under GNU GPL v3 or later
|
||||
#
|
||||
|
||||
import glob
|
||||
import argparse
|
||||
import os
|
||||
import sys
|
||||
import time
|
||||
import subprocess
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
parser = argparse.ArgumentParser()
|
||||
parser.add_argument("-f", "--filename", type=str, default=None, help="CSV Filename to plot.")
|
||||
parser.add_argument("-t", "--title", type=str, default="Correlation Scores", help="Plot Title")
|
||||
parser.add_argument("-s", "--save", type=str, default=None, help="Save plot to file.")
|
||||
parser.add_argument("-v", "--verbose", action='store_true', default=False, help="Show additional debug info.")
|
||||
args = parser.parse_args()
|
||||
|
||||
|
||||
|
||||
# Sonde types to show on the plot.
|
||||
# Lines are shown as dashed when below their threshold value.
|
||||
TYPES = {
|
||||
'RS41':{'color':'C0', 'threshold':0.53},
|
||||
'RS92':{'color':'C1','threshold':0.54},
|
||||
'DFM9':{'color':'C2','threshold':0.62},
|
||||
'M10':{'color':'C3','threshold':0.75},
|
||||
'LMS6':{'color':'C4','threshold':0.70}
|
||||
}
|
||||
|
||||
|
||||
data = np.genfromtxt(args.filename,dtype=float, delimiter=',', names=True,usecols=(0,1,2,3,4,5))
|
||||
|
||||
|
||||
plt.figure(figsize=(10,5))
|
||||
|
||||
for _type in TYPES:
|
||||
# Plot
|
||||
plt.plot(data['snr'], data[_type], color=TYPES[_type]['color'])
|
||||
plt.scatter(data['snr'][ data[_type]>=TYPES[_type]['threshold'] ], data[_type][ data[_type]>=TYPES[_type]['threshold'] ], label=_type, color=TYPES[_type]['color'])
|
||||
|
||||
_max = np.max(data[_type])
|
||||
print("%s Maximum Correlation Score: %.4f" % (_type, _max))
|
||||
|
||||
plt.title("dft_detect Correlation - " + args.title)
|
||||
plt.xlabel("Eb/N0 (db)")
|
||||
plt.ylabel("Correlation Score")
|
||||
plt.legend()
|
||||
plt.grid()
|
||||
|
||||
if args.save != None:
|
||||
plt.savefig(args.save, bbox_inches='tight')
|
||||
|
||||
plt.show()
|
||||
|
|
@ -173,7 +173,6 @@ _sample_fs = 96000.0 # Sample rate of input. Always 96k at the moment.
|
|||
|
||||
# RS41
|
||||
_fm_rate = 15000
|
||||
_decode_timeout = 15
|
||||
# Calculate the necessary conversions
|
||||
_rtlfm_oversampling = 8.0 # Viproz's hacked rtl_fm oversamples by 8x.
|
||||
_shift = -2.0*_fm_rate/_sample_fs # rtl_fm tunes 'up' by rate*2, so we need to shift the signal down by this amount.
|
||||
|
|
@ -188,7 +187,7 @@ else:
|
|||
_resample_command = ""
|
||||
|
||||
_demod_command = "| %s csdr shift_addition_cc %.5f 2>/dev/null | csdr convert_f_u8 |" % (_resample_command, _shift)
|
||||
_demod_command += " timeout %d ./rtl_fm_stdin -M fm -f 401000000 -F9 -s %d 2>/dev/null|" % (_decode_timeout, int(_fm_rate))
|
||||
_demod_command += " ./rtl_fm_stdin -M fm -f 401000000 -F9 -s %d 2>/dev/null|" % (int(_fm_rate))
|
||||
_demod_command += " sox -t raw -r %d -e s -b 16 -c 1 - -r 48000 -b 8 -t wav - lowpass 2600 2>/dev/null |" % int(_fm_rate)
|
||||
|
||||
processing_type['rs41_rtlfm'] = {
|
||||
|
|
@ -205,7 +204,6 @@ processing_type['rs41_rtlfm'] = {
|
|||
|
||||
# RS92
|
||||
_fm_rate = 12000
|
||||
_decode_timeout = 15
|
||||
# Calculate the necessary conversions
|
||||
_rtlfm_oversampling = 8.0 # Viproz's hacked rtl_fm oversamples by 8x.
|
||||
_shift = -2.0*_fm_rate/_sample_fs # rtl_fm tunes 'up' by rate*2, so we need to shift the signal down by this amount.
|
||||
|
|
@ -220,7 +218,7 @@ else:
|
|||
_resample_command = ""
|
||||
|
||||
_demod_command = "| %s csdr shift_addition_cc %.5f 2>/dev/null | csdr convert_f_u8 |" % (_resample_command, _shift)
|
||||
_demod_command += " timeout %d ./rtl_fm_stdin -M fm -f 401000000 -F9 -s %d 2>/dev/null|" % (_decode_timeout, int(_fm_rate))
|
||||
_demod_command += " ./rtl_fm_stdin -M fm -f 401000000 -F9 -s %d 2>/dev/null|" % (int(_fm_rate))
|
||||
_demod_command += " sox -t raw -r %d -e s -b 16 -c 1 - -r 48000 -b 8 -t wav - lowpass 2500 highpass 20 2>/dev/null |" % int(_fm_rate)
|
||||
|
||||
|
||||
|
|
@ -236,7 +234,6 @@ processing_type['rs92_rtlfm'] = {
|
|||
|
||||
# DFM
|
||||
_fm_rate = 20000
|
||||
_decode_timeout = 15
|
||||
# Calculate the necessary conversions
|
||||
_rtlfm_oversampling = 8.0 # Viproz's hacked rtl_fm oversamples by 8x.
|
||||
_shift = -2.0*_fm_rate/_sample_fs # rtl_fm tunes 'up' by rate*2, so we need to shift the signal down by this amount.
|
||||
|
|
@ -251,7 +248,7 @@ else:
|
|||
_resample_command = ""
|
||||
# For some reason the DFM sample breaks type conversion - multiplying it by 0.9 seems to fix it.
|
||||
_demod_command = "| csdr gain_ff 0.90 | %s csdr shift_addition_cc %.5f 2>/dev/null | csdr convert_f_u8 |" % (_resample_command, _shift)
|
||||
_demod_command += " timeout %d ./rtl_fm_stdin -M fm -f 401000000 -F9 -s %d 2>/dev/null|" % (_decode_timeout, int(_fm_rate))
|
||||
_demod_command += " ./rtl_fm_stdin -M fm -f 401000000 -F9 -s %d 2>/dev/null|" % (int(_fm_rate))
|
||||
_demod_command += " sox -t raw -r %d -e s -b 16 -c 1 - -r 48000 -b 8 -t wav - highpass 20 lowpass 2000 2>/dev/null |" % int(_fm_rate)
|
||||
|
||||
processing_type['dfm_rtlfm'] = {
|
||||
|
|
@ -265,7 +262,6 @@ processing_type['dfm_rtlfm'] = {
|
|||
|
||||
# M10
|
||||
_fm_rate = 22000
|
||||
_decode_timeout = 15
|
||||
# Calculate the necessary conversions
|
||||
_rtlfm_oversampling = 8.0 # Viproz's hacked rtl_fm oversamples by 8x.
|
||||
_shift = -2.0*_fm_rate/_sample_fs # rtl_fm tunes 'up' by rate*2, so we need to shift the signal down by this amount.
|
||||
|
|
@ -280,7 +276,7 @@ else:
|
|||
_resample_command = ""
|
||||
|
||||
_demod_command = "| %s csdr shift_addition_cc %.5f 2>/dev/null | csdr convert_f_u8 |" % (_resample_command, _shift)
|
||||
_demod_command += " timeout %d ./rtl_fm_stdin -M fm -f 401000000 -F9 -s %d 2>/dev/null|" % (_decode_timeout, int(_fm_rate))
|
||||
_demod_command += " ./rtl_fm_stdin -M fm -f 401000000 -F9 -s %d 2>/dev/null|" % (int(_fm_rate))
|
||||
_demod_command += " sox -t raw -r %d -e s -b 16 -c 1 - -r 48000 -b 8 -t wav - highpass 20 2>/dev/null |" % int(_fm_rate)
|
||||
|
||||
processing_type['m10_rtlfm'] = {
|
||||
|
|
@ -293,7 +289,7 @@ processing_type['m10_rtlfm'] = {
|
|||
|
||||
# RS_Detect
|
||||
_fm_rate = 22000
|
||||
_decode_timeout = 15
|
||||
#_fm_rate = 15000
|
||||
# Calculate the necessary conversions
|
||||
_rtlfm_oversampling = 8.0 # Viproz's hacked rtl_fm oversamples by 8x.
|
||||
_shift = -2.0*_fm_rate/_sample_fs # rtl_fm tunes 'up' by rate*2, so we need to shift the signal down by this amount.
|
||||
|
|
@ -308,7 +304,7 @@ else:
|
|||
_resample_command = ""
|
||||
|
||||
_demod_command = "| %s csdr shift_addition_cc %.5f 2>/dev/null | csdr convert_f_u8 |" % (_resample_command, _shift)
|
||||
_demod_command += " timeout %d ./rtl_fm_stdin -M fm -f 401000000 -F9 -s %d 2>/dev/null|" % (_decode_timeout, int(_fm_rate))
|
||||
_demod_command += " ./rtl_fm_stdin -M fm -f 401000000 -F9 -s %d 2>/dev/null|" % (int(_fm_rate))
|
||||
_demod_command += " sox -t raw -r %d -e s -b 16 -c 1 - -r 48000 -b 8 -t wav - highpass 20 2>/dev/null |" % int(_fm_rate)
|
||||
|
||||
processing_type['rs_detect_rtlfm'] = {
|
||||
|
|
@ -321,7 +317,7 @@ processing_type['rs_detect_rtlfm'] = {
|
|||
|
||||
# DFT_Detect
|
||||
_fm_rate = 22000
|
||||
_decode_timeout = 15
|
||||
#_fm_rate = 15000
|
||||
# Calculate the necessary conversions
|
||||
_rtlfm_oversampling = 8.0 # Viproz's hacked rtl_fm oversamples by 8x.
|
||||
_shift = -2.0*_fm_rate/_sample_fs # rtl_fm tunes 'up' by rate*2, so we need to shift the signal down by this amount.
|
||||
|
|
@ -336,7 +332,7 @@ else:
|
|||
_resample_command = ""
|
||||
|
||||
_demod_command = "| %s csdr shift_addition_cc %.5f 2>/dev/null | csdr convert_f_u8 |" % (_resample_command, _shift)
|
||||
_demod_command += " timeout %d ./rtl_fm_stdin -M fm -f 401000000 -F9 -s %d 2>/dev/null|" % (_decode_timeout, int(_fm_rate))
|
||||
_demod_command += " ./rtl_fm_stdin -M fm -f 401000000 -F9 -s %d 2>/dev/null|" % (int(_fm_rate))
|
||||
_demod_command += " sox -t raw -r %d -e s -b 16 -c 1 - -r 48000 -b 8 -t wav - highpass 20 2>/dev/null |" % int(_fm_rate)
|
||||
|
||||
processing_type['dft_detect_rtlfm'] = {
|
||||
|
|
@ -348,6 +344,7 @@ processing_type['dft_detect_rtlfm'] = {
|
|||
}
|
||||
|
||||
|
||||
|
||||
def run_analysis(mode, file_mask=None, shift=0.0, verbose=False, log_output = None):
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -17,7 +17,8 @@ static int option_verbose = 0, // ausfuehrliche Anzeige
|
|||
option_inv = 0, // invertiert Signal
|
||||
option_dc = 0,
|
||||
option_silent = 0,
|
||||
wavloaded = 0;
|
||||
wavloaded = 0,
|
||||
option_debug = 0;
|
||||
static int wav_channel = 0; // audio channel: left
|
||||
|
||||
|
||||
|
|
@ -90,11 +91,11 @@ typedef struct {
|
|||
#define idxAB 7
|
||||
#define idxRS 8
|
||||
static rsheader_t rs_hdr[Nrs] = {
|
||||
{ 2500, 0, 0, dfm_header, 1.0, 0.0, 0.65, NULL, "DFM9", 2}, // DFM6: -2 (unsigned)
|
||||
{ 4800, 0, 0, rs41_header, 0.5, 0.0, 0.70, NULL, "RS41", 3},
|
||||
{ 4800, 0, 0, rs92_header, 0.5, 0.0, 0.70, NULL, "RS92", 4},
|
||||
{ 4800, 0, 0, lms6_header, 1.0, 0.0, 0.70, NULL, "LMS6", 8},
|
||||
{ 9600, 0, 0, m10_header, 1.0, 0.0, 0.76, NULL, "M10", 5},
|
||||
{ 2500, 0, 0, dfm_header, 1.0, 0.0, 0.62, NULL, "DFM9", 2}, // DFM6: -2 (unsigned) - was 0.65
|
||||
{ 4800, 0, 0, rs41_header, 0.5, 0.0, 0.53, NULL, "RS41", 3}, // was 0.70
|
||||
{ 4800, 0, 0, rs92_header, 0.5, 0.0, 0.54, NULL, "RS92", 4}, // was 0.70
|
||||
{ 4800, 0, 0, lms6_header, 1.0, 0.0, 0.70, NULL, "LMS6", 8}, //
|
||||
{ 9600, 0, 0, m10_header, 1.0, 0.0, 0.75, NULL, "M10", 5},
|
||||
{ 5800, 0, 0, c34_preheader, 1.5, 0.0, 0.80, NULL, "C34C50", 9}, // C34/C50 2900 Hz tone
|
||||
{ 9600, 0, 0, imet_preamble, 0.5, 0.0, 0.80, NULL, "IMET", 6}, // IMET1AB=7, IMET1RS=8
|
||||
{ 9600, 0, 0, imet1ab_header, 1.0, 0.0, 0.80, NULL, "IMET1AB", 6},
|
||||
|
|
@ -658,7 +659,7 @@ int main(int argc, char **argv) {
|
|||
|
||||
int j;
|
||||
int k, K;
|
||||
float mv[Nrs];
|
||||
float mv[Nrs], mv_max_arr[Nrs];
|
||||
unsigned int mv_pos[Nrs], mv0_pos[Nrs];
|
||||
int mp[Nrs];
|
||||
|
||||
|
|
@ -691,6 +692,10 @@ int main(int argc, char **argv) {
|
|||
else if ( (strcmp(*argv, "-s") == 0) || (strcmp(*argv, "--silent") == 0) ) {
|
||||
option_silent = 1;
|
||||
}
|
||||
else if ( (strcmp(*argv, "--debug") == 0) ) {
|
||||
option_debug = 1;
|
||||
option_silent = 1;
|
||||
}
|
||||
else if ( (strcmp(*argv, "-t") == 0) || (strcmp(*argv, "--time") == 0) ) {
|
||||
++argv;
|
||||
if (*argv) tl = atof(*argv);
|
||||
|
|
@ -736,6 +741,7 @@ int main(int argc, char **argv) {
|
|||
mv[j] = 0;
|
||||
mv_pos[j] = 0;
|
||||
mp[j] = 0;
|
||||
mv_max_arr[j] = 0;
|
||||
}
|
||||
|
||||
k = 0;
|
||||
|
|
@ -852,7 +858,9 @@ int main(int argc, char **argv) {
|
|||
}
|
||||
}
|
||||
}
|
||||
else header_found = 1;
|
||||
else{
|
||||
if (!option_debug) header_found = 1; // Might need to disable this for debug mode.
|
||||
}
|
||||
|
||||
if (header_found) {
|
||||
if (!option_silent) {
|
||||
|
|
@ -866,12 +874,25 @@ int main(int argc, char **argv) {
|
|||
}
|
||||
}
|
||||
|
||||
// Store maximum correlation scores
|
||||
for (j = 0; j < Nrs; j++){
|
||||
if (mv[j] > mv_max_arr[j]){
|
||||
mv_max_arr[j] = mv[j];
|
||||
}
|
||||
}
|
||||
|
||||
if (header_found) break;
|
||||
header_found = 0;
|
||||
for (j = 0; j < Nrs; j++) mv[j] = 0.0;
|
||||
}
|
||||
|
||||
ende:
|
||||
if (option_debug) {
|
||||
//fprintf(stdout, "\nCorrelation scores for all types:\n");
|
||||
for (j=0; j<Nrs; j++){
|
||||
fprintf(stdout, "%s, %.4f, %.4f\n", rs_hdr[j].type, mv_max_arr[j], rs_hdr[j].thres);
|
||||
}
|
||||
}
|
||||
free_buffers();
|
||||
fclose(fp);
|
||||
|
||||
|
|
|
|||
Ładowanie…
Reference in New Issue