Add some benchmarking code

benchmarking/README.md

@@ -0,0 +1,91 @@
+# Wenet Performance Benchmarking
+Some attempts at benchmarking the performance of the Wenet decode chains, so we know if we've broken things in the future.
+
+
+## Setup
+* Have the wenet rx code built. `fsk_demod` and `drs232_ldpc` should exist within wenet/rx/  (../rx/)
+* Have csdr available on path. Doesnt really matter which fork, we just need the convert_u8_f and convert_f_u8 utils.
+
+Make some directories
+```
+mkdir samples
+mkdir generated
+```
+
+We also need numpy available for python3. You could get that via system packages, or create a venv and install with pip.
+
+## Test Samples
+To generate the low SNR files, we need a very high SNR (not overloading though) original sample.
+
+If we have a wenet payload available, we can just dump some samples from rtl_sdr, e.g.:
+
+```
+rtl_sdr -s 921416 -f 443298440 -g 5 - > test_samples.cu8
+```
+
+We need this test sample in float32 format, which we can do using csdr:
+```
+cat test_samples.cu8 | csdr convert_u8_f > samples/wenet_sample_fs921416_float.bin
+```
+
+For 'traditional' Wenet (115177 baud, RS232 framing), a suitable sample set (~95s of received packets) is available here: https://www.dropbox.com/scl/fi/plazem0luo37l2dujbwuo/wenet_sample_fs921416Hz.cu8?rlkey=m4jftwmbazok9ry9kimhpd6kl&dl=0
+(this still needs to be converted to float32 as above).
+
+## Generating Low SNR Samples
+Check generate_lowsnr.py for the list of files to be used as source material for low-snr generation.
+
+Then, run: python generate_lowsnr.py
+
+You should now have a bunch of files in the generated directory.
+
+Note - this can be quite a lot of data! 
+
+Note that the value in dB in the filenames is Eb/N0, so effectively snr-per-bit, normalised for baud rate.
+
+## Running demod tests
+
+Can do a quick check to make sure the highest SNR sample (which should have very good decode) works by running:
+```
+% python test_demod.py -m wenet_rs232_demod --quick
+Command: cat ./generated/wenet_sample_fs921416_float_20.0dB.bin | csdr convert_f_u8 | ../rx/fsk_demod --cu8 -s --stats=100 2 921416 115177 - - 2> stats.txt | ../rx/drs232_ldpc - - 2> /dev/null  |  wc -c
+wenet_sample_fs921416_float_20.0dB.bin, 530688, 11.156
+```
+Output consists of:
+* filename
+* number of bytes received (only packets with valid CRC are output from the decoder)
+* Time taken to run the decode
+
+Our performance metric is the number of bytes received.
+
+We can then go ahead and run the tests using the full set of generated samples:
+
+```
+% python test_demod.py -m wenet_rs232_demod 
+Command: cat ./generated/wenet_sample_fs921416_float_05.0dB.bin | csdr convert_f_u8 | ../rx/fsk_demod --cu8 -s --stats=100 2 921416 115177 - - 2> stats.txt | ../rx/drs232_ldpc - - 2> /dev/null  |  wc -c
+wenet_sample_fs921416_float_05.0dB.bin, 0, 12.870
+wenet_sample_fs921416_float_05.5dB.bin, 0, 13.500
+wenet_sample_fs921416_float_06.0dB.bin, 0, 12.332
+wenet_sample_fs921416_float_06.5dB.bin, 0, 12.375
+wenet_sample_fs921416_float_07.0dB.bin, 0, 13.308
+wenet_sample_fs921416_float_07.5dB.bin, 1280, 19.061
+wenet_sample_fs921416_float_08.0dB.bin, 32512, 20.793
+wenet_sample_fs921416_float_08.5dB.bin, 298240, 12.358
+wenet_sample_fs921416_float_09.0dB.bin, 503040, 12.685
+wenet_sample_fs921416_float_09.5dB.bin, 528128, 13.180
+wenet_sample_fs921416_float_10.0dB.bin, 529920, 14.041
+wenet_sample_fs921416_float_10.5dB.bin, 530176, 12.373
+wenet_sample_fs921416_float_11.0dB.bin, 530432, 15.980
+wenet_sample_fs921416_float_11.5dB.bin, 530176, 12.871
+wenet_sample_fs921416_float_12.0dB.bin, 530432, 13.491
+wenet_sample_fs921416_float_12.5dB.bin, 530432, 12.516
+wenet_sample_fs921416_float_13.0dB.bin, 530432, 12.867
+```
+
+Things to look at:
+* What Eb/N0 the number of received bytes starts to rise. With LDPC FEC it's a fairly quick increase from nothing to complete decodes. In the above case, our 50% packet-error-rate point is around 8.5 dB.
+* How long it takes to run the decode chain. We see a slight increase in runtime around the weak-snr point, when we get a lot of unique-word detections, but where the LDPC decoder runs to maximum iterations (5) without a successful decode.
+
+Currently there are demod chain tests for:
+* `wenet_rs232_demod` - Wenet 'traditional' (v1?), 115177 baud, RS232 framing, with complex u8 samples going into fsk_demod
+* `wenet_rs232_demod_c16` - Same as above, but feeding complex signed-16-bit samoples into fsk_demod (should give the same results).
+

