Thursday, August 6, 2020

How to receive and decode multiple weather sondes with only one RTL-SDR receiver

Note: there's part 2 to this article with the script updated for performance

It could sound strange but I haven't seen a simple solution that would allow me to receive and decode multiple weather sondes at the same time with only one receiver. dxlAPRS software has a channelizer if I'm not mistaken but the software itself is too complex and is designed for things I'm not even interested in. auto-rx on the other hand would require a receiver per frequency. But is it a problem to get me more RTL-SDR dongles to receive as many sondes as I like at the same time? No, getting an RTL-SDR dongle is not a problem, the problem is the antenna - each dongle would require one. So either I would have to install as many antennas as dongles or I would have to use an antenna splitter. And without an amplifier (LNA) a splitter would reduce signal power by a number of its outputs. So both solutions would require additional hardware and investments. But then I had an idea that it should be possible to build a very simple (and very inefficient of course) channelizer even in bash by multiplexing IQ samples baseband stream into multiple decoding processes.

I had this idea for a long time actually but only now I found some time to implement it. The idea is pretty simple - the only thing I need to do is to multiplex the IQ stream from the dongle into multiple decoding processes (slots). There should be a controlling and a frequency scanner processes additionally. As far as I'm going to use the pipe'n'filter pattern I have to preallocate everything statically as a long pipe. This solution allowed me to implement everything even in bash:

#!/bin/bash

. ./defaults.conf

SCAN_BINS=4096
SCAN_AVERAGE_TIMES=100
SCAN_UPDATE_RATE=1
SCAN_UPDATE_RATE_DIV=5 # 5 seconds
SCAN_POWER_THRESHOLD=-69

SCANNER_PORT=5677
DECODER_PORT=5678

MAX_SLOTS=10 # maximum number of rx slots
SLOT_TIMEOUT=10 # i.e 30 seconds *10 = 5 minutes
SLOT_ACTIVATE_TIME=4 # 4 * 5 seconds = 20 seconds min activity

OPTIND=1 #reset index
while getopts "ha:p:f:s:g:p:P:t:" opt; do
  case $opt in
     h)  show_usage $(basename $0); exit 0; ;;
     a)  address="$OPTARG" ;;  # not used atm
     p)  port="$OPTARG" ;;     # not used atm
     f)  TUNER_FREQ="$OPTARG" ;;
     s)  TUNER_SAMPLE_RATE="$OPTARG" ;;
     g)  TUNER_GAIN="$OPTARG" ;;
     P)  DONGLE_PPM="$OPTARG" ;;
     t)  SCAN_POWER_THRESHOLD="$OPTARG" ;;
     \?) exit 1 ;;
     :)  echo "Option -$OPTARG requires an argument" >&2;exit 1 ;;
  esac
done
shift "$((OPTIND-1))"
 
[ ! "$TUNER_FREQ" -eq 0 ] || show_error_exit "Wrong frequency"
[ ! "$TUNER_SAMPLE_RATE" -eq 0 ] || show_error_exit "Wrong sample rate"

DECIMATE=$((TUNER_SAMPLE_RATE/DEMODULATOR_OUTPUT_FREQ))
[ "$((DECIMATE*DEMODULATOR_OUTPUT_FREQ))" -ne "$TUNER_SAMPLE_RATE" ] && show_error_exit "Sample rate should be multiple of $DEMODULATOR_OUTPUT_FREQ"


cleanup()
{
  local children child
  children="$1 $2 $(get_children_pids $1) $(get_children_pids $2)"
  kill $children &>/dev/null;wait $children &>/dev/null
}

