SDR Utilities

I’ve been experimenting with an SDR dongle to see how it can be used as a 1.5GHz scalar network analyser, as a time domain reflectometer, and to estimate digital signal edge speeds. While doing that I’ve developed several general-purpose command-line utilities which capture and post-process power spectrums. These utilities are distributed with a MIT licence, and the code is at the bottom of this note.

  • rtlscan: scans, measures and displays a raw scalar power spectrum, and saves the power spectrum in a CSV file
  • rtlplot: displays one or more raw power spectrums
  • rtldiff: displays one or more calibrated power spectrums, i.e. the difference between multiple raw spectrums and a single reference or calibration spectrum
  • rtltdr: displays the impulse response implicit in a calibrated power spectrum

The graphs can have linear or logarithmic frequency axis, can be zoomed, panned and printed, and have text annotations added. The interaction with the SDR dongle is via Kyle Keen’s rtl_power program, so its capabilities and limitations are reflected in these utilities.

As is traditional, these are an unfinished work; I already know of the next change I would make, but I have not verified that it would be an improvement. Nonetheless, the utilities are usable.

The remainder of this post outlines the general workflow of using the various utilities to make measurements, example parameters for each utility, and prerequisites and installation.

General Workflow

Typically this would be:

  1. sanity check power levels to ensure they are within the dongle’s limited dynamic range: use rtlscan with wide frequency steps and a short integration period. If necessary, insert attenuators or adjust gain and repeat this sanity check
  2. measure the “reference” power spectrum used to “calibrate” the tool: attach the signal source directly to the dongle without the unit-under-test (UUT), and use rtlscan to capture the spectrum in a CSV file
  3. measure the UUT’s response: insert the UUT and use rtlscan to measure the spectrum with exactly the same scan parameters as in step 2
  4. optionally compare spectrums: use rtlplot
  5. optionally display the frequency response of the UUT: use rtldiff based on the spectrums from steps 2 and 3
  6. optionally display the impulse response of the UUT: use rtltdr based on the spectrums from steps 2 and 3

Typically that workflow will create a large number of CSV files, one for each scan. The utilities generate filenames are based on the scan parameters, which enables easy filtering/selection of files using standard shell wildcard “globbing” techniques. Those filenames can be manually changed as and when convenient.


Scan Examples Using rtlscan

rtlscan is invoked with parameters which define a “base” filename plus the parameters of the scan done by rtl_power. The mandatory rtl_power arguments are “-f frequencies” and “-g gain”, whereas the optional parameters are “-i integrationInterval”, “-c cropFactor” which is normally left at 0.5.

The results are stored in a CSV file containing the measured frequency-power tuples, rtl_power’s interpretation of the scan parameters, and a copy of the command-line arguments for future reference. The file’s name is the concatenation of the “base” and the scan parameters, a technique which simplifies selecting multiple spectrums for post-processing by the other utilities. The filename can be changed whenever convenient. (Typography alert: normal shell command options are used: “-” is a single hyphen but ” — ” is actually two hyphens which this font “runs together”, sigh.)

The optional “––gaincheck gain” flag specifies a suitable combination of parameters for the quickly checking the input levels. The optional “––tdrrange range” parameter specifies a suitable combination of parameters for measuring the impulse response. Individual parameters can be overridden by explicitly re-defining them later in the command line. See “Level Check Examples” below for examples of interpreting spectrums to assess whether the power level is too high or too low.

The dongle’s front-end has an LNA, but changing its gain does not change rtl_power’s output. However, the rtl_power output does contain an uncalibrated definition of the gain, and this is used by the utilities to adjust the power shown in graphs. The effect is that there may not be exactly 10dB difference in power levels between scans with, for example, -g 0 and -g 10. The “––ignorerxgain” flag disables the power adjustment, which can occasionally be useful.

Examples when measuring the frequency response of a 2.02m long open-circuit stub transmission line:

  • sanity check scan: “rtlscan.py ––gaincheck 0 sanity” creates file sanity––gaincheck_0_-f_24M_1500M_10M_-1_-c_0.5_-g_0_-i_1.csv
  • reference/calibration scan: “rtlscan.py -f 30M:1500M:.5M -g 0 -i 240 stubReference” creates file stubReference_-f_30M_1500M_.5M_-1_-c_.5_-g_0_-i_240.csv
  • UUT scan: “rtlscan.py -f 30M:1500M:.5M -g 0 -i 60 stub2.02short” creates file stub2.02short_-f_30M_1500M_.5M_-1_-c_.5_-g_0_-i_60.csv

It can be seen that the parameters “-1 -c 0.5” have been automatically inserted when invoking rtl_power, and that more noise would have been removed from the calibration scan by averaging for 240s rather than 60s.


Plot Examples Using rtlplot

