Baseband FM PoC

BBFM.md

@@ -0,0 +1,83 @@
+# Radio Autoencoder - Baseband FM (BBFM)
+
+A version of the Radio Autoencoder (RADE) designed for the baseband FM channel provided by DC coupled and passband FM radios, e.g. land mobile radio (LMR) VHF/UHF use case.
+
+# BBFM ML encoder/decoder
+
+1. First pass training command line:
+    ```
+    python3 ./train_bbfm.py --cuda-visible-devices 0 --sequence-length 400 --batch-size 512 --epochs 100 --lr 0.003 --lr-decay-factor 0.0001 --plot_loss ~/Downloads/tts_speech_16k_speexdsp.f32 model_bbfm_01 --range_EbNo --range_EbNo_start 6 --plot_loss
+    ```
+
+1. Inference (runs encoder and decoder, and outputs symbols `z_hat.f32`):
+    ```
+    ./inference_bbfm.sh model_bbfm_01/checkpoints/checkpoint_epoch_100.pth wav/brian_g8sez.wav - --write_latent z_hat.f32
+    ```
+1. Stand alone decoder, outputs speech from `z_hat.f32` to sound card:
+    ```
+    ./rx_bbfm.sh model_bbfm_01/checkpoints/checkpoint_epoch_100.pth z_hat.f32 -
+    ```
+1. Or save speech out to a wave file:
+    ```
+    ./rx_bbfm.sh model_bbfm_01/checkpoints/checkpoint_epoch_100.pth z_hat.f32 t.wav
+    ```
+
+1. Plot sequence of received symbols:
+    ```
+    octave:4> radae_plots; do_plots_bbfm('z_hat.f32')
+    ```
+
+# Fading channel simulation
+
+HF channel sim (two path Rayleigh) is pretty close to TIA-102.CAAA-E 1.6.33 Faded Channel Simulator. The measured level crossing rate (LCR) seems to meet req (f), for v=60 km/hr, f = 450 MHz, and P=1 when measured over a 10 second sample. We've used Rs=2000 symb/s here, so x-axis of plot is 1 second in time.
+
+![LMR 60](doc/lmr_60.png)
+
+```
+octave:39>  multipath_samples("lmr60",8000, 2000, 1, 10, "h_lmr60.f32")
+Generating Doppler spreading samples...
+fd = 25.000
+path_delay_s = 2.0000e-04
+Nsecplot = 1
+Pav = 1.0366
+P = 1
+LCR_theory = 23.457
+LCR_meas = 24.400
+```
+
+# Single Carrier PSK Modem
+
+A single carrier PSK modem "back end" that connects the ML symbols to the radio.  This particular modem is written in Python, and can work with DC coupled and passband BBFM radios. It uses classical DSP, rather than ML.  Unlike the HF RADE waveform which used OFDM, this modem is single carrier.
+
+1. Run a single test with some plots, Eb/No=4dB, 100ppm sample clock offset, BER should be about 0.01:
+   ```
+   python3 -c "from radae import single_carrier; s=single_carrier(); s.run_test(100,sample_clock_offset_ppm=-100,plots_en=True,EbNodB=4)"
+   ```
+1. Run a suite of tests:
+   ```
+   ctest -V -R bbfm_sc
+   ```
+1. Create a file of BBFM symbols, 80 symbols every 40ms, plays expected output speech:
+    ```
+    ./bbfm_inference.sh model_bbfm_01/checkpoints/checkpoint_epoch_100.pth wav/brian_g8sez.wav - --write_latent z.f32
+    ```
+2. Sanity check of modem, BER test using digital, BPSK symbols, the symbols in z.f32 are replaced with BPSK symbols. `t.int16` is a real valued Fs=9600Hz sample file, that could be played into a FM radio.
+   ```
+   cat z.f32 | python3 sc_tx.py --ber_test > t.int16
+   cat t.int16 | python3 sc_rx.py --ber_test --plots > /dev/null
+   ```
+3. Send the BBFM symbols over the modem, and listen to results:
+    ```
+   cat z.f32 | python3 sc_tx.py > t.int16
+   cat t.int16 | python3 sc_rx.py > z_hat.f32
+   ./bbfm_rx.sh model_bbfm_01/checkpoints/checkpoint_epoch_100.pth z_hat.f32 -
+    ```
+4. Compare MSE of features passed through the system, first with z == z_hat, then with z passed through modem to get z_hat:
+   ```
+     python3 loss.py features_in.f32 features_out.f32
+     loss: 0.033
+     python3 loss.py features_in.f32 features_rx_out.f32
+     loss: 0.035
+   ```
+  This is a really good result, and likely inaudible. The `feature*.f32` files are produced as intermediate outputs form the `bbfm_inference.sh` and `bbfm_rx.sh` scripts.
+