scan_power()
{
  ./csdr convert_u8_f | \
  ./csdr fft_cc $SCAN_BINS $((TUNER_SAMPLE_RATE/(SCAN_UPDATE_RATE*SCAN_AVERAGE_TIMES/SCAN_UPDATE_RATE_DIV))) | \
  ./csdr logaveragepower_cf -70 $SCAN_BINS $SCAN_AVERAGE_TIMES | \
  ./csdr fft_exchange_sides_ff $SCAN_BINS | \
  ./csdr dump_f | tr ' ' '\n' | \
   awk -v f=$TUNER_FREQ -v sr=$TUNER_SAMPLE_RATE -v bins=$SCAN_BINS '
     BEGIN{fstep=sr/bins;fstart=f-sr/2;print fstep;print fstart}
     {printf("%d %.1f\n",fstart+fstep*((NR-1)%bins),$0);
      if(0==(NR%bins)){printf("\n")};fflush()}' | \
   awk -v step=$((TUNER_SAMPLE_RATE/SCAN_BINS)) '
     function abs(x){return (x<0)?-x:x}
     {if(length($1)!=0){if(abs($1-10000*int($1/10000)<step)){print $0}}
      else{print};
      fflush();}' | \
   awk -v thr=$SCAN_POWER_THRESHOLD '{if (length($2)!=0){if(int($2)>thr){print 10000*int(int($1)/10000)" "$2;fflush()}}}'
}

# the line below should come before the m10mod if needed.
#      tee >(c50dft -d1 --ptu --json /dev/stdin > /dev/stderr) | \
decode_sonde()
{
  local bpf3=$(calc_bandpass_param 3000 $DEMODULATOR_OUTPUT_FREQ)
  local bpf10=$(calc_bandpass_param 10000 $DEMODULATOR_OUTPUT_FREQ)
  (
    ./csdr convert_u8_f | \
    ./csdr shift_addition_switchable_cc --fifo "$1" | \
    ./csdr fir_decimate_cc $DECIMATE 0.005 HAMMING | \
    tee >(
      ./csdr bandpass_fir_fft_cc -$bpf10 $bpf10 0.02 | \
      ./csdr fmdemod_quadri_cf | ./csdr limit_ff | ./csdr convert_f_s16 | \
      sox -t raw -esigned-integer -b 16 -r $DEMODULATOR_OUTPUT_FREQ - -b 8 -c 1 -t wav - highpass 10 gain +5 | \
      ./m10mod --ptu --json > /dev/stderr
    ) | \
    ./csdr bandpass_fir_fft_cc -$bpf3 $bpf3 0.02 | \
    ./csdr fmdemod_quadri_cf | ./csdr limit_ff | ./csdr convert_f_s16 | \
    sox -t raw -esigned-integer -b 16 -r $DEMODULATOR_OUTPUT_FREQ - -b 8 -c 1 -t wav - highpass 10 gain +5 | \
    tee >(./dfm09mod --ptu --ecc --json -vv /dev/stdin > /dev/stderr) \
        >(./dfm09mod --ptu --ecc --json -i /dev/stdin > /dev/stderr) \
        >(./rs41mod --ptu --ecc --crc --json -vv /dev/stdin > /dev/stderr) \
        >(./rs92mod -e "$EPHEM_FILE" --crc --ecc --json /dev/stdin > /dev/stderr) | \
    aplay -r 48000 -f S8 -t wav -c 1 -B 500000 &> /dev/null
  ) &>/dev/stdout | while read LINE; do
      tmp_freq="$(echo "$LINE" | awk -F: '/^shift_addition_frequency/{print $2}')"
      [ -n "$tmp_freq" ] && {
        # due to float (un)precision, calculated freq could be +-1, so it is corrected
        freq=$(awk -v f="$tmp_freq" -v sr="$TUNER_SAMPLE_RATE" -v tf="$TUNER_FREQ" 'BEGIN{printf 1000*int((1+tf-int(f*sr))/1000)}')
      }
      echo "$LINE" | grep --line-buffered -E '^{' | jq --unbuffered -rcM '. + {"freq":"'"$freq"'"}' | \
      (flock 200; socat -u - UDP4-DATAGRAM:127.255.255.255:$DECODER_PORT,broadcast,reuseaddr) 200>$MUTEX_LOCK_FILE
    done
}

declare -A actfreq # active frequencies

# create fifos
FIFO_BASE_NAME='/tmp/sonde_scanner_'
for ((fifo=1;fifo<=MAX_SLOTS;fifo++)); do
  fifo_name="${FIFO_BASE_NAME}${fifo}.pipe"
  [ -e "$fifo_name" ] && \rm "$fifo_name"    
  mkfifo "$fifo_name"
  fifos[$fifo]="$fifo_name"
done
 