rtlplot is invoked with the filenames of the spectra to be plotted, plus optional parameters modifying the power and frequency (–ignorerxgain, –logf) axes, plus any annotations.

For some purposes, such as the open-stub transmission line filter where nulls occur at harmonics of the first null, it is conventional to have a linear frequency scale. For other purposes, such as the frequency response of a high/low/band pass filter, it is conventional to use a logarithmic frequency scale. The latter is selected by the ––logf flag.

The plot legend shows the full filename, which is acceptable while testing but would be unacceptably cryptic for published results. Simply renaming the files removes that problem. The scan parameters are not completely lost, because they are still embedded in a comment in the CSV file.

Consider a circuit consisting of a stub transmission line 2.02m long and either open or shorted at the far end. Given the reference and UUT scans processed by

“rtlplot.py -a 302,7.3,”What’s this step?” -a 1070,-3,”Wiggle?” stubR* stub2* “

Two annotations have been added highlighting noteworthy features. In the zoomed plot it can be seen that there is a slightly lower noise level on the reference spectrum, due to a longer integration period, “-i 240” vs “-i 60”.


UUT Response Examples Using rtldiff

rtldiff is the same as rtlplot, except that the first CSV file contains the reference/calibration spectrum, which is subtracted from every other files’ spectrum. Since no interpolation is attempted, the frequency points must be identical in all files. Practically this means all scans must have the same frequency (-f) and crop factors (-c) parameters.

Using the same scans processed by “rtldiff.py stubR* stub2* removes the wiggles and step

rtldiff.py stubR* stub2*

rtldiff.py stubR* stub2*

The reduced depth of the nulls at higher frequencies is implied in the rtlplot graphs, and is due either to inadequate frequency resolution or the dongle’s noise floor. The former can be reduced at the expense of more scanned points (rtlscan’s -f parameter); the latter can be slightly reduced by longer integration times (rtlscan’s -i parameter). Both options increase the scan time.


Level Check Examples

It is obvious when the signal level is below the noise floor, but less obvious when the dongle is overloaded by strong signals. These examples show how to quickly check for overload.

The setup is that a noise source is connected to the dongle with 0dB, 10dB and 20dB attenuator, and the dongle’s gain is set to odB and +10dB. The dongle’s gain is automatically included in the CSV filename, whereas the attenuator’s value is specified in the base filename, either “sanity0dB” or “sanity10dB” or “sanity20dB”. The five combinations are scanned using commands  “rtlscan.py –gaincheck 0 sanity10dB” and obvious permutations. The results are plotted with and without the ––ignorerxgain flag.

Observations:

  • the blue, green and cyan lines are not evenly spaced and have a differing shape: this implies overload and these settings should not be used
  • the red and purple lines correctly spaced ~10dB apart and the same shape: there is no sign of overload and these settings are valid
  • the yellow line is the same as the red and purple lines except below 100MHz: these settings could be used if a 100MHz high-pass filter was inserted. The yellow line ought to be 10dB higher than the purple line; the discrepancy is presumably due to a “sub-optimum” gain calibration figure

The conclusion is to use either:

  • an external 10dB attenuator and -g 0
  • an external 20dB attenuator and -g 0 (or -g 10 above 100MHz)

Prerequisites and Installation

I presume it is possible to use these utilities with other operating systems and SDR hardware, but I haven’t tried anything other than:

  • Xubuntu Linux. I can’t be bothered with Windows any longer, and since Linux exists there is no reason for me to use Apple Macs.
  • a cheap-and-cheerful SDR receiver containing an R802T2 front end and RTL2832U demodulator in a USB dongle

The utilities are written in Python, so Python 2.7 and modules numpy, scipy, matplotlib must be installed.

The SDR is controlled by Kyle Keen’s rtl_power utility. The version supplied with other programs or with Linux distributions is intolerant of errors in the -i parameter. I found it necessary to download the latest source from his github repository. I had no trouble compiling it.


 Code

Unfortunately I can’t upload a zipfile containing the code, so I’ll merely copy and paste it here. You will have to recopy and paste into the named files. Sorry about that.

 

rtlscan.py:

#! /usr/bin/env python
"""
MIT License

Copyright (c) 2016 Thomas Godfrey Gardner

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
"""

#v1.0

import sys
import getopt
from rtlutils import scanToFile, drawGraphFromFiles


