shallow zip-file copy from codec2 e9d726bf20

8 files changed, 851 insertions, 0 deletions

diff --git a/script/gen_phi0 b/script/gen_phi0
new file mode 100755
index 0000000..af01a31
--- /dev/null
+++ b/script/gen_phi0
@@ -0,0 +1,136 @@
+#!/usr/bin/python3
+
+import math
+
+# Save the points in separate lists 
+points_t1 = []  # A list of tuples (x, r)
+points_t2 = []  # A list of tuples (x, r)
+points_t3 = []  # A list of tuples (x, r)
+
+######################################################
+# Generate
+######################################################
+
+adj = 0.50  # Offset ratio for step of 1, estimate
+
+def loop_even(upper, lower, step, ary):
+  result = []
+  for i in range((int)(upper * step), (int)(lower * step), -1):
+    x = i / step
+    z = math.exp(x + (adj / step))
+    r = math.log((z+1)/(z-1))
+    ary.append( (x, r) )
+
+#####
+# Create initial list which will be sorted desending
+
+# Tbl-1 9:5 in linear steps of 0.5
+loop_even(9.999, 4.999, 2, points_t1)
+
+# Tbl-2 5:1 in linear steps of 1/16
+loop_even(4.999, 0.999, 16, points_t2)
+
+# Tbl-3 log steps:  1/(2^0.5), 1/(2^1), 1/(2^1.5), ...
+for i in range(1, 28):
+    e = (2 ** (i/2))
+    x = 1.0 / e
+    z = math.exp(x+ (.19 / e))
+    r = math.log((z+1)/(z-1))
+    points_t3.append( (x, r) )
+
+
+######################################################
+# Output
+######################################################
+
+print("""
+// phi0.c
+//
+// An approximation of the function
+//
+// This file is generated by the gen_phi0 scripts
+// Any changes should be made to that file, not this one
+
+#include <stdint.h>
+
+#define SI16(f) ((int32_t)(f * (1<<16)))
+
+float phi0( float xf ) {
+
+    int32_t x = SI16(xf);
+
+    if (x >= SI16(10.0f)) return(0.0f);
+""", end="")
+
+
+##################################
+# Tbl-1 9:5 in linear steps of 0.5  (9.5, 9.0, 8.5, 8.0, .., 5.0)
+print("""    else {
+      if (x >= SI16(5.0f)) {
+        int i = 19 - (x >> 15);
+        switch (i) {
+""", end="")
+for i in range(len(points_t1)):
+    #assert(points_t1[i][0] == ((18 - i) / 2))
+    print("{}case {}: return({:.9f}f); // ({:.1f})".format(
+        (10*" "), i, points_t1[i][1], points_t1[i][0]))
+print("{}}}".format(( 8*" ")))
+print("{}}}".format(( 6*" ")))
+
+
+##################################
+# Tbl-2 5:1 in linear steps of 1/16  (4-15/16, 4-7/7, ... 1)
+print("""      else {
+        if (x >= SI16(1.0f)) {
+          int i = 79 - (x >> 12);
+          switch (i) {
+""", end="")
+for i in range(len(points_t2)):
+    #assert(points_t2[i][0] == ((18 - i) / 2))
+    print("{}case {}: return({:.9f}f); // ({:.4f})".format(
+        (12*" "), i, points_t2[i][1], points_t2[i][0]))
+print("{}}}".format((10*" ")))
+print("{}}}".format(( 8*" ")))
+
+
+##################################
+# Tbl-3 log steps:  1/(2^0.5), 1/(2^1), 1/(2^1.5), ...
+# Output as a balanced search
+
+def prnt_cmp(x, ind):
+    print("{}if (x > SI16({:.6f}f)) {{".format((" "*ind), x))
+
+def prnt_rtn(r, ind):
+    print("{}return({:.9f}f);".format((" "*ind), r))
+
+def one_level(pts, ind, dft):
+    #print("# One_Level({})".format(pts))
+    mid = (int)(len(pts)/2)
+    lft = pts[:mid]
+    rgt = pts[mid+1:]
+    x = pts[mid][0]
+    r = pts[mid][1]
+    #print("## {}, {}, {}".format(lft, x, rgt))
+    prnt_cmp(x, ind)
+    if (len(lft)): 
+        one_level(lft, ind+2, r)
+    else:
+        prnt_rtn(r, (ind+2))
+    print("{}}} else {{".format((" "*(ind))))
+    if (len(rgt)): 
+        one_level(rgt, ind+2, dft)
+    else:
+        prnt_rtn(dft, (ind+2))
+    print("{}}}".format((" "*(ind))))
+    
+
+# Start recursive process
+print("{}else {{".format((8*" ")))
+indent = 10
+one_level(points_t3, indent, 10)
+print("{}}}".format((8*" ")))  # End of tb1_1
+print("{}}}".format((6*" ")))  # End of tb1_2
+print("{}}}".format((4*" ")))  # End of tb1_3
+print("{}return(10.0f);".format((4*" ")))
+print("}")
+
diff --git a/script/ofdm_stack_use.py b/script/ofdm_stack_use.py
new file mode 100755
index 0000000..9356cf1
--- /dev/null
+++ b/script/ofdm_stack_use.py
@@ -0,0 +1,150 @@
+#! /usr/bin/python3
+
+""" Find stack usage using
+ - compiler generated tables of stack use per function, *.c.su
+ - run time trace output from compiler added enter/exit calls.
+
+Just for ofdm_stack at this point
+
+This script expects to be run in the .../build_linux/unittest directory!
+"""
+
+COMP_DIR   = 'CMakeFiles/ofdm_stack.dir'
+EXE_FILE   = './ofdm_stack'
+
+import sys
+import os
+import argparse
+import pathlib
+import subprocess
+
+##########################
+# Options
+
+## Trace file (name) or default
+## Use existing trace file or run command
+argparser = argparse.ArgumentParser()
+argparser.add_argument('-f', '--trace_file', action='store', default='function_trace.out',
+                    help='Name of trace file, (default is "function_trace.out"')
+argparser.add_argument('-x', '--exec', action='store_true', default=False,
+                    help='Execute program')
+
+args = argparser.parse_args()
+
+
+##########################
+# Checking
+cwd_path = pathlib.Path.cwd()
+# One simple thing we can handle is running from above unittest (in build_linux)
+if ((cwd_path.name != 'unittest') and
+    (pathlib.Path('unittest').exists())):
+    os.chdir('unittest')
+
+# Required files
+assert(pathlib.Path(COMP_DIR).exists())
+assert(pathlib.Path(COMP_DIR + '/ofdm_stack.c.su').exists())
+assert(pathlib.Path(COMP_DIR + '/__/src/ofdm.c.su').exists())
+
+
+##########################
+# If trace file not found, or option set, run command
+if ( not (pathlib.Path(args.trace_file).exists())):
+    print('Trace file "{}" not found, running program'.format(args.trace_file))
+    args.exec = True
+
+if (args.exec):
+    print('Running program: "{}"'.format(EXE_FILE))
+    assert(pathlib.Path(EXE_FILE))
+    result = subprocess.run([EXE_FILE],
+                             stdout=subprocess.PIPE,
+                             stderr=subprocess.STDOUT)
+    if (result.returncode != 0):
+        print('Error: traced program failed! Output:\n{}'.format(result.stdout))
+assert(pathlib.Path(args.trace_file).exists())
+
+
+##########################
+# Data Structures
+su_data = {} # <function> : <stack_size>
+funcs_used = {} # <function> : count
+
+##########################
+# Read compiler generated tables of stack use per function, *.c.su
+#
+#  ofdm_stack.c:184:6:dummy_code	16	static
+#
+
+p = pathlib.Path(COMP_DIR)
+for fpath in p.rglob('*.c.su'):
+    with fpath.open() as f:
+        for line in f.readlines():
+            try:
+                words = line.split()
+                size = int(words[1])
+                words = words[0].split(':')
+                su_data[words[3]] = size
+            except: pass # skip this line if there are errors
+
+##########################
+# Read trace file, convert addresses to names, track stack
+
+max_stack_depth = 0
+cur_stack_depth = 0
+stack = []  # List of tuples of (function names, cur_stack_depth)
+last_func = 'start'
+
+def walk_stack():
+    trace = ''
+    for entry in stack:
+        trace += entry[0] + ' '
+    return(trace)
+
+# Open trace
+with open(args.trace_file, "r") as tf:
+    for line in tf.readlines():
+        #print('Line: "{}"'.format(line.strip()))
+        words = line.split()
+        # Note addr2line needs addr in hex!
+        addr = words[1]
+        if (words[0] == 'e'):
+            # Note: This could be run once with a pipe if needed for faster operation.
+            result = subprocess.run(['addr2line', '-f', addr, '-e', EXE_FILE],
+                                stdout=subprocess.PIPE)
+            result.check_returncode()
+            # function name is first line of stdout
+            if (result.stdout):
+                lines = result.stdout.decode().split('\n')
+                func = lines[0].strip()
+            else: sys.error('unknown function at address {}'.format(addr))
+
+            if (func != "??"):
+
+                # Push last info
+                stack.append((last_func, cur_stack_depth))
+                last_func = func
+
+                # Update
+                cur_stack_depth += su_data[func]
+                #print('func: "{}" = {}'.format(func, cur_stack_depth))
+                if (cur_stack_depth > max_stack_depth):
+                    max_stack_depth = cur_stack_depth
+                    max_stack_trace = walk_stack()
+
+                # Save info
+                if (func in funcs_used):
+                    funcs_used[func] += 1
+                else:
+                    funcs_used[func] = 1
+
+                # end if (func != "??")
+
+            # end if ('e')
+        elif (words[0] == 'x'):
+            # Only pop functions we pushed
+            if (func in funcs_used):
+                # Pop
+                (last_func, cur_stack_depth) = stack.pop()
+                #print('pop:  "{}" = {}'.format(last_func, cur_stack_depth))
+
+print('Max Stack Depth = {}'.format(max_stack_depth))
+print('Max Stack at: {}'.format(max_stack_trace))
diff --git a/script/phi0_plot.py b/script/phi0_plot.py
new file mode 100755
index 0000000..34a436b
--- /dev/null
+++ b/script/phi0_plot.py
@@ -0,0 +1,80 @@
+#!/usr/bin/python3
+
+""" Plot phi0, and (later) anaylze approximations
+
+"""
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+import math
+import itertools
+import sys
+
+def impl_reference(x):
+    if (x< 9.08e-5 ): r = 10
+    else:
+        z = math.exp(x)
+        r = math.log( (z+1) / (z-1) )
+    return(r)
+
+####
+# Table T1
+
+t1_tbl = []  # A list of tuples (x, r)
+t1_adj = 0.50  # Offset ratio for step of 1, estimate
+
+def t1_loop(upper, lower, step):
+  for i in range((int)(upper * step), (int)(lower * step), -1):
+    x = i / step
+    z = math.exp(x + (t1_adj / step))
+    r = math.log((z+1)/(z-1))
+    t1_tbl.append( (x, r) )
+    print(x, r)
+
+# 9:5 step 1
+t1_loop(9, 4.999, 2)
+
+# 5:1 step 1/16
+t1_loop(4.999, 0.999, 16)
+
+# 1/(2^0.5), 1/(2^1), 1/(2^1.5), ...
+for i in range(1, 28):
+    e = (2 ** (i/2))
+    x = 1.0 / e
+    z = math.exp(x + (0.19 / e))
+    r = math.log((z+1)/(z-1))
+    t1_tbl.append( (x, r) )
+
+def impl_t1(x):
+    if (x > 10): return(0)
+    else:
+        for t in t1_tbl:
+            if (x > t[0]): return(t[1])
+    return(10)
+
+
+
+##########################
+# Plot, scanning from 10 to 0
+x_vals = np.logspace(1, -4, 2000, endpoint=False)
+ref_vals = []
+t1_vals = []
+for x in x_vals:
+    ref_vals.append(impl_reference(x))
+    t1_vals.append(impl_t1(x))
+
+# Sum errors
+errsum = errsum2 = 0
+for i in range(len(ref_vals)):
+    error = t1_vals[i] - ref_vals[i]
+    errsum += error
+    errsum2 += error * error
+print("Net error {}".format(errsum))
+print("avg error {}".format(errsum/len(ref_vals)))
+print("rms error {}".format(math.sqrt(errsum2/len(ref_vals))))
+
+plt.xscale('log')
+plt.plot(x_vals, ref_vals, 'g', x_vals, t1_vals, 'r')
+plt.show()
+
diff --git a/script/separate_all.sh b/script/separate_all.sh
new file mode 100755
index 0000000..5fe6d9d
--- /dev/null
+++ b/script/separate_all.sh
@@ -0,0 +1,11 @@
+#!/bin/sh
+sox $1.wav hts1_$1.wav trim 0 3
+sox $1.wav hts2_$1.wav trim 3 3
+sox $1.wav morig_$1.wav trim 6 2
+sox $1.wav forig_$1.wav trim 8 2
+sox $1.wav ve9qrp_$1.wav trim 10 9.5
+sox $1.wav cq_ref_$1.wav trim 20 9
+sox $1.wav kristoff_$1.wav trim 29.5 4
+sox $1.wav vk5qi_$1.wav trim 33.5 13.5
+sox $1.wav vk5dgr_$1.wav trim 47 10
+
diff --git a/script/subsetvq.sh b/script/subsetvq.sh
new file mode 100755
index 0000000..562a4e8
--- /dev/null
+++ b/script/subsetvq.sh
@@ -0,0 +1,148 @@
+#!/bin/bash
+# subsetvq.sh
+# David Rowe August 2021
+#
+# Script to support:
+# 1. Subset VQ training and listening
+# 1. Training Trellis Vector Quantiser for Codec 2 newamp1, supports octave/trellis.m
+# 2. VQ sorting/optimisation experiments, octave/vq_compare.m
+
+TRAIN=~/Downloads/all_speech_8k.sw
+CODEC2_PATH=$HOME/codec2
+PATH=$PATH:$CODEC2_PATH/build_linux/src:$CODEC2_PATH/build_linux/misc
+K=20
+Kst=2
+Ken=16
+
+# train a new VQ and generate quantised training material
+function train() {
+  fullfile=$TRAIN
+  filename=$(basename -- "$fullfile")
+  extension="${filename##*.}"
+  filename="${filename%.*}"
+  
+  c2sim $fullfile --rateK --rateKout ${filename}.f32
+  echo "ratek=load_f32('../build_linux/${filename}.f32',20); vq_700c_eq; ratek_lim=limit_vec(ratek, 0, 40); save_f32('../build_linux/${filename}_lim.f32', ratek_lim); quit" | \
+  octave -p ${CODEC2_PATH}/octave -qf
+  vqtrain ${filename}_lim.f32 $K 4096 vq_stage1.f32 -s 1e-3 --st $Kst --en $Ken
+
+  # VQ the training file
+  cat ${filename}_lim.f32 | vq_mbest --st $Kst --en $Ken -k $K -q vq_stage1.f32 > ${filename}_test.f32
+}
+
+function listen_vq() {
+  vq_fn=$1
+  dec=$2
+  EbNodB=$3
+  fullfile=$4
+  filename=$(basename -- "$fullfile")
+  extension="${filename##*.}"
+  filename="${filename%.*}"
+  
+  fullfile_out=$5
+  do_trellis=$6
+  sox_options='-t raw -e signed-integer -b 16' 
+  sox $fullfile $sox_options - | c2sim - --rateK --rateKout ${filename}.f32
+  
+  echo "ratek=load_f32('../build_linux/${filename}.f32',20); vq_700c_eq; ratek_lim=limit_vec(ratek, 0, 40); save_f32('../build_linux/${filename}_lim.f32', ratek_lim); quit" | \
+  octave -p ${CODEC2_PATH}/octave -qf
+
+  if [ "$do_trellis" -eq 0 ]; then
+     echo "pkg load statistics; vq_compare(action='vq_file', '${vq_fn}', ${dec}, ${EbNodB}, '${filename}_lim.f32', '${filename}_test.f32'); quit" \ |
+     octave -p ${CODEC2_PATH}/octave -qf
+  else
+     echo "pkg load statistics; trellis; vq_file('${vq_fn}', ${dec}, ${EbNodB}, '${filename}_lim.f32', '${filename}_test.f32'); quit" \ |
+     octave -p ${CODEC2_PATH}/octave -qf
+  fi
+  
+  if [ "$fullfile_out" = "aplay" ]; then
+     sox $fullfile $sox_options - | c2sim - --rateK --rateKin ${filename}_test.f32 -o - | aplay -f S16_LE
+  else
+     sox $fullfile $sox_options - | c2sim - --rateK --rateKin ${filename}_test.f32 -o - | sox -t .s16 -r 8000 -c 1 - ${fullfile_out}
+  fi
+     
+}
+
+function print_help {
+    echo
+    echo "Trellis/VQ optimisation support script"
+    echo
+    echo "  usage ./train_trellis.sh [-x] [-t] [-v vq.f32 in.wav out.wav] [-e EbNodB] [-d dec]"
+    echo
+    echo "    -x                         debug mode; trace script execution"
+    echo "    -t                         train VQ and generate a fully quantised version of training vectors"
+    echo "    -v  vq.f32 in.wav out.wav  synthesise an output file out.wav from in.raw, using the VQ vq.f32"
+    echo "    -v  vq.f32 in.wav aplay    synthesise output, play immediately using aplay, using the VQ vq.f32"
+    echo "    -e  EbNodB                 Eb/No in dB for AWGn channel simulation (error insertion)"
+    echo "    -d  dec                    decimation/interpolation rate"
+    echo "    -r                         use trellis decoder"
+    echo
+    exit
+}
+
+# command line arguments to select function
+
+if [ $# -lt 1 ]; then
+    print_help
+fi
+
+do_train=0
+do_vq=0
+do_trellis=0
+EbNodB=100
+dec=1
+POSITIONAL=()
+while [[ $# -gt 0 ]]
+do
+key="$1"
+case $key in
+    -x)
+        set -x	
+        shift
+    ;;
+    -t)
+        do_train=1
+	shift
+    ;;
+    -v)
+        do_vq=1
+	vq_fn="$2"
+	in_wav="$3"
+	out_wav="$4"
+	shift
+	shift
+	shift
+	shift
+    ;;
+    -r)
+        do_trellis=1
+	shift
+    ;;
+    -d)
+	dec="$2"
+	shift
+	shift	
+	;;
+    -e)
+	EbNodB="$2"
+	shift
+	shift	
+	;;
+    -h)
+        print_help	
+    ;;
+    *)
+    POSITIONAL+=("$1") # save it in an array for later
+    shift
+    ;;
+esac
+done
+set -- "${POSITIONAL[@]}" # restore positional parameters
+
+if [ $do_train -eq 1 ]; then
+    train
+fi
+
+if [ $do_vq -eq 1 ]; then
+  listen_vq ${vq_fn} ${dec} ${EbNodB} ${in_wav} ${out_wav} ${do_trellis}
+fi
diff --git a/script/test_2020x.sh b/script/test_2020x.sh
new file mode 100755
index 0000000..541091c
--- /dev/null
+++ b/script/test_2020x.sh
@@ -0,0 +1,166 @@
+#!/bin/bash -x
+# test_2020x.sh
+# David Rowe Feb 2022
+#
+# Script to support testing experimental 2020A and 2020B modes and 700E control.
+
+CODEC2_PATH=$HOME/codec2
+PATH=$PATH:$CODEC2_PATH/build_linux/src:$CODEC2_PATH/build_linux/misc
+FADING_DIR=$CODEC2_PATH/build_linux/unittest
+No_AWGN=-20
+No_AWGN_LOW=-17
+No_Multipath=-25
+serial=0
+compressor_gain=6
+
+# Approximation of Hilbert clipper analog compressor
+function analog_compressor {
+    input_file=$1
+    output_file=$2
+    gain=$3
+    cat $input_file | ch - - 2>/dev/null | \
+    ch - - --No -100 --clip 16384 --gain $gain 2>/dev/null | \
+    # final line prints peak and CPAPR for SSB
+    ch - - --clip 16384 |
+    # manually adjusted to get similar peak levels for SSB and FreeDV
+    sox -t .s16 -r 8000 -c 1 -v 0.85 - -t .s16 $output_file
+}
+
+function run_sim_ssb() {
+    fullfile=$1
+    filename=$(basename -- "$fullfile")
+    extension="${filename##*.}"
+    filename="${filename%.*}"
+    channel=$2
+    No=-100
+    if [ "$channel" == "awgn" ]; then
+        channel_opt=""
+        No=$No_AWGN
+    fi
+    if [ "$channel" == "awgnlow" ]; then
+        channel_opt=""
+        No=$No_AWGN_LOW
+    fi
+    if [ "$channel" == "mpp" ] || [ "$channel" == "mpd" ]; then
+        channel_opt='--'${channel}
+        No=$No_Multipath
+    fi
+    fn=${filename}_ssb_${channel}.wav
+    analog_compressor ${fullfile} ${filename}_ssb.raw ${compressor_gain}
+    tmp=$(mktemp)
+    ch ${filename}_ssb.raw $tmp --No $No ${channel_opt} --fading_dir ${FADING_DIR} 2>t.txt
+    cat $tmp | sox -t .s16 -r 8000 -c 1 - ${fn} trim 0 6
+    snr=$(cat t.txt | grep "SNR3k(dB):" | tr -s ' ' | cut -d' ' -f3)
+    
+    echo "<tr>"
+    echo "<td><a href=\"${fn}\">${serial}</a></td><td>ssb</td><td></td><td></td><td>${channel}</td><td>${snr}</td>"
+    echo "</tr>"
+    serial=$((serial+1))
+}
+
+function run_sim() {
+    fullfile=$1
+    filename=$(basename -- "$fullfile")
+    extension="${filename##*.}"
+    filename="${filename%.*}"
+    mode=$2
+    if [ "$mode" == "700E" ] || [ "$mode" == "700D" ]; then
+        rateHz=8000
+    else
+        rateHz=16000
+    fi
+    clip=$3
+    if [ "$clip" == "clip" ]; then
+        clipflag=1
+        clip_html="yes"
+    else
+        clipflag=0
+        clip_html="no"
+    fi       
+    channel=$4
+    No=-100
+    if [ "$channel" == "awgn" ]; then
+        channel_opt=""
+        No=$No_AWGN
+    fi
+    if [ "$channel" == "awgnlow" ]; then
+        channel_opt=""
+        No=$No_AWGN_LOW
+    fi
+    if [ "$channel" == "mpp" ] || [ "$channel" == "mpd" ]; then
+        channel_opt='--'${channel}
+        No=$No_Multipath
+    fi
+
+    indopt=$5
+    indopt_flag=""
+    indopt_html="no"
+    indopt_str=""
+    if [ "$indopt" == "indopt" ]; then
+        indopt_flag="--indopt 1"
+        indopt_str="_indopt"
+        indopt_html="yes"
+    fi
+    if [ "$indopt" == "no_indopt" ]; then
+        indopt_flag="--indopt 0"
+        indopt_str="_no_indopt"
+    fi
+    
+    fn=${filename}_${mode}_${clip}_${channel}${indopt_str}.wav
+    tmp=$(mktemp)
+    # note we let ch finish to get SNR stats (trim at end of sox causes an early termination)
+    freedv_tx ${mode} ${fullfile} - --clip ${clipflag} ${indopt_flag} | \
+    ch - $tmp --No $No ${channel_opt} --fading_dir ${FADING_DIR} 2>t.txt
+    freedv_rx ${mode} ${indopt_flag} $tmp - | \
+    sox -t .s16 -r ${rateHz} -c 1 - ${fn} trim 0 6
+    snr=$(cat t.txt | grep "SNR3k(dB):" | tr -s ' ' | cut -d' ' -f3)
+
+    echo "<tr>"
+    echo "<td><a href=\"${fn}\">${serial}</a></td><td>${mode}</td><td>${clip_html}</td><td>${indopt_html}</td><td>${channel}</td><td>${snr}</td>"
+    echo "</tr>"
+    serial=$((serial+1))
+}
+
+# convert speech input file to format we need
+SPEECH_IN_16k_WAV=~/Downloads/speech_orig_16k.wav 
+SPEECH_IN_16k_RAW=speech_orig_16k.raw
+SPEECH_IN_8k_RAW=speech_orig_8k.raw
+sox $SPEECH_IN_16k_WAV -t .s16 $SPEECH_IN_16k_RAW
+sox $SPEECH_IN_16k_WAV -t .s16 -r 8000 $SPEECH_IN_8k_RAW
+
+echo "<html><table>"
+echo "<tr><th>Serial</th><th>Mode</th><th>Clip</th><th>index_opt</th><th>Channel</th><th>SNR (dB)</th></tr>"
+
+# run simulations
+
+run_sim_ssb $SPEECH_IN_8k_RAW awgn
+run_sim_ssb $SPEECH_IN_8k_RAW mpp
+run_sim_ssb $SPEECH_IN_8k_RAW mpd
+
+run_sim $SPEECH_IN_16k_RAW 2020 noclip clean
+run_sim $SPEECH_IN_8k_RAW 700E clip clean
+
+run_sim $SPEECH_IN_16k_RAW 2020 noclip awgn
+run_sim $SPEECH_IN_16k_RAW 2020 noclip mpp
+run_sim $SPEECH_IN_16k_RAW 2020 noclip mpd
+run_sim $SPEECH_IN_16k_RAW 2020 clip awgn
+run_sim $SPEECH_IN_16k_RAW 2020 clip mpp
+run_sim $SPEECH_IN_16k_RAW 2020 clip mpd
+
+run_sim $SPEECH_IN_16k_RAW 2020B clip awgn indopt
+run_sim $SPEECH_IN_16k_RAW 2020B clip mpp  indopt
+run_sim $SPEECH_IN_16k_RAW 2020B clip mpp  no_indopt
+run_sim $SPEECH_IN_16k_RAW 2020B clip mpd  indopt
+run_sim $SPEECH_IN_16k_RAW 2020B clip mpd  no_indopt
+
+run_sim $SPEECH_IN_8k_RAW 700E clip awgn
+run_sim $SPEECH_IN_8k_RAW 700E clip mpp
+run_sim $SPEECH_IN_8k_RAW 700E clip mpd
+
+# Low SNR samples
+run_sim_ssb $SPEECH_IN_8k_RAW awgnlow
+run_sim $SPEECH_IN_8k_RAW 700E clip awgnlow
+run_sim $SPEECH_IN_16k_RAW 2020 clip awgnlow
+run_sim $SPEECH_IN_16k_RAW 2020A clip awgnlow indopt
+
+exit
diff --git a/script/train_700c_quant.sh b/script/train_700c_quant.sh
new file mode 100755
index 0000000..4fd0c72
--- /dev/null
+++ b/script/train_700c_quant.sh
@@ -0,0 +1,107 @@
+#!/bin/bash -x
+# train_700C_quant.sh
+# David Rowe May 2019
+#
+# Training a Vector Quantiser (VQ) for Codec 2 700C
+# This is a two stage VQ with 512 entries (9 bits) per stage
+# Also used to support other VQ experiments, such as the effect of the
+# post filter and an experimental equaliser, see octave/vq_700c_eq.m
+
+SRC=~/Downloads/all_speech_8k.sw
+CODEC2_BUILD=/home/david/codec2/build_linux
+K=20
+SAMPLES=~/tmp/c2vec_pass2
+
+# train a new VQ
+function train() {
+  # c2enc can dump "feature vectors" that contain the current VQ input
+  $CODEC2_BUILD/src/c2enc 700C $SRC /dev/null --mlfeat feat.f32
+  # extract VQ input as training data, then train two stage VQ
+  $CODEC2_BUILD/misc/extract -s 1 -e $K -t 41 feat.f32 stage0_in.f32
+  $CODEC2_BUILD/misc/vqtrain stage0_in.f32 $K 512 vq_stage1.f32 -s 1e-3 -r stage1_in.f32
+  $CODEC2_BUILD/misc/vqtrain stage1_in.f32 $K 512 vq_stage2.f32 -s 1e-3
+}
+
+# encode/decode a file with the stock codec2 700C VQ
+function test_a() {
+  b=$(basename "$1" .raw)
+  $CODEC2_BUILD/src/c2enc 700C $1'.raw' - --var | $CODEC2_BUILD/src/c2dec 700C - - | sox -q -t .s16 -c 1 -r 8000 -b 16  - $SAMPLES/$b'_a.wav'
+}
+
+# stock 700C VQ like test_a, but no postfilter
+function test_a_nopf() {
+  b=$(basename "$1" .raw)
+  $CODEC2_BUILD/src/c2enc 700C $1'.raw' - --var | $CODEC2_BUILD/src/c2dec 700C - - --nopf | \
+  sox -q -t .s16 -c 1 -r 8000 -b 16  - $SAMPLES/$b'_a_nopf.wav'
+}
+
+# encode/decode a file with the new VQ we just trained above
+function test_b() {
+  b=$(basename "$1" .raw)
+  $CODEC2_BUILD/src/c2enc 700C $1'.raw' - --loadcb 1 vq_stage1.f32 --loadcb 2 vq_stage2.f32 --var | \
+  $CODEC2_BUILD/src/c2dec 700C - - --loadcb 1 vq_stage1.f32 --loadcb 2 vq_stage2.f32 | sox -q -t .s16 -c 1 -r 8000 -b 16  - $SAMPLES/$b'_b.wav'
+}
+
+# just like b but no postfilter
+function test_b_nopf() {
+  b=$(basename "$1" .raw)
+  $CODEC2_BUILD/src/c2enc 700C $1'.raw' - --loadcb 1 vq_stage1.f32 --loadcb 2 vq_stage2.f32 --var | \
+  $CODEC2_BUILD/src/c2dec 700C - - --loadcb 1 vq_stage1.f32 --loadcb 2 vq_stage2.f32 $2 --nopf | \
+  sox -q -t .s16 -c 1 -r 8000 -b 16  - $SAMPLES/$b'_b_nopf.wav'
+}
+
+# pass an unquantised rate K vector through to the decoder as a control
+function test_uq() {
+  b=$(basename "$1" .raw)
+  $CODEC2_BUILD/src/c2enc 700C $1'.raw' $b'.bin' --mlfeat $b'_feat.f32'
+  $CODEC2_BUILD/misc/extract -s 1 -e 20 -t 41 $b'_feat.f32' $b'_ratek.f32'
+  $CODEC2_BUILD/src/c2dec 700C $b'.bin' - --loadratek $b'_ratek.f32' | sox -q -t .s16 -c 1 -r 8000 -b 16  - $SAMPLES/$b'_uq.wav'
+}
+
+# extract features for use in octave octave/vq_700c_eq.m
+function feat() {
+  RAW_FILES="../raw/hts1a ../raw/hts2a ../raw/vk5qi ../raw/cq_ref ../raw/ve9qrp_10s $HOME/Downloads/ma01_01 $HOME/Downloads/c01_01_8k $HOME/Downloads/cq_freedv_8k $HOME/Downloads/cq_freedv_8k_lfboost $HOME/Downloads/cq_freedv_8k_hfcut "
+  for f in $RAW_FILES
+  do
+    b=$(basename "$f" .raw)
+    $CODEC2_BUILD/src/c2enc 700C $f'.raw' $b'.bin' --mlfeat $b'_feat.f32'
+  done
+}
+
+# generate a bunch of test samples for a listening test
+function listen() {
+  RAW_FILES="../raw/hts1a ../raw/hts2a ../raw/vk5qi ../raw/cq_ref ../raw/ve9qrp_10s $HOME/Downloads/ma01_01 $HOME/Downloads/c01_01_8k $HOME/Downloads/cq_freedv_8k"
+  for f in $RAW_FILES
+  do
+    test_a $f
+    test_a_nopf $f
+    test_b $f
+    test_b_nopf $f
+    test_uq $f
+  done
+}
+
+# Generate a bunch of test samples for VQ equalisation listening tests.  Assumes
+# Octave has generated rate K quantised .f32 files 
+function listen_vq_eq() {
+  FILES="hts1a hts2a vk5qi cq_ref ve9qrp_10s ma01_01 c01_01_8k cq_freedv_8k cq_freedv_8k_lfboost cq_freedv_8k_hfcut"
+  for f in $FILES
+  do
+    # try equaliser wth train_120 VQ
+    $CODEC2_BUILD/src/c2dec 700C $f'.bin' - --loadratek $f'_vq2.f32' | sox -q -t .s16 -c 1 -r 8000 -b 16  - $SAMPLES/$f'_vq2.wav'
+    $CODEC2_BUILD/src/c2dec 700C $f'.bin' - --loadratek $f'_vq2_eq1.f32' | sox -q -t .s16 -c 1 -r 8000 -b 16  - $SAMPLES/$f'_vq2_eq1.wav'
+    $CODEC2_BUILD/src/c2dec 700C $f'.bin' - --loadratek $f'_vq2_eq2.f32' | sox -q -t .s16 -c 1 -r 8000 -b 16  - $SAMPLES/$f'_vq2_eq2.wav'
+    # try equaliser wth train_all_speech VQ
+    #$CODEC2_BUILD/src/c2dec 700C $f'.bin' - --loadratek $f'_vq2_as.f32' | sox -q -t .s16 -c 1 -r 8000 -b 16  - $SAMPLES/$f'_vq2_as.wav'
+    #$CODEC2_BUILD/src/c2dec 700C $f'.bin' - --loadratek $f'_vq2_as_eq.f32' | sox -q -t .s16 -c 1 -r 8000 -b 16  - $SAMPLES/$f'_vq2_as_eq.wav'
+  done
+}
+
+mkdir -p $SAMPLES
+
+# choose which function to run here
+#train
+#listen
+#feat
+listen_vq_eq
+
diff --git a/script/train_sub_quant.sh b/script/train_sub_quant.sh
new file mode 100755
index 0000000..e85b098
--- /dev/null
+++ b/script/train_sub_quant.sh
@@ -0,0 +1,53 @@
+#!/bin/bash -x
+# train_sub_quant.sh
+# David Rowe May 2021
+#
+# Training and testing Vector Quantisers (VQ) for Codec 2 newamp1, in
+# this case training on a subset
+
+TRAIN=~/Downloads/all_speech_8k.sw
+CODEC2_PATH=$HOME/codec2
+PATH=$PATH:$CODEC2_PATH/build_linux/src:$CODEC2_PATH/build_linux/misc
+K=20
+Kst=2
+Ken=16
+
+# train a new VQ
+function train() {
+  fullfile=$TRAIN
+  filename=$(basename -- "$fullfile")
+  extension="${filename##*.}"
+  filename="${filename%.*}"
+  
+  c2sim $fullfile --rateK --rateKout ${filename}.f32
+  echo "ratek=load_f32('../build_linux/${filename}.f32',20); vq_700c_eq; ratek_lim=limit_vec(ratek, 0, 40); save_f32('../build_linux/${filename}_lim.f32', ratek_lim); quit" | \
+  octave -p ${CODEC2_PATH}/octave -qf
+  vqtrain ${filename}_lim.f32 $K 4096 vq_stage1.f32 -s 1e-3 --st $Kst --en $Ken
+}
+
+function listen() {
+  fullfile=$1
+  filename=$(basename -- "$fullfile")
+  extension="${filename##*.}"
+  filename="${filename%.*}"
+  
+  c2sim $fullfile --rateK --rateKout ${filename}.f32
+  echo "ratek=load_f32('../build_linux/${filename}.f32',20); vq_700c_eq; ratek_lim=limit_vec(ratek, 0, 40); save_f32('../build_linux/${filename}_lim.f32', ratek_lim); quit" | \
+  octave -p ${CODEC2_PATH}/octave -qf
+  cat ${filename}_lim.f32 | vq_mbest --st $Kst --en $Ken -k $K -q vq_stage1.f32 > ${filename}_test.f32
+  c2sim $fullfile --rateK --rateKin ${filename}_test.f32 -o - | sox -t .s16 -r 8000 -c 1 - ${filename}_sub.wav
+  c2sim $fullfile --rateK --newamp1vq -o - | sox -t .s16 -r 8000 -c 1 - ${filename}_newamp1.wav
+}
+
+# choose which function to run here
+train
+# these two samples are inside training database
+listen ~/Downloads/fish_8k.sw
+listen ~/Downloads/cap_8k.sw
+# two samples from outside training database
+listen $CODEC2_PATH/raw/big_dog.raw
+listen $CODEC2_PATH/raw/hts2a.raw
+# these two samples are inside training database, but with LPF at 3400 Hz outside of subset
+listen ~/Downloads/fish_8k_lp.sw
+listen ~/Downloads/cap_8k_lp.sw
+