diff --git a/250319_bbfm_lmr60/checkpoints/checkpoint_epoch_100.pth b/250319_bbfm_lmr60/checkpoints/checkpoint_epoch_100.pth
new file mode 100644
index 00000000..37c2d4da
Binary files /dev/null and b/250319_bbfm_lmr60/checkpoints/checkpoint_epoch_100.pth differ
diff --git a/BBFM.md b/BBFM.md
deleted file mode 100644
index c77bef03..00000000
--- a/BBFM.md
+++ /dev/null
@@ -1,83 +0,0 @@
-# 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 from the `bbfm_inference.sh` and `bbfm_rx.sh` scripts.
-
diff --git a/CMakeLists.txt b/CMakeLists.txt
index fc90682e..9792590a 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -586,3 +586,10 @@ add_test(NAME bbfm_sc_bpf_loss
                        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")
 
+# BBFM streaming rx
+add_test(NAME bbfm_stream
+        COMMAND sh -c "cd ${CMAKE_SOURCE_DIR}; \
+                       ./bbfm_inference.sh 250319_bbfm_lmr60/checkpoints/checkpoint_epoch_100.pth wav/brian_g8sez.wav /dev/null --write_latent z_hat.f32; \
+                       cat z_hat.f32 | python3 bbfm_rx_stream.py 250319_bbfm_lmr60/checkpoints/checkpoint_epoch_100.pth > features_rx_stream.f32; \
+                       python3 loss.py features_in.f32 features_out.f32 --features_hat2 features_rx_stream.f32 --compare")
+                       set_tests_properties(bbfm_stream PROPERTIES PASS_REGULAR_EXPRESSION "PASS")
diff --git a/README.md b/README.md
index 023b9b3a..76d40b7e 100644
--- a/README.md
+++ b/README.md
@@ -1,74 +1,45 @@
-# Radio Autoencoder - RADAE
+# RADE for Land Mobile Radio: A Neural Codec for Transmission of Speech over Baseband FM Radio Channels
 
-A hybrid Machine Learning/DSP system for sending speech over HF radio channels.  The RADAE encoder takes vocoder features such as pitch, short term spectrum and voicing and generates a sequence of analog PSK symbols.  These are placed on OFDM carriers and sent over the HF radio channel. The decoder takes the received PSK symbols and produces vocoder features that can be sent to a speech decoder.  The system is trained to minimise the distortion of the vocoder features over HF channels.  Compared to classical DSP approaches to digital voice, the key innovation is jointly performing transformation, prediction, quantisation, channel coding, and modulation in the ML network.
+A version of the Radio Autoencoder (RADE) designed for the baseband FM channel for land mobile radio (LMR).
+This branch contains source code to support the paper: 
 
-## Scope 
+D. Rowe, T. Bece, [RADE for Land Mobile Radio: A Neural Codec for Transmission of Speech over Baseband FM Radio Channels](URL TBD)
 
-This repo is intended to support the authors experimental work, with just enough information for the advanced experimenter to reproduce aspects of the work.  The focus is on waveform development, not software configuration.  It is not intended to be a polished distribution for general use or to work across multiple Linux distros and operating systems - that will come later.  Unless otherwise stated, the code is this repo is intended to run only on Ubuntu Linux 22 on a non-virtual machine.
+The RADE source code is released under the two-clause BSD license.
 
-# Quickstart
-
-1. Installation section
-1. Inference section
-1. If you would like to transmit/receive test files over HF radios: Over the Air/Over the Cable (OTA/OTC)
-
-# Attributions and License
-
-This software was derived from RDOVAE Python source (Github xiph/opus.git opus-ng branch opus/dnn/torch/rdovae):
-
-J.-M. Valin, J. Büthe, A. Mustafa, [Low-Bitrate Redundancy Coding of Speech Using a Rate-Distortion-Optimized Variational Autoencoder](https://jmvalin.ca/papers/valin_dred.pdf), *Proc. ICASSP*, arXiv:2212.04453, 2023. ([blog post](https://www.amazon.science/blog/neural-encoding-enables-more-efficient-recovery-of-lost-audio-packets))
-
-The RDOVAE derived Python source code is released under the two-clause BSD license.
+For acronyms RADAE and RADE can be used interchangeably.  RADAE was the original name of the project.
 
 # Files
 
 | File | Description |
 | --- | --- |
-| radae/radae.py | ML model and channel simulation |
+| radae/bbfm.py | RADE BBFM ML model and channel simulation |
 | radae/dataset.py | loading data for training |
-| train.py | trains models |
-| inference.py | Testing models, injecting channel impairments, simulate modem sync |
-| rx.py | Stand alone receiver, use inference.py as transmitter |
-| inference.sh | helper script for inference.py |
-| rx.sh | helper script for rx.py |
-| ofdm_sync.sh | generates curves to evaluate classical DSP sync performance |
-| evaluate.sh | script to compare radae to SSB, generates plots and speech samples |
-| evaluate_loop.sh | script to run evaluate.sh over a range of SNRs and channels |
+| radae_base.py | Shared ML code between models |
+| train_bbfm.py | trains models |
+| bbfm_inference.py | Testing models, injecting channel impairments |
+| bbfm_rx.py | Stand alone receiver, use bbfm_inference.py as transmitter |
+| bbfm_inference.sh | helper script for bbfm_inference.py |
+| bbfm_rx.sh | helper script for bbfm_rx.py |
+| bbfm_rx_stream.py | streaming receiver |
 | doppler_spread.m | Octave script to generate Doppler spreading samples |
 | load_f32.m | Octave script to load .f32 samples |
-| multipath_samples.m | Octave script to generate multipath magnitude sample over a time/freq grid |
-| plot_specgram.m | Plots sepctrogram of radae modem signals |
+| multipath_samples.m | Octave script to generate multipath channel simulation samples |
 | radae_plots.m | Helper Octave script to generate various plots |
-| radio_ae.[tex,pdf] | Latex documenation |
-| ota_test.sh | Script to automate Over The Air (OTA) testing |
-| Radio Autoencoder Waveform Design.ods | Working for OFDM waveform, including pilot and cyclic prefix overheads |
-| compare_models.sh | Builds loss versus Eq/No curves for models to objectively compare |
-| est_snr.py | Prototype pilot sequence based SNR estimator - doesn't work for multipath |
 | test folder | Helper scripts for ctests |
 | loss.py | Tool to calculate mean loss between two feature files, a useful objective measure |
-| ml_pilot.py | Training low PAPR pilot sequence |
-| stateful_decoder.[py,sh] | Inference test that compares stateful to vanilla decoder |
-| stateful_encoder.[py,sh] | Inference test that compares stateful to vanilla encoder |
-| 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
 
+This code is designed to run on Ubuntu 22.
+
 ## Packages
 
 sox, python3, python3-matplotlib and python3-tqdm, octave, octave-signal, cmake.  Pytorch should be installed using the instructions from the [pytorch](https://pytorch.org/get-started/locally/) web site. 
 
 ## codec2-dev
 
-Supplies some utilities used for `ota_test.sh` and `evaluate.sh`
+Supplies some utilities used for the automated ctests.
 ```
 cd ~
 git clone https://github.com/drowe67/codec2-dev.git
@@ -76,7 +47,6 @@ cd codec2-dev
 mkdir build_linux
 cd build_linux
 cmake -DUNITTEST=1 ..
-make ch mksine tlininterp
 ```
 
 ## RADAE
@@ -91,510 +61,192 @@ cd build
 cmake ..
 make
 ```
+# BBFM ML encoder/decoder
 
-### Building on Windows
-
-While most of RADAE is in Python, there is a `lpcnet_demo` application
-that is required to be compiled. To do this for Windows, you can run
-something like the following from a Linux machine:
-
-```
-wget https://github.com/mstorsjo/llvm-mingw/releases/download/20240619/llvm-mingw-20240619-ucrt-ubuntu-20.04-x86_64.tar.xz
-tar xzf llvm-mingw-20240619-ucrt-ubuntu-20.04-x86_64.tar.xz
-export PATH=`pwd`/llvm-mingw-20240619-ucrt-ubuntu-20.04-x86_64/bin:$PATH
-export RADAE_PATH=`pwd`/radae
-cd $RADAE_PATH
-mkdir build_windows
-cd build_windows
-cmake -DCMAKE_TOOLCHAIN_FILE=$RADAE_PATH/cross-compile/mingw-llvm-x86_64.cmake ..
-make
-```
-
-(Replace `x86_64` in `mingw-llvm-x86_64.cmake` with either `i686` or `aarch64` for 32-bit x86 or 64-bit ARM, respectively.)
-
-Once done, `lpcnet_demo.exe` will be inside the `src` folder inside `build_windows`.
-
-#### Limitations
-
-* ctests are untested and likely do not work without additional changes.
-* Visual Studio is not supported, only MinGW.
-* Generating a Windows installer is currently not supported. `lpcnet_demo.exe` is intended to be included with other applications built with MinGW (such as freedv-gui).
-
-# Automated Tests
-
-The `cmake/ctest` framework is being used as a build and test framework. The command lines in `CmakeLists.txt` are a good source of examples, if you are interested in running the code in this repo. The ctests are a work in progress and may not pass on all systems (see Scope above).
-
-To run the cests:
-```
-cd radae/build
-ctest
-```
-To list tests `ctest -N`, to run just one test `ctest -R inference_model5`, to run in verbose mode `ctest -V -R inference_model5`.  You can change the paths to `codec2-dev` on the `cmake` command line:
-```
-cmake -DCODEC2_DEV=~/tmp/codec2-dev ..
-```
-A lot of the tests generate a float IQ sample file.  You can listen to this file with: 
-```
-cat rx.f32 | python3 f32toint16.py --real --scale 8192 | play -t .s16 -r 8000 -c 1 - bandpass 300 2000
-```
-The scaling `--scale` is required as the low SNRs mean the noise peak amplitude can clip 16 bit samples if not carefully scaled.
-
-# Inference
-
-`inference.py` is used for inference, which has been wrapped up in a helper script `inference.sh`.  Inference runs by default on the CPU, but will run on the GPU with the `--cuda-visible-devices 0` option.
-
-1. Generate `out.wav` at the default Eb/No = 100 dB:
-   ```
-   ./inference.sh model05/checkpoints/checkpoint_epoch_100.pth wav/brian_g8sez.wav out.wav
-   ```
-
-1. Play output sample to your default `aplay` sound device at BPSK Eb/No = 3dB:
-   ```
-   ./inference.sh model05/checkpoints/checkpoint_epoch_100.pth wav/brian_g8sez.wav - --EbNodB 3
-   ```
+Quickstart - ignore the first two steps (training), and use pre-trained `250319_bbfm_lmr60` model.
 
-1. Vanilla LPCNet-fargan (ie no analog VAE) for comparison:
+1. Generate 10 hours of fading samples for training:
    ```
-   ./inference.sh model05/checkpoints/checkpoint_epoch_100.pth wav/brian_g8sez.wav - --passthru
+   octave:67> multipath_samples("lmr60",8000, 2000, 1, 10*60*60, "h_lmr60_train.f32")
    ```
-
-1. Multipath demo at approx 0dB B=3000 Hz SNR. First generate multipath channel samples using GNU Octave (only need to be generated once): 
-   ```
-   octave:85> Rs=50; Nc=20; multipath_samples("mpp", Rs, Rs, Nc, 60, "h_mpp.f32")
-   $ ./inference.sh model05/checkpoints/checkpoint_epoch_100.pth ~/LPCNet/wav/all.wav tmp.wav --EbNodB 3 --write_latent z_hat.f32 --h_file h_mpp.f32
-   ```
-   Then use Octave to plot scatter diagram using z_hat latents from channel:
-   ```
-   octave:91> radae_plots; do_plots('z_hat.f32') 
-   ```
-
-# Multipath rate Fs
-
-1. Baseline no noise simulation on Multipath Poor MPP channel:
-   ```
-   octave:85> Fs=8000; Rs=50; Nc=20; multipath_samples("mpp", Fs, Rs, Nc, 60, "h_mpp.f32","g_mpp.f32")
-   ./inference.sh model05/checkpoints/checkpoint_epoch_100.pth wav/peter.wav /dev/null --rate_Fs --write_latent z.f32 --write_rx rx.f32 --pilots --pilot_eq --eq_ls --ber_test --EbNo 100 --g_file g_mpp.f32 --cp 0.004
-   octave:87> radae_plots; do_plots('z.f32','rx.f32')
-   ```
-
-1. Multipath Disturbed (MPD) demo:
-   ```
-   ./inference.sh model17/checkpoints/checkpoint_epoch_100.pth wav/brian_g8sez.wav brian_g8sez_mpd_snr3dB.wav --rate_Fs --pilots --pilot_eq --eq_ls --cp 0.004 --bottleneck 3 --EbNodB 6 --g_file g_mpd.f32 --write_rx rx.f32
-   ```
-   Optional plots (e.g. spectrogram):
-   ```
-   octave:40> radae_plots; do_plots('z.f32','rx.f32')
-   ```
-   Listen to "off air" signal.
-   ```
-   cat rx.f32 | python3 f32toint16.py --real --scale 8192 | sox -t .s16 -r 8000 -c 1 - brian_g8sez_mpd_snr3dB_rx.wav sinc 0.3-2.7k
-   ```
-
-# Evaluate script
-
-Automates joint simulation of SSB and RADAE, generates wave files and spectrograms.  Can adjust noise for equal C/No or equal P/No.
-
-1. Run `peter.wav`` through RADAE and SSB using Eb/No=16dB set point, calibrating SSB noise such that P/No is the same for both samples
-   ```
-   ./evaluate.sh model17/checkpoints/checkpoint_epoch_100.pth wav/peter.wav 240521 16 --bottleneck 3 -d --peak
-   ```
-1. As above but MPP channel using Eb/No=6dB set point
-   ```
-   ./evaluate.sh model17/checkpoints/checkpoint_epoch_100.pth wav/peter.wav 240521 6 --bottleneck 3 -d --peak --g_file g_mpp.f32
-   ```
-1. With dim=40 mixed rate model, and we run Eb/No 3dB higher than a dim=80 model.  `g_nc20_mpp.f32` is a time domain fading file so works for Nc=20 and Nc=10.
-   ```
-   ./evaluate.sh model18/checkpoints/checkpoint_epoch_100.pth wav/peter.wav 240524 9 --bottleneck 3 -d --peak --latent_dim 40 --g_file g_mpp.f32
-   ```
-
-# Simulation of Seperate Tx and Rx
-
-We separate the system into a transmitter `inference.py` and stand alone receiver `rx.py`.  These examples test the OFDM waveform, including pilot symbol insertion, cyclic prefix, least squares phase EQ, and coarse magnitude EQ.
-
-BER tests are useful to calibrate the system, and measure loss from classical DSP based synchronisation.  We have a good theoretical model for the expected BER on AWGN and multipath channels.
-
-1. First generate no-noise reference symnbols for BER measurement `z_100dB.f32`:
-   ```
-   ./inference.sh model05/checkpoints/checkpoint_epoch_100.pth wav/peter.wav /dev/null --rate_Fs --pilots --write_latent z_100dB.f32 --write_rx rx_100dB.f32 --EbNodB 100 --cp 0.004 --pilot_eq --eq_ls --ber_test
-   ```
-   `rx_100dB.f32` is the rate Fs IQ sample file, the actual modem signal we could send over the air.
-
-1. The file `z_100dB.f32` can then be used as a reference to measure BER at the receiver, e.g. using the no noise `rx_100dB.f32` sample:
-   ```
-   ./rx.sh model05/checkpoints/checkpoint_epoch_100.pth rx_100dB.f32 /dev/null --pilots --pilot_eq --cp 0.004 --plots --time_offset -16 --coarse_mag --ber_test z_100dB.f32
-   ```
-
-1. An AWGN channel at Eb/No = 0dB, first generate `rx_0dB.f32`:
-   ```
-   ./inference.sh model05/checkpoints/checkpoint_epoch_100.pth wav/peter.wav /dev/null --rate_Fs --pilots --write_rx rx_0dB.f32 --EbNodB 0 --cp 0.004 --pilot_eq --eq_ls --ber_test
-   ```
-   Then demodulate
+   Note 10 hours << 205 hours in the speech training dataset.  The same fading data is therefore repeated 205/10 times in each training epoch. I think 10 hours or so might be the max I can generate due to memory limitations in the current Octave code (TBC). It should be enough, based on argument for dataset length with similar models used for HF fading (ITU-R F.1487) which suggests a test length of 3000*(1/Doppler Spread Hz), which for 60 km/hr is 3000/25 = 120 seconds.
+   
+1. If you wish to perform training, a serious NVIDIA GPU is required - the author used a RTX4090. Training with fading (multipath):
    ```
-   ./rx.sh model05/checkpoints/checkpoint_epoch_100.pth rx_0dB.f32 /dev/null --pilots --pilot_eq --cp 0.004 --plots --time_offset -16 --coarse_mag --ber_test z_100dB.f32
+   ./lpcnet_demo -features training_input.pcm training_features_file.f32
+   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 training_features_file.f32 250319_bbfm_lmr60 --RdBm -100 --plot_loss --range_RdBm --h_file h_lmr60_train.f32
    ```
-   This will give a BER of around 0.094, compared to 0.079 theoretical for BPSK, a loss of about 1dB due to non-ideal synchronisation.
 
-1. Compare to the BER with the ideal phase estimate (pilot based EQ disabled):
-   ```
-   ./inference.sh model05/checkpoints/checkpoint_epoch_100.pth wav/peter.wav t.wav --pilots --rate_Fs --EbNodB 0 --ber_test
+1. Inference (runs encoder and decoder, plays result to sound card, and outputs symbols `z_hat.f32`):
    ```
-   Which is exactly the theoretical 0.078. Note the ideal BER for AWGN is given by `BER = 0.5*erfc(sqrt(Eb/No))`, where Eb/No is the linear Eb/No.
-
-1. Lets introduce frequency, and magnitude (gain) offsets typical of real radio channels:
+   ./bbfm_inference.sh 250319_bbfm_lmr60/checkpoints/checkpoint_epoch_100.pth wav/brian_g8sez.wav - --write_latent z_hat.f32
    ```
-   /inference.sh model05/checkpoints/checkpoint_epoch_100.pth wav/peter.wav /dev/null --rate_Fs --pilots --write_rx rx_0dB.f32 --EbNodB 0 --cp 0.004 --pilot_eq --eq_ls --ber_test --freq_offset 2 --gain 0.1
+1. Inference (-120dBm, fading, per sample Rx levels written to r.f32 for plotting):
    ```
-   ./rx.sh model05/checkpoints/checkpoint_epoch_100.pth rx_0dB.f32 /dev/null --pilots --pilot_eq --cp 0.004 --plots --time_offset -16 --coarse_mag --ber_test z_100dB.f32
+   octave:67> multipath_samples("lmr60", 8000, 2000, 1, 60, "h_lmr60.f32")
+   ./bbfm_inference.sh 250319_bbfm_lmr60/checkpoints/checkpoint_epoch_100.pth wav/brian_g8sez.wav - --h_file h_lmr60.f32 --RdBm -120 --write_RdBm r.f32
    ```
-   The acquisition system detected and corrected the 2Hz frequency offset, and the -20dB (--gain 0.1) magntitude offset, and the resulting BER was about the same at 0.095.
+1. Stand alone decoder, outputs speech from `z_hat.f32` to sound card:
+    ```
+    ./bbfm_rx.sh 250319_bbfm_lmr60/checkpoints/checkpoint_epoch_100.pth z_hat.f32 -
+    ```
+1. Or save speech out to a wave file:
+    ```
+    ./bbfm_rx.sh 250319_bbfm_lmr60/checkpoints/checkpoint_epoch_100.pth z_hat.f32 t.wav
+    ```
 
-1. Typical HF multipath channels evolve at around 1 Hz, so it's a good idea to use longer samples to get a meaningful average.  First generate the reference `z_100dB.f32` file for `all.wav`:
-   ```
-   ./inference.sh model05/checkpoints/checkpoint_epoch_100.pth wav/all.wav /dev/null --rate_Fs --pilots --write_latent z_100dB.f32 --write_rx rx_100dB.f32 --EbNodB 100 --cp 0.004 --pilot_eq --eq_ls --ber_test
-   ```
-   Check that it's doing sensible things with no noise (BER=0):
-   ```
-   ./rx.sh model05/checkpoints/checkpoint_epoch_100.pth rx_100dB.f32 /dev/null --pilots --pilot_eq --cp 0.004 --plots --time_offset -16 --coarse_mag --ber_test z_100dB.f32
-   ```
-   Lets generate a MPP sample `rx_100dB_mpp.f32` with no noise, and check the scatter diagram is a nice cross shape with BER close to 0:
-   ```
-   ./inference.sh model05/checkpoints/checkpoint_epoch_100.pth wav/all.wav /dev/null --rate_Fs --pilots --write_rx rx_100dB_mpp.f32 --EbNodB 100 --cp 0.004 --pilot_eq --eq_ls --ber_test --g_file g_mpp.f32
-   ```
-   ```
-   ./rx.sh model05/checkpoints/checkpoint_epoch_100.pth rx_100dB_mpp.f32 /dev/null --pilots --pilot_eq --cp 0.004 --plots --time_offset -16 --coarse_mag --ber_test z_100dB.f32
-   ```
-   You can listen to the modulated OFDM waveform over the MPP channel simulation with:
-   ```
-   play -r 8k -e float -b 32 -c 2 rx_100dB_mpp.f32 sinc 0.3-2.7k
-   ```
-   Lets add some channel impairments:
-   ```
-   ./inference.sh model05/checkpoints/checkpoint_epoch_100.pth wav/all.wav /dev/null --rate_Fs --pilots --write_rx rx_0dB_mpp.f32 --EbNodB 0 --cp 0.004 --pilot_eq --eq_ls --ber_test --g_file g_mpp.f32
-   ```
+1. Streaming decoder, reads a stream of z_hat float[80] vectors and synthesises decoded speech:
    ```
-   ./rx.sh model05/checkpoints/checkpoint_epoch_100.pth rx_0dB_mpp.f32 /dev/null --pilots --pilot_eq --cp 0.004 --plots --time_offset -16 --coarse_mag --ber_test z_100dB.f32
+   cat z_hat.f32 | python3 bbfm_rx_stream.py 250319_bbfm_lmr60/checkpoints/checkpoint_epoch_100.pth | build/src/lpcnet_demo -fargan-synthesis - - | aplay -f S16_LE -r 16000
    ```
-   Which gives us a BER of 0.172, about 1.5dB from the ideal Rayleigh multipath channel BER of 0.15 (which the MPP model approximates).
 
-# Over the Air/Over the Cable (OTA/OTC)
+1. Plot sequence of received symbols:
+    ```
+    octave:4> radae_plots; do_plots_bbfm('z_hat.f32')
+    ```
 
-1. Example of `ota_test.sh` script. `ota_test.sh -x` generates `tx.wav` which contains the simulated SSB and radae modem signals ready to run through a HF radio.  We add noise to create `rx.wav`, then use `ota_test.sh -r` to generate the demodulated audio files `rx_ssb.wav` and `rx_radae.wav`:
-   ```
-   ./ota_test.sh wav/peter.wav -x 
-   ~/codec2-dev/build_linux/src/ch tx.wav - --No -20 | sox -t .s16 -r 8000 -c 1 - rx.wav
-   ./ota_test.sh -r rx.wav
-   aplay rx_ssb.wav rx_radae.wav
-   ```
+# Faded (multipath) channel simulation
 
-1. Testing OTA over HF channels. Using my IC7200 as the Tx station:
-   ```
-   ./ota_test.sh wav/david_vk5dgr.wav -g 6 -t -d -f 14236
-   ```
-   The `-g 6` is the SSB compressor gain (default 6 so in this case optional); this can be adjusted by experiment, e.g. listening to the `tx.wav` file, and looking for signs of a compressed waveform on Audacity.  To receive the signal I tune into a convenient KiwiSDR, and manually start recording when my radio starts transmitting.  I stop recording when I hear the transmission end.  This will result in a wave file being downloaded.  It's a good idea to trim any excess off the start and end of the rx wave file. It can be decoded with:
-   ```
-   ./ota_test.sh -d -r ~/Downloads/kiwisdr_usb.wav
-   ```
-   The output will be a `~/Downloads/kiwisdr_usb_radae.wav` and `~/Downloads/kiwisdr_usb_ssb.wav`, which you can listen to and compare, `~/Downloads/kiwisdr_usb_spec.png` is the spectrogram.  The C/No will be estimated and displayed but this is unreliable at present for non-AWGN channels.  The `ota_test.sh` script is capable of automatically recording from KiwiSDRs, however this code hasn't been debugged yet.
+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.
 
-# Streaming
+![LMR 60](doc/lmr_60.png)
 
-`radae_rx.py` is s streaming receiver that accepts IQ samples on stdin, and outputs z vectors on stdout.  To listen to an example decode:
-```
-./inference.sh model17/checkpoints/checkpoint_epoch_100.pth wav/brian_g8sez.wav /dev/null --rate_Fs --pilots --pilot_eq --eq_ls --cp 0.004 --bottleneck 3 --write_rx rx.f32
-cat rx.f32 | python3 radae_rx.py model17/checkpoints/checkpoint_epoch_100.pth | build/src/lpcnet_demo -fargan-synthesis - - | aplay -f S16_LE -r 16000
-```
-To run just the core streaming decoder:
 ```
-cat rx.f32 | python3 radae_rx.py model17/checkpoints/checkpoint_epoch_100.pth > features_rx_out.f32
+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
 ```
-Full RADAE Streaming Rx in real time (fom off air audio samples to speaker):
-```
-cd ~/radae
-cat rx.f32 | python3 radae_rx.py model17/checkpoints/checkpoint_epoch_100.pth -v 1 | ./build/src/lpcnet_demo -fargan-synthesis - - | aplay -f S16_LE -r 16000
-```
-Ctest that measures % CPU used:
-```
-cd ~/radae/build
-ctest -V -R radae_rx_fargan
-<snip>
-run time:  6.41 duration:  9.82 percent CPU: 65.26
-```
-
-## Profiling example
 
-Total run time for 50 second `all/wav`:
+You can also generate fading samples with other speeds, e.g. 10 seconds of 120 km/hr:
 ```
-ctest -R radae_rx_dfdt
-time cat rx.f32 | python3 radae_rx.py model17/checkpoints/checkpoint_epoch_100.pth -v 0 --no_stdout
-```
-Profiling each function, using shorter wave file:
-```
-ctest -V -R radae_rx_basic
-cat rx.f32 | python3 -m cProfile -s time radae_rx.py model17/checkpoints/checkpoint_epoch_100.pth -v 0 --no_stdout | more
+octave:93> multipath_samples("lmr120", 8000, 2000, 1, 10, "h_lmr120_train.f32")
 ```
 
-# Training
+# Analog FM simulation
 
-This section is optional - pre-trained models that run on a standard laptop CPU are available for experimenting with RADAE. If you wish to perform training, a serious NVIDIA GPU is required - the author used a RTX4090.
+Analog FM simulation using same linearised FM model as we use for ML training/simulation:
 
-1. Generate a training features file using your speech training database `training_input.pcm`, we used 200 hours of speech from open source databases:
+1. Generate some Fs=8kHz LMR 60 samples in Octave:
    ```
-   ./lpcnet_demo -features training_input.pcm training_features_file.f32
+   multipath_samples("lmr60",8000, 8000, 1, 60, "h_lmr60_Fs_8000Hz.f32")
    ```
-   
-1. Vanilla fixed Eb/No:
+1. AWGN sim, play output to sound card (default RdBm = -100dBm)
    ```
-   python3 ./train.py --cuda-visible-devices 0 --sequence-length 400 --batch-size 512 --epochs 100 --lr 0.003 --lr-decay-factor 0.0001 --plot_loss training_features_file.f32 model_dir_name
+   ./bbfm_analog.sh wav/brian_g8sez.wav -
    ```
 
-1. Rate Rs with multipath, over range of Eb/No:
+1. AWGN sim, play output to sound card, lower RdBm
    ```
-   python3 ./train.py --cuda-visible-devices 0 --sequence-length 400 --batch-size 512 --epochs 100 --lr 0.003 --lr-decay-factor 0.0001 ~/Downloads/tts_speech_16k_speexdsp.f32 model05 --mp_file h_mpp.f32 --range_EbNo --plot_loss
+   ./bbfm_analog.sh wav/brian_g8sez.wav - --RdBm -110
    ```
 
-1. Rate Fs with simulated PA:
+1. AWGN sim, save output to wave file:
    ```
-   python3 ./train.py --cuda-visible-devices 0 --sequence-length 400 --batch-size 512 --epochs 100 --lr 0.003 --lr-decay-factor 0.0001 ~/Downloads/tts_speech_16k_speexdsp.f32 model06 --plot_loss --rate_Fs --range_EbNo
+   ./bbfm_analog.sh wav/brian_g8sez.wav brian_g8sez_analog_100dBm_awgn.wav
    ```
 
-1. Rate Fs with phase and freq offsets:
+1. Fading sim using LMR 60 km/hr model, save output to wave file:
    ```
-   python3 ./train.py --cuda-visible-devices 0 --sequence-length 400 --batch-size 512 --epochs 100 --lr 0.003 --lr-decay-factor 0.0001 ~/Downloads/tts_speech_16k_speexdsp.f32 model07 --range_EbNo --plot_loss --rate_Fs --freq_rand
+   ./bbfm_analog.sh wav/brian_g8sez.wav brian_g8sez_analog_100dBm_lmr60.wav --h_file h_lmr60_Fs_8000Hz.f32
    ```
 
-1. Generate `loss` versus Eb/No curves for a model:
-   ```
-   python3 ./train.py --cuda-visible-devices 0 --sequence-length 400 --batch-size 512 --epochs 1 --lr 0.003 --lr-decay-factor 0.0001 ~/Downloads/tts_speech_16k_speexdsp.f32 tmp --range_EbNo --plot_EqNo model05 --rate_Fs --initial-checkpoint model05/checkpoints/checkpoint_epoch_100.pth
-   ```
-   This runs another training epoch (results not used so saved in `tmp` folder), but the results won't change much as the network has converged.  While training across a range of Eb/No, it gathers stats on `loss` against Eq/No, and plots them in a PNG and dumps a text file.  The text output is useful for plotting curves from different training runs together. TODO: reconcile original Eb/No simulation parameter with latest model
-   that also trains constellations which means symbol energy Eq is a better parameter.
+# ASR Tests
 
-   Octave can be used to plot several loss curves together:
+1. Create LMR 60 sample at Fs=8000 Hz for the analog FM simulation, and Rs=2000 Hz for RADE. Assume Librispeech samples max length 10 seconds each, we want to use around 500 samples:
    ```
-   octave:120> radae_plots; loss_EqNo_plot("loss_models",'model05_loss_EbNodB.txt','m5_Rs_mp','model07_loss_EbNodB.txt','m7_Fs_offets','model08_loss_EbNodB.txt','m8_Fs')
+   octave:47> multipath_samples("lmr60", 8000, 8000, 1, 10*500, "h_lmr60_Fs_8000Hz.f32")
+   octave:48> multipath_samples("lmr60", 8000, 2000, 1, 10*500, "h_lmr60_Rs_2000Hz.f32")
    ```
 
-1. (May 2024) Training dim=80 mixed rate PAPR optimised model.  Note we need the Nc=20 version of the multipath H matrix `h_nc20_train_mpp.f32` as fading is aplies at rate Rs.  Bottleneck 3 is a tanh() on the magnitude of the complex rate Fs time domain samples.  The SNR ends up about 3dB higher, as discussed in the mixed rate/noise calibration section of the Latext doc: 
+1. Run a single RdBm analog FM ASR test with fading, using 10 samples:
    ```
-   python3 ./train.py --cuda-visible-devices 0 --sequence-length 400 --batch-size 512 --epochs 100 --lr 0.003 --lr-decay-factor 0.0001 ~/Downloads/tts_speech_16k_speexdsp.f32 model19 --bottleneck 3 --h_file h_nc20_train_mpp.f32 --range_EbNo --plot_loss
+   ./asr_test.sh fm -n 10 --RdBm -110 --h_file h_lmr60_Fs_8000Hz.f32
    ```
 
-1. (May 2024) Training dim=40 mixed rate PAPR optimised model.  Note we need the Nc=10 version of the multipath H matrix `h_nc10_train_mpp.f32`.  Bottleneck 3 is a tanh() on the magnitude of the complex rate Fs time domain samples.  We bump the range of Eb/Nos trained over by 3dB `--range_EbNo_start -3` as a 10 carrier waveform will have 3dB more energy per symbol.
+1. Locate simulation output files (e.g. for manual listening):
    ```
-   python3 ./train.py --cuda-visible-devices 0 --sequence-length 400 --batch-size 512 --epochs 100 --lr 0.003 --lr-decay-factor 0.0001 ~/Downloads/tts_speech_16k_speexdsp.f32 model18 --latent-dim 40 --bottleneck 3 --h_file h_nc10_train_mpp.f32 --range_EbNo_start -3 --range_EbNo --plot_loss
+   find /home/david/.cache/LibriSpeech/test-other/ -name *.flac
    ```
 
-1. (Aug 2024) Training model with auxillary/embedded data at 25 bits/s:
-    ```
-    python3 train.py --cuda-visible-devices 0 --sequence-length 400 --batch-size 512 --epochs 100 --lr 0.003 --lr-decay-factor 0.0001 ~/Downloads/tts_speech_16k_speexdsp.f32 model19_check3 --bottleneck 3 --h_file h_nc20_train_mpp.f32 --range_EbNo --plot_loss --auxdata
-    ```
-
-# Models & samples
-
-A log of models trained by the author.
-
-| Model | Description | Train at | Samples |
-| ---- | ---- | ---- | ---- |
-| model01 | trained at Eb/No 0 dB | Rs | - |
-| model02 | trained at Eb/No 10 dB | Rs | - |
-| model03 | --range_EbNo -2 ... 13 dB, modified sqrt loss |Rs | - |
-| model04 | --range_EbNo -2 ... 13 dB, orginal loss, noise might be 3dB less after calibration | Rs | - | 
-| model05 | --range_EbNo, --mp_file h_mpp.f32, sounds good on MPP and AWGN at a range of SNR - no pops | Rs | 240221_m5_Rs_mp |
-| model06 | --range_EbNo, --rate_Fs, trained on AWGN with PA model, PAPR about 1dB, OK at a range of Eb/No | Fs |240223_m6_Fs_papr |
-| model07 | --range_EbNo, -6 ... 14, --rate_Fs, AWGN freq, phase, gain offsets, some degredation at 0dB | Fs | 240301_m7_Fs_offets | Fs |
-| model08 | --range_EbNo, -6 ... 14, --rate_Fs, AWGN no offsets (vanilla rate Fs), similar to model 05 | Fs | 240301_m8_Fs | 
-| model05 | practical OFDM with 4ms CP and pilots, increased Rs', MPP, 4dB sync loss, speech dropping in and out | Fs | 240319_m5_Fs_mp | 
-| model05 | practical OFDM with 120ms modem frame, 4ms CP and pilots, reduced Rs', 1dB improvement, mooneer sample | Fs | 240320_m5_Fs | 
-| model05 | HF OTA tests of up to 2000km, including weak signal, EMI, NVIS | Fs | 240326_ota_hf | 
-| model09 | repeat of model05 as a sanity check, similar loss v Eq/No | Rs | |
-| model10 | First attempt at dim=40 1D bottleneck, AWGN, poor audio quality, relatively high loss v Eq/No, square constellation | Rs | |
-| model11 | dim=40 with 2D bottleneck #1, AWGN, good audio quality, similar to model05, good loss v Eq/No, circular constellation | Rs | |
-| model12 | dim=40 with 2D bottleneck #2, AWGN, a more sensible --range_EbNo_start 0, improved loss v Eq/No, circular constellation | Rs | |
-| model13 | dim=80 with 2D bottleneck 3 on rate Fs, AWGN, 0.5db PAPR, loss > m5, a few dB poorer at low SNR, very similar at high SNR | Fs | |
-| model14 | dim=80 with 2D bottleneck 3 on rate Fs, 10 hour --h_file h_nc20_test.f32 --range_EbNo_start 0, 0.7dB PAPR, "accident" as it introduces phase distortion with no EQ, but does a reasonable job (however speech quality < m5), handles phase and small freq offsets with no pilots, worth exploring further | Fs | |
-| model15 | repeat of model05/09 with 250 hour --h_file h_nc20_train_mpp.f32, after refactoring dataloader, loss v epoch curve v close to model09, ep 100 loss 0.150 | Rs | |
-| model16 | repeat of model05/09 with 10 hour --h_file h_nc20_test.f32, testing short h file, ep 100 loss 0.149 | Rs | |
-| model17 | `--bottleneck 3 --h_file h_nc20_train_mpp.f32` mixed rate Rs with time domain bottelneck 3, ep 100 loss 0.112 | Rs | 240601_m17 |
-| model18 | `--latent-dim 40 --bottleneck 3 --h_file h_nc10_train_mpp.f32 --range_EbNo_start -3` like model17 but dim 40, ep 100 loss 0.123 | Rs | 240601_m18 |
-| model05_auxdata | model05 (rate Rs h_nc20_train_mpp.f32) with --auxdata 100 bits/s see PR#13 | Rs | - |
-| model05_auxdata25 | model05 (rate Rs h_nc20_train_mpp.f32) with --auxdata 25 bits/s see PR#13 | Rs | - |
-| model19 | like model17 but with 25 bits/s auxdata, ep 100 loss 0.124 | Fs | - |
-| model19_check3 | model19 but loss function weighting for data symbols reduced fom 1/18 to 0.5/18, which reduced vocoder feature loss with just a small impact on BER.  Loss at various op points and channels very close to model17 | Fs | - |
-| model20 | model19_check3 but loss function weighting for pitch and corr doubled, attempt to improve rick samples.  Didn't help. |  Fs | - |
-| model21 | very based fixed EbNodB model like model02 but at 20dB at epoch 30, produced good quality speech on rick, loss 0.02, but not a practical solution |  Rs | - |
-| model22 | very based fixed EbNodB model like model02 but at 10dB at epoch 30, rick sample starting to get tonal/pitch shift artefact, in between, loss 0.045 |  Rs | - |
-
-Note the samples are generated with `evaluate.sh`, which runs inference at rate Fs. even if (e.g model 05), trained at rate Rs.
-
-# Specifications
-
-Using model19_check3 waveform:
-
-| Parameter | Value | Comment |
-| --- | --- | --- |
-| Audio Bandwith | 100-8000 Hz | |
-| RF Bandwidth | 1500 Hz (-6dB) | |
-| Tx Peak Average Power Ratio | < 1dB | |
-| Threshold SNR | -3dB | AWGN channel, 3000 Hz noise bandwidth |
-| Threshold C/No | 32 dBHz | AWGN channel |
-| Threshold SNR | 0dB | MPP channel (1Hz Doppler, 2ms path delay), 3000 Hz noise bandwidth |
-| Threshold C/No | 35 dBHz | MPP channel |
-| Frame size | 120ms | algorithmic latency |
-| Modulation | OFDM | discrete time, continously valued symbols |
-| Vocoder | FARGAN | low complexity ML vocoder |
-| Total Payload Symbol rate | 2000 Hz | payload data symbols, all carriers combined |
-| Number of Carriers | 30 | |
-| Per Carrier Symbol rate | 50 Hz | |
-| Cyclic Prefix | 4ms | |
-| Worst case channel | MPD: 2Hz Doppler spread, 4ms path delay | Two path Watterson model |
-| Mean acquisition time | < 1.5s | 0dB SNR MPP channel | 
-| Acquisition frequency range | +/- 50 Hz | |
-| Acquisition co-channel interference tolerence | -3dBC | Interfering sine wave level, <2s mean acquisition time |
-| Auxilary text channel | 25 bits/s | currently used for sync |
-| SNR measurement | No | |
-| Tx and Rx sample clock offset | 200ppm | e.g. Tx sample clock 8000 Hz, Rx sample clock 8001 Hz |
-
-# Web based Stored File Processing
-
-This section contains some notes on setting up a web server to run `ota_test.sh`.  The idea is to make it easier for non-Linux users to contribute to the stored file test program.  The general idea is a CGI script interfaces to `ota_test.sh` to perform the Tx and Rx processing.  We configure the web server so that the HTML forms and CGI scripts run in `~/public_html`.  The notes below are for Apache on Ubuntu 22. 
-
-1. The Python packages need to be available system wide , so `www-data` can use them: 
-   ```
-   sudo pip3 install torch torchvision torchaudio --index-url https://download.pytorch.org/whl/cu118
-   sudo -u www-data python3 -c "import torch"
-   ```
-
-   ```
-   sudo pip3 install matplotlib
-   sudo -u www-data python3 -c "import matplotlib"
-   ```
-   The presence of the packages can be checked by mimicing the www-data user (the last line in each step above should not fail if all is well).
-
-1. Configure Apache for CGI and serving pages from our `~/public_html` dir.
-   ```
-   sudo a2enmod cgid
-   sudo a2enmod userdir
-   sudo systemctl restart apache2
-   ```
-   We want html and cgi to run out of ~/public_html, so permissions have to be `755` and `www-data` has to be added to the users group.
-   ```
-   mkdir ~/public_html
-   chmod 755 public_html
-   sudo usermod -a -G <username> www-data
-   ```
-   To let CGI scripts run from ~/public_html I placed this in my `/etc/apache2/apache2.conf`:
-   ```
-   <Directory "/home/<username>/public_html">
-      Options +ExecCGI
-      AddHandler cgi-script .cgi
-   </Directory>  
-   ```   
-   Then restart apache as above.
-
-1. Create sym links to HTML/CGI scripts in `radae` repo, this allows the script to be part of the RADAE repo:
-   ```
-   cd ~/public_html
-   ln -s ~/radae/public_html/tx_form.html tx_form.html
-   ln -s ~/radae/public_html/tx_process.cgi tx_process.cgi
+1. Top level ASR test script to generate ASR results across a range of RdBm:
+   ```
+   ./asr_test_top.sh bbfm -n 100
    ``` 
-
-1. Edit the path to `CODEC2_DEV` in `radae/public_html/tx_process.cgi`:
+   Eyeball results:
    ```
-   my_env["CODEC2_DEV"] = "/home/<username>/codec2-dev"
+   ls 250610_asr*
+   cat 250610_asr_lmr60_bbfm.txt
+   ...
    ```
-
-1. Note that files created when the CGI process run (e.g. `/tmp/input.wav`) get put in a sandbox rather than directly in `/tmp`.  This is a systemd security feature.  You can find the files with:
+1. Plot curves, save to .png:
    ```
-   sudo find /tmp -name input.wav | xargs sudo ls -ld
-   -rw-r--r-- 1 www-data www-data 3918458 Aug 15 15:28 /tmp/systemd-private-2fcf85ad243b4da08d79d2e27e0375af-apache2.service-vDE2Dg/tmp/input.wav
+   octave:47> radae_plots; plot_wer_bbfm("250610","250610_bbfm_wer.png")
    ```
 
-1. Apache error log, good for viewing `ota_test.sh` progress and spotting any issues:
-   ```
-   tail -f /var/log/apache2/error.log
-   ```
+# Single Carrier PSK Modem
 
-# Real Time PTT
+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.
 
-WIP notes
-
-## Real Time decode using KiwiSDR
-
-1. Install pulse audio null module
-   ```
-   pactl load-module module-null-sink sink_name=vsink
-   ```
-1. Start your web browser, and open a tab to a KiwiSDR.
-1. Open `pavucontrol`, *Playback* tab, send web browser sound to NULL module.  Audio from web browser should go silent.
-1. We take the audio from the null device monitor output for the input to the RADAE Rx:
-   ```
-   parec --device=vsink.monitor --rate=8000 --channels=1 | python3 int16tof32.py --zeropad | python3 radae_rx.py model19_check3/checkpoints/checkpoint_epoch_100.pth -v 2 --auxdata | ./build/src/lpcnet_demo -fargan-synthesis - - | aplay -f S16_LE -r 16000
+1. Run a single test with some plots, Eb/No=4dB, 100ppm sample clock offset, BER should be about 0.01:
    ```
-1. Try transmitting a RADAE signal:
+   python3 -c "from radae import single_carrier; s=single_carrier(); s.run_test(100,sample_clock_offset_ppm=-100,plots_en=True,EbNodB=4)"
    ```
-   ./ota_test.sh -t radae_test.raw -d -f 7175
+1. Run a suite of tests:
    ```
-   Where `radae_test.raw` is a RADAE-only sample (i.e. without the chirp and SSB, copied from a temp file generated by `ota_test.sh -x`). If you can't open the SSB radio playback device to radio try closing `pavucontrol`.
-1. Other useful pulse audio commands:
+   ctest -V -R bbfm_sc
    ```
-   pactl list sinks short
-   pactl list sources short
-   pactl list modules
-   ```
-
-## Real Time Tx from mic to SSB radio
-
-Work in progress notes, needs a clean up once this settles down.
-
-1. Install null module as above.  Using Settings redirect system sounds to null so default analog sound output is free.
-
-1. Test headset mic to audio:
+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.
    ```
-   parec --device=17 --rate=16000 --channels=1 --latency=1024 | pacat --device=9 --rate=16000 --channels=1 --latency=1024
+   cat z.f32 | python3 sc_tx.py --ber_test > t.int16
+   cat t.int16 | python3 sc_rx.py --ber_test --plots > /dev/null
    ```
-   However this is unreliable, doesn't always start.
-
-1. Input from headset mic, save to file.
+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:
    ```
-   arecord --device "plughw:CARD=LX3000,DEV=0" -f S16_LE -c 1 -r 16000 | ./build/src/lpcnet_demo -features - - | python3 radae_tx.py model19_check3/checkpointscheckpoint_epoch_100.pth --auxdata | python3 f32toint16.py --real --scale 8192 > t.raw
+     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 from the `bbfm_inference.sh` and `bbfm_rx.sh` scripts.
 
-1. Test decode with:
+5. Playing samples over a USB sounds card connected to a radio, note selection of sample rate:
    ```
-   cat t.raw | python3 int16tof32.py --zeropad | python3 radae_rx.py model19_check3/checkpoints/checkpoint_epoch_100.pth -v 2 --auxdata | ./build/src/lpcnet_demo -fargan-synthesis - - | aplay -f S16_LE -r 16000
+   aplay --device="plughw:CARD=Audio,DEV=0" -r 9600 -f S16_LE t1.int16
    ```
 
-1. Real time transmit:
-   ```
-   arecord --device "plughw:CARD=LX3000,DEV=0" -f S16_LE -c 1 -r 16000 | ./build/src/lpcnet_demo -features - - | python3 radae_tx.py model19_check3/checkpoints/checkpoint_epoch_100.pth --auxdata | python3 f32toint16.py --real --scale 8192 | aplay -f S16_LE --device "plughw:CARD=CODEC,DEV=0
-   ```
-   I keyed radio manually.  I recorded the transmission on a local SDR, then decoded with:
+6. Feeding samples from an off air wave file captured from a Rx to demod. Note `sc_xx` tools default to a centre freq of 1500Hz
    ```
-   sox ~/Downloads/sdr.ironstonerange.com_2024-08-19T22_03_13Z_7185.00_lsb.wav -t .s16 -r 8000 -c 1 - | python3 int16tof32.py --zeropad | python3 radae_rx.py model19_check3/checkpoints/checkpoint_epoch_100.pth -v 2 --auxdata | ./build/src/lpcnet_demo -fargan-synthesis - - | aplay -f S16_LE -r 16000
+   sox ~/Desktop/sc-ber-003.wav -t .s16 -r 9600 -c 1 - highpass 100 | python3 sc_rx.py --plots > z_hat.f32
    ```
-# C Port of Core Encoder and Decoder
 
-The model weights can be compiled in or loaded at init-time from a binary blob.  The actual model is hard coded in `rade_enc.c` and `rade_dec.c`, and can't be easily changed.
+# Automated Tests
 
-To compile-in the weights:
-1. Export weights:
-   ```
-   cd radae
-   python3 export_rade_weights.py model19_check3/checkpoints/checkpoint_epoch_100.pth src
-   ```
-1. We need to make some manual changes to the weight files to support changing input dimension at run time.  In `rade_enc_dat.c`, the first call to `linear_init()` should look like:
-   ```
-   int init_radeenc(RADEEnc *model, const WeightArray *arrays, int input_dim) {
-     if (linear_init(&model->enc_dense1, arrays, "enc_dense1_bias", NULL, NULL,"enc_dense1_weights_float", NULL, NULL, NULL, input_dim, 64)) return 1;
-   ```
-   e.g. the fixed input dimension (84 for `model19_check3`, 80 for earlier models without auxdata) should be changed to the `input_dim` variable. This allows us to enable/disable `auxdata` at init time, without changing the C code for the model.
-1. Also make manual changes to support `output_dim` in `rade_dec_dat.c`, `init_radedec()`.
-3. Build C code.
-4. Run ctests.
+The `cmake/ctest` framework is being used as a build and test framework. The command lines in `CMakeLists.txt` are a good source of examples, if you are interested in running the code in this repo.
 
-To export the compiled in weights to a binary blob:
+To run the cests:
 ```
 cd radae/build
-./src/write_rade_weights ../bin/model05.bin
+ctest
 ```
-These can then be loaded at init-time, see examples in `src/test_rand_enc.c` and `src/test_rand_dec.c`.
\ No newline at end of file
+To list tests `ctest -N`, to run just one test `ctest -R inference_model5`, to run in verbose mode `ctest -V -R inference_model5`.  You can change the path to `codec2-dev` on the `cmake` comma
+nd line:
+```
+cmake -DCODEC2_DEV=~/tmp/codec2-dev ..
+```
+
diff --git a/analog_bbfm.sh b/analog_bbfm.sh
deleted file mode 100755
index 2250a1f0..00000000
--- a/analog_bbfm.sh
+++ /dev/null
@@ -1,44 +0,0 @@
-#!/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
diff --git a/asr_test.sh b/asr_test.sh
new file mode 100755
index 00000000..6f2b831f
--- /dev/null
+++ b/asr_test.sh
@@ -0,0 +1,367 @@
+#!/usr/bin/env bash
+# asr_test.sh
+#
+# Automatic Speech Recognition (ASR) testing for the  Radio Autoencoder. This script
+# takes the samples from a clean dataset (e.g. Librispeech test-clean), and generates
+# a dataset with channel simulations (RADE, SSB etc) applied.
+
+CODEC2_DEV=${CODEC2_DEV:-${HOME}/codec2-dev}
+PATH=${PATH}:${CODEC2_DEV}/build_linux/src:${CODEC2_DEV}/build_linux/misc:${PWD}/build/src
+
+which ch >/dev/null || { printf "\n**** Can't find ch - check CODEC2_PATH **** \n\n"; exit 1; }
+
+source utils.sh
+
+function print_help {
+    echo
+    echo "Automated Speech Recognition (ASR) dataset processing for Radio Autoencoder testing"
+    echo
+    echo "  usage ./asr_test.sh ssb|rade|700D|fargan|4kHz|fm|bbfm [test option below]"
+    echo "  usage ./asr_test.sh ssb --No -30"
+    echo "  usage ./asr_test.sh rade --EbNodB 10"
+    echo "  usage ./asr_test.sh fm --RdBm -120"
+    echo "  usage ./asr_test.sh bbfm --RdBm -120"
+    echo
+    echo "    --EbNodB EbNodB           inference.py simulation noise level (experiment to get desired SNR)"
+    echo "    --No NodB                 ch channel simulation No value (experiment to get desired SNR)"
+    echo "    --RdBm RdBm               fm/bbfm simulation received power level"
+    echo "    -n numSamples             number of dataset samples to process (default all)"
+    echo "    --results resultsFile     name of results file (default results.txt)"
+    echo "    -d                        verbose debug information"
+    exit
+}
+
+n_samples=0
+No=-100
+EbNodB=100
+RdBm=-100
+setpoint_rms=2048
+setpoint_rms_fm=4096
+comp_gain=6
+results=asr_results.txt
+inference_args=""
+ch_args=""
+sil=0.5
+
+POSITIONAL=()
+while [[ $# -gt 0 ]]
+do
+key="$1"
+case $key in
+    --EbNodB)
+        EbNodB="$2"
+        shift
+        shift
+    ;;
+    --g_file)
+        g_file="$2"
+        if [ ! -f $2 ]; then
+            echo "can't find $2"
+            exit 1
+        fi
+        inference_args="${inference_args} --g_file ${2}"
+        cp ${2} fast_fading_samples.float	
+        ch_args="${ch_args} --fading_dir . --mpp --gain 0.5"
+        shift
+        shift
+    ;;
+    --h_file)
+        h_file="--h_file $2"
+        if [ ! -f $2 ]; then
+            echo "can't find $2"
+            exit 1
+        fi
+        shift
+        shift
+    ;;
+    --No)
+        No="$2"
+        shift
+        shift
+    ;;
+    --RdBm)
+        RdBm="$2"
+        echo $RdBm
+        shift
+        shift
+    ;;
+    --sil)
+        sil="$2"
+        shift
+        shift
+    ;;
+    --results)
+        results="$2"
+        shift
+        shift
+    ;;
+    -n)
+        n_samples="$2"
+        shift
+        shift
+    ;;
+    -d)
+        set -x;
+        shift
+    ;;
+    -h)
+        print_help	
+    ;;
+    *)
+    POSITIONAL+=("$1") # save it in an array for later
+    shift
+    ;;
+esac
+done
+set -- "${POSITIONAL[@]}" # restore positional parameters
+
+if [ $# -lt 1 ]; then
+    print_help
+fi
+mode=$1
+
+source=~/.cache/LibriSpeech/test-clean
+if [ ! -d $source ]; then
+  echo "cant find Librispeech source directory" $source
+  exit 1 
+fi
+# results must be written to a directory known by Librispeech package (can't be any name)
+dest=~/.cache/LibriSpeech/test-other
+rm -Rf $dest
+
+# cp translation files to new dataset directory
+function cp_translation_files {
+    pushd $source > /dev/null; trans=$(find . -name '*.txt'); popd  > /dev/null
+    for f in $trans
+    do
+        d=$(dirname $f)
+        mkdir -p ${dest}/${d}
+        cp ${source}/${f} ${dest}/${f}
+    done
+}
+
+function print_mean_text_file {
+    file_name=$1
+    python3 - <<END
+import numpy as np
+s=np.loadtxt("${file_name}")
+print(f"{np.mean(s):5.2f}", end='')
+END
+}
+
+# process audio files and place in new dataset directory
+function process {
+    pushd $source > /dev/null; flac=$(find . -name '*.flac'); popd > /dev/null
+    if [ $n_samples -ne 0 ]; then
+        flac=$(echo "$flac" | shuf --random-source=<(yes 42) | head -n $n_samples)
+    fi
+
+    n=$(echo "$flac" | wc -l)
+    printf "Processing %d samples in dataset\n" $n
+
+    in=in.raw
+    comp=comp.raw
+    ch_log=ch_log.txt
+    rade_log=rade_log.txt
+    snr_log=snr_log.txt
+    RdBm_log=RdBm_log.txt
+    asr_log=asr.txt
+    rm -f ${snr_log}
+    rm -f ${RdBm_log}
+    CNo_log=CNo_log.txt
+    rm -f ${CNo_log}
+    sox -n -r 16000 -c 1 /tmp/silence.wav trim 0.0 ${sil}
+
+    if [ $mode == "ssb" ] || [ $mode == "4kHz" ]; then
+        
+        fading_adv=0
+        for f in $flac
+        do
+            d=$(dirname $f)
+            mkdir -p ${dest}/${d}
+
+            if [ $mode == "ssb" ]; then
+                sox ${source}/${f} -t .s16 -r 8000 ${in}
+                # AGC and Hilbert compression
+                set_rms ${in} $setpoint_rms
+                analog_compressor  ${in} ${comp} ${comp_gain} 2>/dev/null
+                ch ${comp} - --No ${No} ${ch_args} --fading_adv ${fading_adv} 2>${ch_log} | sox -t .s16 -r 8000 -c 1 - -r 16000 ${dest}/${f}
+                grep "Fading file finished" $ch_log
+                if [ $? -eq 0 ]; then
+                    echo "Error - fading file too short after" $fading_adv " seconds"
+                    exit 1
+                fi
+                snr=$(cat $ch_log | grep "SNR3k" | tr -s ' ' | cut -d' ' -f3)
+                CNo=$(cat $ch_log | grep "SNR3k" | tr -s ' ' | cut -d' ' -f5)
+                echo $snr >> ${snr_log}
+                echo $CNo >> ${CNo_log}
+
+                # advance through fading simulation file
+                dur=$(sox --info -D ${source}/${f})
+                fading_adv=$(python3 -c "print(${fading_adv} + ${dur})")
+            else
+              # $mode == "4kHz" (4kHz bandwidth, representing ideal Fs=8kHz vocoder)
+              sox ${source}/${f} -r 8000 -t .s16 -c 1 - | sox -r 8000 -t .s16 -c 1 - -r 16000 ${dest}/${f}
+            fi
+        done
+        if [ $mode == "ssb" ]; then
+          SNR_mean=$(print_mean_text_file ${snr_log})
+          CNo_mean=$(print_mean_text_file ${CNo_log})
+        fi
+    fi
+    
+    if [ $mode == "700D" ]; then
+        
+        fading_adv=0
+        for f in $flac
+        do
+            d=$(dirname $f)
+            mkdir -p ${dest}/${d}
+
+            # silence either side of sample to allow time for acquisition and latency
+            sox /tmp/silence.wav /tmp/silence.wav ${source}/${f} /tmp/silence.wav -t .s16 -r 8000 ${in}
+
+            # trim start to remove acquisition noise
+            freedv_tx 700D ${in} - | \
+            ch - - --No ${No} ${ch_args} --fading_adv ${fading_adv} 2>${ch_log} | \
+            freedv_rx 700D - out.raw 2>/dev/null
+            cat out.raw | sox -t .s16 -r 8000 -c 1 - -r 16000 ${dest}/${f} trim 0.5
+            # error check
+            grep "Fading file finished" $ch_log
+            if [ $? -eq 0 ]; then
+                echo "Error - fading file too short after" $fading_adv " seconds"
+                exit 1
+            fi
+            snr=$(cat $ch_log | grep "SNR3k" | tr -s ' ' | cut -d' ' -f3)
+            CNo=$(cat $ch_log | grep "SNR3k" | tr -s ' ' | cut -d' ' -f5)
+            echo $snr >> ${snr_log}
+            echo $CNo >> ${CNo_log}
+
+            # advance through fading simulation file
+            dur=$(sox --info -D ${source}/${f})
+            fading_adv=$(python3 -c "print(${fading_adv} + ${dur})")
+
+        done
+        SNR_mean=$(print_mean_text_file ${snr_log})
+        CNo_mean=$(print_mean_text_file ${CNo_log})
+    fi
+
+    if [ $mode == "rade" ] || [ $mode == "fargan" ] || [ $mode == "bbfm" ] || [ $mode == "3200" ]; then
+        # find length of each file
+        duration_log=""
+        flac_full=""
+        pushd $source > /dev/null;
+        for f in $flac
+        do
+          duration_log+=$(sox --info -D ${f})
+          duration_log+=" "
+          flac_full+="${source}/${f} /tmp/silence.wav "
+        done
+        popd > /dev/null; 
+        
+        # cat samples into one long input file, insert 500ms at end of sample to allow for processing at output
+        sox $flac_full -t .s16 ${in}
+
+        # process all samples as one file to save time
+
+        if [ $mode == "rade" ]; then
+            ./inference.sh model19_check3/checkpoints/checkpoint_epoch_100.pth ${in} out.wav \
+            --rate_Fs --pilots --pilot_eq --eq_ls --cp 0.004 --bottleneck 3 --auxdata  --time_offset -16 \
+            --EbNodB $EbNodB ${inference_args} | tee ${rade_log}
+            grep "Multipath Doppler spread file too short" $rade_log
+            if [ $? -eq 0 ]; then
+                echo "Error - fading file too short"
+                exit 1
+            fi
+
+            SNR_mean=$(cat $rade_log | grep "Measured" | tr -s ' ' | cut -d' ' -f4)
+            CNo_mean=$(cat $rade_log | grep "Measured" | tr -s ' ' | cut -d' ' -f3)
+        fi
+
+        if [ $mode == "fargan" ]; then
+            lpcnet_demo -features ${in} - | lpcnet_demo -fargan-synthesis - - | sox -t .s16 -r 16000 -c 1 - out.wav
+        fi
+
+        # Codece 2 3200 as a control
+        if [ $mode == "3200" ]; then
+            sox -t.s16 -r 16000 -c 1 ${in} -t .s16 -r 8000 - | c2enc 3200 - - | c2dec 3200 - - | sox -t .s16 -r 8000 -c 1 - -r 16000 out.wav
+        fi
+
+        if [ $mode == "bbfm" ]; then
+            ./bbfm_inference.sh 250319_bbfm_lmr60/checkpoints/checkpoint_epoch_100.pth ${in} out.wav --RdBm $RdBm $h_file
+            if [ $? -ne 0 ]; then
+                exit 1
+            fi
+        fi
+
+        # extract individual output files
+        duration_array=( ${duration_log} )
+        i=0
+        st=0
+        for f in $flac
+        do
+          dur=${duration_array[i]}
+          dur=$(python3 -c "print($dur + ${sil})")
+          #printf "%4d %s %5.2f %5.2f\n" $i $f $st $dur
+          ((i++))
+          if [ $i -eq ${#duration_array[@]} ]; then
+            sox out.wav ${dest}/${f} trim $st
+          else
+            sox out.wav ${dest}/${f} trim $st $dur
+          fi
+          st=$(python3 -c "print($st + $dur)")
+        done
+    fi
+
+    # test mode that just copies files
+    if [ $mode == "clean" ]; then
+        for f in $flac
+        do
+          cp ${source}/${f} ${dest}/${f}
+        done        
+    fi
+
+    if [ $mode == "fm" ]; then
+        
+        fading_adv=0
+        for f in $flac
+        do
+            d=$(dirname $f)
+            mkdir -p ${dest}/${d}
+
+            sox ${source}/${f} -t .s16 -r 8000 ${in}
+            # AGC
+            set_rms ${in} $setpoint_rms_fm
+            sox -t .s16 -r 8000 -c 1 ${in} in.wav
+            ./bbfm_analog.sh in.wav out.wav --RdBm $RdBm $h_file --fading_adv $fading_adv
+            if [ $? -eq 1 ]; then
+                exit 1
+            fi
+            sox out.wav -r 16000 -c 1 ${dest}/${f}
+            echo $RdBm >> ${RdBm_log}
+            echo $CNo >> ${CNo_log}
+
+            # advance through fading simulation file
+            dur=$(sox --info -D ${source}/${f})
+            fading_adv=$(python3 -c "print(${fading_adv} + ${dur})")
+        done
+
+    fi
+
+    python3 asr_wer.py test-other -n $n_samples --model turbo | tee > $asr_log
+    wer=$(tail -n1 $asr_log | tr -s ' ' | cut -d' ' -f2)
+    if [ $mode == "ssb" ] || [ $mode == "rade" ] || [ $mode == "700D" ]; then
+      printf "%-6s %5.2f %5.2f %5.2f\n" $mode $SNR_mean $CNo_mean $wer | tee -a $results
+    fi
+    if [ $mode == "clean" ] || [ $mode == "fargan" ] || [ $mode == "4kHz" ] || [ $mode == "3200" ]; then
+      printf "%-6s %5.2f\n" $mode $wer | tee -a $results
+    fi
+    if [ $mode == "fm" ] || [ $mode == "bbfm" ]; then
+      printf "%-6s %5.2f %5.2f\n" $mode $RdBm $wer | tee -a $results
+    fi
+
+}
+
+cp_translation_files
+process
+
diff --git a/asr_test_top.sh b/asr_test_top.sh
new file mode 100755
index 00000000..a9e6a667
--- /dev/null
+++ b/asr_test_top.sh
@@ -0,0 +1,165 @@
+#!/usr/bin/env bash
+# asr_test_awgn.sh
+#
+# Top level ASR test script for AWGN and MPP channels
+
+n=500
+
+function print_help {
+    echo
+    echo "Generate curve data for ASR Radio Autoencoder testing"
+    echo
+    echo "  usage ./asr_test_top.sh rade|bbfm  [test option below]"
+    echo "  usage ./asr_test.sh rade"
+    echo "  usage ./asr_test.sh bbfm -n 100"
+    echo
+    echo "    -n numSamples             number of dataset samples to process (default 500)"
+    echo "    --results resultsFile     name of results file (default results.txt)"
+    echo "    -d                        verbose debug information"
+    exit
+}
+
+function ssb {
+    local results_file=$1
+    No_range=$2
+    for No in $No_range
+    do
+        ./asr_test.sh ssb --No $No -n $n --results ${results_file} $3
+    done
+    cat ${results_file} | grep ssb | sed -e "s/ssb//" > tmp.txt 
+    mv tmp.txt ${results_file}
+}
+
+function rade {
+    local results_file=$1
+    EbNodB_range=$2
+    for EbNodB in $EbNodB_range
+    do
+        ./asr_test.sh rade --EbNodB $EbNodB -n $n --results ${results_file} $3
+    done
+    cat ${results_file} | grep rade | sed -e "s/rade//" > tmp.txt
+    mv tmp.txt ${results_file}
+}
+
+function freedv_700D {
+    local results_file=$1
+    No_range=$2
+    for No in $No_range
+    do
+        ./asr_test.sh 700D --No $No -n $n --results ${results_file} $3
+    done
+    cat ${results_file} | grep 700D | sed -e "s/700D//" > tmp.txt 
+    mv tmp.txt ${results_file}
+}
+
+function fm {
+    local results_file=$1
+    RdBm_range=$2
+    for RdBm in $RdBm_range
+    do
+        ./asr_test.sh fm --RdBm $RdBm -n $n --results ${results_file} $3
+    done
+    cat ${results_file} | grep fm | sed -e "s/fm//" > tmp.txt
+    mv tmp.txt ${results_file}
+}
+
+function bbfm {
+    local results_file=$1
+    RdBm_range=$2
+    for RdBm in $RdBm_range
+    do
+        ./asr_test.sh bbfm --RdBm $RdBm -n $n --results ${results_file} $3
+    done
+    cat ${results_file} | grep bbfm | sed -e "s/bbfm//" > tmp.txt
+    mv tmp.txt ${results_file}
+}
+
+# Curves for 2024 HF RADE paper
+function rade_hf_top {
+    freedv_700D ${results_file}_awgn_700D.txt  "-100 -30 -26 -23 -20 -17 -15 -13"
+    freedv_700D ${results_file}_mpp_700D.txt   "-100 -39 -36 -33 -30 -27" "--g_file g_mpp.f32"
+    #freedv_700D ${results_file}_awgn_700D.txt  "-100 -38 -35 -32 -29 -26 -23 -20 -17"
+    #freedv_700D ${results_file}_mpp_700D.txt   "-100 -44 -39 -36 -33 -30 -27" "--g_file g_mpp.f32"
+    exit 0
+
+    # run the controls
+    controls_file=${results_file}_controls.txt
+    rm -f ${controls_file}
+    ./asr_test.sh clean -n $n --results ${controls_file}
+    ./asr_test.sh fargan -n $n --results ${controls_file}
+    ./asr_test.sh 4kHz -n $n --results ${controls_file}
+    ./asr_test.sh ssb -n $n --results ${controls_file}
+    ./asr_test.sh rade -n $n --results ${controls_file}
+    # strip off all but last column for Octave plotting
+    cat ${controls_file} | awk '{print $NF}' > ${results_file}_c.txt
+
+
+    ssb  ${results_file}_awgn_ssb.txt  "-100 -38 -35 -32 -29 -26 -23 -20 -17"
+    rade ${results_file}_awgn_rade.txt "100 15 10 5 2.5 0 -2.5"
+    ssb  ${results_file}_mpp_ssb.txt   "-100 -44 -39 -36 -33 -30 -27" "--g_file g_mpp.f32"
+    rade ${results_file}_mpp_rade.txt  "100 15 10 5 2.5 0" "--g_file g_mpp.f32"
+}
+
+# Curves for 2025 LMR paper
+function bbfm_top {
+    # run the single point controls
+    controls_file=${results_file}_controls.txt
+    rm -f ${controls_file}
+    ./asr_test.sh clean -n $n --results ${controls_file}
+    # use --sil 0 as FARGAN WER was high without it, could hear an artefact in the
+    # intersample silence section at the end of some output files.  OK to remove
+    # intersample silence as very little processing delay with FARGAN alone
+    ./asr_test.sh fargan -n $n --results ${controls_file} --sil 0
+    ./asr_test.sh 3200 -n $n --results ${controls_file}
+    # strip off all but last column for Octave plotting
+    cat ${controls_file} | awk '{print $NF}' > ${results_file}_c.txt
+
+    fm ${results_file}_awgn_fm.txt "-100 -110 -115 -118 -120 -122 -125"
+    bbfm ${results_file}_awgn_bbfm.txt "-100 -110 -120 -125 -126 -127"
+    fm ${results_file}_lmr60_fm.txt "-100 -105 -110 -113 -115 -117 -120" "--h_file h_lmr60_Fs_8000Hz.f32"
+    bbfm ${results_file}_lmr60_bbfm.txt "-100 -105 -110 -115 -120 -122 -123 -125 -126" "--h_file h_lmr60_Rs_2000Hz.f32"
+}
+
+POSITIONAL=()
+while [[ $# -gt 0 ]]
+do
+key="$1"
+case $key in
+    -n)
+        n="$2"
+        shift
+        shift
+    ;;
+    -d)
+        set -x;
+        shift
+    ;;
+    -h)
+        print_help	
+    ;;
+    *)
+    POSITIONAL+=("$1") # save it in an array for later
+    shift
+    ;;
+esac
+done
+set -- "${POSITIONAL[@]}" # restore positional parameters
+
+if [ $# -lt 1 ]; then
+    print_help
+fi
+
+mode=$1
+if [ $mode == 'rade' ]; then
+  results_file=241221_asr
+  rade_hf_top
+  exit 0
+fi
+if [ $mode == 'bbfm' ]; then
+  results_file=250807_asr
+  bbfm_top
+  exit 0
+fi
+
+print_help
+
diff --git a/asr_wer.py b/asr_wer.py
new file mode 100644
index 00000000..68b441a9
--- /dev/null
+++ b/asr_wer.py
@@ -0,0 +1,91 @@
+# coding: utf-8
+
+# derived from: https://github.com/openai/whisper/blob/main/notebooks/LibriSpeech.ipynb
+
+import os,argparse
+import numpy as np
+import torch
+import pandas as pd
+import whisper
+import torchaudio
+from tqdm.notebook import tqdm
+
+DEVICE = "cuda" if torch.cuda.is_available() else "cpu"
+
+
+class LibriSpeech(torch.utils.data.Dataset):
+    """
+    A simple class to wrap LibriSpeech and trim/pad the audio to 30 seconds.
+    It will drop the last few seconds of a very small portion of the utterances.
+    """
+    def __init__(self, n_mels, split="test-clean", device=DEVICE):
+        self.dataset = torchaudio.datasets.LIBRISPEECH(
+            root=os.path.expanduser("~/.cache"),
+            url=split,
+            download=True,
+        )
+        self.device = device
+        self.n_mels = n_mels
+        print(n_mels)
+
+    def __len__(self):
+        return len(self.dataset)
+
+    def __getitem__(self, item):
+        audio, sample_rate, text, _, _, _ = self.dataset[item]
+        assert sample_rate == 16000
+        audio = whisper.pad_or_trim(audio.flatten()).to(self.device)
+        mel = whisper.log_mel_spectrogram(audio,n_mels=self.n_mels)
+        
+        return (mel, text)
+
+
+parser = argparse.ArgumentParser()
+parser.add_argument('test_name', type=str, help='Librispeech dataset name (e.g. test-clean)')
+parser.add_argument('-n', type=str, help='Number of dataset entries to use (default all of them)')
+parser.add_argument('--model', default='base.en',type=str, help='Whisper model')
+args = parser.parse_args()
+
+model = whisper.load_model(args.model)
+print(
+    f"Model is {'multilingual' if model.is_multilingual else 'English-only'} "
+    f"and has {sum(np.prod(p.shape) for p in model.parameters()):,} parameters."
+)
+# predict without timestamps for short-form transcription
+options = whisper.DecodingOptions(language="en", without_timestamps=True)
+
+dataset = LibriSpeech(model.dims.n_mels, args.test_name)
+if args.n:
+    dataset = torch.utils.data.Subset(dataset,list(range(0,int(args.n))))
+print("dataset length:", dataset.__len__())
+loader = torch.utils.data.DataLoader(dataset, batch_size=16)
+
+
+hypotheses = []
+references = []
+
+for mels, texts in loader:
+    results = model.decode(mels, options)
+    hypotheses.extend([result.text for result in results])
+    references.extend(texts)
+
+data = pd.DataFrame(dict(hypothesis=hypotheses, reference=references))
+
+
+# # Calculating the word error rate
+# 
+# Now, we use our English normalizer implementation to standardize the transcription and calculate the WER.
+
+import jiwer
+from whisper.normalizers import EnglishTextNormalizer
+
+normalizer = EnglishTextNormalizer()
+
+data["hypothesis_clean"] = [normalizer(text) for text in data["hypothesis"]]
+data["reference_clean"] = [normalizer(text) for text in data["reference"]]
+print(data)
+
+wer = jiwer.wer(list(data["reference_clean"]), list(data["hypothesis_clean"]))
+
+print(f"WER: {wer * 100:.2f} %")
+
diff --git a/bbfm_analog.py b/bbfm_analog.py
new file mode 100644
index 00000000..ab14f43d
--- /dev/null
+++ b/bbfm_analog.py
@@ -0,0 +1,128 @@
+"""
+  Baseband simulation of analog FM using (11) from paper.
+
+  Fs=8000 Hz int16 speech samples on stdin, output samples on stdout.
+
+/* Copyright (c) 2025 David Rowe */
+   
+/*
+   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 sys,struct
+import argparse
+import numpy as np
+import math as m
+
+parser = argparse.ArgumentParser()
+
+parser.add_argument('--RdBm', type=float, default=-100, help='Receive level set point in dBm')
+parser.add_argument('--h_file', type=str, default="", help='Path to rate Fs fading channel magnitude samples, rate Fs time steps by Nc=1 carriers .f32 format')
+parser.add_argument('--fading_advance', type=float, default=0, help='Where to start sampling fading samples in seconds (default 0)')
+parser.add_argument('-v', action='store_true', help='Verbose debug info')
+args = parser.parse_args()
+RdBm = args.RdBm
+
+Am = 16384  # peak input int16 level, corresponds to max deviation f_d 
+A = 1       # normalised peak input level assumed in SNR expression
+x_bar = 0.5 # for sine wave with peak A=1
+
+k = 1.38E-23; T=274; NFdB = 5
+Fs = 8000
+fd_Hz = 2500
+fm_Hz = 3000
+beta = fd_Hz/fm_Hz
+Gfm = 10*m.log10(3*(beta**2)*x_bar/(1E3*k*T*fm_Hz)) - NFdB
+TdBm = 12 - Gfm
+
+print(f"Fs: {Fs:5.2f} Deviation: {fd_Hz} Hz  Max Modn freq: {fm_Hz} Hz Beta: {beta:3.2f}", file=sys.stderr)
+print(f"x_bar: {x_bar:5.2f} Gfm: {Gfm:5.2f} dB TdB: {TdBm:5.2f} dB  RdBm: {RdBm:5.2f}", file=sys.stderr)
+
+# user supplied rate Rs multipath model, sequence of H magnitude samples
+if len(args.h_file):
+  H = np.fromfile(args.h_file, dtype=np.float32)
+  fading_index = int(args.fading_advance*Fs)
+  print(f"fading_adv: {args.fading_advance:f} offset (samples): {fading_index:d}",  file=sys.stderr)
+
+# average noise and signal power
+n2_sum = 0.0
+x2_sum = 0.0
+n_sum = 0
+n_clipped = 0
+sigma = 0.0
+
+while True:
+    buffer = sys.stdin.buffer.read(struct.calcsize("h"))
+    if len(buffer) != struct.calcsize("h"):
+      break
+    x = np.frombuffer(buffer,np.int16).astype(np.float32)[0]
+
+    RdBm_dash = RdBm
+    if args.h_file:
+      if fading_index > len(H)-1:
+        print(f"ERROR; h_file too short for sample! Quitting", file=sys.stderr)
+        sys.exit(1)
+      RdBm_dash = 20*m.log10(H[fading_index]) + RdBm_dash
+      fading_index += 1
+
+    if RdBm_dash > TdBm:
+      SNRdB = RdBm_dash + Gfm
+    else:
+      SNRdB = 3*RdBm_dash + Gfm - 2*TdBm
+
+    # Work out sigma of the noise generator.  Eq (11) is the SNR in noise bandwidth f_m Hz.
+    # We want to simulate at Fs Hz. The noise generator spreads noise uniformly across Fs/2 Hz.
+    # After generation of the noise at any sample rate, the power in f_m Hz should be the same
+    # as given by (11).  This can be achieved by keeping the noise density constant across sample
+    # rate changes.  
+    SNR = 10**(SNRdB/10)   # Linear SNR in fm Hz
+    N_fm = x_bar/SNR       # Noise power in fm Hz
+    No   = N_fm/fm_Hz      # noise density, we want to preserve this with change in sample rate
+    N_Fs2 = No*Fs/2        # total noise power in Fs/2 Hz
+    sigma = N_Fs2 **0.5
+    sigma *= Am            # map A=1 to int16 peak
+    if args.v:
+        print(f"SNRdB: {SNRdB:5.2f}", file=sys.stderr)
+    n = sigma*np.random.randn()
+
+    n2_sum += n*n
+    x2_sum += x*x
+    n_sum += 1
+    x += n
+    if x > 32767.0:
+        x = 32767.0
+        n_clipped += 1
+    if x < -32776.0:
+        x = -32776.0
+        n_clipped += 1
+
+    x = np.float32(x).astype(np.int16)
+    sys.stdout.buffer.write(x) 
+
+SNRdB_ = 10.0*m.log10(x2_sum/n2_sum)
+x_bar_ = x2_sum/(n_sum*Am*Am)
+percent_clipped = 100 * n_clipped/n_sum
+print(f"SNRdB setpoint: {RdBm+Gfm:5.2f} SNRdB_ measured: {SNRdB_:5.2f} SNRdB_ - SNRdB: {SNRdB_-SNRdB:5.2f}", file=sys.stderr)
+print(f"x_bar measured: {x_bar_:5.2f} %clipped {percent_clipped:5.2f} sigma: {sigma:5.2f}", file=sys.stderr)
diff --git a/bbfm_analog.sh b/bbfm_analog.sh
new file mode 100755
index 00000000..d6561ebf
--- /dev/null
+++ b/bbfm_analog.sh
@@ -0,0 +1,51 @@
+#!/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:${CODEC2_DEV}/build_linux/misc
+
+which ch >/dev/null || { printf "\n**** Can't find ch - check CODEC2_PATH **** \n\n"; exit 1; }
+
+source utils.sh
+
+if [ $# -lt 2 ]; then
+    echo "usage (write output to file):"
+    echo "  ./analog_bbfm.sh in.wav out.wav [--RdBm XX --h_file YY]"
+    echo "usage (play output with aplay):"
+    echo "  ./analog_bbfm.sh in.wav - [--RdBm XX --h_file YY]"
+    exit 1
+fi
+
+if [ ! -f $1 ]; then
+    echo "can't find $1"
+    exit 1
+fi
+
+input_speech=$1
+output_speech=$2
+
+# eat first 2 args before passing rest to bbfm_analog.py in $@
+shift; shift;
+
+#tmp_comp=$(mktemp)
+tmp_comp=comp.raw
+tmp_out=$(mktemp)
+tmp_fm=$(mktemp)
+log=$(mktemp)
+
+# input wav -> BPF 300-3100 -> pre-emp -> Hilbert compressor -> de-emp -> linearised BBFM sim -> BPF 300-3100 -> output wav
+
+sox ${input_speech} -t .s16 -r 8000 -c 1 - | ch - - --ssbfilt 2 | pre - - | ch - - --clip 8192 --gain 2 --ssbfilt 2 | de - $tmp_comp
+cat $tmp_comp | python3 bbfm_analog.py $@ 2> >(tee $log >&2) | ch - - --ssbfilt 2 | sox -t .s16 -r 8000 -c 1 - -t .s16 -r 8000 ${tmp_out}
+grep "h_file too short" $log
+if [ $? -eq 0 ]; then
+    echo "Error - fading file too short"
+    exit 1
+fi
+if [ $output_speech == "-" ]; then
+    aplay ${tmp_out} -f S16_LE 2>/dev/null
+elif [ $output_speech != "/dev/null" ]; then
+    sox -t .s16 -r 8000 -c 1 ${tmp_out} -r 8000 ${output_speech}
+fi
diff --git a/bbfm_inference.py b/bbfm_inference.py
index 3b494cf0..f2272ed0 100644
--- a/bbfm_inference.py
+++ b/bbfm_inference.py
@@ -46,10 +46,10 @@
 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('--RdBm', type=float, default=-100, help='Receive level set point in dBm')
 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('--write_RdBm', type=str, default="", help='path to output file of RdBm 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()
 
@@ -66,8 +66,8 @@
 num_features = 20
 num_used_features = 20
 
-# load model from a checkpoint file
-model = BBFM(num_features, latent_dim, args.CNRdB)
+model = BBFM(num_features, latent_dim, args.RdBm)
+# load model weights from a checkpoint file
 checkpoint = torch.load(args.model_name, map_location='cpu', weights_only=True)
 model.load_state_dict(checkpoint['state_dict'], strict=False)
 checkpoint['state_dict'] = model.state_dict()
@@ -81,10 +81,10 @@
 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)
+# default AWGN channel (H=1)
 H = torch.ones((1,num_timesteps_at_rate_Rs,Nc))
 
 # user supplied rate Rs multipath model, sequence of H magnitude samples
@@ -93,7 +93,7 @@
    print(H.shape, num_timesteps_at_rate_Rs)
    if H.shape[1] < num_timesteps_at_rate_Rs:
       print("Multipath H file too short")
-      quit()
+      exit(1)
    H = H[:,:num_timesteps_at_rate_Rs,:]
    H = torch.tensor(H)
 
@@ -130,13 +130,13 @@
       else:
          print("FAIL")
 
-   # write output symbols (latent vectors)
+   # optionally 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)
+   # optionally write RdBm after fading
+   if len(args.write_RdBm):
+      RdBm = output["RdBm"].cpu().detach().numpy().flatten().astype('float32')
+      RdBm.tofile(args.write_RdBm)
       
diff --git a/bbfm_inference.sh b/bbfm_inference.sh
index 04bd3ae1..809b8cd1 100755
--- a/bbfm_inference.sh
+++ b/bbfm_inference.sh
@@ -32,7 +32,7 @@ features_out=features_out.f32
 shift; shift; shift
 
 lpcnet_demo -features ${input_speech} ${features_in}
-python3 ./bbfm_inference.py ${model} ${features_in} ${features_out} "$@"
+python3 ./bbfm_inference.py ${model} ${features_in} ${features_out} "$@" 2> >(tee $bbfm_log >&2)
 if [ $? -ne 0 ]; then
   exit 1
 fi
diff --git a/bbfm_rx.py b/bbfm_rx.py
index c55db0d8..71e20dec 100644
--- a/bbfm_rx.py
+++ b/bbfm_rx.py
@@ -44,6 +44,7 @@
 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('--stateful', action='store_true', help='Use stateful decoder')
 args = parser.parse_args()
 
 # set visible devices
@@ -60,17 +61,17 @@
 num_used_features = 20
 
 # load model from a checkpoint file
-model = BBFM(num_features, latent_dim, CNRdB=100)
+model = BBFM(num_features, latent_dim, RdBm=-100, stateful_decoder=args.stateful)
 checkpoint = torch.load(args.model_name, map_location='cpu', weights_only=True)
 model.load_state_dict(checkpoint['state_dict'], strict=False)
 checkpoint['state_dict'] = model.state_dict()
+if args.stateful:
+   model.core_decoder_statefull_load_state_dict()
 
 # dataloader
 z_hat = np.reshape(np.fromfile(args.z_hat, dtype=np.float32), (1, -1, args.latent_dim))
-nb_frames_rounded = model.num_10ms_times_steps_rounded_to_modem_frames(z_hat.shape[1])
-z_hat = z_hat[:,:nb_frames_rounded,:]
 z_hat = torch.tensor(z_hat)
-print(f"Processing: {nb_frames_rounded} modem frames")
+print(f"Processing: {z_hat.shape[1]} modem frames")
 
 if __name__ == '__main__':
 
diff --git a/bbfm_rx.sh b/bbfm_rx.sh
index 94bb7a26..a0202d3a 100755
--- a/bbfm_rx.sh
+++ b/bbfm_rx.sh
@@ -1,6 +1,6 @@
 #!/bin/bash
 #
-# The usual wrapper around rx_bbfm.py
+# The usual wrapper around bbfm_rx.py
 
 OPUS=build/src
 PATH=${PATH}:${OPUS}
diff --git a/bbfm_rx_stream.py b/bbfm_rx_stream.py
new file mode 100644
index 00000000..f4dc29d1
--- /dev/null
+++ b/bbfm_rx_stream.py
@@ -0,0 +1,75 @@
+"""
+  BBFM stand alone streaming Rx:
+    float z_hat[80] on stdin
+    float features[36] on stdout
+
+/* Copyright (c) 2024 David Rowe */
+   
+/*
+   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,sys, struct, argparse
+import numpy as np
+import torch
+from radae import BBFM
+
+parser = argparse.ArgumentParser()
+
+parser.add_argument('model_name', type=str, help='path to model in .pth format')
+parser.add_argument('--latent-dim', type=int, help="number of symbols produces by encoder, default: 80", default=80)
+args = parser.parse_args()
+
+# make sure we don't use a GPU
+os.environ['CUDA_VISIBLE_DEVICES'] = ""
+device = 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, RdBm=-100, stateful_decoder=True)
+checkpoint = torch.load(args.model_name, map_location='cpu', weights_only=True)
+model.load_state_dict(checkpoint['state_dict'], strict=False)
+checkpoint['state_dict'] = model.state_dict()
+model.core_decoder_statefull_load_state_dict()
+
+if __name__ == '__main__':
+
+    while True:
+        buffer = sys.stdin.buffer.read(args.latent_dim*struct.calcsize("f"))
+        if len(buffer) != args.latent_dim*struct.calcsize("f"):
+            break
+        z_hat = np.reshape(np.frombuffer(buffer,np.float32),(1,1,args.latent_dim))
+        z_hat = torch.tensor(z_hat)
+        z_hat = z_hat.to(device)
+        features_hat = model.receiver(z_hat)
+        features_hat = torch.cat([features_hat, torch.zeros_like(features_hat)[:,:,:16]], dim=-1)
+        features_hat = features_hat.cpu().detach().numpy().flatten().astype('float32')
+        sys.stdout.buffer.write(features_hat)
diff --git a/compare_models_bbfm.sh b/compare_models_bbfm.sh
new file mode 100755
index 00000000..1cb910bf
--- /dev/null
+++ b/compare_models_bbfm.sh
@@ -0,0 +1,49 @@
+#!/bin/bash -x
+#
+# Compare BBFM models by plotting loss v Eq/No (and C/NO and SNR) curves generated by train_bbfm.py
+
+# Run through training dataset on each trained model to build loss versus Eq/No curve
+function run_model() {
+    model=$1
+    dim=$2
+    epoch=$3
+    chan=$4
+    shift
+    shift
+    shift
+    shift
+    python3 ./train_bbfm.py --cuda-visible-devices 0 --sequence-length 400 --batch-size 512 --latent-dim ${dim} \
+        --epochs 1 --lr 0.003 --lr-decay-factor 0.0001 \
+        ~/Downloads/tts_speech_16k_speexdsp.f32 tmp \
+        --range_RdBm --plot_R ${model}_${chan} --initial-checkpoint ${model}/checkpoints/checkpoint_epoch_${epoch}.pth $@
+}
+
+#run_model model_bbfm_01 80 100 awgn
+#run_model 250319_bbfm 80 100 awgn
+#run_model 250319_bbfm_lmr60 80 100 awgn
+#run_model model_bbfm_01 80 100 lmr60 --h_file h_lmr60_train.f32
+#run_model 250319_bbfm 80 100 lmr60 --h_file h_lmr60_train.f32
+#run_model 250319_bbfm_lmr60 80 100 lmr60 --h_file h_lmr60_train.f32
+#run_model 250319_bbfm_lmr60 80 100 lmr30 --h_file h_lmr30_train.f32
+#run_model 250319_bbfm_lmr60 80 100 lmr120 --h_file h_lmr120_train.f32
+
+plot="250402"
+
+if [ $plot == "250402" ]; then
+  model_list='250319_bbfm_lmr60_awgn 250319_bbfm_lmr60_lmr60 250319_bbfm_lmr60_lmr30 250319_bbfm_lmr60_lmr120'
+  declare -a model_legend=("AWGN" "lmr60" "lmr30" "lmr120")
+fi
+
+# Generate the plots in PNG and EPS form, file names have suffix of ${plot}
+
+loss_R=""
+i=0;
+for model in $model_list
+  do
+    loss_R="${loss_R},'${model}_loss_RdBm.txt','${model_legend[i]}'"
+    ((i++))
+  done
+echo "radae_plots; loss_RdBm_plot('${plot}_loss_RdBm_models',''${loss_R}); quit" | octave-cli -qf  # PNG
+echo "radae_plots; loss_RdBm_plot('','${plot}_loss_R_models'${loss_R}); quit" | octave-cli -qf     # EPS
+
+
diff --git a/multipath_samples.m b/multipath_samples.m
index 208f7113..21b0b571 100644
--- a/multipath_samples.m
+++ b/multipath_samples.m
@@ -14,9 +14,14 @@ function multipath_samples(ch, Fs, Rs, Nc, Nseconds, H_fn, G_fn="")
         dopplerSpreadHz = 1.0; path_delay_s = 2E-3;
     elseif strcmp(ch,"mpd")
         dopplerSpreadHz = 2.0; path_delay_s = 4E-3;
-    elseif strcmp(ch,"lmr60")
-        % 60 km/hr, 450 MHz
-        fd = 450E6*(60*1E3/3600/3E8)
+    elseif strncmp(ch,"lmr",3)
+        % extract speed in km/hr from last few digits
+        v_km_hr = str2num(substr(ch,4,length(ch)-3))
+        % convert to m/s
+        v =  v_km_hr*1E3/3600
+        % freq in Hz
+        f = 450E6
+        fd = f*(v/3E8)
         dopplerSpreadHz = 2*fd;
         path_delay_s = 200E-6
     else
@@ -55,11 +60,15 @@ function multipath_samples(ch, Fs, Rs, Nc, Nseconds, H_fn, G_fn="")
         p2 = H(n+1,1).^2;
         if p1 < P && p2 > P
           LC++;
-          LC_log = [LC_log n];
+          if n < Nsecplot*Rs
+            % prevent repeated memory allocations for large samples, just enough to plot
+            LC_log = [LC_log n];
+          end
         end
       end
       LCR_meas = LC/Nseconds
-      subplot(211); hold on; stem(LC_log,sqrt(P)*ones(length(LC_log))); hold off; axis([0 Nsecplot*Rs 0 3]);
+      # Plot zero crossings on top of |H|
+      subplot(211); hold on; stem(LC_log,sqrt(P)*ones(length(LC_log)),'r'); hold off; axis([0 Nsecplot*Rs 0 3]);
     end
     printf("H file size is Nseconds*Rs*Nc*(4 bytes/sample) = %d*%d*%d*4 = %d bytes\n", Nseconds,Rs,Nc,Nseconds*Rs*Nc*4)
     f=fopen(H_fn,"wb");
diff --git a/radae/bbfm.py b/radae/bbfm.py
index efaea8a0..bf6d906f 100644
--- a/radae/bbfm.py
+++ b/radae/bbfm.py
@@ -43,20 +43,22 @@ class BBFM(nn.Module):
     def __init__(self,
                  feature_dim,
                  latent_dim,
-                 CNRdB,
-                 fd_Hz=5000,
-                 fm_Hz=3000,
-                 stateful_decoder = False
+                 RdBm,
+                 fd_Hz=1800,
+                 fm_Hz=2880,
+                 stateful_decoder = False,
+                 range_RdBm = False
                 ):
 
         super(BBFM, self).__init__()
 
         self.feature_dim = feature_dim
         self.latent_dim  = latent_dim
-        self.CNRdB = CNRdB
+        self.RdBm = RdBm
         self.fd_Hz = fd_Hz
         self.fm_Hz = fm_Hz
         self.stateful_decoder = stateful_decoder
+        self.range_RdBm = range_RdBm
 
         # TODO: nn.DataParallel() shouldn't be needed
         self.core_encoder =  nn.DataParallel(radae_base.CoreEncoder(feature_dim, latent_dim, bottleneck=1))
@@ -75,11 +77,14 @@ def __init__(self,
         self.Rz = 1/self.Tz
         self.Rb =  latent_dim/self.Tz                  # payload data BPSK symbol rate (symbols/s or Hz)
 
-        self.beta = self.fd_Hz/self.fm_Hz              # deviation
-        self.BWfm = 2*(self.fd_Hz + self.fm_Hz)        # BW estimate using Carsons rule
-        self.Gfm = 10*m.log10(3*(self.beta**2)*(self.beta+1))
-
-        print(f"Rb: {self.Rb:5.2f} Deviation: {self.fd_Hz}Hz Max Modn freq: {self.fm_Hz}Hz Beta: {self.beta:3.2f}", file=sys.stderr)
+        x_bar = 1                                      # average power of modulating symbols wrt peak deviation          
+        k = 1.38E-23; T=274; NFdB = 5
+        self.beta = self.fd_Hz/self.fm_Hz
+        self.Gfm = 10*m.log10(3*(self.beta**2)*x_bar/(1E3*k*T*self.fm_Hz)) - NFdB
+        self.TdBm = 12 - self.Gfm
+ 
+        print(f"Rb: {self.Rb:5.2f} Deviation: {self.fd_Hz} Hz  Max Modn freq: {self.fm_Hz} Hz Beta: {self.beta:3.2f}", file=sys.stderr)
+        print(f"x_bar: {x_bar:5.2f} Gfm: {self.Gfm:5.2f} dB TdB: {self.TdBm:5.2f} dB  RdBm: {self.RdBm:5.2f}", file=sys.stderr)
 
     # Stateful decoder wasn't present during training, so we need to load weights from existing decoder
     def core_decoder_statefull_load_state_dict(self):
@@ -134,7 +139,6 @@ def key_transformation(old_key):
     # stand alone receiver, takes received symbols z and returns features f
     def receiver(self, z_hat):
         if self.stateful_decoder:
-            print("stateful!", file=sys.stderr)
             features_hat = torch.empty(1,0,self.feature_dim)
             for i in range(z_hat.shape[1]):
                 features_hat = torch.cat([features_hat, self.core_decoder_statefull(z_hat[:,i:i+1,:])],dim=1)
@@ -171,21 +175,30 @@ def forward(self, features, H):
         z_shape = z.shape
         z_hat = torch.reshape(z,(num_batches,num_timesteps_at_rate_Rs,1))
         
-        # determine FM demod SNR using piecewise approximation implemented with relus to be torch-friendly
-        # note SNR is a vector, 1 sample for symbol as SNR evolves with H
-        CNRdB = 20*torch.log10(H) + self.CNRdB
-        print(H.shape,CNRdB.shape)
-        SNRdB_relu = torch.relu(CNRdB-12) + 12 + self.Gfm
-        SNRdB_relu += -torch.relu(-(CNRdB-12))*(1 + self.Gfm/3)
-        SNR = 10**(SNRdB_relu/10)
-
+        # training time option to use a range of R, one value per batch
+        if self.range_RdBm:
+            RdBm = self.RdBm - 20*torch.rand(num_batches,1,1,device=features.device)
+        else:        
+            RdBm = self.RdBm*torch.ones(num_batches,1,1,device=features.device)        
+        RdBm_ = RdBm.reshape(num_batches)
+
+        # determine FM demod SNR using piecewise approximation expressed as sum of
+        # heaviside step functions for efficient implementation during training.
+        # (avoids a for loop with per-sample if-then-else)
+        # Note SNR is a vector, 1 sample per symbol as SNR evolves with H
+        values = torch.zeros(1, device=H.device) 
+        RdBm = 20*torch.log10(H) + RdBm
+        SNRdB = (RdBm+self.Gfm)*torch.heaviside(RdBm-self.TdBm, values) \
+              + (3*RdBm+self.Gfm-2*self.TdBm)*torch.heaviside(-RdBm+self.TdBm, values)
+        SNR = 10**(SNRdB/10)
         # note sigma is a vector, noise power evolves across each symbol with H
+        # TODO: should be A/sqrt(SNR); correct small error due to noise bandwidth f_m used 
+        # for SNR expression v R_s for noise generation below
         sigma = 1/(SNR**0.5)
         n = sigma*torch.randn_like(z_hat)
         z_hat = torch.clamp(z_hat + n, min=-1.0,max=1.0)
                         
         z_hat = torch.reshape(z_hat,z_shape)
-        #print(z.shape, z_hat.shape)
         
         features_hat = self.core_decoder(z_hat)
         
@@ -193,5 +206,7 @@ def forward(self, features, H):
             "features_hat" : features_hat,
             "z_hat"  : z_hat,
             "sigma"  : sigma,
-            "CNRdB"  : CNRdB
-        }
+            "SNRdB"  : SNRdB,
+            "RdBm_"   : RdBm_, # setpoint R for entire sequence
+            "RdBm"   : RdBm    # per sample R after fading added
+       }
diff --git a/radae_plots.m b/radae_plots.m
index 39258622..239e7056 100644
--- a/radae_plots.m
+++ b/radae_plots.m
@@ -72,20 +72,26 @@ function do_plots(z_fn='l.f32',rx_fn='', png_fn='', epslatex='')
     end
 endfunction
 
-function do_plots_bbfm(z1_fn, z2_fn="", png_fn='')
+function do_plots_bbfm(z1_fn, z2_fn='', png_fn='', epslatex='')
+    if length(epslatex)
+      [textfontsize linewidth] = set_fonts(20);
+    end
     z1=load_f32(z1_fn,1);
     figure(1); clf; 
-    stem(z1(1:40),'g');
+    stem(z1(1:80),'g');
+    axis([0 80 -1.2 1.2]);
     if length(z2_fn) 
       z2=load_f32(z2_fn,1);
       hold on;
       stem(z2(1:40),'r');
       hold off;
     end
-    title('Rx Symbols');
     if length(png_fn)
       print("-dpng",sprintf("%s.png",png_fn));
     end
+    if length(epslatex)
+      print_eps_restore(sprintf("%s.eps",epslatex),"-S300,200",textfontsize,linewidth);
+    end
 endfunction
 
 
@@ -176,6 +182,29 @@ function loss_CNo_plot(png_fn, Rs, B, varargin)
     end
 endfunction
 
+% Plots loss v R curves from text files dumped by train_bbfm.py, pass in pairs from *loss_RdBm.txt,legend
+function loss_RdBm_plot(png_fn, epslatex, varargin)
+    if length(epslatex)
+        [textfontsize linewidth] = set_fonts(20);
+    end
+    figure(1); clf; hold on;
+    i = 1;
+    while i <= length(varargin)
+        fn = varargin{i};
+        data = load(fn);
+        i++; leg = varargin{i}; leg = strrep (leg, "_", "-")
+        plot(data(:,1),data(:,2),sprintf("+-;%s;",leg))
+        i++;
+    end
+    hold off; grid; xlabel('R (dBm)'); ylabel('loss'); legend('boxoff');
+    if length(png_fn)
+        print("-dpng",png_fn);
+    end
+    if length(epslatex)
+        print_eps_restore(epslatex,"-S300,200",textfontsize,linewidth);
+    end
+endfunction
+
 % usage:
 %   radae_plots; ofdm_sync_plots("","ofdm_sync.txt","go-;genie;","ofdm_sync_pilot_eq.txt","r+-;mean6;","ofdm_sync_pilot_eq_f2.txt","bx-;mean6 2 Hz;","ofdm_sync_pilot_eq_g0.1.txt","gx-;mean6 gain 0.1;","ofdm_sync_pilot_eq_ls.txt","ro-;LS;","ofdm_sync_pilot_eq_ls_f2.txt","bo-;LS 2 Hz;")
 
@@ -383,42 +412,66 @@ function test_rayleigh(epslatex="")
   y(find(x<0)) = 0;
 end
 
-% Plot SNR v CNR for FM demod model
-function plot_SNR_CNR(epslatex="")
+function y = heaviside(x)
+ y = x>0;
+end
+
+% Plot SNR v R for FM demod model
+function bbfm_plot_SNR_R(epslatex="")
     if length(epslatex)
-        [textfontsize linewidth] = set_fonts();
+        [textfontsize linewidth] = set_fonts(20);
     end
-    figure(1); clf; hold on;
-    fd=5000; fm=3000; 
-    beta= fd/fm;
-    Gfm=10*log10(3*(beta^2)*(beta+1))
-    BWfm = 2*(fd+fm);
+
+    fd=2500; fm=3000; A = 1; k=1.38E-23; T=274; NFdB=5;
+    beta = fd/fm;
+    x_bar = A^2/2;
+    Gfm=10*log10(3*(beta^2)*x_bar/(1E3*k*T*fm)) - NFdB;
+    TdBm = 12 - Gfm;
+    printf("fd: %6.0f fm: %6.0f Beta: %f A: %5.2f x_bar: %5.2f Gfm: %5.2f dB TdBm: %5.2f\n", fd, fm, beta, A, x_bar, Gfm, TdBm);
 
     % vanilla implementation of curve
-    CNRdB=0:20;
-    for i=1:length(CNRdB)
-      if CNRdB(i) >= 12
-        SNRdB(i) = CNRdB(i) + Gfm;
+    RdBm=-130:-105;
+    for i=1:length(RdBm)
+      if RdBm(i) >= TdBm
+        SNRdB(i) = RdBm(i) + Gfm;
       else
-        SNRdB(i) = (1+Gfm/3)*CNRdB(i) - 3*Gfm;
+        SNRdB(i) = 3*RdBm(i) + Gfm - 2*TdBm;
       end
     end
 
-    % implementation using relus (suitable for PyTorch)
-    SNRdB_relu = relu(CNRdB-12) + 12 + Gfm;
-    SNRdB_relu += -relu(-(CNRdB-12))*(1+Gfm/3);
-
-    plot(CNRdB,SNRdB,'g;FM;'); 
-    plot(CNRdB,SNRdB_relu,'r+;FM relu;'); 
-    SSBdB = CNRdB + 10*log10(BWfm) - 10*log10(fm);
-    plot(CNRdB,SSBdB,'b;SSB;'); 
-    axis([min(CNRdB) max(CNRdB) 10 30]);
-    hold off; grid('minor'); xlabel('CNR (dB)'); ylabel('SNR (dB)'); legend('boxoff'); legend('location','northwest');
+    % implementation using common ML toolkit non-linearity rather than if/then for efficiency in training
+    SNRdB_heaviside = (RdBm+Gfm).*heaviside(RdBm-TdBm) + (3*RdBm+Gfm-2*TdBm).*heaviside(-RdBm+TdBm);
+  
+    figure(1); clf; hold on;
+    plot(RdBm,SNRdB,'g+-');
+    if length(epslatex) == 0
+      hold on; plot(RdBm,SNRdB_heaviside,'bx'); hold off;
+    end
+    grid('minor'); xlabel('R (dBm)'); ylabel('SNR (dB)'); legend('off');
     if length(epslatex)
-        print_eps_restore(epslatex,"-S300,300",textfontsize,linewidth);
+        print_eps_restore(epslatex,"-S300,200",textfontsize,linewidth);
     end
 endfunction
 
+% test expression derived from Carslon (17)
+function bbfm_carlson()
+    fd=2500; fm=3000; 
+    beta = fd/fm;
+    Sx = 0.5;
+    Bt = 9E3;
+    CNR_dB = 0:20; CNR = 10.^(CNR_dB/10);
+    SNR1 = 3*(beta^2)*Sx*CNR;
+    num = 3*(beta^2)*Sx*CNR;
+    denom = (1+(12*beta/pi)*CNR.*exp(-fm*CNR/Bt));
+    SNR3 = num./denom;
+    SNR1_dB = 10*log10(SNR1);
+    SNR3_dB = 10*log10(SNR3);
+    figure(1); clf; hold on;
+    plot(CNR_dB,SNR1_dB,'b;Eq (1);');
+    plot(CNR_dB,SNR3_dB,'r;Eq (4);');
+    hold off; grid;
+endfunction
+
 % test handling of single sample per symbol phase jumps
 function test_phase_est
   theta = 0:0.01:2*pi;
@@ -483,3 +536,43 @@ function plot_sample_spec(wav_fn,png_spec_fn="")
     print("-dpng",png_spec_fn,"-S800,600");
   end
 end
+
+function plot_wer_bbfm(prefix_fn, png_fn="", epslatex="")
+  fm_awgn_fn = sprintf("%s_asr_awgn_fm.txt",prefix_fn);
+  rade_awgn_fn = sprintf("%s_asr_awgn_bbfm.txt",prefix_fn);
+  fm_lmr60_fn = sprintf("%s_asr_lmr60_fm.txt",prefix_fn);
+  rade_lmr60_fn = sprintf("%s_asr_lmr60_bbfm.txt",prefix_fn);
+  controls_fn = sprintf("%s_asr_c.txt",prefix_fn);
+
+  fm_awgn = load(fm_awgn_fn);
+  rade_awgn = load(rade_awgn_fn);
+  fm_lmr60 = load(fm_lmr60_fn);
+  rade_lmr60 = load(rade_lmr60_fn);
+  c = load(controls_fn);
+
+  if length(epslatex)
+    [textfontsize linewidth] = set_fonts(30);
+  end
+
+  # WER v RdBm plot
+  figure(1); clf;
+  plot(fm_awgn(:,1),fm_awgn(:,2),'b+-;FM AWGN;');
+  hold on;
+  plot(rade_awgn(:,1),rade_awgn(:,2),'g+-;RADE AWGN;');
+  plot(fm_lmr60(:,1),fm_lmr60(:,2),'bo--;FM LMR60;');
+  plot(rade_lmr60(:,1),rade_lmr60(:,2),'go--;RADE LMR60;');
+  xmax=-100; xmin=-130; 
+  plot([xmin xmax],[c(3) c(3)],'r-;Codec 2 3200;')
+  plot([xmin xmax],[c(2) c(2)],'m-;FARGAN;')
+  plot([xmin xmax],[c(1) c(1)],'c-;clean;')
+  hold off;
+  axis([xmin,xmax,0,40]); grid; ylabel('WER \%'); xlabel("R (dBm)");
+  legend('boxoff');
+
+  if length(png_fn)
+    print("-dpng",png_fn,"-S800,600");
+  end
+  if length(epslatex)
+    print_eps_restore(epslatex,"-S250,250",textfontsize,linewidth);
+  end  
+end
\ No newline at end of file
diff --git a/train_bbfm.py b/train_bbfm.py
index d1b794c9..a697e862 100644
--- a/train_bbfm.py
+++ b/train_bbfm.py
@@ -49,7 +49,8 @@
 parser.add_argument('output', type=str, help='path to output folder')
 parser.add_argument('--cuda-visible-devices', type=str, help="comma separates list of cuda visible device indices, default: ''", default="")
 parser.add_argument('--latent-dim', type=int, help="number of symbols produced by encoder, default: 80", default=80)
-parser.add_argument('--CNRdB', type=float, default=0, help='FM demod input CNR in dB')
+parser.add_argument('--RdBm', type=float, default=-100.0, help='Receive level set point in dBm (default -120)')
+parser.add_argument('--range_RdBm',  action='store_true', help='Sweep receive level during training')
 parser.add_argument('--h_file', type=str, default="", help='path to rate Rs multipath file, rate Rs time steps by 1 carriers .f32 format')
 
 training_group = parser.add_argument_group(title="training parameters")
@@ -61,8 +62,7 @@
 
 training_group.add_argument('--initial-checkpoint', type=str, help='initial checkpoint to start training from, default: None', default=None)
 training_group.add_argument('--plot_loss', action='store_true', help='plot loss versus epoch as we train')
-training_group.add_argument('--plot_EqNo', type=str, default="", help='plot loss versus Eq/No for final epoch')
-training_group.add_argument('--auxdata', action='store_true', help='inject auxillary data symbol')
+training_group.add_argument('--plot_R', type=str, default="", help='plot loss versus RdBm for final epoch, arg is suffix')
 
 args = parser.parse_args()
 
@@ -103,15 +103,13 @@
 latent_dim = args.latent_dim
 
 num_features = 20
-if args.auxdata:
-    num_features += 1
 
 # training data
 feature_file = args.features
 
 # model
-checkpoint['model_args'] = (num_features, latent_dim, args.CNRdB)
-model = BBFM(num_features, latent_dim, args.CNRdB)
+checkpoint['model_args'] = (num_features, latent_dim, args.RdBm)
+model = BBFM(num_features, latent_dim, args.RdBm, range_RdBm=args.range_RdBm)
 
 if type(args.initial_checkpoint) != type(None):
     print(f"Loading from checkpoint: {args.initial_checkpoint}")
@@ -142,11 +140,72 @@
     # push model to device
     model.to(device)
 
+    # -----------------------------------------------------------------------------------------------
+    # run through dataset once with current model but training disabled, to gather loss v R stats
+    # -----------------------------------------------------------------------------------------------
+    if len(args.plot_R):
+        # TODO: move this to a function
+        print("Measuring loss v RdBm over training set with training disabled")
+        model.eval()
+        R_loss = np.zeros((dataloader.__len__()*batch_size,2))
+        running_total_loss      = 0
+        previous_total_loss     = 0
+        current_loss            = 0.
+
+        with torch.no_grad():
+            with tqdm.tqdm(dataloader, unit='batch') as tepoch:
+                for i, (features,H,G) in enumerate(tepoch):
+                    features = features.to(device)
+                    H = H.to(device)
+                    output = model(features,H)
+                    loss_by_batch = distortion_loss(features[..., :20], output["features_hat"][..., :20])
+                    total_loss = torch.mean(loss_by_batch)
+                    
+                    # collect running stats, R and loss for each sequence in batch
+                    R_loss[i*batch_size:(i+1)*batch_size,0] = output["RdBm_"].cpu().detach().numpy()
+                    R_loss[i*batch_size:(i+1)*batch_size,1] = loss_by_batch.cpu().detach().numpy()                       
+
+                    running_total_loss += float(total_loss.detach().cpu())
+                    
+                    if (i + 1) % log_interval == 0:
+                        current_loss = (running_total_loss - previous_total_loss) / log_interval
+                        tepoch.set_postfix(
+                            current_loss=current_loss,
+                            total_loss=running_total_loss / (i + 1),
+                        )
+                        previous_total_loss = running_total_loss
+
+        # Plot loss against EqNodB for final epoch, using log of loss and Eq/No for
+        # each sequence.  We group losses into 1dB R bins, kind of like a histogram.
+        R_min = int(np.ceil(np.min(R_loss[:,0])))
+        R_max = int(np.ceil(np.max(R_loss[:,0])))
+        R_mean_loss = np.zeros((R_max-R_min,2))
+        # group the losses from training into 1dB wide bins, and find mean for that bin
+        r = np.arange(R_min,R_max)
+        for i in np.arange(len(r)):
+            R = r[i]
+            x = np.where(np.abs(R_loss[:,0] - R) < 0.5)
+            R_mean_loss[i,0] = R
+            R_mean_loss[i,1] = np.mean(R_loss[x,1])
+        plt.figure(2)
+        plt.plot(R_mean_loss[:,0],R_mean_loss[:,1],'b+-')
+        plt.grid()
+        plt.xlabel('R (dBm)')
+        plt.ylabel('Loss')
+        plt.show(block=False)
+        plt.savefig(args.plot_R + '_loss_RdBm.png')
+
+        np.savetxt(args.plot_R + '_loss_RdBm' + '.txt', R_mean_loss)
+        quit()
+
+    # ---------------------------------------------------------------------------------------------
+    # Regular training loop
+    # ---------------------------------------------------------------------------------------------
+
     if args.plot_loss:
         plt.figure(1)
         loss_epoch=np.zeros((args.epochs+1))
     
-    # Main training loop
     for epoch in range(1, epochs + 1):
 
         print(f"training epoch {epoch}...",file=sys.stderr)