def helpArgs(msg):
    if msg != "":
        print "**** " + msg + " ****"
    print "rtlscan [-c crop_factor] [-g tuner_gain] [-f frequencies] [-i integration_interval] [-h --help]"
    print "      [--gaincheck gain] [--tdrrange r] [--logf] [--ignorerxgain] filenameRoot"
    print "  Invoke rtl_power to do one scan, display the results, and save its raw output for processing "
    print "  by other utilities. The file's name is filenameRoot_options.csv, where options is the"
    print "  unix-friendly concatenation of the rtl_power options with ' ' and ':' characters replaced by '_'."
    print "  -c crop_factor            from rtl_power, range 0 to 1 (default 0.5)"
    print "  -g tuner_gain             from rtl_power"
    print "  -i integration_interval   from rtlPower (default 20)"
    print "  -f frequencies            from rtl_power, lower:upper:bin_size"
    print "  --gaincheck gain          -c .5 -g gain -i 1 -f 24M:1500M:10M, unless parameters overriden"
    print "  --tdrrange range          -c .5 -g 10 -i 20 -f 30M:1500M:X, unless parameters overriden, where X "
    print "                            puts distance r at the graph's centre, assuming vf=0.66"
    print "  --logf                    display graph with logarithmic frequency axis"
    print "  --ignorerxgain            don't include presumed rx gain in displayed power"
    raise Exception("")  # exit fast, but don't re-note the cause


def parseArgs():
    opts = ""
    args = ""
    try:
        opts, args = getopt.getopt(sys.argv[1:], "c:g:i:f:h", ["help", "tdrrange=", "gaincheck=", "logf", "ignorerxgain"])
    except getopt.GetoptError as err:
        helpArgs(str(err))  # will print something like "option -a not recognized"
    if len(args) < 1:
        helpArgs("missing filenameRoot")
    # set default/undefined values, where they make sense
    cropArg = ".5"  # a good compromise
    gainArg = "undefined"  # don't want to leave that to chance
    integrationArg = "10"  # the rtl_power default
    freqArg = "undefined"  # must be specified
    logf = False
    ignoreRxGain = False
    compositeArgument = ""  # might be e.g. --tdrrange 5
    for o, a in opts:
        if o == "-c":
            cropArg = a
        elif o == "-g":
            gainArg = a
        elif o == "-i":
            integrationArg = a
        elif o == "-f":
            freqArg = a
        elif o == "--tdrrange":
            compositeArgument = "--tdrrange " + a
            rM = float(a)
            fHighHz = 1500e6
            velocityFactor = 0.66
            binLengthM = 1.0 / fHighHz * (3e8 * velocityFactor) / 2
            numBins = 2 * (int(rM / binLengthM) + 1)  # round up, *2 to get in centre
            bin_sizeHz = int(fHighHz / numBins / 2.0)
            cropArg = "0.5"
            gainArg = "10"
            integrationArg = "20"
            freqArg = "30M:1500M:" + str(bin_sizeHz)
        elif o == "--gaincheck":
            compositeArgument = "--gaincheck " + a
            cropArg = "0.5"
            gainArg = a
            integrationArg = "1"
            freqArg = "24M:1500M:10M"
        elif o in ("-h", "--help"):
            helpArgs("")
        elif o in "--logf":
            logf = True
        elif o in "--ignorerxgain":
            ignoreRxGain = True
        else:
            helpArgs("option " + o + " not recognised")
    rtl_powerArgs = ["-f " + freqArg]
    rtl_powerArgs = rtl_powerArgs + ["-1"]  # exit after 1 sweep
    rtl_powerArgs = rtl_powerArgs + ["-c " + cropArg]
    rtl_powerArgs = rtl_powerArgs + ["-g " + gainArg]
    rtl_powerArgs = rtl_powerArgs + ["-i " + integrationArg]
    if cropArg == "undefined" or gainArg == "undefined" or integrationArg == "undefined" or freqArg == "undefined":
        msg = " ".join(["undefined parameter:"] + rtl_powerArgs)
        helpArgs(msg)
    # create filename, replace chars that are painful in unix filemanes, ensure not .csv
    filenameRoot = args[0] + compositeArgument + '_' + '_'.join(rtl_powerArgs)
    filenameRoot = filenameRoot.replace(' ', '_').replace(':', '_')
    return [filenameRoot, rtl_powerArgs, logf, ignoreRxGain]


def doScan():
    # noinspection PyBroadException
    try:
        [filenameRoot, rtlpowerArgs, logf, ignoreRxGain] = parseArgs()
        scanToFile(filenameRoot, rtlpowerArgs)
        drawGraphFromFiles([filenameRoot], "Power", "rtlscan.png", ignoreRxGain, logf, [])
    except:
        # catch everything
        if sys.exc_info()[1] != "":
            print "**** ", sys.exc_info()[1], " ****"
        sys.exit(1)


doScan()

 

rtlplot.py:

#! /usr/bin/env python
"""
MIT License

Copyright (c) 2016 Thomas Godfrey Gardner

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
"""

#v1.0

import sys
import getopt
from rtlutils import drawGraphFromFiles, toFilenameRoot


def helpArgs(msg):
    if msg != "":
        print "**** " + msg + " ****"
    print("rtlplot [-a|--annotation x,y,text]* [--logf] [--ignorerxgain] filename*")
    print("  draw spectrum of one or more filenames, with or without .csv")
    print("  --logf             a logarithmic frequency axis")
    print("  --ignorerxgain     don't include presumed rx gain in displayed power")
    print "  -a, --annotation   write 'text' at (x,y)"
    raise Exception("")  # exit fast, but don't re-note the cause