(socat -u UDP-RECVFROM:$SCANNER_COM_PORT,fork,reuseaddr - | while read LINE; do
  case "$LINE" in
    TIMER30)
       for freq in "${!actfreq[@]}"; do 
         actfreq[$freq]=$((actfreq[$freq]-1))
         [ "${actfreq[$freq]}" -gt "$SLOT_TIMEOUT" ] && actfreq[$freq]=$SLOT_TIMEOUT
         if [ "${actfreq[$freq]}" -eq 0 ]; then
           # deactivate slot
           for slot in "${!slots[@]}"; do
             [ "$freq" = "${slots[$slot]}" ] && {
               echo "1.0" > "${fifos[$slot]}"
               unset slots[$slot]
               break
             }
           done
           unset actfreq[$freq]
         elif [ "${actfreq[$freq]}" -ge $SLOT_ACTIVATE_TIME ]; then
           [[ "${slots[@]}" =~ "$freq" ]] || {
             # activate slot
             for ((slot=1;slot<=MAX_SLOTS;slot++)); do
               [ -z "${slots[$slot]}" ] && {
                 slots[$slot]="$freq"
                 calc_bandpass_param "$((TUNER_FREQ-freq))" "$TUNER_SAMPLE_RATE" > "${fifos[$slot]}"
                 break
               }
             done
           }
         fi
       done
       ;;
    *) freq="${LINE% *}"
       [ -z "${_FREQ_BLACK_LIST[$freq]}" ] && {
         [ -n "$freq" ] && actfreq[$freq]=$((actfreq[$freq]+1))           
       }
       ;;
  esac
done) &
pid1=$!

(while sleep 30; do echo "TIMER30" | socat -u - UDP4-DATAGRAM:127.255.255.255:$SCANNER_PORT,broadcast,reuseaddr; done) &
pid2=$!

trap "cleanup $pid1 $pid2" EXIT INT TERM

# This is a magic line that kick-starts the decode_sonde processes
# I really don't understand why it is needed because the way I
# implemented the shift_addition_switchable_cc should make it work
# without this line but it doesn't. Some kind of weired magic
(sleep 2;eval "$(printf 'echo '1.0' > %q;' "${fifos[@]}")") &
 
rtl_sdr -p $DONGLE_PPM -f $TUNER_FREQ -g $TUNER_GAIN -s $TUNER_SAMPLE_RATE - | \
eval "tee $(printf '>(decode_sonde %q) ' "${fifos[@]}")" | \
scan_power | socat -u - UDP4-DATAGRAM:127.255.255.255:$SCANNER_PORT,broadcast,reuseaddr

There's some magic hidden behind the syntax so I'll try to explain how it works.

Here is the basic structure of the script:


IQ samples are multiplexed into a number of decode_sonde processes and into the scan_power process. scan_power process outputs integrated results of the FFT for the whole band once in 5 seconds. The bandwidth is defined by the TUNER_SAMPLE_RATE parameter. scan_power outputs only frequencies divisible by 10000 and only when the power level is greater than SCAN_POWER_THRESHOLD. scan_power communicates with the slot_activate/deactivate_logic over UDP on the localhost. There's an additional process though that simply notifies the Slot_activate/deactivate_logic process each 30 seconds (I forgot to put it into the diagram above). If a frequency is active at least SLOT_ACTIVATE_TIME+1 times (20 seconds) within 30 seconds period then a slot is activated for that frequency and one of the decode_sonde processes is assigned to the slot and all decoders of that particular decode_sonde process will also become active. If that frequency is not active then the slot will be deactivated within 30*(SLOT_ACTIVATE_TIME+1) minimum to 30*SLOT_TIMEOUT maximum seconds. Simply put it's just a hysteresis to avoid activating slots too early and deactivate them too soon. Maximum number of slots are limited by MAX_SLOTS constant - this is the theoretical maximum number of sondes that can be received simultaneously. Theoretical because some slots could be occupied by noise signals and because the number of maximum decode_sonde processes is limited rather by the CPU and not by that constant.

slot_activate/deactivate_logic process communicates with decoder processes over named pipes. When it wants a slot to be activated it writes the shift rate into that pipe so that the decoder process knows how should it shift the baseband input signal in the frequency domain.