CMakeLists.txt

@@ -168,13 +168,13 @@ add_test(NAME chirp_mpp
 # Low SNR ota_test.sh, with chirp measurement, AWGN
 add_test(NAME ota_test_awgn
         COMMAND sh -c "cd ${CMAKE_SOURCE_DIR}; \
-                       ./test/ota_test_cal.sh ~/codec2-dev/build_linux/ -21 0.4")
+                       ./test/ota_test_cal.sh ${CODEC2_DEV_BUILD_DIR} -21 0.4")
 
 # Low SNR ota_test.sh, with chirp measurement, MPP, high loss threshold as we only care about gross errors,
 # like stuck in false sync
 add_test(NAME ota_test_mpp
         COMMAND sh -c "cd ${CMAKE_SOURCE_DIR}; \
-                       ./test/ota_test_cal.sh ~/codec2-dev/build_linux/ -24 0.4 --mpp --freq -25")
+                       ./test/ota_test_cal.sh ${CODEC2_DEV_BUILD_DIR} -24 0.4 --mpp --freq -25")
 
 
 # Acquisition tests ------------------------------------------------------------------------------------
@@ -432,6 +432,7 @@ add_test(NAME radae_rx_aux_mpp
                        python3 loss.py features_in.f32 features_rx_out.f32 --loss 0.3 --clip_start 300")
                        set_tests_properties(radae_rx_aux_mpp PROPERTIES PASS_REGULAR_EXPRESSION "PASS")
 
+
 # Embedding Python in C callable library ------------------------------------------------------------------------------------------------------
 
 # Test Embedded version of streaming Tx, running with Python __main__
@@ -471,3 +472,40 @@ add_test(NAME radae_rx_embed_c
                        cat rx.f32 | PYTHONPATH='.' ${CMAKE_CURRENT_BINARY_DIR}/src/radae_rx > features_out.f32;
                        python3 loss.py features_in.f32 features_out.f32 --loss_test 0.15 --acq_time_test 0.5")
                        set_tests_properties(radae_rx_embed_c PROPERTIES PASS_REGULAR_EXPRESSION "PASS")
+
+
+# BBFM -----------------------------------------------------------------------------------------------
+
+# single carrier modem internal (inside single_carrier class) tests
+add_test(NAME bbfm_sc_internal
+        COMMAND sh -c "cd ${CMAKE_SOURCE_DIR}; \
+                       python3 -c 'from radae import single_carrier,single_carrier_tests; single_carrier_tests()'")
+                       set_tests_properties(bbfm_sc_internal PROPERTIES PASS_REGULAR_EXPRESSION "ALL PASS")
+
+# single carrier modem stand alone tx/rx, using BPSK symbol/BER test mode
+add_test(NAME bbfm_sc_ber
+        COMMAND sh -c "cd ${CMAKE_SOURCE_DIR}; \
+                       ./bbfm_inference.sh model_bbfm_01/checkpoints/checkpoint_epoch_100.pth wav/brian_g8sez.wav /dev/null --write_latent z.f32; \
+                       cat z.f32 | python3 sc_tx.py --ber_test > t.int16; \
+                       cat t.int16 | python3 sc_rx.py --ber_test --target_ber 0 > /dev/null")
+                       set_tests_properties(bbfm_sc_ber PROPERTIES PASS_REGULAR_EXPRESSION "PASS")
+
+# single carrier modem stand alone tx/rx, using continuously valued BBFM symbols, compare loss in decoded features
+add_test(NAME bbfm_sc_loss
+        COMMAND sh -c "cd ${CMAKE_SOURCE_DIR}; \
+                       ./bbfm_inference.sh model_bbfm_01/checkpoints/checkpoint_epoch_100.pth wav/brian_g8sez.wav /dev/null --write_latent z.f32; \
+                       cat z.f32 | python3 sc_tx.py | python3 sc_rx.py > z_hat.f32; \
+                       ./bbfm_rx.sh model_bbfm_01/checkpoints/checkpoint_epoch_100.pth z_hat.f32 /dev/null; \
+                       python3 loss.py features_in.f32 features_out.f32 --features_hat2 features_rx_out.f32 --compare")
+                       set_tests_properties(bbfm_sc_loss PROPERTIES PASS_REGULAR_EXPRESSION "PASS")
+
+# single carrier modem stand alone tx/rx with external 300-2700Hz band pass filter and no noise, measure loss.  In practice, SNR will be
+# quite high, so channel distortions other than noise may dominate
+add_test(NAME bbfm_sc_bpf_loss
+        COMMAND sh -c "cd ${CMAKE_SOURCE_DIR}; \
+                       ./bbfm_inference.sh model_bbfm_01/checkpoints/checkpoint_epoch_100.pth wav/brian_g8sez.wav /dev/null --write_latent z.f32; \
+                       cat z.f32 | python3 sc_tx.py | ${CODEC2_DEV_BUILD_DIR}/src/ch - - | python3 sc_rx.py > z_hat.f32; \
+                       ./bbfm_rx.sh model_bbfm_01/checkpoints/checkpoint_epoch_100.pth z_hat.f32 /dev/null; \
+                       python3 loss.py features_in.f32 features_out.f32 --features_hat2 features_rx_out.f32 --compare")
+                       set_tests_properties(bbfm_sc_bpf_loss PROPERTIES PASS_REGULAR_EXPRESSION "PASS")
+

README.md

@@ -52,6 +52,13 @@ The RDOVAE derived Python source code is released under the two-clause BSD licen
 | radae_tx.[py,sh] | streaming RADAE encoder and helper script |
 | radae_rx.[py,sh] | streaming RADAE decoder and helper script |
 | resource_est.py | WIP estimate CPU/memory resources |
+| radae_base.py | Shared ML code between models |
+| radae/bbfm.py | Baseband FM PyTorch model |
+| train_bbfm.py | Training for BBFM model |
+| inference_bbfm.py | Baseband FM inference |
+| inference_bbfm.sh | helper script for infereence_bbfm.sh |
+| fm.m | Octave analog FM mod/demod simulation |
+| analog_bbfm.sh | helper script for analog FM simulation |
 
 # Installation
 

analog_bbfm.sh

@@ -0,0 +1,44 @@
+#!/bin/bash -x
+#
+# Analog FM simulation, for comparison to ML BBFM
+
+CODEC2_DEV=${CODEC2_DEV:-${HOME}/codec2-dev}
+OPUS=build/src
+PATH=${PATH}:${OPUS}:${CODEC2_DEV}/build_linux/src
+gain=6
+
+which ch >/dev/null || { printf "\n**** Can't find ch - check CODEC2_PATH **** \n\n"; exit 1; }
+
+source utils.sh
+
+if [ $# -lt 3 ]; then
+    echo "usage (write output to file):"
+    echo "  ./analog_bbfm.sh in.wav out.wav CNRdB"
+    echo "usage (play output with aplay):"
+    echo "  ./analog_bbfm.sh in.wav - CNRdB"
+    exit 1
+fi
+
+if [ ! -f $1 ]; then
+    echo "can't find $1"
+    exit 1
+fi
+
+input_speech=$1
+output_speech=$2
+CNRdB=$3
+
+tmp_in=$(mktemp)
+tmp_out=$(mktemp)
+tmp_fm=$(mktemp)
+
+# We use hilbert clipper in ch util for speech compressor. Octave FM simulation uses 48 kHz sample rate.
+# input wav -> 300-3100Hz Fs=8kHz -> ch compressor -> 300-3100Hz Fs=48kHz -> FM mod/demod
+sox ${input_speech} -t .s16 -r 8000 -c 1 - sinc 0.3-3.1k | ch - - --clip 16384 --gain $gain 2>/dev/null | sox -t .s16 -r 8000 -c 1 - -t .s16 -r 48000 ${tmp_in} sinc 0.3-3.1k
+echo "fm; pkg load signal; fm_mod_file('${tmp_fm}','${tmp_in}',${CNRdB}); fm_demod_file('${tmp_out}','${tmp_fm}'); quit;" | octave-cli -qf
+
+if [ $output_speech == "-" ]; then
+    aplay ${tmp_out} -r 48000 -f S16_LE 2>/dev/null
+elif [ $output_speech != "/dev/null" ]; then
+    sox -t .s16 -r 48000 -c 1 ${tmp_out} -r 8000 ${output_speech}
+fi

bbfm_bpf.m

@@ -0,0 +1,30 @@
+% Octave script to explore BPF of basband symbols
+
+function bbfm_bpf()
+    pkg load signal;
+    Fs = 8000; T = 1/Fs; Rs = 2000; M = Fs/Rs; Nsym = 6; alpha = 0.25;
+    rn = gen_rn_coeffs(alpha, T, Rs, Nsym, M);
+    bpf = fir2(100,[0 250 350 3000 3100 4000]/(Fs/2),[0.001 0.001 1 1 0.001 0.001]);
+
+    figure(1); clf; 
+    subplot(211);
+    [h,w] = freqz(rn); plot(w*Fs/(2*pi), 20*log10(abs(h))); grid('minor'); ylabel('RRC');
+    subplot(212);
+    [h,w] = freqz(bpf); plot(w*Fs/(2*pi), 20*log10(abs(h))); grid; ylabel('BPF');
+
+    Nsymb = 1000;
+    tx_symb = 1 - 2*(rand(Nsymb,1)>0.5);
+    tx_pad = zeros(1,M*Nsymb);
+    tx_pad(1:M:end) = tx_symb;
+    tx = filter(rn,1,tx_pad);
+    tx = filter(bpf,1,tx)
+    rx = filter(rn,1,tx);
+    rx_symb = rx(1:M:end);
+    figure(2); clf; 
+    subplot(211); stem(tx_symb(1:100)); ylabel('Tx symbols');
+    subplot(212); stem(rx_symb(1:100)); ylabel('Rx Symbols');
+
+    figure(3); clf;
+    plot(20*log10(abs(fft(tx)(1:length(tx)/2))))
+end
+

bbfm_inference.py

@@ -0,0 +1,142 @@
+"""
+/* Copyright (c) 2024 modifications for radio autoencoder project
+   by David Rowe */
+
+/* Copyright (c) 2022 Amazon
+   Written by Jan Buethe */
+/*
+   Redistribution and use in source and binary forms, with or without
+   modification, are permitted provided that the following conditions
+   are met:
+
+   - Redistributions of source code must retain the above copyright
+   notice, this list of conditions and the following disclaimer.
+
+   - Redistributions in binary form must reproduce the above copyright
+   notice, this list of conditions and the following disclaimer in the
+   documentation and/or other materials provided with the distribution.
+
+   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+   ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+   A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
+   OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+"""
+
+import os
+import argparse
+
+import numpy as np
+import torch
+
+from radae import BBFM, distortion_loss
+
+parser = argparse.ArgumentParser()
+
+parser.add_argument('model_name', type=str, help='path to model in .pth format')
+parser.add_argument('features', type=str, help='path to input feature file in .f32 format')
+parser.add_argument('features_hat', type=str, help='path to output feature file in .f32 format')
+parser.add_argument('--latent-dim', type=int, help="number of symbols produces by encoder, default: 80", default=80)
+parser.add_argument('--cuda-visible-devices', type=str, help="set to 0 to run using GPU rather than CPU", default="")
+parser.add_argument('--write_latent', type=str, default="", help='path to output file of latent vectors z[latent_dim] in .f32 format')
+parser.add_argument('--CNRdB', type=float, default=100, help='FM demod input CNR in dB')
+parser.add_argument('--passthru', action='store_true', help='copy features in to feature out, bypassing ML network')
+parser.add_argument('--h_file', type=str, default="", help='path to rate Rs fading channel magnitude samples, rate Rs time steps by Nc=1 carriers .f32 format')
+parser.add_argument('--write_CNRdB', type=str, default="", help='path to output file of CNRdB per sample after fading in .f32 format')
+parser.add_argument('--loss_test', type=float, default=0.0, help='compare loss to arg, print PASS/FAIL')
+args = parser.parse_args()
+
+# set visible devices
+os.environ['CUDA_VISIBLE_DEVICES'] = args.cuda_visible_devices
+
+# device
+device = torch.device("cuda") if torch.cuda.is_available() else torch.device("cpu")
+
+latent_dim = args.latent_dim
+
+# not exposed
+nb_total_features = 36
+num_features = 20
+num_used_features = 20
+
+# load model from a checkpoint file
+model = BBFM(num_features, latent_dim, args.CNRdB)
+checkpoint = torch.load(args.model_name, map_location='cpu')
+model.load_state_dict(checkpoint['state_dict'], strict=False)
+checkpoint['state_dict'] = model.state_dict()
+
+# load features from file
+feature_file = args.features
+features_in = np.reshape(np.fromfile(feature_file, dtype=np.float32), (1, -1, nb_total_features))
+nb_features_rounded = model.num_10ms_times_steps_rounded_to_modem_frames(features_in.shape[1])
+features = features_in[:,:nb_features_rounded,:]
+features = features[:, :, :num_used_features]
+features = torch.tensor(features)
+print(f"Processing: {nb_features_rounded} feature vectors")
+
+# default rate Rb multipath model H=1
+Rb = model.Rb
+Nc = 1
+num_timesteps_at_rate_Rs = model.num_timesteps_at_rate_Rs(nb_features_rounded)
+H = torch.ones((1,num_timesteps_at_rate_Rs,Nc))
+
+# user supplied rate Rs multipath model, sequence of H magnitude samples
+if args.h_file:
+   H = np.reshape(np.fromfile(args.h_file, dtype=np.float32), (1, -1, Nc))
+   print(H.shape, num_timesteps_at_rate_Rs)
+   if H.shape[1] < num_timesteps_at_rate_Rs:
+      print("Multipath H file too short")
+      quit()
+   H = H[:,:num_timesteps_at_rate_Rs,:]
+   H = torch.tensor(H)
+
+if __name__ == '__main__':
+
+   if args.passthru:
+      features_hat = features_in.flatten()
+      features_hat.tofile(args.features_hat)
+      quit()
+
+   # push model to device and run test
+   model.to(device)
+   features = features.to(device)
+   H = H.to(device)
+   output = model(features,H)
+
+   # Lets check actual SNR at output of FM demod
+   tx_sym = output["z_hat"].cpu().detach().numpy()
+   S = np.mean(np.abs(tx_sym)**2)
+   N = np.mean(output["sigma"].cpu().detach().numpy()**2)
+   SNRdB_meas = 10*np.log10(S/N)
+   print(f"SNRdB Measured: {SNRdB_meas:6.2f}")
+
+   features_hat = output["features_hat"][:,:,:num_used_features]
+   features_hat = torch.cat([features_hat, torch.zeros_like(features_hat)[:,:,:16]], dim=-1)
+   features_hat = features_hat.cpu().detach().numpy().flatten().astype('float32')
+   features_hat.tofile(args.features_hat)
+
+   loss = distortion_loss(features,output['features_hat']).cpu().detach().numpy()[0]
+   print(f"loss: {loss:5.3f}")
+   if args.loss_test > 0.0:
+      if loss < args.loss_test:
+         print("PASS")
+      else:
+         print("FAIL")
+
+   # write output symbols (latent vectors)
+   if len(args.write_latent):
+      z_hat = output["z_hat"].cpu().detach().numpy().flatten().astype('float32')
+      z_hat.tofile(args.write_latent)
+
+   # write CNRdB after fading
+   if len(args.write_CNRdB):
+      CNRdB = output["CNRdB"].cpu().detach().numpy().flatten().astype('float32')
+      CNRdB.tofile(args.write_CNRdB)
+