def parseArgs():
    if len(sys.argv) <= 1:
        helpArgs("missing filenames")
    opts, filenames = getopt.getopt(sys.argv[1:], "ha:", ["logf", "help", "ignorerxgain", "annotations="])
    logf = False
    ignoreRxGain = False
    annotations = []
    for o, a in opts:
        if o in "--logf":
            logf = True
        elif o in ("--help", "-h"):
            helpArgs("")
        elif o in "--ignorerxgain":
            ignoreRxGain = True
        elif o in ('--annotation', '-a'):
            # noinspection PyBroadException
            try:
                xs, ys, txt = a.split(',')
                x = float(xs)
                y = float(ys)
                annotations.append((x, y, txt))
            except:
                helpArgs("annotation '" + a + "' should be number,number,text")
        else:
            helpArgs("option " + o + " not recognised")
    filenamesRoot = map(toFilenameRoot, filenames)
    return [filenamesRoot, logf, ignoreRxGain, annotations]


def doPlot():
    # noinspection PyBroadException
    try:
        filenamesRoot, logf, ignoreRxGain, annotations = parseArgs()
        drawGraphFromFiles(filenamesRoot, 'Power', "rtlplot.png", ignoreRxGain, logf, annotations)
    except:
        # catch everything
        if sys.exc_info()[1] != "":
            print "**** ", sys.exc_info()[1], " ****"
        sys.exit(1)


doPlot()

 

rtldiff.py:

#! /usr/bin/env python
"""
MIT License

Copyright (c) 2016 Thomas Godfrey Gardner

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
"""

#v1.0

import sys
import getopt
from rtlutils import parseRtlOutputFile, calculateSpectrum, plotLine, initGraph, \
    displayGraph, standardSpectrumRubric, toFilenameRoot


def helpArgs(msg):
    if msg != "":
        print "**** " + msg + " ****"
    print("rtldiff [-a|annotation x,y,text]* [--logf] [--ignorerxgain] referenceFilename filename*")
    print("  draw one or more filename-referenceFilename spectrums, with or without .csv")
    print("  --logf: a logarithmic frequency axis")
    print("  --ignorerxgain: don't include presumed rx gain in displayed power")
    print "  -a, --annotation write 'text' at (x,y)"
    raise Exception("")  # exit fast, but don't re-note the cause


def parseArgs():
    if len(sys.argv) <= 2:
        helpArgs("missing filenames")
    annotations = []
    opts, filenames = getopt.getopt(sys.argv[1:], "ha:", ["logf", "help", "ignorerxgain", "annotations="])
    logf = False
    ignoreRxGain = False
    for o, a in opts:
        if o in "--logf":
            logf = True
        elif o in ("--help", "-h"):
            helpArgs("")
        elif o in "--ignorerxgain":
            ignoreRxGain = True
        elif o in ('--annotation', '-a'):
            # noinspection PyBroadException
            try:
                xs, ys, txt = a.split(',')
                x = float(xs)
                y = float(ys)
                annotations.append((x, y, txt))
            except:
                helpArgs("annotation '" + a + "' should be number,number,text")
        else:
            helpArgs("option " + o + " not recognised")
    filenamesRoot = map(toFilenameRoot, filenames)
    return [filenamesRoot, logf, ignoreRxGain, annotations]


def displayDiff():
    # noinspection PyBroadException,PyBroadException
    try:
        [filenamesRoot, logf, ignoreRxGain, annotations] = parseArgs()
        [refFreqs, refDbms, refFftBinSize, _] = calculateSpectrum(parseRtlOutputFile(filenamesRoot[0]), ignoreRxGain)

        initGraph("Relative Power", "MHz", "dB", annotations)
        fftBinSize = 0
        for filenameRoot in filenamesRoot[1:]:
            # print filenameRoot
            [freqs, dbms, fftBinSize, _] = calculateSpectrum(parseRtlOutputFile(filenameRoot), ignoreRxGain)
            # check files compatible with reference
            if freqs[0] != refFreqs[0]:
                raise Exception("Start frequency: {!s}/{!s} in {!s}".format(refFreqs[0], freqs[0], filenameRoot))
            if len(freqs) != len(refFreqs):
                raise Exception("Number of points: {!s}/{!s} in {!s}".format(len(refFreqs), len(freqs), filenameRoot))
            if fftBinSize != refFftBinSize:
                raise Exception("FFT Bin Size: {!s}/{!s} in {!s}".format(refFftBinSize, fftBinSize, filenameRoot))
            deltaDbms = []
            for [dbm, refDbm] in zip(dbms, refDbms):
                deltaDbms.append(dbm - refDbm)
            plotLine([refFreqs, deltaDbms], filenameRoot)
        displayGraph("Reference: " + filenamesRoot[0] + ", " + standardSpectrumRubric(fftBinSize, ignoreRxGain),
                     len(filenamesRoot) - 1, "rtldiff.png", logf)
    except:
        # catch everything
        if sys.exc_info()[1] != "":
            print "**** ", sys.exc_info()[1], " ****"
        sys.exit(1)