To make this concept work I had to patch the csdr framework by adding a new switchable shift_addition_cc. I simply added the following into the csdr.c
 
    if(!strcmp(argv[1],"shift_addition_switchable_cc"))
    {
        bigbufs=1;

        float starting_phase=0;
        float rate = 1.0f;

        int fd = init_fifo(argc,argv);
        if (!fd)
        {
            if(argc<=2) return badsyntax("need required parameter (rate)");
            sscanf(argv[2],"%g",&rate);
        }

        if(!sendbufsize(initialize_buffers())) return -2;
 
        for(;;)
        {
            shift_addition_data_t data=shift_addition_init(rate);
            errhead();
            if (0==(int)abs(rate)) {
              fprintf(stderr,"\nshift_addition_frequency:%.8f\n",rate);
            }
            int remain, current_size;
            float* ibufptr;
            float* obufptr;
            for(;;)
            {
                FEOF_CHECK;
                if(!FREAD_C) break;
                if (0==(int)abs(rate))
                {
                  remain=the_bufsize;
                  ibufptr=input_buffer;
                  obufptr=output_buffer;
                  while(remain)
                  {
                      current_size=(remain>1024)?1024:remain;
                      starting_phase=shift_addition_cc((complexf*)ibufptr, (complexf*)obufptr, current_size, data, starting_phase);
                      ibufptr+=current_size*2;
                      obufptr+=current_size*2;
                      remain-=current_size;
                  }
                  FWRITE_C;
                }
                if(read_fifo_ctl(fd,"%g\n",&rate)) break;
                TRY_YIELD;
            }
        }
        return 0;
    }

It work the same way as the shift_addition_cc block except when the rate is set above or equal 1.0 it would only consume the input without sending anything to its output. So if I want a slot to be activated I'd write the shift ratio into the fifo file the block is reading. And if I want the slot to be deactivating then I'd write 1.0 into the fifo effectively blocking any process behind the shift_addition_switchable_cc block. Without this all decoders would be activated just after the start and this would overload my CPU. Well on my machine I can receive and decode 5 to 6 sondes at the same time so if I would set MAX_SLOTS to 6 I wouldn't need to patch the csdr but my CPU would be 100% busy constantly. And with this block its only 30% without activated slots.

There's an additional trick I'm using in the shift_addition_switchable_cc block. Normally the decode_sonde process doesn't know to which frequency it is tuned to. To provide this information in the decoder JSON output the shift_addition_switchable_cc writes the shift ratio value into the process's stderr stream and decode_sonde process would convert it to the frequency because it knows the sampling rate and the tuner frequency.

It makes it clear from the code that I start all the decoders at the same time for a single frequency. This isnt' very efficient of course but doesn't require any detection logic which simplifies the things a lot. Also I have separated sondes into two groups - sondes like M10 that uses 20 kHz bandwidth and the sondes that only require 6 kHz. This improves decoders sensitivity.

The decoder output is then broadcasted on another UDP port locally and can easily be changed to be broadcasted on my local network. Additionally instead of the rtl_sdr I could connect to a rtl_tcp daemon using my client to get the IQ samples from a remote computer.

There's a tricky line in the script that starts all the decode_sonde processes:
eval "tee $(printf '>(decode_sonde %q) ' "${fifos[@]}")" | \

This line expands into this:
tee >(decode_sonde "${fifos[1]}") \
    >(decode_sonde "${fifos[2]}") \
    >(decode_sonde "${fifos[3]}") \
    >(decode_sonde "${fifos[4]}") \
    >(decode_sonde "${fifos[5]}") \
    >(decode_sonde "${fifos[6]}") \
    >(decode_sonde "${fifos[7]}") \
    >(decode_sonde "${fifos[8]}") \
    >(decode_sonde "${fifos[9]}") \
    >(decode_sonde "${fifos[10]}")

How to use the script

Note: this script is proven to work on Intel n5000 processor and Ubuntu 20.04. On my HummingBoard Pro which is equipped with the NXP i.MX6 Quad core Arm Cortex A9 CPU the script would not work even if I set MAX_SLOTS to 1. It is simply not powerful enough.


I usually start the script with the following line:
./receivemultisonde.sh -f 403003000 -s 2400000 -P 35 -g 40