benchmarking/generate_lowsnr.py

@@ -0,0 +1,124 @@
+#!/usr/bin/env python
+#
+#   Generate Noisy Sonde Samples, with a calibrated Eb/No
+#
+#   Run from ./scripts/ with
+#   $ python generate_lowsnr.py
+#
+#   The generated files will end up in the 'generated' directory.
+#
+#   Copyright (C) 2018  Mark Jessop <vk5qi@rfhead.net>
+#   Released under GNU GPL v3 or later
+#
+import numpy as np
+import os
+
+
+# Where to find the samples files.
+# These are all expected to be 96khz float (dtype='c8') files.
+SAMPLE_DIR = "./samples"
+
+# Directory to output generated files
+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(5,15.0,0.5)
+
+# Normalise the samples to +/- 1.0! 
+# If we don't do this, bad things can happen later down the track...
+NORMALISE = True
+
+# List of samples
+# [filename, baud_date, threshold, sample_rate]
+# filename = string, without path
+# baud_rate = integer
+# threshold = threshold for calculating variance. Deterimined by taking 20*np.log10(np.abs(data)) and looking for packets.
+# sample_rate = input file sample rate.
+
+SAMPLES = [
+    ['wenet_sample_fs921416_float.bin', 115200, -100, 921416], # No threshold set, as signal is continuous.
+]
+
+
+
+def load_sample(filename):
+    _filename = os.path.join(SAMPLE_DIR, filename)
+    return np.fromfile(_filename, dtype='c8')
+
+
+def save_sample(data, filename):
+    _filename = os.path.join(GENERATED_DIR, filename)
+    # We have to make sure to convert to complex64..
+    data.astype(dtype='c8').tofile(_filename)
+
+    # TODO: Allow saving as complex s16 - see view solution here: https://stackoverflow.com/questions/47086134/how-to-convert-a-numpy-complex-array-to-a-two-element-float-array
+
+
+
+def calculate_variance(data, threshold=-100.0):
+    # Calculate the variance of a set of radiosonde samples.
+    # Optionally use a threshold to limit the sample the variance 
+    # is calculated over to ones that actually have sonde packets in them.
+
+    _data_log = 20*np.log10(np.abs(data))
+
+    return np.var(data[_data_log>threshold])
+
+
+def add_noise(data, variance, baud_rate, ebno, fs=96000,  bitspersymbol=1.0):
+    # Add calibrated noise to a sample.
+
+    # Calculate Eb/No in linear units.
+    _ebno = 10.0**((ebno)/10.0)
+
+    # Calculate the noise variance we need to add
+    _noise_variance = variance*fs/(baud_rate*_ebno*bitspersymbol)
+
+    # Generate complex random samples
+    _rand_i = np.sqrt(_noise_variance/2.0)*np.random.randn(len(data))
+    _rand_q = np.sqrt(_noise_variance/2.0)*np.random.randn(len(data))
+
+    _noisy = (data + (_rand_i + 1j*_rand_q))
+
+    if NORMALISE:
+        print("Normalised to 1.0")
+        return _noisy/np.max(np.abs(_noisy))
+    else:
+        return _noisy
+
+
+
+
+if __name__ == '__main__':
+
+    for _sample in SAMPLES:
+        # Extract the stuff we need from the entry.
+        _source = _sample[0]
+        _baud_rate = _sample[1]
+        _threshold = _sample[2]
+        _fs = _sample[3]
+
+        print("Generating samples for: %s" % _source)
+
+        # Read in source file.
+        _data = load_sample(_source)
+
+        # Calculate variance
+        _var = calculate_variance(_data, _threshold)
+        print("Calculated Variance: %.5f" % _var)
+
+        # Now loop through the ebno's and generate the output.
+        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] + "_%04.1fdB"%ebno + ".bin"
+
+            save_sample(_data_noise, _out_file)
+            print("Saved file: %s" % _out_file)
+
+
+
+
+
+

benchmarking/test_demod.py