displayDiff()

 

rtltdr.py:

#! /usr/bin/env python
"""
MIT License

Copyright (c) 2016 Thomas Godfrey Gardner

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
"""

#v1.0

import sys
import getopt
import numpy as np
from rtlutils import parseRtlOutputFile, calculateSpectrum, plotLine, \
    initGraph, displayGraph, standardTimeDomainRubric, dbToPower, \
    lfExtendLine2, plotMarkers
from functools import partial


def helpArgs(msg):
    if msg != "":
        print "**** " + msg + " ****"
    print "rtltdr [-a|annotation x,y,text]* [--vf v] [-s sf] [--hann] [--rect]"
    print "       [--kaiser beta] [--hamming] [--bartlett] referenceFilename filename"
    print "  Draw TDR of one or more (filename-referenceFilename), with or without .csv"
    print "  Ticks indicate the bin boundaries"
    print "  --vf vf           velocity factor, typically 0.66 for cables/PCBs or "
    print "                     0.75 for TV cables (default 0.66)"
    print "  --sf sf            smoothing factor a.k.a. interpolation factor when displaying "
    print "                     TDR graphs, >=1 (default 4)"
    print "  --hann             Hanning window function (default). Best for discontinuity amplitude"
    print "  --rect             rectangular window function. Best for distance resolution"
    print "  --kaiser beta      Kaiser window function"
    print "  --hamming          Hamming window function"
    print "  --bartlett         Bartlett window function"
    print "  -a, --annotation   write 'text' at (x,y)"
    raise Exception("")  # exit fast, but don't re-note the cause


def toFilenameRoot(s):
    # strip off trailing .csv, if any
    if s.endswith(".csv"):
        return s[0:len(s) - 4]
    return s


def parseArgs():
    opts = ""
    args = ""
    try:
        opts, args = getopt.getopt(sys.argv[1:], "ha:", ["help", "vf=", "sf=", "hann", "rect", "kaiser=", "hamming", "bartlett", "annotation="])
    except getopt.GetoptError as err:
        helpArgs(str(err))  # will print something like "option -a not recognized"
    if len(sys.argv) <= 2:
        helpArgs("missing filenames")
    # set default/undefined values, where they make sense
    velocityFactor = 0.66
    smoothingFactor = 4
    annotations = []
    windowingFunction = (np.hanning, "Hann")
    for o, a in opts:
        if o in ("--help", "-h"):
            helpArgs("")
        elif o in "--vf":
            # noinspection PyBroadException
            try:
                velocityFactor = float(a)
                if velocityFactor <= 0 or velocityFactor > 1:
                    helpArgs("velocity factor must be >0 and <=1") except: helpArgs("velocity factor must be >0 and <=1")
        elif o in "--sf":
            # noinspection PyBroadException
            try:
                smoothingFactor = int(a)
                if a < 1: helpArgs("smoothing factor must be an integer >=1")
            except:
                helpArgs("smoothing factor must be an integer >=1")
        elif o in "--hann":
            windowingFunction = (np.hanning, "Hann")
        elif o in "--rect":
            windowingFunction = (partial(np.kaiser, beta=0), "Rectangular")
        elif o in "--hamming":
            windowingFunction = (partial(np.hamming), "Hamming")
        elif o in "--bartlett":
            windowingFunction = (partial(np.bartlett), "Bartlett")
        elif o in "--kaiser":
            try:
                beta = float(a)
                if beta < 1: helpArgs("beta >= 0")
                windowingFunction = (partial(np.kaiser, beta=beta), "Kaiser beta=" + str(beta))
            except:
                helpArgs("beta >= 0")
        elif o in ('--annotation', '-a'):
            # noinspection PyBroadException
            try:
                xs, ys, txt = a.split(',')
                x = float(xs)
                y = float(ys)
                annotations.append((x, y, txt))
            except:
                helpArgs("annotation '" + a + "' should be number,number,text")
        else:
            helpArgs("option " + o + " not recognised")
    filenamesRoot = map(toFilenameRoot, args)
    return [filenamesRoot, velocityFactor, smoothingFactor, windowingFunction, annotations]