I tune to a frequency band where the most sondes are around my place but I don't tune to a multiple of 10000 Hz

Sampling rate is set to 2400000. It's tempting to set it to a higher value to get more baseband but it won't work. Sampling at higher frequencies would result in samples loss and this will make decoders unhappy. Simply don't go above 2400000. There's another HW bug which makes it impossible to set the sampling rate below 225001 and between 300000 and 900000 Hz.

Then the scanner threshold depends on the gain value. If the gain is changed then the scanner threshold needs to be adjusted. It's always a compromise between decoder sensitivity and the noise tolerance. If the threshold is too low then a lot of slots will be consumed by parasitic signals. The higher the threshold the less sensitive the scanner will be. Usually I'd choose a new gain value, set a lower threshold and then I'd listen the scanner on its port to check reported power levels. Then I'd adjust the threshold until nothing is reported assuming that no sondes are visible. Frequencies with parasitic signals can be blacklisted in the defaults.conf .

Note: currently the script assumes that all the decoders and the csdr binary are in the same folder. Maybe later I'll make it configurable. 

Note2: csdr is compiled as a binary and a lib. As I never do 'make install' I have a little script to help csdr binary find its lib (check csdr in the scripts)


Problem with the RS92 decoder


There's a little problem with the RS92 sondes. Unlike other sonde types RS92 decoder requires either an almanac or ephemeris information files be downloaded from the internet. Without this information the GPS coordinates cannot' be correctly decoded. As far as I canno't restart the rs92mod decoder process without restarting the script I can't actually decode these kind of sondes. These sondes are extremely rare now so it can be just skipped completely but I still hope to catch one so I patched the rs92mod code to make it reread the ephemeris or almanac file upon receiving the SIGUSR1 signal.

So I start another cron job that get the file from the internet regularly using getephemeris.sh or getsemalmanac.sh scripts and then it simply sends SIGUSR1 signal to all rs92mod processes:

pgrep rs92mod | xargs -I{} kill -SIGUSR1 {}

Note: check if gnss-sdr project or ublox_ephemeris would allow to generate RINEX files to make system work completely offline


How to display multiple sondes


The next problem I faced as everything started to work is how to display multiple sondes. In my previous article about radiosondes hunting I explained how to use Foxtrotgps program to show position of a single sonde. But this program cannot be used to display multiple targets at the same time. So I patched the auto-rx to simply listen on a UDP port without it interacting with the rtl-sdr hardware. Here for example I'm receiving 3 sondes simultaneously and showing them using auto-rx:




What else could be done


Currently I only scan signal with 10000 Hz steps because almost all sondes are started on frequencies divisible by 10 kHz. But some sondes aren't. My current solution is to change bandpass filter from 3kHz to 5kHz but later I could try to reduce the step to make tuning more precise without triggering a lot of false positives.

Another thing is that I haven't tested the script with M10 type of sondes that use wider frequency band. I don't know how the scanner will work with such wide band signals. There simply no such sondes flying around at my place to test.

Additionally it might be necessary to define left and right limits to the whole bandwidth defined by the sampling rate. This is needed to avoid aliased signal appearing at the edges and also the signal level drops at the edges

Now, to cover the whole frequency range which is 6 MHz I would need 3 dongles (well in my place even two) and using say two computers and two dongles I could easily combine two scripts output and then feed the auto-rx with it.


Note: all the scripts can be found here 

EDIT 8.Aug.2020
I've changed the script in the repository so that it also outputs power scanner results in JSON format per UDP on SCANNER_OUT_PORT

Additionally I've added a script plotpowerjson.sh to plot the signal power:

There're two sondes visible here and some parasitic signals. Only one parasitic signal is close to the 10 kHz boundary and it can be seen on the power map. Then this particular signal is blacklisted in my defaults.conf

gnuplotblock.sh script is described here

EDIT 9.Aug.2020
The script described in the article is renamed to legacy.receivemultisonde.sh because a new concept emerged since yesterday. The new concept uses new iq_server/iq_client code from Zilog80 and allows dynamic allocation and disposal of decode_sonde processes. I'll describe it shortly in a new article

1 comment:

Anonymous said...

Hi,

würde das ganze auch auf einem Raspberry4 laufen?