@@ -0,0 +1,148 @@
+#!/usr/bin/env python
+#
+#   Run a set of files through a processing and decode chain, and handle the output.
+#
+#   Copyright (C) 2018  Mark Jessop <vk5qi@rfhead.net>
+#   Released under GNU GPL v3 or later
+#
+#   Refer to the README.md in this directory for instructions on use.
+#
+import glob
+import argparse
+import os
+import sys
+import time
+import traceback
+import subprocess
+
+
+# Dictionary of available processing types.
+
+processing_type = {
+
+    # Wenet, RS232 modulation
+    # Convert to u8 using csdr, then pipe into fsk_demod, then drs232_ldpc.
+    # Count bytes at the output as a metric of performance.
+    'wenet_rs232_demod': {
+        'demod': '| csdr convert_f_u8 | ../rx/fsk_demod --cu8 -s --stats=100 2 921416 115177 - - 2> stats.txt | ',
+        'decode': '../rx/drs232_ldpc - - 2> /dev/null ',
+        "post_process" : " |  wc -c", #
+        'files' : "./generated/wenet_sample_fs921416*.bin"
+    },
+    'wenet_rs232_demod_c16': {
+        'demod': '| csdr convert_f_s16 | ../rx/fsk_demod --cs16 -s --stats=100 2 921416 115177 - - 2> stats.txt | ',
+        'decode': '../rx/drs232_ldpc - - 2> /dev/null ',
+        "post_process" : " |  wc -c", #
+        'files' : "./generated/wenet_sample_fs921416*.bin"
+    },
+}
+
+def run_analysis(mode, file_mask=None, shift=0.0, verbose=False, log_output = None, dry_run = False, quick=False, show=False):
+
+
+    _mode = processing_type[mode]
+
+    # If we are not supplied with a file mask, use the defaults.
+    if file_mask is None:
+        file_mask = _mode['files']
+
+    # Get the list of files.
+    _file_list = glob.glob(file_mask)
+    if len(_file_list) == 0:
+        print("No files found matching supplied path.")
+        return
+
+    # Sort the list of files.
+    _file_list.sort()
+
+    # If we are only running a quick test, just process the last file in the list.
+    if quick:
+        _file_list = [_file_list[-1]]
+
+    _first = True
+
+    # Calculate the frequency offset to apply, if defined.
+    _shiftcmd = "| csdr shift_addition_cc %.5f 2>/dev/null" % (shift/96000.0)
+
+    if log_output is not None:
+        _log = open(log_output,'w')
+
+    # Iterate over the files in the supplied list.
+    for _file in _file_list:
+
+        # Generate the command to run.
+        _cmd = "cat %s "%_file 
+
+        # Add in an optional frequency error if supplied.
+        if shift != 0.0:
+            _cmd += _shiftcmd
+
+        # Add on the rest of the demodulation and decoding commands.
+        _cmd += _mode['demod'] + _mode['decode'] 
+        
+        if args.show:
+            _cmd += " | head -n 10"
+        else:
+            _cmd += _mode['post_process']
+
+
+        if _first or dry_run:
+            print("Command: %s" % _cmd)
+            _first = False
+
+        if dry_run:
+            continue
+
+        # Run the command.
+        try:
+            _start = time.time()
+            _output = subprocess.check_output(_cmd, shell=True, stderr=None)
+            _output = _output.decode()
+        except:
+            #traceback.print_exc()
+            _output = "error"
+
+        _runtime = time.time() - _start
+
+        _result = "%s, %s, %.3f" % (os.path.basename(_file), _output.strip(), _runtime)
+
+        print(_result)
+        if log_output is not None:
+            _log.write(_result + '\n')
+
+        if verbose:
+            print("Runtime: %.1d" % _runtime)
+
+    if log_output is not None:
+        _log.close()
+
+
+
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser()
+    parser.add_argument("-m", "--mode", type=str, default="rs41_fsk_demod_soft", help="Operation mode.")
+    parser.add_argument("-f", "--files", type=str, default=None, help="Glob-path to files to run over.")
+    parser.add_argument("-v", "--verbose", action='store_true', default=False, help="Show additional debug info.")
+    parser.add_argument("-d", "--dry-run", action='store_true', default=False, help="Show additional debug info.")
+    parser.add_argument("--shift", type=float, default=0.0, help="Shift the signal-under test by x Hz. Default is 0.")
+    parser.add_argument("--batch", action='store_true', default=False, help="Run all tests, write results to results directory.")
+    parser.add_argument("--quick", action='store_true', default=False, help="Only process the last sample file in the list (usually the strongest). Useful for checking the demodulators are still working.")
+    parser.add_argument("--show", action='store_true', default=False, help="Show the first few lines of output, instead of running the post-processing step.")
+    args = parser.parse_args()
+
+    # Check the mode is valid.
+    if args.mode not in processing_type:
+        print("Error - invalid operating mode.")
+        print("Valid Modes: %s" % str(processing_type.keys()))
+        sys.exit(1)
+
+
+    batch_modes = []
+
+    if args.batch:
+        for _mode in batch_modes:
+            _log_name = "./results/" + _mode + ".txt"
+            run_analysis(_mode, file_mask=None, shift=args.shift, verbose=args.verbose, log_output=_log_name, dry_run=args.dry_run, quick=args.quick, show=args.show)
+    else:
+        run_analysis(args.mode, args.files, shift=args.shift, verbose=args.verbose, dry_run=args.dry_run, quick=args.quick, show=args.show)