def createBinMarkers(numBins, binLengthM, irPeak):
    bx = [0.0, 0.0, 0.0, 0.0] * numBins
    by = [-0.0001, 0.01 * irPeak, 0.01 * irPeak, -0.0001] * numBins
    for i in range(0, numBins):
        bx[4 * i + 0] = i * binLengthM - 0.0011 * binLengthM
        bx[4 * i + 1] = i * binLengthM - 0.0010 * binLengthM
        bx[4 * i + 2] = i * binLengthM + 0.0010 * binLengthM
        bx[4 * i + 3] = i * binLengthM + 0.0011 * binLengthM
    bx[0] = 0.0  # not < 0
    bx[1] = 0.0  # not < 0
    return bx, by


# read a spectrum and check it is compatible with the reference spectrum
def readSpectrumAndCheckCompatibility(filenameRoot, ignoreRxGain, refFftBinSize, refFreqs):
    [freqs, powersDB, fftBinSize, _] = calculateSpectrum(parseRtlOutputFile(filenameRoot), ignoreRxGain)
    # check files compatible with reference
    if refFreqs[0] != freqs[0]:
        raise Exception("Incompatible start frequency: {!s} vs {!s} in {!s}".format(refFreqs[0], freqs[0], filenameRoot))
    if len(freqs) != len(refFreqs):
        raise Exception(
            "Incompatible number of points: {!s} vs {!s} in {!s}".format(len(refFreqs), len(freqs), filenameRoot))
    if fftBinSize != refFftBinSize:
        raise Exception("Incompatible FFT Bin Size: {!s} vs {!s} in {!s}".format(refFftBinSize, fftBinSize, filenameRoot))
    return powersDB


def calculateImpulseResponse(powersDB, refPowersDB, refFreqs, fftBinSize, interpolationFactor, windowingFunction):
    # calculate spectrum relative to the reference spectrum
    deltaPowersDB = []
    for [powerDB, refPowerDB] in zip(powersDB, refPowersDB):
        deltaPowersDB.append(powerDB - refPowerDB)

    # pad down to zero frequency, and convert from logarithmic scale to linear scale
    zeroPaddedPowersDB = lfExtendLine2(refFreqs, deltaPowersDB, fftBinSize)
    powers = map(dbToPower, zeroPaddedPowersDB)

    # apply window
    windowFunc = windowingFunction[0]
    window = windowFunc(len(powers))
    windowedPowers = np.multiply(window, powers)

    # interpolate in the time domain by padding with 0 in the frequency domains
    paddedPowers = [0.0] * interpolationFactor * len(windowedPowers)
    for i in range(0, len(windowedPowers) - 1):
        paddedPowers[i] = windowedPowers[i]

    # convert to time domain
    # remove values irrelevant because of the lack of phase information in the
    # spectrum and due to any reflections in the TDR
    timeDomain = np.fft.ifft(paddedPowers)
    n = 2
    for x in range(0, n*interpolationFactor+1, 1):  # zero out the first n bins to get rid of crap
        timeDomain[x] = 0.0
    timeDomain = np.absolute(timeDomain)
    impulseResponse = timeDomain[:len(timeDomain) / 2]
    return impulseResponse


def displayTdr():
    # noinspection PyBroadException
    try:
        [filenamesRoot, velocityFactor, smoothingFactor, windowingFunction, annotations] = parseArgs()
        ignoreRxGain = False
        logf = False

        # read the reference spectrum
        [refFreqs, refPowersDB, refFftBinSize, _] = calculateSpectrum(parseRtlOutputFile(filenamesRoot[0]), ignoreRxGain)

        # calculate x axis' parameters
        fHighHz = (refFreqs[-1] + (refFreqs[1] - refFreqs[0]))  # not n-1
        numBins = int(fHighHz / (refFreqs[1] - refFreqs[0])) / 2  # 2 since only half of ifft is used
        binLengthM = 1.0 / (fHighHz * 1000000) * (3e8 * velocityFactor) / 2

        initGraph("Impulse Response", "Distance (m), v=" + str(velocityFactor) + "c", "Amplitude", annotations)

        irPeak = 0.0  # for scaling the markers
        for filenameRoot in filenamesRoot[1:]:
            # read next spectrum from file
            powersDB = readSpectrumAndCheckCompatibility(filenameRoot, ignoreRxGain, refFftBinSize, refFreqs)
            # convert to time domain
            impulseResponse = calculateImpulseResponse(powersDB, refPowersDB, refFreqs, refFftBinSize, smoothingFactor, windowingFunction)
            irPeak = max(irPeak, max(impulseResponse))
            # calculate the xaxis
            distance = [d * binLengthM / smoothingFactor for d in range(0, len(impulseResponse))]
            # plot the line
            plotLine([distance, impulseResponse], filenameRoot)

        # plot markers showing the bin length
        bx, by = createBinMarkers(numBins, binLengthM, irPeak)
        plotMarkers([bx, by])

        graphHeader = "Reference: " + filenamesRoot[0] + ", Window: " + windowingFunction[1] + ", " + standardTimeDomainRubric(binLengthM)
        displayGraph(graphHeader, len(filenamesRoot) - 1, "rtltdr.png", logf)
    except:
        # catch everything
        if sys.exc_info()[1] != "":
            print "**** ", sys.exc_info()[1], " ****"
        sys.exit(1)

displayTdr()

 

rtlutils.py:

"""
MIT License

Copyright (c) 2016 Thomas Godfrey Gardner

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
"""

#v1.0

import sys
import subprocess
from collections import defaultdict
from pylab import frange, plt, draw, gca, connect


class DeltaCursor:
    """
    Enable measurements of delta x and delta y
    When button 2 is pressed and dragged the original position is marked
    and the delta position is displayed
    """

    def __init__(self, ax):
        self.ax = ax
        self.lx = ax.axhline(color='k', alpha=0.5, linestyle="dashed")  # the horiz line
        self.ly = ax.axvline(color='k', alpha=0.5, linestyle="dashed")  # the vert line

        # text location in axes coords
        self.txt = ax.text(0.01, 0.95, '', transform=ax.transAxes)
        self.xdown = None
        self.ydown = None
        self.lx.set_ydata(-1000000000)
        self.ly.set_xdata(-1000000000)

    def mouse_move(self, event):
        if not event.inaxes:
            return
        x, y = event.xdata, event.ydata

        if self.xdown is not None:
            xdelta = x - self.xdown
            ydelta = y - self.ydown
            self.txt.set_text('[%2.3f, %2.3f]' % (xdelta, ydelta))
        draw()

    def button_press(self, event):
        if event.button != 3:
            return
        x, y = event.xdata, event.ydata
        self.xdown = x
        self.ydown = y
        # update the line positions
        self.lx.set_ydata(y)
        self.ly.set_xdata(x)
        draw()

    def button_release(self, event):
        if event.button != 3:
            return
        self.txt.set_text("")
        self.xdown = None
        self.ydown = None
        # move lines out of the way
        self.lx.set_ydata(-1000000000)
        self.ly.set_xdata(-1000000000)
        draw()


def toFilenameRoot(s):
    # strip off trailing .csv, if any
    if s.endswith(".csv"):
        return s[0:len(s) - 4]
    return s


def fftBinSizeToHz(fftBinSize):
    return int((fftBinSize.split()[3])[:-2])  # assume string showing 00000Hz


def dbToPower(db):
    return 10 ** (db / 10)


def binLength(binHz, numBins, vFactor):
    return 1.0 / binHz / numBins * (3e8 * vFactor) / 4


def scanToFile(filenameRoot, rtlpowerArgs):
    filenameCsv = filenameRoot + ".csv"
    print("rtl_power output directed to file: " + filenameCsv)
    print(["rtl_power"] + rtlpowerArgs)

    subP = subprocess.Popen(["rtl_power"] + rtlpowerArgs, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
    stdout, stderr = subP.communicate()

    isDisplaySamples = True
    with open(filenameCsv, "w") as rtlOutputFile:
        # record the command line and the derived conditions
        rtlOutputFile.write("# rtl_power {} > {}\n".format(rtlpowerArgs, filenameCsv))
        for ln in stderr.splitlines():
            rtlOutputFile.write("# {!s}\n".format(ln))
            sys.stdout.write("rtl_power> {!s}\n".format(ln))
        # record the data
        for ln in stdout.splitlines():
            if isDisplaySamples:
                isDisplaySamples = False
                sys.stdout.write("rtl_power> samples/bin: {!s}\n".format(ln.strip().split(', ')[5]))
            rtlOutputFile.write("{!s}\n".format(ln))
    rtlOutputFile.close()


def parseRtlOutputFile(filenameRoot):
    filename = filenameRoot + ".csv"
    sums = defaultdict(float)
    counts = defaultdict(int)
    fftBinSize = ""
    rxGain = 0.0
    with open(filename, "r") as rtlOutputFile:
        for line in rtlOutputFile:
            if line.startswith("#"):
                # ignore comments except for any that we want to show on the graph
                comment = line[1:].strip()
                if comment.startswith("FFT bin size"):
                    fftBinSize = comment
                gainStr = "Tuner gain set to "
                if comment.startswith(gainStr):
                    comment = comment[len(gainStr):-3]
                    rxGain = float(comment)
            else:
                # process the data
                line = line.strip().split(', ')
                fLow = int(line[2])
                fHigh = int(line[3])
                fStep = float(line[4])
                weight = int(line[5])
                dbm = [float(d) for d in line[6:]]
                for freq, d in zip(frange(fLow, fHigh, fStep), dbm):
                    sums[freq] += d * weight
                    counts[freq] += weight
        rtlOutputFile.close()
    return [sums, counts, fftBinSize, rxGain]


def calculateSpectrum(resultsFromFile, ignoreRxGain):
    [sums, counts, fftBinSize, rxGain] = resultsFromFile
    if ignoreRxGain:
        rxGain = 0.0
    ave = defaultdict(float)
    for freq in sums:
        ave[freq] = sums[freq] / counts[freq]
    freqs = []  # for the plot
    powersDB = []  # for the plot
    for freq in sorted(ave):
        freqs.append(freq / 1e6)  # in MHz
        powersDB.append(ave[freq] - rxGain)
    return [freqs, powersDB, fftBinSize, rxGain]


# In order to do the ifft, must have a power for each frequency bin down to f=zero.
# Somewhat arbitrary, but pad a line down to zero frequency, with the lowest value.
# This would be a good place to have a standard windowing function
def lfExtendLine(freqs, dbms, fftBinSize):
    binSize = int((fftBinSize.split()[3])[:-2])  # assume string showing 00000Hz
    extendFreqs = [binSize * freq / 1e6 for freq in range(int(1e6 * freqs[0] / binSize))]
    extendDbms = []
    for freq in range(int(1e6 * freqs[0] / binSize)):
        extendDbms.append(dbms[0])
    return [extendFreqs + freqs, extendDbms + dbms]


# In order to do the ifft, must have a power for each frequency bin down to f=zero.
# The line is padded down to zero frequency with the lowest available value.
# That is somewhat arbitrary, but in most cases will be unimportant since the
# low-frequency powers will be significantly reduced by any sane windowing function
def lfExtendLine2(freqs, dbms, fftBinSize):
    binSize = int((fftBinSize.split()[3])[:-2])  # assume string showing 00000Hz
    extendDbms = []
    for freq in range(int(1e6 * freqs[0] / binSize)):
        extendDbms.append(dbms[0])
    return extendDbms + dbms


def plotLine(points, lineLabel):
    [freqs, dbms] = points
    plt.plot(freqs, dbms, label=lineLabel)


def plotMarkers(points):
    [x, y] = points
    plt.plot(x, y, linewidth=0.2, color='#000000', alpha=0.5, linestyle="solid")


def plotLineFromFile(filenameRoot, ignoreRxGain):
    [freqs, dbms, fftBinSize, rxGain] = calculateSpectrum(parseRtlOutputFile(filenameRoot), ignoreRxGain)
    plotLine([freqs, dbms], filenameRoot)
    return [fftBinSize, rxGain]


def setHeader(header):
    plt.title(header, fontsize="x-small")


def standardSpectrumRubric(fftBinSize, ignoreRxGain):
    if ignoreRxGain:
        header = "RxGain Ignored, " + fftBinSize
    else:
        header = "RxGain Included, " + fftBinSize
    return header


def standardTimeDomainRubric(lengthM):
    return "Resolution: %dmm" % int(1000 * lengthM)


def drawGraphFromFiles(filenamesRoot, titleString, pngFilename, ignoreRxGain, logf, annotations):
    initGraph(titleString, "MHz", "dB", annotations)
    fftBinSize = 0
    for filenameRoot in filenamesRoot:
        [fftBinSize, _] = plotLineFromFile(filenameRoot, ignoreRxGain)
    displayGraph(standardSpectrumRubric(fftBinSize, ignoreRxGain), len(filenamesRoot), pngFilename, logf)


def initGraph(titleString, xLabel, yLabel, annotations):
    """
    Everything that can be setup before the lines are known
    :param titleString:
    :param xLabel:
    :param yLabel:
    :param annotations:
    :return:
    """
    plt.figure(figsize=(10, 6), dpi=96)
    plt.suptitle(titleString, fontsize="large")
    plt.xlabel(xLabel)
    plt.ylabel(yLabel)
    plt.grid(which='major', axis='y', linestyle=':', color='gray')
    plt.grid(which='major', axis='x', linestyle=':', color='gray')
    for (x, y, text) in annotations:
        plt.gca().annotate(text, xy=(x, y), xycoords='data')
        plt.plot(x, y, 'x', color='black')
    global cursor
    cursor = DeltaCursor(gca())
    connect('motion_notify_event', cursor.mouse_move)
    connect('button_press_event', cursor.button_press)
    connect('button_release_event', cursor.button_release)


def displayGraph(header, numLegendCols, filenamePng, logf):
    """
    Now that everything is known, finish off drawing the graph
    :param header:
    :param numLegendCols:
    :param filenamePng:
    :param logf:
    :return:
    """
    setHeader(header)
    if logf:
        plt.xscale('log', subsx=[2, 3, 4, 5, 6, 7, 8, 9])
    legendString = plt.legend(loc="lower center", ncol=min(numLegendCols, 2), prop={'size': 'x-small'})
    legendString.get_frame().set_alpha(0.5)
    with open(filenamePng, "w") as pngFile:
        plt.savefig(filenamePng, bbox_inches='tight')
    pngFile.close()
    plt.show()
Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s