From ac7c48b4dee99d4c772f133d70d8d1b38262fcd2 Mon Sep 17 00:00:00 2001
From: Author Name <david@rowetel.com>
Date: Fri, 7 Jul 2023 12:20:59 +0930
Subject: shallow zip-file copy from codec2 e9d726bf20

---
 misc/16_8_short.c       |  47 ++++++
 misc/CMakeLists.txt     |  35 +++++
 misc/de.c               |  59 +++++++
 misc/extract.c          | 117 ++++++++++++++
 misc/ge_train.c         | 306 +++++++++++++++++++++++++++++++++++++
 misc/genlsp.c           | 183 ++++++++++++++++++++++
 misc/mksine.c           |  54 +++++++
 misc/pre.c              |  59 +++++++
 misc/raw2h.c            |  45 ++++++
 misc/speexnoisesup.c    |  58 +++++++
 misc/tcodec2.c          | 220 ++++++++++++++++++++++++++
 misc/tdec.c             | 141 +++++++++++++++++
 misc/timpulse.c         | 128 ++++++++++++++++
 misc/tinterp.c          | 151 ++++++++++++++++++
 misc/tlininterp.c       | 153 +++++++++++++++++++
 misc/tnlp.c             | 164 ++++++++++++++++++++
 misc/tprede.c           |  53 +++++++
 misc/tquant.c           | 214 ++++++++++++++++++++++++++
 misc/tsrc.c             | 109 +++++++++++++
 misc/vq_binary_switch.c | 296 +++++++++++++++++++++++++++++++++++
 misc/vq_mbest.c         | 302 ++++++++++++++++++++++++++++++++++++
 misc/vqtrain.c          | 399 ++++++++++++++++++++++++++++++++++++++++++++++++
 22 files changed, 3293 insertions(+)
 create mode 100644 misc/16_8_short.c
 create mode 100644 misc/CMakeLists.txt
 create mode 100644 misc/de.c
 create mode 100644 misc/extract.c
 create mode 100644 misc/ge_train.c
 create mode 100644 misc/genlsp.c
 create mode 100644 misc/mksine.c
 create mode 100644 misc/pre.c
 create mode 100644 misc/raw2h.c
 create mode 100644 misc/speexnoisesup.c
 create mode 100644 misc/tcodec2.c
 create mode 100644 misc/tdec.c
 create mode 100644 misc/timpulse.c
 create mode 100644 misc/tinterp.c
 create mode 100644 misc/tlininterp.c
 create mode 100644 misc/tnlp.c
 create mode 100644 misc/tprede.c
 create mode 100644 misc/tquant.c
 create mode 100644 misc/tsrc.c
 create mode 100644 misc/vq_binary_switch.c
 create mode 100644 misc/vq_mbest.c
 create mode 100644 misc/vqtrain.c

(limited to 'misc')

diff --git a/misc/16_8_short.c b/misc/16_8_short.c
new file mode 100644
index 0000000..a23722c
--- /dev/null
+++ b/misc/16_8_short.c
@@ -0,0 +1,47 @@
+/*
+   16_8_short.c
+   David Rowe
+   October 2018
+
+   Utility for resampling raw files from 16 to 8 kHz.
+*/
+
+#include <assert.h>
+#include <math.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include "codec2_fdmdv.h"
+
+#define N8                        160 /* procssing buffer size at 8 kHz */
+#define N16             (N8*FDMDV_OS)
+
+int main(int argc, char *argv[]) {
+    short in16k_short[FDMDV_OS_TAPS_16K + N16];
+    FILE *f16;
+    short out8k_short[N16];
+    FILE *f8;
+    int i;
+
+    if (argc != 3) {
+        fprintf(stderr, "usage: %s 16kHz.raw 8kHz.raw\n", argv[0]);
+        exit(1);
+    }
+    f16 = fopen(argv[1], "rb");
+    assert(f16 != NULL);
+    f8 = fopen(argv[2], "wb");
+    assert(f8 != NULL);
+
+    /* clear filter memories */
+
+    for(i=0; i<FDMDV_OS_TAPS_16K; i++)
+	in16k_short[i] = 0;
+
+    while(fread(&in16k_short[FDMDV_OS_TAPS_16K], sizeof(short), N16, f16) == N16) {
+	fdmdv_16_to_8_short(out8k_short, &in16k_short[FDMDV_OS_TAPS_16K], N8);
+	fwrite(out8k_short, sizeof(short), N8, f8);
+    }
+
+    fclose(f16);
+    fclose(f8);
+    return 0;
+}
diff --git a/misc/CMakeLists.txt b/misc/CMakeLists.txt
new file mode 100644
index 0000000..55728b4
--- /dev/null
+++ b/misc/CMakeLists.txt
@@ -0,0 +1,35 @@
+add_definitions(-DFLOATING_POINT -DVAR_ARRAYS)
+include_directories(../src)
+
+add_executable(mksine mksine.c)
+target_link_libraries(mksine m)
+
+add_executable(16_8_short 16_8_short.c ../src/fdmdv.c ../src/kiss_fft.c)
+target_link_libraries(16_8_short codec2)
+
+add_executable(extract extract.c)
+
+add_executable(vqtrain vqtrain.c)
+target_link_libraries(vqtrain m)
+
+add_executable(raw2h raw2h.c)
+target_link_libraries(raw2h codec2)
+
+add_executable(tnlp tnlp.c)
+target_link_libraries(tnlp codec2)
+
+add_executable(tlininterp tlininterp.c)
+
+add_executable(tdec tdec.c)
+
+add_executable(timpulse timpulse.c)
+target_link_libraries(timpulse m)
+
+add_executable(vq_mbest vq_mbest.c)
+target_link_libraries(vq_mbest codec2)
+
+add_executable(vq_binary_switch vq_binary_switch.c)
+target_link_libraries(vq_binary_switch m)
+
+add_executable(pre pre.c)
+target_link_libraries(pre codec2)
diff --git a/misc/de.c b/misc/de.c
new file mode 100644
index 0000000..9726743
--- /dev/null
+++ b/misc/de.c
@@ -0,0 +1,59 @@
+/*
+   de.c
+   David Rowe
+   Sep 26 2012
+
+   Takes audio from a file, de-emphasises, and sends to output file.
+*/
+
+#include <assert.h>
+#include <errno.h>
+#include <math.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include "lpc.h"
+
+#define N 80
+
+int main(int argc, char *argv[]) {
+    FILE  *fin, *fout;
+    short  buf[N];
+    float  Sn[N], Sn_de[N];
+    float  de_mem = 0.0;
+    int    i;
+
+    if (argc != 3) {
+	printf("usage: de InputRawSpeechFile OutputRawSpeechFile\n");
+	printf("e.g    de input.raw output.raw");
+	exit(1);
+    }
+
+    if (strcmp(argv[1], "-")  == 0) fin = stdin;
+    else if ( (fin = fopen(argv[1],"rb")) == NULL ) {
+	fprintf(stderr, "Error opening input speech file: %s: %s.\n",
+         argv[1], strerror(errno));
+	exit(1);
+    }
+
+    if (strcmp(argv[2], "-") == 0) fout = stdout;
+    else if ( (fout = fopen(argv[2],"wb")) == NULL ) {
+	fprintf(stderr, "Error opening output speech file: %s: %s.\n",
+         argv[2], strerror(errno));
+	exit(1);
+    }
+
+    while(fread(buf, sizeof(short), N, fin) == N) {
+	for(i=0; i<N; i++)
+	    Sn[i] = buf[i];
+	de_emp(Sn_de, Sn, &de_mem, N);
+	for(i=0; i<N; i++)
+	    buf[i] = Sn_de[i];
+	fwrite(buf, sizeof(short), N, fout);
+    }
+
+    fclose(fin);
+    fclose(fout);
+
+    return 0;
+}
diff --git a/misc/extract.c b/misc/extract.c
new file mode 100644
index 0000000..b4f38db
--- /dev/null
+++ b/misc/extract.c
@@ -0,0 +1,117 @@
+/*
+  extract.c
+  david Rowe Jan 2019
+
+  Extracts sub sets of vectors from .f32 files, used for LPCNet VQ experiments.
+*/
+
+#include <assert.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <getopt.h>
+
+#define NB_FEATURES 55 /* number of cols per row */
+
+int main(int argc, char *argv[]) {
+    FILE *fin, *fout;
+    int st = 0;
+    int en = 17;
+    int stride = NB_FEATURES;
+    float gain = 1.0;
+    int frame_delay = 1;
+    float pred = 0.0;
+    int removemean = 0;
+    float lower = -1E32;
+    
+    static struct option long_options[] = {
+        {"startcol",   required_argument, 0, 's'},
+        {"endcol",     required_argument, 0, 'e'},
+        {"stride",     required_argument, 0, 't'},
+        {"gain",       required_argument, 0, 'g'},
+        {"pred",       required_argument, 0, 'p'},
+        {"delay",      required_argument, 0, 'd'},
+        {"removemean", no_argument, 0, 'm'},
+        {"lower",      required_argument, 0, 'l'},
+        {0, 0, 0, 0}
+    };
+
+    int opt_index = 0;
+    int c;
+    
+    while ((c = getopt_long (argc, argv, "s:e:t:g:p:d:ml:", long_options, &opt_index)) != -1) {
+        switch (c) {
+        case 's':
+            st = atoi(optarg);
+            break;
+        case 'e':
+            en = atoi(optarg);
+            break;
+        case 't':
+            stride = atoi(optarg);
+            break;
+        case 'g':
+            gain = atof(optarg);
+            break;
+        case 'p':
+            pred = atof(optarg);
+            break;
+        case 'd':
+            frame_delay = atoi(optarg);
+            break;
+        case 'm':
+            removemean = 1;
+            break;
+        case 'l':
+            lower = atof(optarg);
+            break;
+        default:
+        helpmsg:
+            fprintf(stderr, "usage: %s  -s startCol -e endCol [-t strideCol -g gain -p predCoeff -d framesDelay --removemean --lower] input.f32 output.f32\n", argv[0]);
+            exit(1);
+        }
+    }
+    if ( (argc - optind) < 2) {
+        fprintf(stderr, "Too few arguments\n");
+        goto helpmsg;
+    }
+ 
+    fin = fopen(argv[optind],"rb"); assert(fin != NULL);
+    fout = fopen(argv[optind+1],"wb"); assert(fout != NULL);
+    printf("extracting from %d to %d inclusive (stride %d) ... gain = %f pred = %f frame_delay = %d\n",
+           st, en, stride, gain, pred, frame_delay);
+   
+    float features[stride], features_prev[frame_delay][stride], delta[stride];
+    int i,f,wr=0;
+    
+    for (f=0; f<frame_delay; f++)
+        for(i=0; i<stride; i++)
+            features_prev[f][i] = 0.0;
+
+    while((fread(features, sizeof(float), stride, fin) == stride)) {
+	float mean = 0.0;
+	for(i=st; i<=en; i++)
+	    mean += features[i];
+	mean /= (en-st+1);
+	if (removemean) {
+	    for(i=st; i<=en; i++)
+		features[i] -= mean;
+	}
+	for(i=st; i<=en; i++) {
+	    delta[i] = gain*(features[i] - pred*features_prev[frame_delay-1][i]);
+	}
+	if (mean > lower) {
+	    fwrite(&delta[st], sizeof(float), en-st+1, fout);
+	    wr++;
+	}
+	for (f=frame_delay-1; f>0; f--)
+	    for(i=0; i<stride; i++)
+		features_prev[f][i] = features_prev[f-1][i];
+	for(i=0; i<stride; i++)
+	    features_prev[0][i] = features[i];
+    }
+
+    fclose(fin); fclose(fout);
+    fprintf(stderr, "%d extracted\n", wr);
+    return 0;
+}
+
diff --git a/misc/ge_train.c b/misc/ge_train.c
new file mode 100644
index 0000000..db786fc
--- /dev/null
+++ b/misc/ge_train.c
@@ -0,0 +1,306 @@
+/*
+ ge_train.c
+ Jean Marc Valin Feb 2012
+
+ Joint pitch and energy VQ training program
+
+ usage:
+
+   cat GE | ./ge_train 2 1000000 8 > quantized
+
+ The first column is the log2 of the pitch compared to the lowest freq,
+ so log2(wo/pi*4000/50) where wo is the frequency your patch outputs. The
+ second column is the energy in dB, so 10*log10(1e-4+E)
+*/
+
+/*
+  Copyright (C) 2012 Jean-Marc Valin
+
+  All rights reserved.
+
+  This program is free software; you can redistribute it and/or modify
+  it under the terms of the GNU Lesser General Public License version 2.1, as
+  published by the Free Software Foundation.  This program is
+  distributed in the hope that it will be useful, but WITHOUT ANY
+  WARRANTY; without even the implied warranty of MERCHANTABILITY or
+  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public
+  License for more details.
+
+  You should have received a copy of the GNU Lesser General Public License
+  along with this program; if not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <valgrind/memcheck.h>
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+
+#define MIN(a,b) ((a)<(b)?(a):(b))
+//#define COEF 0.0
+
+static float COEF[2] = {0.8, 0.9};
+//static float COEF[2] = {0.0, 0.};
+
+#define MAX_ENTRIES 16384
+
+void compute_weights2(const float *x, const float *xp, float *w, int ndim)
+{
+  w[0] = 30;
+  w[1] = 1;
+  if (x[1]<0)
+  {
+     w[0] *= .6;
+     w[1] *= .3;
+  }
+  if (x[1]<-10)
+  {
+     w[0] *= .3;
+     w[1] *= .3;
+  }
+  /* Higher weight if pitch is stable */
+  if (fabs(x[0]-xp[0])<.2)
+  {
+     w[0] *= 2;
+     w[1] *= 1.5;
+  } else if (fabs(x[0]-xp[0])>.5) /* Lower if not stable */
+  {
+     w[0] *= .5;
+  }
+
+  /* Lower weight for low energy */
+  if (x[1] < xp[1]-10)
+  {
+     w[1] *= .5;
+  }
+  if (x[1] < xp[1]-20)
+  {
+     w[1] *= .5;
+  }
+
+  //w[0] = 30;
+  //w[1] = 1;
+
+  /* Square the weights because it's applied on the squared error */
+  w[0] *= w[0];
+  w[1] *= w[1];
+
+}
+
+int find_nearest_weighted(const float *codebook, int nb_entries, float *x, const float *w, int ndim)
+{
+  int i, j;
+  float min_dist = 1e15;
+  int nearest = 0;
+
+  for (i=0;i<nb_entries;i++)
+  {
+    float dist=0;
+    for (j=0;j<ndim;j++)
+      dist += w[j]*(x[j]-codebook[i*ndim+j])*(x[j]-codebook[i*ndim+j]);
+    if (dist<min_dist)
+    {
+      min_dist = dist;
+      nearest = i;
+    }
+  }
+  return nearest;
+}
+
+int quantize_ge(const float *x, const float *codebook1, int nb_entries, float *xq, int ndim)
+{
+  int i, n1;
+  float err[ndim];
+  float w[ndim];
+
+  compute_weights2(x, xq, w, ndim);
+
+  for (i=0;i<ndim;i++)
+    err[i] = x[i]-COEF[i]*xq[i];
+  n1 = find_nearest_weighted(codebook1, nb_entries, err, w, ndim);
+
+  for (i=0;i<ndim;i++)
+  {
+    xq[i] = COEF[i]*xq[i] + codebook1[ndim*n1+i];
+    err[i] -= codebook1[ndim*n1+i];
+  }
+  return 0;
+}
+
+void split(float *codebook, int nb_entries, int ndim)
+{
+  int i,j;
+  for (i=0;i<nb_entries;i++)
+  {
+    for (j=0;j<ndim;j++)
+    {
+      float delta = .01*(rand()/(float)RAND_MAX-.5);
+      codebook[i*ndim+j] += delta;
+      codebook[(i+nb_entries)*ndim+j] = codebook[i*ndim+j] - delta;
+    }
+  }
+}
+
+
+void update_weighted(float *data, float *weight, int nb_vectors, float *codebook, int nb_entries, int ndim)
+{
+  int i,j;
+  float count[MAX_ENTRIES][ndim];
+  int nearest[nb_vectors];
+
+  //fprintf(stderr, "weighted: %d %d\n", nb_entries, ndim);
+  for (i=0;i<nb_entries;i++)
+    for (j=0;j<ndim;j++)
+      count[i][j] = 0;
+
+  for (i=0;i<nb_vectors;i++)
+  {
+    nearest[i] = find_nearest_weighted(codebook, nb_entries, data+i*ndim, weight+i*ndim, ndim);
+  }
+  for (i=0;i<nb_entries*ndim;i++)
+    codebook[i] = 0;
+
+  for (i=0;i<nb_vectors;i++)
+  {
+    int n = nearest[i];
+    for (j=0;j<ndim;j++)
+    {
+      float w = sqrt(weight[i*ndim+j]);
+      count[n][j]+=w;
+      codebook[n*ndim+j] += w*data[i*ndim+j];
+    }
+  }
+
+  //float w2=0;
+  for (i=0;i<nb_entries;i++)
+  {
+    for (j=0;j<ndim;j++)
+      codebook[i*ndim+j] *= (1./count[i][j]);
+    //w2 += (count[i]/(float)nb_vectors)*(count[i]/(float)nb_vectors);
+  }
+  //fprintf(stderr, "%f / %d\n", 1./w2, nb_entries);
+}
+
+void vq_train_weighted(float *data, float *weight, int nb_vectors, float *codebook, int nb_entries, int ndim)
+{
+  int i, j, e;
+  e = 1;
+  for (j=0;j<ndim;j++)
+    codebook[j] = 0;
+  for (i=0;i<nb_vectors;i++)
+    for (j=0;j<ndim;j++)
+      codebook[j] += data[i*ndim+j];
+  for (j=0;j<ndim;j++)
+    codebook[j] *= (1./nb_vectors);
+
+
+  while (e< nb_entries)
+  {
+#if 1
+    split(codebook, e, ndim);
+    e<<=1;
+#else
+    split1(codebook, e, data, nb_vectors, ndim);
+    e++;
+#endif
+    fprintf(stderr, "%d\n", e);
+    for (j=0;j<10;j++)
+      update_weighted(data, weight, nb_vectors, codebook, e, ndim);
+  }
+}
+
+
+int main(int argc, char **argv)
+{
+  int i,j;
+  int nb_vectors, nb_entries, ndim;
+  float *data, *pred, *codebook, *codebook2, *codebook3;
+  float *weight, *weight2, *weight3;
+  float *delta;
+  double err[2] = {0, 0};
+  double werr[2] = {0, 0};
+  double wsum[2] = {0, 0};
+
+  ndim = atoi(argv[1]);
+  nb_vectors = atoi(argv[2]);
+  nb_entries = 1<<atoi(argv[3]);
+
+  data = malloc(nb_vectors*ndim*sizeof(*data));
+  weight = malloc(nb_vectors*ndim*sizeof(*weight));
+  weight2 = malloc(nb_vectors*ndim*sizeof(*weight2));
+  weight3 = malloc(nb_vectors*ndim*sizeof(*weight3));
+  pred = malloc(nb_vectors*ndim*sizeof(*pred));
+  codebook = malloc(nb_entries*ndim*sizeof(*codebook));
+  codebook2 = malloc(nb_entries*ndim*sizeof(*codebook2));
+  codebook3 = malloc(nb_entries*ndim*sizeof(*codebook3));
+
+  for (i=0;i<nb_vectors;i++)
+  {
+    if (feof(stdin))
+      break;
+    for (j=0;j<ndim;j++)
+    {
+      scanf("%f ", &data[i*ndim+j]);
+    }
+  }
+  nb_vectors = i;
+  VALGRIND_CHECK_MEM_IS_DEFINED(data, nb_entries*ndim);
+
+  for (i=0;i<nb_vectors;i++)
+  {
+    if (i==0)
+       compute_weights2(data+i*ndim, data+i*ndim, weight+i*ndim, ndim);
+    else
+       compute_weights2(data+i*ndim, data+(i-1)*ndim, weight+i*ndim, ndim);
+  }
+  for (i=0;i<ndim;i++)
+    pred[i] = data[i];
+  for (i=1;i<nb_vectors;i++)
+  {
+    for (j=0;j<ndim;j++)
+      pred[i*ndim+j] = data[i*ndim+j] - COEF[j]*data[(i-1)*ndim+j];
+  }
+
+  VALGRIND_CHECK_MEM_IS_DEFINED(pred, nb_entries*ndim);
+  vq_train_weighted(pred, weight, nb_vectors, codebook, nb_entries, ndim);
+  printf("%d %d\n", ndim, nb_entries);
+  for (i=0;i<nb_entries;i++)
+  {
+   for (j=0;j<ndim;j++)
+    {
+      printf("%f ", codebook[i*ndim+j]);
+    }
+    printf("\n");
+  }
+
+  delta = malloc(nb_vectors*ndim*sizeof(*data));
+  float xq[2] = {0,0};
+  for (i=0;i<nb_vectors;i++)
+  {
+    //int nearest = find_nearest_weighted(codebook, nb_entries, &pred[i*ndim], &weight[i*ndim], ndim);
+    quantize_ge(&data[i*ndim], codebook, nb_entries, xq, ndim);
+    //printf("%f %f\n", xq[0], xq[1]);
+    for (j=0;j<ndim;j++)
+    {
+      delta[i*ndim+j] = xq[j]-data[i*ndim+j];
+      err[j] += (delta[i*ndim+j])*(delta[i*ndim+j]);
+      werr[j] += weight[i*ndim+j]*(delta[i*ndim+j])*(delta[i*ndim+j]);
+      wsum[j] += weight[i*ndim+j];
+      //delta[i*ndim+j] = pred[i*ndim+j] - codebook[nearest*ndim+j];
+      //printf("%f ", delta[i*ndim+j]);
+      //err[j] += (delta[i*ndim+j])*(delta[i*ndim+j]);
+    }
+    //printf("\n");
+  }
+  fprintf(stderr, "GE RMS error: %f %f\n", sqrt(err[0]/nb_vectors), sqrt(err[1]/nb_vectors));
+  fprintf(stderr, "Weighted GE error: %f %f\n", sqrt(werr[0]/wsum[0]), sqrt(werr[1]/wsum[1]));
+
+  free(codebook);
+  free(codebook2);
+  free(codebook3);
+  free(weight);
+  free(weight2);
+  free(weight3);
+  free(delta);
+  return 0;
+}
diff --git a/misc/genlsp.c b/misc/genlsp.c
new file mode 100644
index 0000000..536651f
--- /dev/null
+++ b/misc/genlsp.c
@@ -0,0 +1,183 @@
+/*--------------------------------------------------------------------------*\
+
+	FILE........: genlsp.c
+	AUTHOR......: David Rowe
+	DATE CREATED: 23/2/95
+
+	This program generates a text file of LSP vectors from an input
+	speech file.
+
+\*--------------------------------------------------------------------------*/
+
+/*
+  Copyright (C) 2009 David Rowe
+
+  All rights reserved.
+
+  This program is free software; you can redistribute it and/or modify
+  it under the terms of the GNU Lesser General Public License version 2.1, as
+  published by the Free Software Foundation.  This program is
+  distributed in the hope that it will be useful, but WITHOUT ANY
+  WARRANTY; without even the implied warranty of MERCHANTABILITY or
+  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public
+  License for more details.
+
+  You should have received a copy of the GNU Lesser General Public License
+  along with this program; if not, see <http://www.gnu.org/licenses/>.
+*/
+
+#define P 	12	/* LP order					*/
+#define LSP_DELTA1 0.01 /* grid spacing for LSP root searches */
+#define NW	279	/* frame size in samples 			*/
+#define	N  	80 	/* frame to frame shift				*/
+#define THRESH	40.0	/* threshold energy/sample for frame inclusion 	*/
+#define PI         3.141592654	/* mathematical constant                */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include "lpc.h"	/* LPC analysis functions 			*/
+#include "lsp.h"	/* LSP encode/decode functions 			*/
+
+int switch_present(sw,argc,argv)
+  char sw[];     /* switch in string form */
+  int argc;      /* number of command line arguments */
+  char *argv[];  /* array of command line arguments in string form */
+{
+  int i;       /* loop variable */
+
+  for(i=1; i<argc; i++)
+    if (!strcmp(sw,argv[i]))
+      return(i);
+
+  return 0;
+}
+
+int main(int argc, char *argv[]) {
+    FILE   *fspc;	/* input file ptr for test database		*/
+    FILE   *flsp;	/* output text file of LSPs 			*/
+    short  buf[N];	/* input frame of speech samples 		*/
+    float  Sn[NW];	/* float input speech samples 			*/
+    float  ak[P+1];	/* LPCs for current frame 			*/
+    float  lsp[P];	/* LSPs for current frame 			*/
+    float  lsp_prev[P];	/* LSPs for previous frame 			*/
+    float  E;		/* frame energy 				*/
+    long   f;		/* number of frames                             */
+    long   af;		/* number frames with "active" speech 		*/
+    float  Eres;	/* LPC residual energy 				*/
+    int    i;
+    int    roots;
+    int    unstables;
+    int    lspd, log, lspdt;
+    float  diff;
+
+    /* Initialise ------------------------------------------------------*/
+
+    if (argc < 3) {
+	printf("usage: %s RawFile LSPTextFile [--lspd] [--log] [--lspdt] \n", argv[0]);
+	exit(1);
+    }
+
+    /* Open files */
+
+    fspc = fopen(argv[1],"rb");
+    if (fspc == NULL) {
+	printf("Error opening input SPC file: %s",argv[1]);
+	exit(1);
+    }
+
+    flsp = fopen(argv[2],"wt");
+    if (flsp == NULL) {
+	printf("Error opening output LSP file: %s",argv[2]);
+	exit(1);
+    }
+
+    lspd = switch_present("--lspd", argc, argv);
+    log = switch_present("--log", argc, argv);
+    lspdt = switch_present("--lspdt", argc, argv);
+
+    for(i=0; i<NW; i++)
+	Sn[i] = 0.0;
+
+    /* Read SPC file, and determine aks[] for each frame ------------------*/
+
+    f = af = 0;
+    unstables = 0;
+    while(fread(buf,sizeof(short),N,fspc) == N) {
+
+	for(i=0; i<NW-N; i++)
+	    Sn[i] = Sn[i+N];
+	E = 0.0;
+	for(i=0; i<N; i++) {
+	    Sn[i+NW-N] = buf[i];
+	    E += Sn[i]*Sn[i];
+	}
+
+	E = 0.0;
+	for(i=0; i<NW; i++) {
+	    E += Sn[i]*Sn[i];
+	}
+	E = 10.0*log10(E/NW);
+
+	/* If energy high enough, include this frame */
+
+	f++;
+	if (E > THRESH) {
+	    af++;
+	    printf("Active Frame: %ld  unstables: %d\n",af, unstables);
+
+	    find_aks(Sn, ak, NW, P, &Eres);
+	    roots = lpc_to_lsp(ak, P , lsp, 5, LSP_DELTA1);
+	    if (roots == P) {
+		if (lspd) {
+		    if (log) {
+			fprintf(flsp,"%f ",log10(lsp[0]));
+			for(i=1; i<P; i++) {
+			    diff = lsp[i]-lsp[i-1];
+			    if (diff < (PI/4000.0)*25.0) diff = (PI/4000.0)*25.0;
+			    fprintf(flsp,"%f ",log10(diff));
+			}
+		    }
+		    else {
+			fprintf(flsp,"%f ",lsp[0]);
+			for(i=1; i<P; i++)
+			    fprintf(flsp,"%f ",lsp[i]-lsp[i-1]);
+		    }
+
+		    fprintf(flsp,"\n");
+
+		}
+		else if (lspdt) {
+		    for(i=0; i<P; i++)
+			fprintf(flsp,"%f ",lsp[i]-lsp_prev[i]);
+		    fprintf(flsp,"\n");
+
+		}
+		else {
+		    if (log) {
+			for(i=0; i<P; i++)
+			    fprintf(flsp,"%f ",log10(lsp[i]));
+			fprintf(flsp,"\n");
+		    }
+		    else {
+			for(i=0; i<P; i++)
+			    fprintf(flsp,"%f ",lsp[i]);
+			fprintf(flsp,"\n");
+		    }
+
+		}
+		memcpy(lsp_prev, lsp, sizeof(lsp));
+	    }
+	    else
+		unstables++;
+	}
+    }
+
+    printf("%3.2f %% active frames\n", 100.0*(float)af/f);
+    fclose(fspc);
+    fclose(flsp);
+
+    return 0;
+}
+
diff --git a/misc/mksine.c b/misc/mksine.c
new file mode 100644
index 0000000..fb5f822
--- /dev/null
+++ b/misc/mksine.c
@@ -0,0 +1,54 @@
+/*
+  mksine.c
+  David Rowe
+  10 Nov 2010
+
+  Creates a file of sine wave samples.
+*/
+
+#include <assert.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <math.h>
+#include <errno.h>
+
+#define TWO_PI     6.283185307
+#define FS         8000.0
+
+int main(int argc, char *argv[]) {
+    FILE *f;
+    int   i,n;
+    float freq, length;
+    short *buf;
+    float amp = 1E4;
+
+    if (argc < 4) {
+	      printf("usage: %s outputFile frequencyHz lengthSecs [PeakAmp]\n", argv[0]);
+	      exit(1);
+    }
+
+    if (strcmp(argv[1], "-") == 0) {
+        f = stdout;
+    } else if ( (f = fopen(argv[1],"wb")) == NULL ) {
+	      fprintf(stderr, "Error opening output file: %s: %s.\n", argv[3], strerror(errno));
+	      exit(1);
+    }
+    freq = atof(argv[2]);
+    length = atof(argv[3]);
+    if (argc == 5) amp = atof(argv[4]);
+
+    n = length*FS;
+    buf = (short*)malloc(sizeof(short)*n);
+    assert(buf != NULL);
+
+    for(i=0; i<n; i++)
+	      buf[i] = amp*cos(freq*i*(TWO_PI/FS));
+
+    fwrite(buf, sizeof(short), n, f);
+
+    fclose(f);
+    free(buf);
+
+    return 0;
+}
diff --git a/misc/pre.c b/misc/pre.c
new file mode 100644
index 0000000..8016c3f
--- /dev/null
+++ b/misc/pre.c
@@ -0,0 +1,59 @@
+/*
+   pre.c
+   David Rowe
+   Sep 26 2012
+
+   Takes audio from a file, pre-emphasises, and sends to output file.
+*/
+
+#include <assert.h>
+#include <math.h>
+#include <errno.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include "lpc.h"
+
+#define N 80
+
+int main(int argc, char*argv[]) {
+    FILE  *fin, *fout;
+    short  buf[N];
+    float  Sn[N], Sn_pre[N];
+    float  pre_mem = 0.0;
+    int    i;
+
+    if (argc != 3) {
+	printf("usage: pre InputRawSpeechFile OutputRawSpeechFile\n");
+	printf("e.g    pre input.raw output.raw\n");
+	exit(1);
+    }
+
+    if (strcmp(argv[1], "-")  == 0) fin = stdin;
+    else if ( (fin = fopen(argv[1],"rb")) == NULL ) {
+	fprintf(stderr, "Error opening input speech file: %s: %s.\n",
+         argv[1], strerror(errno));
+	exit(1);
+    }
+
+    if (strcmp(argv[2], "-") == 0) fout = stdout;
+    else if ( (fout = fopen(argv[2],"wb")) == NULL ) {
+	fprintf(stderr, "Error opening output speech file: %s: %s.\n",
+         argv[2], strerror(errno));
+	exit(1);
+    }
+
+    while(fread(buf, sizeof(short), N, fin) == N) {
+	for(i=0; i<N; i++)
+	    Sn[i] = buf[i];
+	pre_emp(Sn_pre, Sn, &pre_mem, N);
+	for(i=0; i<N; i++)
+	    buf[i] = Sn_pre[i];
+	fwrite(buf, sizeof(short), N, fout);
+    }
+
+    fclose(fin);
+    fclose(fout);
+
+    return 0;
+}
diff --git a/misc/raw2h.c b/misc/raw2h.c
new file mode 100644
index 0000000..cb17f7c
--- /dev/null
+++ b/misc/raw2h.c
@@ -0,0 +1,45 @@
+/*
+  raw2h.c
+  David Rowe
+  10 April 2013
+
+  Converts a raw sound file to a C header file.  Used for generating arrays to
+  test Codec2 on embedded systems without disk I/O.
+*/
+
+#include <assert.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+
+int main(int argc, char *argv[]) {
+    FILE *fraw, *fheader;
+    int   i, samples, ret;
+    short sam;
+
+    if (argc != 5) {
+	printf("usage: %s inputRawFile outputHeaderFile arrayName samples\n", argv[0]);
+	exit(1);
+    }
+
+    fraw = fopen(argv[1] ,"rb");
+    assert(fraw != NULL);
+    fheader = fopen(argv[2],"wt");
+    assert(fheader != NULL);
+    samples = atoi(argv[4]);
+
+    fprintf(fheader, "short %s[] = {\n", argv[3]);
+    for(i=0; i<samples-1; i++) {
+	ret = fread(&sam, sizeof(short), 1, fraw);
+        assert(ret == 1);
+        fprintf(fheader, "%d,\n", sam);
+    }
+    ret = fread(&sam, sizeof(short), 1, fraw);
+    assert(ret == 1);
+    fprintf(fheader, "%d\n};\n", sam);
+
+    fclose(fraw);
+    fclose(fheader);
+
+    return 0;
+}
diff --git a/misc/speexnoisesup.c b/misc/speexnoisesup.c
new file mode 100644
index 0000000..7328235
--- /dev/null
+++ b/misc/speexnoisesup.c
@@ -0,0 +1,58 @@
+/*----------------------------------------------------------------------------*\
+
+  FILE....: speexnoisesup.c
+  AUTHOR..: David Rowe
+  CREATED.: Sun 22 June 2014
+
+  File I/O based test program for Speex pre-processor, used for
+  initial testing of Speech noise suppression.
+
+\*----------------------------------------------------------------------------*/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <stdint.h>
+#include <speex/speex_preprocess.h>
+
+#define    N          80
+#define    FS       8000
+
+int main(int argc, char *argv[]) {
+    FILE *fin, *fout;
+    short buf[N];
+    SpeexPreprocessState *st;
+
+    if (argc < 2) {
+	    printf("usage: %s InFile OutFile\n", argv[0]);
+	    exit(0);
+    }
+
+    if (strcmp(argv[1], "-")  == 0) fin = stdin;
+    else if ( (fin = fopen(argv[1],"rb")) == NULL ) {
+        fprintf(stderr, "Error opening %s\n", argv[1]);
+        exit(1);
+    }
+
+    if (strcmp(argv[2], "-") == 0) fout = stdout;
+    else if ((fout = fopen(argv[2],"wb")) == NULL) {
+	    fprintf(stderr, "Error opening %s\n", argv[2]);
+	    exit(1);
+    }
+
+    st = speex_preprocess_state_init(N, FS);
+
+    while(fread(buf, sizeof(short), N, fin) == N) {
+        speex_preprocess_run(st, buf);
+	    fwrite(buf, sizeof(short), N, fout);
+        if (fout == stdout) fflush(stdout);
+    }
+
+    speex_preprocess_state_destroy(st);
+
+    fclose(fin);
+    fclose(fout);
+
+    return 0;
+}
diff --git a/misc/tcodec2.c b/misc/tcodec2.c
new file mode 100644
index 0000000..ecb81d7
--- /dev/null
+++ b/misc/tcodec2.c
@@ -0,0 +1,220 @@
+/*---------------------------------------------------------------------------*\
+
+  FILE........: tcodec2.c
+  AUTHOR......: David Rowe
+  DATE CREATED: 24/8/10
+
+  Test program for codec2.c functions.
+
+\*---------------------------------------------------------------------------*/
+
+/*
+  Copyright (C) 2010 David Rowe
+
+  All rights reserved.
+
+  This program is free software; you can redistribute it and/or modify
+  it under the terms of the GNU Lesser General Public License version 2.1, as
+  published by the Free Software Foundation.  This program is
+  distributed in the hope that it will be useful, but WITHOUT ANY
+  WARRANTY; without even the implied warranty of MERCHANTABILITY or
+  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public
+  License for more details.
+
+  You should have received a copy of the GNU Lesser General Public License
+  along with this program; if not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <assert.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <math.h>
+#include "defines.h"
+#include "comp.h"
+#include "codec2.h"
+#include "quantise.h"
+#include "interp.h"
+
+/* CODEC2 struct copies from codec2.c to help with testing */
+
+struct CODEC2 {
+    int    mode;
+    float  w[M];	        /* time domain hamming window                */
+    COMP   W[FFT_ENC];	        /* DFT of w[]                                */
+    float  Pn[2*N];	        /* trapezoidal synthesis window              */
+    float  Sn[M];               /* input speech                              */
+    float  hpf_states[2];       /* high pass filter states                   */
+    void  *nlp;                 /* pitch predictor states                    */
+    float  Sn_[2*N];	        /* synthesised output speech                 */
+    float  ex_phase;            /* excitation model phase track              */
+    float  bg_est;              /* background noise estimate for post filter */
+    float  prev_Wo;             /* previous frame's pitch estimate           */
+    MODEL  prev_model;          /* previous frame's model parameters         */
+    float  prev_lsps_[LPC_ORD]; /* previous frame's LSPs                     */
+    float  prev_energy;         /* previous frame's LPC energy               */
+};
+
+void analyse_one_frame(struct CODEC2 *c2, MODEL *model, short speech[]);
+void synthesise_one_frame(struct CODEC2 *c2, short speech[], MODEL *model, float ak[]);
+
+int test1()
+{
+    FILE   *fin, *fout;
+    short   buf[N];
+    struct CODEC2 *c2;
+    MODEL   model;
+    float   ak[LPC_ORD+1];
+    float   lsps[LPC_ORD];
+
+    c2 = codec2_create(CODEC2_MODE_2400);
+
+    fin = fopen("../raw/hts1a.raw", "rb");
+    assert(fin != NULL);
+    fout = fopen("hts1a_test.raw", "wb");
+    assert(fout != NULL);
+
+    while(fread(buf, sizeof(short), N, fin) == N) {
+	analyse_one_frame(c2, &model, buf);
+	speech_to_uq_lsps(lsps, ak, c2->Sn, c2->w, LPC_ORD);
+	synthesise_one_frame(c2, buf, &model, ak);
+	fwrite(buf, sizeof(short), N, fout);
+    }
+
+    codec2_destroy(c2);
+
+    fclose(fin);
+    fclose(fout);
+
+    return 0;
+}
+
+int test2()
+{
+    FILE   *fin, *fout;
+    short   buf[2*N];
+    struct  CODEC2 *c2;
+    MODEL   model, model_interp;
+    float   ak[LPC_ORD+1];
+    int     voiced1, voiced2;
+    int     lsp_indexes[LPC_ORD];
+    int     energy_index;
+    int     Wo_index;
+    char    *bits;
+    int     nbit;
+    int     i;
+    float   lsps[LPC_ORD];
+    float   e;
+
+    c2 = codec2_create(CODEC2_MODE_2400);
+    bits = (char*)malloc(codec2_bits_per_frame(c2));
+    assert(bits != NULL);
+    fin = fopen("../raw/hts1a.raw", "rb");
+    assert(fin != NULL);
+    fout = fopen("hts1a_test.raw", "wb");
+    assert(fout != NULL);
+
+    while(fread(buf, sizeof(short), 2*N, fin) == 2*N) {
+	/* first 10ms analysis frame - we just want voicing */
+
+	analyse_one_frame(c2, &model, buf);
+	voiced1 = model.voiced;
+
+	/* second 10ms analysis frame */
+
+	analyse_one_frame(c2, &model, &buf[N]);
+	voiced2 = model.voiced;
+
+	Wo_index = encode_Wo(model.Wo);
+	e = speech_to_uq_lsps(lsps, ak, c2->Sn, c2->w, LPC_ORD);
+	encode_lsps_scalar(lsp_indexes, lsps, LPC_ORD);
+	energy_index = encode_energy(e);
+	nbit = 0;
+	pack((unsigned char*)bits, (unsigned *)&nbit, Wo_index, WO_BITS);
+	for(i=0; i<LPC_ORD; i++) {
+	    pack((unsigned char*)bits, (unsigned *)&nbit, lsp_indexes[i], lsp_bits(i));
+	}
+	pack((unsigned char*)bits, (unsigned *)&nbit, energy_index, E_BITS);
+	pack((unsigned char*)bits, (unsigned *)&nbit, voiced1, 1);
+	pack((unsigned char*)bits, (unsigned *)&nbit, voiced2, 1);
+
+	nbit = 0;
+	Wo_index = unpack((unsigned char*)bits, (unsigned *)&nbit, WO_BITS);
+	for(i=0; i<LPC_ORD; i++) {
+	    lsp_indexes[i] = unpack((unsigned char*)bits, (unsigned *)&nbit, lsp_bits(i));
+	}
+	energy_index = unpack((unsigned char*)bits, (unsigned *)&nbit, E_BITS);
+	voiced1 = unpack((unsigned char*)bits, (unsigned *)&nbit, 1);
+	voiced2 = unpack((unsigned char*)bits, (unsigned *)&nbit, 1);
+
+	model.Wo = decode_Wo(Wo_index);
+	model.L = PI/model.Wo;
+	decode_amplitudes(&model,
+			  ak,
+			  lsp_indexes,
+			  energy_index,
+			  lsps,
+			  &e);
+
+	model.voiced = voiced2;
+	model_interp.voiced = voiced1;
+	interpolate(&model_interp, &c2->prev_model, &model);
+
+	synthesise_one_frame(c2,  buf,     &model_interp, ak);
+	synthesise_one_frame(c2, &buf[N],  &model, ak);
+
+	memcpy(&c2->prev_model, &model, sizeof(MODEL));
+	fwrite(buf, sizeof(short), 2*N, fout);
+    }
+
+    free(bits);
+    codec2_destroy(c2);
+
+    fclose(fin);
+    fclose(fout);
+
+    return 0;
+}
+
+int test3()
+{
+    FILE   *fin, *fout, *fbits;
+    short   buf1[2*N];
+    short   buf2[2*N];
+    char   *bits;
+    struct CODEC2 *c2;
+
+    c2 = codec2_create(CODEC2_MODE_2400);
+    int numBits  = codec2_bits_per_frame(c2);
+    int numBytes = (numBits+7)>>3;
+
+    bits = (char*)malloc(numBytes);
+
+    fin = fopen("../raw/hts1a.raw", "rb");
+    assert(fin != NULL);
+    fout = fopen("hts1a_test.raw", "wb");
+    assert(fout != NULL);
+    fbits = fopen("hts1a_test3.bit", "wb");
+    assert(fout != NULL);
+
+    while(fread(buf1, sizeof(short), 2*N, fin) == 2*N) {
+	codec2_encode(c2, (void*)bits, buf1);
+	fwrite(bits, sizeof(char), numBytes, fbits);
+	codec2_decode(c2, buf2, (void*)bits);
+	fwrite(buf2, sizeof(short), numBytes, fout);
+    }
+
+    free(bits);
+    codec2_destroy(c2);
+
+    fclose(fin);
+    fclose(fout);
+    fclose(fbits);
+
+    return 0;
+}
+
+int main() {
+    test3();
+    return 0;
+}
diff --git a/misc/tdec.c b/misc/tdec.c
new file mode 100644
index 0000000..501edc6
--- /dev/null
+++ b/misc/tdec.c
@@ -0,0 +1,141 @@
+/*
+   tdec.c
+   David Rowe
+   Jan 2017
+
+   Trivial non filtered decimator for high ration sample rate conversion.
+
+   build: gcc tdec.c -o tdec -Wall -O2
+
+*/
+
+#include <assert.h>
+#include <getopt.h>
+#include <string.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <stdint.h>
+
+#define SIGNED_16BIT   0
+#define SIGNED_8BIT    1
+#define UNSIGNED_8BIT  2
+
+void freq_shift_complex_buf(short buf[], int n, int lo_i[], int lo_q[]);
+
+void display_help(void) {
+    fprintf(stderr, "\nusage: tdec inputRawFile OutputRawFile DecimationRatio [-c]\n");
+    fprintf(stderr, "\nUse - for stdin/stdout\n\n");
+    fprintf(stderr, "-c complex signed 16 bit input and output\n");
+    fprintf(stderr, "-d complex signed 8 bit input (e.g. HackRF), complex signed 16 bit output\n");
+    fprintf(stderr, "-e complex unsigned 8 bit input (e.g. RTL-SDR), complex signed 16 bit output\n");
+    fprintf(stderr, "-f -Fs/4 freq shift\n\n");
+}
+
+int main(int argc, char *argv[]) {
+    FILE       *fin, *fout;
+    short       dec;
+    int         lo_i[3], lo_q[3];
+
+    if (argc < 3) {
+	display_help();
+	exit(1);
+    }
+
+    if (strcmp(argv[1], "-") == 0) 
+        fin = stdin;
+    else
+        fin = fopen(argv[1], "rb");
+    assert(fin != NULL);
+
+    if (strcmp(argv[2], "-") == 0) 
+        fout = stdout;
+    else
+        fout = fopen(argv[2], "wb");
+    assert(fout != NULL);
+
+    dec = atoi(argv[3]);
+
+    int channels = 1;
+    int freq_shift = 0;
+    lo_i[0] = -1; lo_i[1] =  0;
+    lo_q[0] =  0; lo_q[1] = -1;
+    int opt;
+    int format = SIGNED_16BIT;
+    while ((opt = getopt(argc, argv, "cdef")) != -1) {
+        switch (opt) {
+        case 'c': channels = 2; break;
+        case 'd': channels = 2; format = SIGNED_8BIT; break;
+        case 'e': channels = 2; format = UNSIGNED_8BIT; break;
+        case 'f': freq_shift = 1; break;
+        default:
+            display_help();
+            exit(1);
+        }
+    }
+
+    if (format == SIGNED_16BIT) {
+        short buf[dec*channels];
+        while(fread(buf, sizeof(short)*channels, dec, fin) == dec) {
+            if (freq_shift)            
+                freq_shift_complex_buf(buf, dec*channels, lo_i, lo_q);
+            fwrite(buf, sizeof(short), channels, fout);
+        }
+    }
+    else {
+        uint8_t inbuf[dec*channels];
+        short   outbuf[dec*channels];
+        short   sam, i;
+        
+        while(fread(inbuf, sizeof(uint8_t)*channels, dec, fin) == dec) {
+            for (i=0; i<dec*channels; i++) {
+                assert((format == SIGNED_8BIT) || (format == UNSIGNED_8BIT));
+                if (format == SIGNED_8BIT) {
+                    sam = (short)inbuf[i];
+                    sam <<= 8;
+                } else {
+                    sam = (short)inbuf[i] - 127;
+                    sam <<= 8;
+                }
+                outbuf[i] = sam;
+            }
+            if (freq_shift)
+                freq_shift_complex_buf(outbuf, dec*channels, lo_i, lo_q);
+            fwrite(outbuf, sizeof(short), channels, fout);
+        }
+
+    }
+
+    fclose(fout);
+    fclose(fin);
+
+    return 0;
+}
+                
+
+void freq_shift_complex_buf(short buf[], int n, int lo_i[], int lo_q[]) {
+    int i;
+
+    for (i=0; i<n; i+=2) {
+        /* update local osc recursion */
+
+        lo_i[2] = -lo_i[0]; lo_q[2] = -lo_q[0];
+
+        /* freq shift down input samples */
+
+        int a = buf[i];
+        int b = buf[i+1];
+        int c = lo_i[2];
+        int d = -lo_q[2];  /* conj LO as down shift */
+
+        buf[i]   = a*c - b*d;
+        buf[i+1] = b*c + a*d;
+
+        //fprintf(stderr, "% d % d % 5d % 5d\n", lo_i[2], lo_q[2], buf[i], buf[i+1]);
+
+        /* update LO memory */
+
+        lo_i[0] = lo_i[1]; lo_i[1] = lo_i[2]; 
+        lo_q[0] = lo_q[1]; lo_q[1] = lo_q[2]; 
+    }
+}
+
diff --git a/misc/timpulse.c b/misc/timpulse.c
new file mode 100644
index 0000000..8227c0d
--- /dev/null
+++ b/misc/timpulse.c
@@ -0,0 +1,128 @@
+/*
+  timpulse.c
+  David Rowe Dec 2019
+
+  Generate a synthetic speech signal from a sum of sinusoids.  Generates a known test
+  signals for phaseNN and ampNN projects.
+*/
+
+#include <assert.h>
+#include <math.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <getopt.h>
+
+#define FS 8000
+
+int main(int argc, char *argv[]) {
+    short buf[FS] = {0};
+    float f0 = 60.0;
+    float n0 = 0.0;
+    int   Nsecs = 1;
+    int   randf0 = 0;
+    int   filter = 0;
+    int   rande = 0;
+    
+    int o = 0;
+    int opt_idx = 0;
+    while( o != -1 ) {
+        static struct option long_opts[] = {
+            {"help",   no_argument,       0, 'h'},
+            {"n0",     required_argument, 0, 'n'},
+            {"f0",     required_argument, 0, 'f'},
+            {"secs",   required_argument, 0, 's'},
+            {"randf0", no_argument, 0, 'r'},
+            {"rande",  required_argument, 0, 'e'},
+            {"filter", no_argument, 0, 'i'},
+            {0, 0, 0, 0}
+        };
+        
+        o = getopt_long(argc,argv,"hn:f:s:r",long_opts,&opt_idx);
+        
+        switch(o) {
+        case 'n':
+	    n0 = atof(optarg);
+            break;
+        case 'f':
+            f0 = atof(optarg);
+	    break;
+        case 's':
+            Nsecs = atoi(optarg);
+	    break;
+        case 'r':
+            randf0 = 1;
+	    break;
+        case 'i':
+            filter = 1;
+	    break;
+        case 'e':
+            rande = atoi(optarg);
+	    break;
+        case '?':
+        case 'h':
+	    fprintf(stderr,
+		    "usage: %s\n"
+		    "[--f0 f0Hz]          fixed F0\n" 
+                    "[--n0 samples]       time offset\n" 
+                    "[--secs Nsecs]       number of seconds to generate\n"
+	            "[--randf0]           choose a random F0 every second\n"
+	            "[--rande Ndiscrete]  choose a random frame energy every second, Ndiscrete values\n"
+		    "\n", argv[0]);
+	    exit(1);      
+	break;
+        }
+    }
+
+    int t = 0;
+    float A = 100.0;
+    
+    /* optionally filter with 2nd order system */
+    float alpha = 0.25*M_PI, gamma=0.99;
+    float a[2] = {-2.0*gamma*cos(alpha), gamma*gamma};
+    float mem[2] = {0};
+    
+    for (int j=0; j<Nsecs; j++) {
+	if (rande) {
+	    float AdB_min = 20.0*log10(100.0);
+	    float AdB_step = 6.0;
+	    float num_values = rande;
+
+	    // discrete RV between 0..1
+	    float r = (float)rand()/RAND_MAX;
+	    r = floor(r*num_values);
+	    
+	    float AdB = AdB_min + r*AdB_step;
+	    A = pow(10.0,AdB/20.0);
+	    fprintf(stderr, "r: %f AdB: %f A: %f\n", r, AdB, A);
+	}
+	if (randf0) {
+	    float pitch_period = FS/400.0 + (FS/80.0 - FS/400.0)*rand()/RAND_MAX;
+	    f0 = (float)FS/pitch_period;
+	    //fprintf(stderr, "P: %f f0: %f\n", pitch_period, f0);
+	}
+	float Wo = 2.0*M_PI*f0/FS;
+	int L = M_PI/Wo;
+	float e = 0.0;
+	for(int i=0; i<FS; i++) {
+	    buf[i] = 0;
+	    // 1/sqrt(L) term makes power constant across Wo
+	    for(int m=1; m<L; m++)
+		buf[i] += (A/sqrt(L))*cos(m*Wo*(t + n0));
+	    e += pow(buf[i], 2.0);
+	    t++;
+	}
+	//fprintf(stderr, "e (dB): %f\n", 10*log10(e));
+	if (filter) {
+	    for(int i=0; i<FS; i++) {
+		float x = (float)buf[i];
+		float y = (x - mem[0]*a[0] - mem[1]*a[1]);
+		mem[1] = mem[0]; mem[0] = y;
+		buf[i] = (short)y;
+	    }
+	}
+    
+	fwrite(buf, sizeof(short), FS, stdout);
+    }
+    
+    return 0;
+}
diff --git a/misc/tinterp.c b/misc/tinterp.c
new file mode 100644
index 0000000..34f2eff
--- /dev/null
+++ b/misc/tinterp.c
@@ -0,0 +1,151 @@
+/*---------------------------------------------------------------------------*\
+
+  FILE........: tinterp.c
+  AUTHOR......: David Rowe
+  DATE CREATED: 22/8/10
+
+  Tests interpolation functions.
+
+\*---------------------------------------------------------------------------*/
+
+/*
+  Copyright (C) 2010 David Rowe
+
+  All rights reserved.
+
+  This program is free software; you can redistribute it and/or modify
+  it under the terms of the GNU Lesser General Public License version 2.1, as
+  published by the Free Software Foundation.  This program is
+  distributed in the hope that it will be useful, but WITHOUT ANY
+  WARRANTY; without even the implied warranty of MERCHANTABILITY or
+  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public
+  License for more details.
+
+  You should have received a copy of the GNU Lesser General Public License
+  along with this program; if not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <assert.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <math.h>
+#include <ctype.h>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <unistd.h>
+
+#include "defines.h"
+#include "sine.h"
+#include "interp.h"
+
+void make_amp(MODEL *model, float f0, float cdB, float mdBHz)
+{
+    int   i;
+    float mdBrad = mdBHz*FS/TWO_PI;
+
+    model->Wo = f0*TWO_PI/FS;
+    model->L  = PI/model->Wo;
+    for(i=0; i<=model->L; i++)
+	model->A[i] = pow(10.0,(cdB + (float)i*model->Wo*mdBrad)/20.0);
+    model->voiced = 1;
+}
+
+void write_amp(char file[], MODEL *model)
+{
+    FILE  *f;
+    int    i;
+
+    f = fopen(file,"wt");
+    for(i=1; i<=model->L; i++)
+	fprintf(f, "%f\t%f\n", model->Wo*i, model->A[i]);
+    fclose(f);
+}
+
+const char *get_next_float(const char *s, float *num)
+{
+    const char *p = s;
+    char  tmp[MAX_STR];
+
+    while(*p && !isspace(*p))
+	p++;
+    assert((p-s) < (int)(sizeof(tmp)-1));
+    memcpy(tmp, s, p-s);
+    tmp[p-s] = 0;
+    *num = atof(tmp);
+
+    return p+1;
+}
+
+const char *get_next_int(const char *s, int *num)
+{
+    const char *p = s;
+    char  tmp[MAX_STR];
+
+    while(*p && !isspace(*p))
+	p++;
+    assert((p-s) < (int)(sizeof(tmp)-1));
+    memcpy(tmp, s, p-s);
+    tmp[p-s] = 0;
+    *num = atoi(tmp);
+
+    return p+1;
+}
+
+void load_amp(MODEL *model, const char * file, int frame)
+{
+    FILE *f;
+    int   i;
+    char  s[1024];
+    const char *ps;
+
+    f = fopen(file,"rt");
+    assert(f);
+
+    for(i=0; i<frame; i++)
+	ps = fgets(s, 1023, f);
+
+    /// can frame ever be 0? what if fgets fails?
+    ps = s;
+    ps = get_next_float(ps, &model->Wo);
+    ps = get_next_int(ps, &model->L);
+    for(i=1; i<=model->L; i++)
+	ps = get_next_float(ps, &model->A[i]);
+
+    fclose(f);
+}
+
+void load_or_make_amp(MODEL *model,
+                      const char * filename, int frame,
+                      float f0, float cdB, float mdBHz)
+{
+    struct stat buf;
+    int rc = stat(filename, &buf);
+    if (rc || !S_ISREG(buf.st_mode) || ((buf.st_mode & S_IRUSR) != S_IRUSR))
+    {
+        make_amp(model, f0, cdB, mdBHz);
+    }
+    else
+    {
+        load_amp(model, filename, frame);
+    }
+}
+int main() {
+    MODEL  prev, next, interp;
+
+    load_or_make_amp(&prev,
+                     "../src/hts1a_model.txt", 32,
+                     50.0, 60.0, 6E-3);
+    load_or_make_amp(&next,
+                     "../src/hts1a_model.txt", 34,
+                     50.0, 40.0, 6E-3);
+
+    interp.voiced = 1;
+    interpolate(&interp, &prev, &next);
+
+    write_amp("tinterp_prev.txt", &prev);
+    write_amp("tinterp_interp.txt", &interp);
+    write_amp("tinterp_next.txt", &next);
+
+    return 0;
+}
diff --git a/misc/tlininterp.c b/misc/tlininterp.c
new file mode 100644
index 0000000..7db6791
--- /dev/null
+++ b/misc/tlininterp.c
@@ -0,0 +1,153 @@
+/*
+   tlininterp.c
+   David Rowe
+   Jan 2017
+
+   Fast linear interpolator for high oversampling rates.  Upsample
+   with a decent filter first such that the signal is "low pass" wrt
+   to the input sample rate.
+
+   build: gcc tlininterp.c -o tlininterp -Wall -O2
+
+*/
+
+#include <assert.h>
+#include <getopt.h>
+#include <string.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <stdint.h>
+
+#define NBUF         1000
+#define SIGNED_16BIT    0
+#define SIGNED_8BIT     1
+
+void display_help(void) {
+    fprintf(stderr, "\nusage: tlininterp inputRawFile OutputRawFile OverSampleRatio [-c]\n");
+    fprintf(stderr, "\nUse - for stdin/stdout\n\n");
+    fprintf(stderr, "-c complex signed 16 bit input and output\n");
+    fprintf(stderr, "-d complex signed 16 bit input, complex signed 8 bit output\n");
+    fprintf(stderr, "-f +Fs/4 freq shift\n\n");
+}
+
+int main(int argc, char *argv[]) {
+    FILE       *fin, *fout;
+    short       left[2], right[2], out[2*NBUF], i;
+    float       oversample, t;
+    int8_t      out_s8[2*NBUF];
+    int         lo_i[3], lo_q[3];
+
+    if (argc < 3) {
+	display_help();
+	exit(1);
+    }
+
+    if (strcmp(argv[1], "-") == 0) 
+        fin = stdin;
+    else
+        fin = fopen(argv[1], "rb");
+    assert(fin != NULL);
+
+    if (strcmp(argv[2], "-") == 0) 
+        fout = stdout;
+    else
+        fout = fopen(argv[2], "wb");
+    assert(fout != NULL);
+
+    oversample = atof(argv[3]);
+    if (oversample <= 1) {
+ 	display_help();
+	exit(1);
+    }
+       
+    int channels = 1;
+    int freq_shift = 0;
+    lo_i[0] = -1; lo_i[1] =  0;
+    lo_q[0] =  0; lo_q[1] = -1;
+    int format = SIGNED_16BIT;
+    int opt;
+    while ((opt = getopt(argc, argv, "cdf")) != -1) {
+        switch (opt) {
+        case 'c': channels = 2; break;
+        case 'd': channels = 2; format = SIGNED_8BIT; break;
+        case 'f': freq_shift = 1; break;
+        default:
+            display_help();
+            exit(1);
+        }
+    }
+
+    for (i=0; i<channels; i++)
+        left[i] = 0;
+    t = 0.0;
+    int j = 0;
+    while(fread(&right, sizeof(short)*channels, 1, fin) == 1) {
+        while (t < 1.0) {
+
+            for (i=0; i<channels; i++) {
+                out[2*j+i] = (1.0 - t)*left[i] + t*right[i];
+            }
+
+            if (freq_shift) {
+
+                /* update local osc recursion */
+
+                lo_i[2] = -lo_i[0]; lo_q[2] = -lo_q[0];
+                
+                /* complex mixer to up-shift complex samples */
+
+                int a = out[2*j];
+                int b = out[2*j+1];
+                int c = lo_i[2];
+                int d = lo_q[2];
+
+                out[2*j]   = a*c - b*d;
+                out[2*j+1] = b*c + a*d;
+
+                //fprintf(stderr, "% d % d % 5d % 5d\n", lo_i[2], lo_q[2], out[0], out[1]);
+
+                /* update LO memory */
+
+                lo_i[0] = lo_i[1]; lo_i[1] = lo_i[2]; 
+                lo_q[0] = lo_q[1]; lo_q[1] = lo_q[2]; 
+            }
+
+            /* once we have enough samples write to disk */
+
+            j++;
+            if (j == NBUF) {
+                if (format == SIGNED_16BIT) {
+                    fwrite(&out, sizeof(short)*channels, NBUF, fout);
+                } else {
+                    for (i=0; i<channels*NBUF; i++) {
+                        out_s8[i] = out[i] >> 8;
+                    }
+                    fwrite(&out_s8, sizeof(int8_t)*channels, NBUF, fout);
+                }
+                j = 0;
+            }
+
+            t += 1.0/oversample;
+        }
+        
+        t -= 1.0;
+        for (i=0; i<channels; i++)
+            left[i] = right[i];
+    }
+
+    /* write remaining samples to disk */
+
+    if (format == SIGNED_16BIT) {
+        fwrite(&out, sizeof(short), j, fout);
+    } else {
+        for (i=0; i<j; i++) {
+            out_s8[i] = out[i] >> 8;
+        }
+        fwrite(&out_s8, sizeof(int8_t)*channels, j, fout);
+    }
+
+    fclose(fout);
+    fclose(fin);
+
+    return 0;
+}
diff --git a/misc/tnlp.c b/misc/tnlp.c
new file mode 100644
index 0000000..35e2ea4
--- /dev/null
+++ b/misc/tnlp.c
@@ -0,0 +1,164 @@
+/*---------------------------------------------------------------------------*\
+
+  FILE........: tnlp.c
+  AUTHOR......: David Rowe
+  DATE CREATED: 23/3/93
+
+  Test program for non linear pitch estimation functions.
+
+\*---------------------------------------------------------------------------*/
+
+/*
+  Copyright (C) 2009 David Rowe
+
+  All rights reserved.
+
+  This program is free software; you can redistribute it and/or modify
+  it under the terms of the GNU Lesser General Public License version 2.1, as
+  published by the Free Software Foundation.  This program is
+  distributed in the hope that it will be useful, but WITHOUT ANY
+  WARRANTY; without even the implied warranty of MERCHANTABILITY or
+  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public
+  License for more details.
+
+  You should have received a copy of the GNU Lesser General Public License
+  along with this program; if not, see <http://www.gnu.org/licenses/>.
+*/
+
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <math.h>
+
+#include "defines.h"
+#include "dump.h"
+#include "sine.h"
+#include "nlp.h"
+#include "kiss_fft.h"
+
+int   frames;
+
+/*---------------------------------------------------------------------------*\
+
+   switch_present()
+
+   Searches the command line arguments for a "switch".  If the switch is
+   found, returns the command line argument where it ws found, else returns
+   NULL.
+
+\*---------------------------------------------------------------------------*/
+
+int switch_present(sw,argc,argv)
+  char sw[];     /* switch in string form */
+  int argc;      /* number of command line arguments */
+  char *argv[];  /* array of command line arguments in string form */
+{
+  int i;       /* loop variable */
+
+  for(i=1; i<argc; i++)
+    if (!strcmp(sw,argv[i]))
+      return(i);
+
+  return 0;
+}
+
+/*---------------------------------------------------------------------------*\
+
+                                    MAIN
+
+\*---------------------------------------------------------------------------*/
+
+int main(int argc, char *argv[])
+{
+    if (argc < 3) {
+	printf("\nusage: tnlp InputRawSpeechFile Outputf0PitchTextFile "
+	       "[--dump DumpFile] [--Fs SampleRateHz]\n");
+        exit(1);
+    }
+    
+    int Fs = 8000;
+    if (switch_present("--Fs",argc,argv)) {
+        Fs = atoi(argv[argc+1]);
+    }
+
+    C2CONST c2const = c2const_create(Fs, N_S);
+    int   n = c2const.n_samp;
+    int   m = c2const.m_pitch;
+    FILE *fin,*fout;
+    short buf[n];
+    float Sn[m];	        /* float input speech samples */
+    kiss_fft_cfg  fft_fwd_cfg;
+    COMP  Sw[FFT_ENC];	        /* DFT of Sn[] */
+    float w[m];	                /* time domain hamming window */
+    float W[FFT_ENC];	        /* DFT of w[] */
+    float pitch_samples;
+    int   i;
+    float f0, prev_f0;
+    void  *nlp_states;
+    #ifdef DUMP
+    int   dump;
+    #endif
+
+    /* Input file */
+
+    if ((fin = fopen(argv[1],"rb")) == NULL) {
+      printf("Error opening input speech file: %s\n",argv[1]);
+      exit(1);
+    }
+
+    /* Output file */
+
+    if ((fout = fopen(argv[2],"wt")) == NULL) {
+      printf("Error opening output text file: %s\n",argv[2]);
+      exit(1);
+    }
+
+    #ifdef DUMP
+    dump = switch_present("--dump",argc,argv);
+    if (dump)
+      dump_on(argv[dump+1]);
+    #else
+    /// TODO
+    /// #warning "Compile with -DDUMP if you expect to dump anything."
+    #endif
+
+    for(i=0; i<m; i++) {
+      Sn[i] = 0.0;
+    }
+
+    nlp_states = nlp_create(&c2const);
+    fft_fwd_cfg = kiss_fft_alloc(FFT_ENC, 0, NULL, NULL);
+    make_analysis_window(&c2const, fft_fwd_cfg, w, W);
+
+    frames = 0;
+    prev_f0 = 1/P_MAX_S;
+    while(fread(buf, sizeof(short), n, fin)) {
+      /* Update input speech buffers */
+
+      for(i=0; i<m-n; i++)
+        Sn[i] = Sn[i+n];
+      for(i=0; i<n; i++)
+        Sn[i+m-n] = buf[i];
+      dft_speech(&c2const, fft_fwd_cfg, Sw, Sn, w);
+      #ifdef DUMP
+      dump_Sn(m, Sn); dump_Sw(Sw);
+      #endif
+
+      f0 = nlp(nlp_states, Sn, n, &pitch_samples, Sw, W, &prev_f0);
+
+      fprintf(stderr,"%d %f %f\n", frames++, f0, pitch_samples);
+      fprintf(fout,"%f %f\n", f0, pitch_samples);
+    }
+
+    fclose(fin);
+    fclose(fout);
+    #ifdef DUMP
+    if (dump) dump_off();
+    #endif
+    nlp_destroy(nlp_states);
+
+    return 0;
+}
+
+
diff --git a/misc/tprede.c b/misc/tprede.c
new file mode 100644
index 0000000..4d3d09c
--- /dev/null
+++ b/misc/tprede.c
@@ -0,0 +1,53 @@
+/*
+   tpre_de.c
+   David Rowe
+   Sep 24 2012
+
+   Unit test to generate the combined impulse response of pre & de-emphasis filters.
+
+     pl("../unittest/out48.raw",1,3000)
+     pl("../unittest/out8.raw",1,3000)
+
+   Listening to it also shows up anything nasty:
+
+     $ play -s -2 -r 48000 out48.raw
+     $ play -s -2 -r 8000 out8.raw
+
+  */
+
+#include <assert.h>
+#include <math.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include "lpc.h"
+
+#define N                        10
+#define F                        10
+
+int main() {
+    FILE  *fprede;
+    float  Sn[N], Sn_pre[N], Sn_de[N];
+    float  pre_mem = 0.0, de_mem = 0.0;
+    int    i, f;
+
+    fprede = fopen("prede.txt", "wt");
+    assert(fprede != NULL);
+
+    for(i=0; i<N; i++)
+	Sn[i] = 0.0;
+
+    Sn[0]= 1.0;
+
+    for(f=0; f<F; f++) {
+	pre_emp(Sn_pre, Sn, &pre_mem, N);
+	de_emp(Sn_de, Sn_pre, &de_mem, N);
+	for(i=0; i<N; i++) {
+	    fprintf(fprede, "%f\n", Sn_de[i]);
+	}
+	Sn[0] = 0.0;
+    }
+
+    fclose(fprede);
+
+    return 0;
+}
diff --git a/misc/tquant.c b/misc/tquant.c
new file mode 100644
index 0000000..2492123
--- /dev/null
+++ b/misc/tquant.c
@@ -0,0 +1,214 @@
+/*---------------------------------------------------------------------------*\
+
+  FILE........: tquant.c
+  AUTHOR......: David Rowe
+  DATE CREATED: 22/8/10
+
+  Generates quantisation curves for plotting on Octave.
+
+\*---------------------------------------------------------------------------*/
+
+/*
+  Copyright (C) 2010 David Rowe
+
+  All rights reserved.
+
+  This program is free software; you can redistribute it and/or modify
+  it under the terms of the GNU Lesser General Public License version 2.1, as
+  published by the Free Software Foundation.  This program is
+  distributed in the hope that it will be useful, but WITHOUT ANY
+  WARRANTY; without even the implied warranty of MERCHANTABILITY or
+  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public
+  License for more details.
+
+  You should have received a copy of the GNU Lesser General Public License
+  along with this program; if not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <assert.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <math.h>
+
+#include "defines.h"
+#include "dump.h"
+#include "quantise.h"
+
+int test_Wo_quant();
+int test_lsp_quant();
+int test_lsp(int lsp_number, int levels, float max_error_hz);
+int test_energy_quant(int levels, float max_error_dB);
+
+int main() {
+    quantise_init();
+    test_Wo_quant();
+    test_lsp_quant();
+    test_energy_quant(E_LEVELS, 0.5*(E_MAX_DB - E_MIN_DB)/E_LEVELS);
+
+    return 0;
+}
+
+int test_lsp_quant() {
+    test_lsp( 1, 16,  12.5);
+    test_lsp( 2, 16,  12.5);
+    test_lsp( 3, 16,  25);
+    test_lsp( 4, 16,  50);
+    test_lsp( 5, 16,  50);
+    test_lsp( 6, 16,  50);
+    test_lsp( 7, 16,  50);
+    test_lsp( 8,  8,  50);
+    test_lsp( 9,  8,  50);
+    test_lsp(10,  4, 100);
+
+    return 0;
+}
+
+int test_energy_quant(int levels, float max_error_dB) {
+    FILE  *fe;
+    float  e,e_dec, error, low_e, high_e;
+    int    index, index_in, index_out, i;
+
+    /* check 1:1 match between input and output levels */
+
+    for(i=0; i<levels; i++) {
+	index_in = i;
+	e = decode_energy(index_in, E_BITS);
+	index_out = encode_energy(e, E_BITS);
+	if (index_in != index_out) {
+	    printf("edB: %f index_in: %d index_out: %d\n",
+		   10.0*log10(e), index_in, index_out);
+	    exit(0);
+	}
+    }
+
+    /* check error over range of quantiser */
+
+    low_e = decode_energy(0, E_BITS);
+    high_e = decode_energy(levels-1, E_BITS);
+    fe = fopen("energy_err.txt", "wt");
+
+    for(e=low_e; e<high_e; e +=(high_e-low_e)/1000.0) {
+	index = encode_energy(e, E_BITS);
+	e_dec = decode_energy(index, E_BITS);
+	error = 10.0*log10(e) - 10.0*log10(e_dec);
+	fprintf(fe, "%f\n", error);
+	if (fabs(error) > max_error_dB) {
+	    printf("error: %f %f\n", error, max_error_dB);
+	    exit(0);
+	}
+    }
+
+    fclose(fe);
+    return 0;
+}
+
+int test_lsp(int lsp_number, int levels, float max_error_hz) {
+    float lsp[LPC_ORD];
+    int   indexes_in[LPC_ORD];
+    int   indexes_out[LPC_ORD];
+    int   indexes[LPC_ORD];
+    int   i;
+    float lowf, highf, f, error;
+    char  s[MAX_STR];
+    FILE *flsp;
+    float max_error_rads;
+
+    lsp_number--;
+    max_error_rads = max_error_hz*TWO_PI/FS;
+
+    for(i=0; i<LPC_ORD; i++)
+	indexes_in[i] = 0;
+
+    for(i=0; i<levels; i++) {
+	indexes_in[lsp_number] = i;
+	decode_lsps_scalar(lsp, indexes_in, LPC_ORD);
+	encode_lsps_scalar(indexes_out, lsp,LPC_ORD);
+	if (indexes_in[lsp_number] != indexes_out[lsp_number]) {
+	    printf("freq: %f index_in: %d index_out: %d\n",
+		   lsp[lsp_number]+1, indexes_in[lsp_number],
+		   indexes_out[lsp_number]);
+	    exit(0);
+	}
+    }
+
+    for(i=0; i<LPC_ORD; i++)
+	indexes[i] = 0;
+    indexes[lsp_number] = 0;
+    decode_lsps_scalar(lsp, indexes, LPC_ORD);
+    lowf = lsp[lsp_number];
+    indexes[lsp_number] = levels - 1;
+    decode_lsps_scalar(lsp, indexes, LPC_ORD);
+    highf = lsp[lsp_number];
+    sprintf(s,"lsp%d_err.txt", lsp_number+1);
+    flsp = fopen(s, "wt");
+
+    for(f=lowf; f<highf; f +=(highf-lowf)/1000.0) {
+	lsp[lsp_number] = f;
+	encode_lsps_scalar(indexes, lsp, LPC_ORD);
+	decode_lsps_scalar(lsp, indexes, LPC_ORD);
+	error = f - lsp[lsp_number];
+	fprintf(flsp, "%f\n", error);
+	if (fabs(error) > max_error_rads) {
+	    printf("%d error: %f %f\n", lsp_number+1, error, max_error_rads);
+	    exit(0);
+	}
+    }
+
+    fclose(flsp);
+
+    printf("OK\n");
+
+    return 0;
+}
+
+int test_Wo_quant() {
+    int    c;
+    FILE  *f;
+    float  Wo,Wo_dec, error, step_size;
+    int    index, index_in, index_out;
+
+    /* output Wo quant curve for plotting */
+
+    f = fopen("quant_pitch.txt","wt");
+
+    for(Wo=0.9*(TWO_PI/P_MAX); Wo<=1.1*(TWO_PI/P_MIN); Wo += 0.001) {
+	index = encode_Wo(Wo, WO_BITS);
+	fprintf(f, "%f %d\n", Wo, index);
+    }
+
+    fclose(f);
+
+    /* check for all Wo codes we get 1:1 match between encoder
+       and decoder Wo levels */
+
+    for(c=0; c<WO_LEVELS; c++) {
+	index_in = c;
+	Wo = decode_Wo(index_in, WO_BITS);
+        index_out = encode_Wo(Wo, WO_BITS);
+	if (index_in != index_out)
+	    printf("  Wo %f index_in %d index_out %d\n", Wo,
+		   index_in, index_out);
+    }
+
+    /* measure quantisation error stats and compare to expected.  Also
+       plot histogram of error file to check. */
+
+    f = fopen("quant_pitch_err.txt","wt");
+    step_size = ((TWO_PI/P_MIN) - (TWO_PI/P_MAX))/WO_LEVELS;
+
+    for(Wo=TWO_PI/P_MAX; Wo<0.99*TWO_PI/P_MIN; Wo += 0.0001) {
+	index = encode_Wo(Wo, WO_BITS);
+	Wo_dec = decode_Wo(index, WO_BITS);
+	error = Wo - Wo_dec;
+	if (fabs(error) > (step_size/2.0)) {
+	    printf("error: %f  step_size/2: %f\n", error, step_size/2.0);
+	    exit(0);
+	}
+	fprintf(f,"%f\n",error);
+    }
+    printf("OK\n");
+
+    fclose(f);
+    return 0;
+}
diff --git a/misc/tsrc.c b/misc/tsrc.c
new file mode 100644
index 0000000..6791b51
--- /dev/null
+++ b/misc/tsrc.c
@@ -0,0 +1,109 @@
+/*
+   tsrc.c
+   David Rowe
+   Sat Nov 3 2012
+
+   Unit test for libresample code.
+
+   build: gcc tsrc.c -o tsrc -lm -lsamplerate
+
+  */
+
+#include <assert.h>
+#include <math.h>
+#include <string.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <samplerate.h>
+#include <unistd.h>
+
+#define N    10000                   /* processing buffer size */
+
+void display_help(void) {
+    fprintf(stderr, "\nusage: tsrc inputRawFile OutputRawFile OutSampleRatio [-l] [-c]\n");
+    fprintf(stderr, "\nUse - for stdin/stdout\n\n");
+    fprintf(stderr, "-l fast linear resampler\n");
+    fprintf(stderr, "-c complex (two channel) resampling\n\n");
+}
+
+int main(int argc, char *argv[]) {
+    FILE       *fin, *fout;
+    short       in_short[N], out_short[N];
+    float       in[N], out[N];
+    SRC_STATE  *src;
+    SRC_DATA    data;
+    int         error, nin, nremaining, i;
+
+    if (argc < 3) {
+	display_help();
+	exit(1);
+    }
+
+    if (strcmp(argv[1], "-") == 0) 
+        fin = stdin;
+    else
+        fin = fopen(argv[1], "rb");
+    assert(fin != NULL);
+
+    if (strcmp(argv[2], "-") == 0) 
+        fout = stdout;
+    else
+        fout = fopen(argv[2], "wb");
+    assert(fout != NULL);
+
+    data.data_in = in;
+    data.data_out = out;
+    data.end_of_input = 0;
+    data.src_ratio = atof(argv[3]);
+
+    int channels = 1;
+    int resampler = SRC_SINC_FASTEST;
+    int opt;
+    while ((opt = getopt(argc, argv, "lc")) != -1) {
+        switch (opt) {
+        case 'l': resampler = SRC_LINEAR; break;
+        case 'c': channels = 2; break;
+        default:
+            display_help();
+            exit(1);
+        }
+    }
+
+    data.input_frames = N/channels;
+    data.output_frames = N/channels;
+
+    src = src_new(resampler, channels, &error);
+    assert(src != NULL);
+
+    int total_in = 0;
+    int total_out = 0;
+
+    nin = data.input_frames;
+    nremaining = 0;
+    while(fread(&in_short[nremaining*channels], sizeof(short)*channels, nin, fin) == nin) {
+	src_short_to_float_array(in_short, in, N);
+	error = src_process(src, &data);
+        assert(error == 0);
+	src_float_to_short_array(out, out_short, data.output_frames_gen*channels);
+
+	fwrite(out_short, sizeof(short), data.output_frames_gen*channels, fout);
+        if (fout == stdout) fflush(stdout);
+
+        nremaining = data.input_frames - data.input_frames_used;
+        nin = data.input_frames_used;
+	//fprintf(stderr, "input frames: %d output_frames %d nremaining: %d\n", 
+        //        (int)data.input_frames_used, (int)data.output_frames_gen, nremaining);
+        for(i=0; i<nremaining*channels; i++)
+            in_short[i] = in_short[i+nin*channels];
+
+        total_in  += data.input_frames_used;
+        total_out += data.output_frames_gen;
+    }
+
+    //fprintf(stderr, "total_in: %d total_out: %d\n", total_in, total_out);
+
+    fclose(fout);
+    fclose(fin);
+
+    return 0;
+}
diff --git a/misc/vq_binary_switch.c b/misc/vq_binary_switch.c
new file mode 100644
index 0000000..4270c3d
--- /dev/null
+++ b/misc/vq_binary_switch.c
@@ -0,0 +1,296 @@
+/*
+  vq_binary_switch.c
+  David Rowe Dec 2021
+
+  C implementation of [1], that re-arranges VQ indexes so they are robust to single
+  bit errors.
+
+  [1] Pseudo Gray Coding, Zeger & Gersho 1990 
+*/
+
+#include <assert.h>
+#include <getopt.h>
+#include <math.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <limits.h>
+#include "mbest.h"
+
+#define MAX_DIM     20
+#define MAX_ENTRIES 4096
+
+// equation (33) of [1], total cost of all hamming distance 1 vectors of vq index k
+float cost_of_distance_one(float *vq, int n, int dim, float *prob, int k, int st, int en, int verbose) {
+  int log2N = log2(n);
+  float c = 0.0;
+  for (int b=0; b<log2N; b++) {
+    unsigned int index_neighbour = k ^ (1<<b);
+    float dist = 0.0;
+    for(int i=st; i<=en; i++)
+      dist += pow(vq[k*dim+i] - vq[index_neighbour*dim+i], 2.0);
+    c += prob[k]*dist;
+    if (verbose)
+      printf("k: %d b: %d index_neighbour: %d dist: %f prob: %f c: %f \n", k, b, index_neighbour, dist, prob[k], c);
+  }
+  return c;
+}
+
+// equation (39) of [1]
+float distortion_of_current_mapping(float *vq, int n, int dim, float *prob, int st, int en) {
+  float d = 0.0;
+  for(int k=0; k<n; k++)
+    d += cost_of_distance_one(vq, n, dim, prob, k, st, en, 0);
+  return d;
+}
+ 
+// we sort the cost array c[], returning the indexes of sorted elements
+float c[MAX_ENTRIES];
+
+/* Note how the compare function compares the values of the
+ * array to be sorted. The passed value to this function
+ * by `qsort' are actually the `idx' array elements.
+ */
+int compare_increase (const void * a, const void * b) {
+    int aa = *((int *) a), bb = *((int *) b);
+    if (c[aa] < c[bb]) {
+        return 1;
+    } else if (c[aa] == c[bb]) {
+        return 0;
+    } else {
+        return -1;
+     }
+}
+
+void sort_c(int *idx, const size_t n) {
+    for (size_t i=0; i<n; i++) idx[i] = i;
+    qsort(idx, n, sizeof(int), compare_increase);
+}
+
+void swap(float *vq, int dim, float *prob, int index1, int index2) {
+  float tmp[dim];
+  for(int i=0; i<dim; i++) tmp[i] = vq[index1*dim+i];
+  for(int i=0; i<dim; i++) vq[index1*dim+i] = vq[index2*dim+i];
+  for(int i=0; i<dim; i++) vq[index2*dim+i] = tmp[i];
+
+  tmp[0] = prob[index1];
+  prob[index1] = prob[index2];
+  prob[index2] = tmp[0];
+}
+
+int main(int argc, char *argv[]) {
+    float vq[MAX_DIM*MAX_ENTRIES];
+    int   dim = MAX_DIM;
+    int   max_iter = INT_MAX;
+    int   st = -1;
+    int   en = -1;
+    int   verbose = 0;
+    int   n = 0;
+    int   fast_en = 0;
+    char  prob_fn[80]="";
+    
+    int o = 0; int opt_idx = 0;
+    while (o != -1) {
+       static struct option long_opts[] = {
+           {"prob",    required_argument, 0, 'p'},
+           {"st",      required_argument, 0, 't'},
+           {"en",      required_argument, 0, 'e'},
+	   {0, 0, 0, 0}
+        };
+        o = getopt_long(argc,argv,"hd:m:vt:e:n:fp:",long_opts,&opt_idx);
+        switch (o) {
+	case 'd':
+	    dim = atoi(optarg);
+	    assert(dim <= MAX_DIM);
+	    break; 
+        case 'm':
+            max_iter = atoi(optarg);
+            break;
+        case 't':
+            st = atoi(optarg);
+            break;
+        case 'e':
+            en = atoi(optarg);
+            break;
+        case 'f':
+	    fast_en = 1;
+            break;
+        case 'n':
+            n = atoi(optarg);
+            break;
+        case 'p':
+	    strcpy(prob_fn,optarg);
+            break;	    
+        case 'v':
+            verbose = 1;
+            break;
+	help:
+	    fprintf(stderr, "\n");
+            fprintf(stderr, "usage: %s -d dimension [-m max_iterations -v --st Kst --en Ken -n nVQ] vq_in.f32 vq_out.f32\n", argv[0]);
+	    fprintf(stderr, "\n");
+            fprintf(stderr, "-n nVQ           Run with just the first nVQ entries of the VQ\n");
+            fprintf(stderr, "--st Kst         Start vector element for error calculation (default 0)\n");
+            fprintf(stderr, "--en Ken         End vector element for error calculation (default K-1)\n");
+            fprintf(stderr, "--prob probFile  f32 file of probabilities for each VQ element (default 1.0)\n");
+            fprintf(stderr, "-v               verbose\n");
+            exit(1);
+        }
+    }
+
+    int dx = optind;
+    if ((argc - dx) < 2) {
+        fprintf(stderr, "Too few arguments\n");
+	goto help;
+    }
+    if (dim == 0) goto help;
+
+    /* default to measuring error on entire vector */
+    if (st == -1) st = 0; 
+    if (en == -1) en = dim-1;
+
+    /* load VQ quantiser file --------------------*/
+
+    fprintf(stderr, "loading %s ... ", argv[dx]);
+    FILE *fq=fopen(argv[dx], "rb");
+    if (fq == NULL) {
+      fprintf(stderr, "Couldn't open: %s\n", argv[dx]);
+      exit(1);
+    }
+
+    if (n==0) {
+      /* count how many entries m of dimension k are in this VQ file */
+      float dummy[dim];
+      while (fread(dummy, sizeof(float), dim, fq) == (size_t)dim)
+	n++;
+      assert(n <= MAX_ENTRIES);
+      fprintf(stderr, "%d entries of vectors width %d\n", n, dim);
+
+      rewind(fq);
+    }
+
+    /* load VQ into memory */
+    int nrd = fread(vq, sizeof(float), n*dim, fq);
+    assert(nrd == n*dim);
+    fclose(fq);
+   
+    /* set probability of each vector to 1.0 as default */
+    float prob[n];
+    for(int l=0; l<n; l++) prob[l] = 1.0;
+    if (strlen(prob_fn)) {
+      fprintf(stderr, "Reading probability file: %s\n", prob_fn);
+      FILE *fp = fopen(prob_fn,"rb");
+      assert(fp != NULL);
+      int nrd = fread(prob, sizeof(float), n, fp);
+      assert(nrd == n);
+      fclose(fp);
+      float sum = 0.0;
+      for(int l=0; l<n; l++) sum += prob[l];
+      fprintf(stderr, "sum = %f\n", sum);
+    }
+    
+    int iteration = 0;
+    int i = 0;
+    int finished = 0;
+    int switches = 0;
+    int log2N = log2(n);
+    float distortion0 = distortion_of_current_mapping(vq, n, dim, prob, st, en);
+    fprintf(stderr, "distortion0: %f\n", distortion0);
+
+    while(!finished) {
+
+      // generate a list A(i) of which vectors have the largest cost of bit errors
+      for(int k=0; k<n; k++) {
+	c[k] = cost_of_distance_one(vq, n, dim, prob, k, st, en, verbose);
+      }      
+      int A[n];
+      sort_c(A, n);
+      
+      // Try switching each vector with A(i)
+      float best_delta = 0; int best_j = 0;
+      for(int j=1; j<n; j++) {
+	float distortion1, distortion2, delta = 0.0;
+	
+	// we can't switch with ourself
+	if (j != A[i]) {
+	  if (fast_en) {
+	    // subtract just those contributions to delta that will change
+	    delta -= cost_of_distance_one(vq, n, dim, prob, A[i], st, en, verbose);
+	    delta -= cost_of_distance_one(vq, n, dim, prob, j, st, en, verbose);
+	    for (int b=0; b<log2N; b++) {
+	      unsigned int index_neighbour;
+	      index_neighbour = A[i] ^ (1<<b);
+	      if ((index_neighbour != j) && (index_neighbour != A[i]))
+		   delta -= cost_of_distance_one(vq, n, dim, prob, index_neighbour, st, en, verbose);
+	      index_neighbour = j ^ (1<<b);
+	      if ((index_neighbour != j) && (index_neighbour != A[i]))
+		   delta -= cost_of_distance_one(vq, n, dim, prob, index_neighbour, st, en, verbose);
+	    }
+	  }
+	  else
+	    distortion1 = distortion_of_current_mapping(vq, n, dim, prob, st, en);
+
+	  // switch vq entries A(i) and j
+	  swap(vq, dim, prob, A[i], j);
+
+	  if (fast_en) {
+	    // add just those contributions to delta that will change
+	    delta += cost_of_distance_one(vq, n, dim, prob, A[i], st, en, verbose);
+	    delta += cost_of_distance_one(vq, n, dim, prob, j, st, en, verbose);
+	    for (int b=0; b<log2N; b++) {
+	      unsigned int index_neighbour;
+	      index_neighbour = A[i] ^ (1<<b);
+	      if ((index_neighbour != j) && (index_neighbour != A[i]))
+		   delta += cost_of_distance_one(vq, n, dim, prob, index_neighbour, st, en, verbose);
+	      index_neighbour = j ^ (1<<b);
+	      if ((index_neighbour != j) && (index_neighbour != A[i]))
+		   delta += cost_of_distance_one(vq, n, dim, prob, index_neighbour, st, en, verbose);
+	    }
+	  }
+	  else {
+	    distortion2 = distortion_of_current_mapping(vq, n, dim, prob, st, en);
+	    delta = distortion2 - distortion1;
+	  }
+
+	  if (delta < 0.0) {
+	    if (fabs(delta) > best_delta) {
+	      best_delta = fabs(delta);
+	      best_j = j;
+	    }
+	  }
+	  // unswitch
+	  swap(vq, dim, prob, A[i], j);
+	}
+      } //next j
+
+      // printf("best_delta: %f best_j: %d\n", best_delta, best_j);
+      if (best_delta == 0.0) {
+	// Hmm, no improvement, lets try the next vector in the sorted cost list
+	if (i == n-1) finished = 1; else i++;
+      } else {
+	// OK keep the switch that minimised the distortion
+	swap(vq, dim, prob, A[i], best_j);
+	switches++;
+
+	// save results
+	FILE *fq=fopen(argv[dx+1], "wb");
+	if (fq == NULL) {
+	    fprintf(stderr, "Couldn't open: %s\n", argv[dx+1]);
+	    exit(1);
+	}
+	int nwr = fwrite(vq, sizeof(float), n*dim, fq);
+	assert(nwr == n*dim);
+	fclose(fq);
+	
+	// set up for next iteration
+	iteration++;
+	float distortion = distortion_of_current_mapping(vq, n, dim, prob, st, en);
+	fprintf(stderr, "it: %3d dist: %f %3.2f i: %3d sw: %3d\n", iteration, distortion,
+	                distortion/distortion0, i, switches);
+	if (iteration >= max_iter) finished = 1;
+	i = 0;
+      }
+    }
+    
+    return 0;
+}
+
diff --git a/misc/vq_mbest.c b/misc/vq_mbest.c
new file mode 100644
index 0000000..a247f61
--- /dev/null
+++ b/misc/vq_mbest.c
@@ -0,0 +1,302 @@
+/*
+  vq_mbest.c
+  David Rowe Dec 2019
+
+  Utility to perform a mbest VQ search on vectors from stdin, sending
+  quantised vectors to stdout.
+*/
+
+#include <assert.h>
+#include <getopt.h>
+#include <math.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <limits.h>
+#include "mbest.h"
+
+#define MAX_K       20
+#define MAX_ENTRIES 4096
+#define MAX_STAGES  5
+
+void quant_mbest(float vec_out[],
+                 int   indexes[],
+                 float vec_in[],
+                 int   num_stages,
+                 float vqw[], float vq[], 
+                 int   m[], int k,
+                 int   mbest_survivors);
+
+int verbose = 0;
+
+int main(int argc, char *argv[]) {
+    float vq[MAX_STAGES*MAX_K*MAX_ENTRIES];
+    float vqw[MAX_STAGES*MAX_K*MAX_ENTRIES];
+    int   m[MAX_STAGES];
+    int   k=0, mbest_survivors=1, num_stages=0;
+    char  fnames[256], fn[256], *comma, *p;
+    FILE *fq;
+    float lower = -1E32;
+    int   st = -1;
+    int   en = -1;
+    int   num = INT_MAX;
+    int   output_vec_usage = 0;
+    
+    int o = 0; int opt_idx = 0;
+    while (o != -1) {
+       static struct option long_opts[] = {
+	   {"k",         required_argument, 0, 'k'},
+	   {"quant",     required_argument, 0, 'q'},
+	   {"mbest",     required_argument, 0, 'm'},
+	   {"lower",     required_argument, 0, 'l'},
+	   {"verbose",   required_argument, 0, 'v'},
+           {"st",        required_argument, 0, 't'},
+           {"en",        required_argument, 0, 'e'},
+           {"num",       required_argument, 0, 'n'},
+           {"vec_usage", no_argument, 0, 'u'},
+	   {0, 0, 0, 0}
+        };
+
+        o = getopt_long(argc,argv,"hk:q:m:vt:e:n:u",long_opts,&opt_idx);
+        switch (o) {
+	case 'k':
+	    k = atoi(optarg);
+	    assert(k <= MAX_K);
+	    break; 
+	case 'q':
+	    /* load up list of comma delimited file names */
+            strcpy(fnames, optarg);
+            p = fnames;
+            num_stages = 0;
+            do {
+                assert(num_stages < MAX_STAGES);
+                strcpy(fn, p);
+                comma = strchr(fn, ',');
+                if (comma) {
+                    *comma = 0;
+                    p = comma+1;
+                }
+                /* load quantiser file */
+                fprintf(stderr, "stage: %d loading %s ... ", num_stages, fn);
+                fq=fopen(fn, "rb");
+                if (fq == NULL) {
+                    fprintf(stderr, "Couldn't open: %s\n", fn);
+                    exit(1);
+                }
+                /* count how many entries m of dimension k are in this VQ file */
+                m[num_stages] = 0;
+		float dummy[k];
+                while (fread(dummy, sizeof(float), k, fq) == (size_t)k)
+                    m[num_stages]++;
+                assert(m[num_stages] <= MAX_ENTRIES);
+                fprintf(stderr, "%d entries of vectors width %d\n", m[num_stages], k);
+                /* now load VQ into memory */
+                rewind(fq);                       
+                int rd = fread(&vq[num_stages*k*MAX_ENTRIES], sizeof(float), m[num_stages]*k, fq);
+                assert(rd == m[num_stages]*k);
+                num_stages++;
+                fclose(fq);
+            } while(comma);
+            break;
+        case 'm':
+            mbest_survivors = atoi(optarg);
+            fprintf(stderr, "mbest_survivors = %d\n",  mbest_survivors);
+            break;
+        case 'n':
+	    num = atoi(optarg);
+            break;
+        case 'l':
+            lower = atof(optarg);
+            break;
+        case 't':
+            st = atoi(optarg);
+            break;
+        case 'e':
+            en = atoi(optarg);
+            break;
+        case 'u':
+            output_vec_usage = 1;
+            break;
+	case 'v':
+            verbose = 1;
+            break;
+	help:
+	    fprintf(stderr, "\n");
+            fprintf(stderr, "usage: %s -k dimension -q vq1.f32,vq2.f32,.... [Options]\n", argv[0]);
+	    fprintf(stderr, "\n");
+	    fprintf(stderr, "input vectors on stdin, output quantised vectors on stdout\n");
+	    fprintf(stderr, "\n");
+	    fprintf(stderr, "--lower lowermeanLimit   Only count vectors with average above this level in distortion calculations\n");
+	    fprintf(stderr, "--mbest N                number of survivors at each stage, set to 0 for standard VQ search\n");
+            fprintf(stderr, "--st    Kst              start vector element for error calculation (default 0)\n");
+            fprintf(stderr, "--en    Ken              end vector element for error calculation (default K-1)\n");
+	    fprintf(stderr, "--num   numToProcess     number of vectors to quantise (default to EOF)\n");
+	    fprintf(stderr, "--vec_usage              Output a record of how many times each vector is used\n");
+ 	    fprintf(stderr, "-v                       Verbose\n");
+            exit(1);
+        }
+    }
+
+    if ((num_stages == 0) || (k == 0))
+	goto help;
+
+    /* default to measuring error on entire vector */
+    if (st == -1) st = 0; 
+    if (en == -1) en = k-1;
+
+    float w[k];
+    for(int i=0; i<st; i++)
+        w[i] = 0.0;
+    for(int i=st; i<=en; i++)
+        w[i] = 1.0;
+    for(int i=en+1; i<k; i++)
+        w[i] = 0.0;
+
+    /* apply weighting to codebook (rather than in search) */
+    memcpy(vqw, vq, sizeof(vq));
+    for(int s=0; s<num_stages; s++) {
+        mbest_precompute_weight(&vqw[s*k*MAX_ENTRIES], w, k, m[s]);
+    }
+    
+    int indexes[num_stages], nvecs = 0; int vec_usage[m[0]];
+    for(int i=0; i<m[0]; i++) vec_usage[i] = 0;
+    float target[k], quantised[k];
+    float sqe = 0.0;
+    while(fread(&target, sizeof(float), k, stdin) && (nvecs < num)) {
+	for(int i=0; i<k; i++)
+            target[i] *= w[i];
+	int dont_count = 0;
+	/* optional clamping to lower limit or mean */
+	float mean = 0.0;
+	for(int i=0; i<k; i++)
+	    mean += target[i];
+	mean /= k;         
+	float difference = mean - lower;
+	if (difference < 0.0) {
+	    /* bring target up to lower clamping limit */
+	    for(int i=0; i<k; i++)
+		target[i] += -difference;
+	    dont_count = 1;
+	}
+	quant_mbest(quantised, indexes, target, num_stages, vqw, vq, m, k, mbest_survivors);
+	if (dont_count == 0) {
+	    for(int i=st; i<=en; i++)
+		sqe += pow(target[i]-quantised[i], 2.0);
+	}
+	fwrite(&quantised, sizeof(float), k, stdout);
+	nvecs++;
+	// count number f time each vector is used (just for first stage)
+	vec_usage[indexes[0]]++;
+    }
+
+    fprintf(stderr, "MSE: %4.2f\n", sqe/(nvecs*(en-st+1)));
+
+    if (output_vec_usage) {
+      for(int i=0; i<m[0]; i++)
+	fprintf(stderr, "%d\n", vec_usage[i]);
+    }
+
+    return 0;
+}
+
+// print vector debug function
+
+void pv(char s[], float v[], int k) {
+    int i;
+    if (verbose) {
+        fprintf(stderr, "%s",s);
+        for(i=0; i<k; i++)
+            fprintf(stderr, "%4.2f ", v[i]);
+        fprintf(stderr, "\n");
+    }
+}
+
+// mbest algorithm version, backported from LPCNet/src
+
+void quant_mbest(float vec_out[],
+                 int   indexes[],
+                 float vec_in[],
+                 int   num_stages,
+                 float vqw[], float vq[],
+                 int   m[], int k,
+                 int   mbest_survivors)
+{
+    float err[k], se1;
+    int i,j,s,s1,ind;
+    
+    struct MBEST *mbest_stage[num_stages];
+    int index[num_stages];
+    float target[k];
+    
+    for(i=0; i<num_stages; i++) {
+        mbest_stage[i] = mbest_create(mbest_survivors);
+        index[i] = 0;
+    }
+
+    se1 = 0.0;
+    for(i=0; i<k; i++) {
+        err[i] = vec_in[i];
+        se1 += err[i]*err[i];
+    }
+    se1 /= k;
+                
+    /* now quantise err[] using multi-stage mbest search, preserving
+       mbest_survivors at each stage */
+    
+    mbest_search(vqw, err, k, m[0], mbest_stage[0], index);
+    if (verbose) mbest_print("Stage 1:", mbest_stage[0]);
+    
+    for(s=1; s<num_stages; s++) {
+
+        /* for each candidate in previous stage, try to find best vector in next stage */
+        for (j=0; j<mbest_survivors; j++) {
+            /* indexes that lead us this far */
+            for(s1=0; s1<s; s1++) {
+                index[s1+1] = mbest_stage[s-1]->list[j].index[s1];
+            }
+            /* target is residual err[] vector given path to this candidate */
+            for(i=0; i<k; i++)
+                target[i] = err[i];
+            for(s1=0; s1<s; s1++) {
+                ind = index[s-s1];
+                if (verbose) fprintf(stderr, "   s: %d s1: %d s-s1: %d ind: %d\n", s,s1,s-s1,ind);
+                for(i=0; i<k; i++) {
+                    target[i] -= vqw[s1*k*MAX_ENTRIES+ind*k+i];
+                }
+            }
+            pv("   target: ", target, k);
+            mbest_search(&vqw[s*k*MAX_ENTRIES], target, k, m[s], mbest_stage[s], index);
+        }
+        char str[80]; sprintf(str,"Stage %d:", s+1);
+        if (verbose) mbest_print(str, mbest_stage[s]);
+    }
+
+    for(s=0; s<num_stages; s++) {
+        indexes[s] = mbest_stage[num_stages-1]->list[0].index[num_stages-1-s];
+    }
+
+    /* OK put it all back together using best survivor */
+    for(i=0; i<k; i++)
+        vec_out[i] = 0.0;
+    for(s=0; s<num_stages; s++) {
+        int ind = indexes[s];
+        float se2 = 0.0;
+        for(i=0; i<k; i++) {
+            err[i] -= vqw[s*k*MAX_ENTRIES+ind*k+i];
+            vec_out[i] += vq[s*k*MAX_ENTRIES+ind*k+i];
+            se2 += err[i]*err[i];
+        }
+        se2 /= k;
+        pv("    err: ", err, k);
+        if (verbose) fprintf(stderr, "    se2: %f\n", se2);
+    }
+    pv("  vec_out: ",vec_out, k);
+
+    pv("\n  vec_in: ", vec_in, k);
+    pv("  vec_out: ", vec_out, k);
+    pv("    err: ", err, k);
+    if (verbose) fprintf(stderr, "    se1: %f\n", se1);
+
+    for(i=0; i<num_stages; i++)
+        mbest_destroy(mbest_stage[i]);
+}
diff --git a/misc/vqtrain.c b/misc/vqtrain.c
new file mode 100644
index 0000000..59e88f1
--- /dev/null
+++ b/misc/vqtrain.c
@@ -0,0 +1,399 @@
+/*--------------------------------------------------------------------------*\
+
+	FILE........: VQTRAIN.C
+	AUTHOR......: David Rowe
+	DATE CREATED: 23/2/95
+
+	This program trains vector quantisers using K dimensional Lloyd-Max
+	method.
+
+\*--------------------------------------------------------------------------*/
+
+/*
+  Copyright (C) 2009 David Rowe
+
+  All rights reserved.
+
+  This program is free software; you can redistribute it and/or modify
+  it under the terms of the GNU Lesser General Public License version 2.1, as
+  published by the Free Software Foundation.  This program is
+  distributed in the hope that it will be useful, but WITHOUT ANY
+  WARRANTY; without even the implied warranty of MERCHANTABILITY or
+  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public
+  License for more details.
+
+  You should have received a copy of the GNU Lesser General Public License
+  along with this program; if not, see <http://www.gnu.org/licenses/>.
+*/
+
+/*-----------------------------------------------------------------------*\
+
+				INCLUDES
+
+\*-----------------------------------------------------------------------*/
+
+#include <assert.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <ctype.h>
+#include <getopt.h>
+#include <time.h>
+
+/*-----------------------------------------------------------------------*\
+
+				DEFINES
+
+\*-----------------------------------------------------------------------*/
+
+#define	DELTAQ 	0.005		/* quitting distortion			*/
+#define	MAX_STR	80		/* maximum string length		*/
+
+/*-----------------------------------------------------------------------*\
+
+			FUNCTION PROTOTYPES
+
+\*-----------------------------------------------------------------------*/
+
+void zero(float v[], int k);
+void acc(float v1[], float v2[], int k);
+void norm(float v[], int k, long n);
+long quantise(float cb[], float vec[], int k, int m, int st, int en, float *beste, float *se);
+
+/*-----------------------------------------------------------------------* \
+
+				MAIN
+
+\*-----------------------------------------------------------------------*/
+
+int main(int argc, char *argv[]) {
+    long   k,m;		/* dimension and codebook size			*/
+    float  *vec;	/* current vector 				*/
+    float  *cb;		/* vector codebook				*/
+    float  *cent;	/* centroids for each codebook entry		*/
+    long   *n;		/* number of vectors in this interval		*/
+    long   J;		/* number of vectors in training set		*/
+    long   i,j;
+    long   ind;	     	/* index of current vector			*/
+    float  e;           /* squared error for current vector             */
+    float  se;		/* squared error for this iteration		*/
+    float  var,var_1;	/* current and previous iterations distortion	*/
+    float  delta;	/* improvement in distortion 			*/
+    long   noutliers[3];/* number of vectors quantisers with > 3*sd     */
+    FILE   *ftrain;	/* file containing training set			*/
+    FILE   *fvq;	/* file containing vector quantiser		*/
+    int     ret;
+    float   deltaq_stop = DELTAQ;
+    FILE   *fres = NULL;
+    int     st = -1;
+    int     en = -1;
+    int     init_rand = 0;
+    
+    int o = 0;
+    int opt_idx = 0;
+    while( o != -1 ) {
+        static struct option long_opts[] = {
+            {"help",     no_argument,       0, 'h'},
+            {"residual", required_argument, 0, 'r'},
+            {"stop",     required_argument, 0, 's'},
+            {"st",       required_argument, 0, 't'},
+            {"en",       required_argument, 0, 'e'},
+            {"rand",     no_argument,       0, 'i'},
+            {0, 0, 0, 0}
+        };
+        
+        o = getopt_long(argc,argv,"hr:s:t:e:",long_opts,&opt_idx);
+        
+        switch(o) {
+        case 'r':
+            fres = fopen(optarg,"wb"); assert(fres != NULL);
+            //fprintf(stderr, "writing res to : %s \n", optarg);
+            break;
+        case 's':
+            deltaq_stop = atof(optarg);
+            //fprintf(stderr, "deltaq_stop :%f\n", deltaq_stop);
+            break;
+        case 't':
+            st = atoi(optarg);
+            break;
+        case 'e':
+            en = atoi(optarg);
+            break;
+        case 'i':
+            init_rand = 1;
+            break;
+        case 'h':
+        case '?':
+            goto helpmsg;
+            break;
+        }
+    }
+    int dx = optind;
+
+    if ((argc - dx) < 4) {
+        fprintf(stderr, "Too few arguments\n");
+    helpmsg:
+        fprintf(stderr, "usage: %s [Options] TrainFile.f32 K(dimension) M(codebook size) VQFile.f32\n", argv[0]);
+        fprintf(stderr, "  -r --residual VQResidualErrorFile.f32\n");
+        fprintf(stderr, "  -s --stop StopDelta\n");
+        fprintf(stderr, "     --st   Kst        start vector element for error calculation (default 0)\n");
+        fprintf(stderr, "     --en   Ken        end vector element for error calculation (default K-1)\n");
+        fprintf(stderr, "     --rand            use random sampling for initial VQ population\n");
+        exit(1);
+    }
+
+    /* Open training file */
+    ftrain = fopen(argv[dx],"rb");
+    if (ftrain == NULL) {
+	printf("Error opening training database file: %s\n",argv[dx]);
+	exit(1);
+    }
+
+    /* determine k and m, and allocate arrays */
+    k = atol(argv[dx+1]);
+    m = atol(argv[dx+2]);
+
+    /* default to measuring error on entire vector */
+    if (st == -1) st = 0; 
+    if (en == -1) en = k-1;
+
+    printf("vector dimension K=%ld  codebook size M=%ld ", k, m);
+    vec = (float*)malloc(sizeof(float)*k);
+    cb = (float*)malloc(sizeof(float)*k*m);
+    cent = (float*)malloc(sizeof(float)*k*m);
+    n = (long*)malloc(sizeof(long)*m);
+    if (vec == NULL || cb == NULL || cent == NULL || n == NULL) {
+	printf("Error in malloc.\n");
+	exit(1);
+    }
+
+    /* determine size of training set */
+    J = 0; zero(cent, k);
+    while(fread(vec, sizeof(float), k, ftrain) == (size_t)k) {
+        J++;
+        acc(cent, vec, k);
+    }
+    printf("J=%ld vectors in training set\n", J);
+
+    /* Lets measure 0 bit VQ (i.e. mean of training set) as starting point */   
+    norm(cent, k, J);
+    memcpy(cb, cent, k*sizeof(float));
+    se = 0.0;
+    rewind(ftrain);
+    for(i=0; i<J; i++) {
+        ret = fread(vec, sizeof(float), k, ftrain);
+        assert(ret == k);
+        quantise(cb, vec, k, 1, st, en, &e, &se);
+    }
+    var = se/(J*(en-st+1));
+    printf("\r  It: 0, var: %f sd: %f\n", var, sqrt(var));
+
+    /* set up initial codebook state from samples of training set */
+    if (init_rand) srand(time(NULL));
+    for(i=0; i<m; i++) {
+        if (init_rand)
+            j = J*(float)rand()/RAND_MAX;
+        else
+            j = i*(J/m);
+        fseek(ftrain, j*k*sizeof(float), SEEK_SET);
+        ret = fread(&cb[i*k], sizeof(float), k, ftrain);
+        assert(ret == k);
+    }
+
+    /* main loop */
+    j = 1;
+    do {
+	var_1 = var;
+
+	/* zero centroids */
+	for(i=0; i<m; i++) {
+	    zero(&cent[i*k], k);
+	    n[i] = 0;
+	}
+
+	/* quantise training set */
+	se = 0.0; noutliers[0] = noutliers[1] = noutliers[2] = 0;
+	rewind(ftrain);
+	for(i=0; i<J; i++) {
+	    ret = fread(vec, sizeof(float), k, ftrain);
+            assert(ret == k);
+	    ind = quantise(cb, vec, k, m, st, en, &e, &se);
+	    n[ind]++;
+	    acc(&cent[ind*k], vec, k);
+            if (sqrt(e/(en-st+1)) > 1.0) noutliers[0]++;
+            if (sqrt(e/(en-st+1)) > 2.0) noutliers[1]++;
+            if (sqrt(e/(en-st+1)) > 3.0) noutliers[2]++;
+	}
+	var = se/(J*(en-st+1));
+	delta = (var_1-var)/var;
+        int n_min = J;
+        int n_max = 0;
+	for(i=0; i<m; i++) {
+            if (n[i] < n_min) n_min = n[i];
+            if (n[i] > n_max) n_max = n[i];
+        }
+	printf("\r  It: %2ld, var: %5f sd: %f outliers > 1/2/3 dB = %3.2f/%3.2f/%3.2f Delta = %5.4f %d %d\n", j, var, sqrt(var),
+               (float)noutliers[0]/J, (float)noutliers[1]/J, (float)noutliers[2]/J, delta, n_min, n_max);
+	j++;
+
+	/* determine new codebook from centroids */
+	if (delta > deltaq_stop)
+	    for(i=0; i<m; i++) {
+		if (n[i] != 0) {
+		    norm(&cent[i*k], k, n[i]);
+		    memcpy(&cb[i*k], &cent[i*k], k*sizeof(float));
+		}
+	    }
+
+    } while (delta > deltaq_stop);
+
+    /* save VQ to disk */
+    fvq = fopen(argv[dx+3],"wt");
+    if (fvq == NULL) {
+	printf("Error opening VQ file: %s\n",argv[dx+3]);
+	exit(1);
+    }
+
+    fwrite(cb, sizeof(float), m*k, fvq);
+    
+    /* optionally output residual error for next stage VQ */
+    if (fres != NULL) {
+        float res[k];
+	rewind(ftrain);
+	for(j=0; j<J; j++) {
+	    ret = fread(vec, sizeof(float), k, ftrain);
+	    ind = quantise(cb, vec, k, m, st, en, &e, &se);
+            for(i=0; i<k; i++) {
+                res[i] = vec[i] - cb[k*ind+i];
+            }
+            fwrite(res, sizeof(float), k, fres);
+	}
+        fclose(fres);
+    }
+    
+    fclose(fvq);
+    fclose(ftrain);
+    free(vec);
+    free(n);
+
+    return 0;
+}
+
+/*-----------------------------------------------------------------------*\
+
+				FUNCTIONS
+
+\*-----------------------------------------------------------------------*/
+
+/*-----------------------------------------------------------------------*\
+
+	FUNCTION....: zero()
+
+	AUTHOR......: David Rowe
+	DATE CREATED: 23/2/95
+
+	Zeros a vector of length k.
+
+\*-----------------------------------------------------------------------*/
+
+void zero(float v[], int k)
+/*  float  v[];		ptr to start of vector		*/
+/*  int    k;		lngth of vector			*/
+{
+    int	i;
+
+    for(i=0; i<k; i++)
+	v[i] = 0.0;
+}
+
+/*---------------------------------------------------------------------------*\
+
+	FUNCTION....: acc()
+
+	AUTHOR......: David Rowe
+	DATE CREATED: 23/2/95
+
+	Adds k dimensional vectors v1 to v2 and stores the result back in v1.
+
+\*---------------------------------------------------------------------------*/
+
+void acc(float v1[], float v2[], int k)
+/*  float  v1[];	ptr to start of vector to accumulate	*/
+/*  float  v2[];	ptr to start of vector to add		*/
+/*  int	   k;		dimension of vectors			*/
+{
+    int	   i;
+
+    for(i=0; i<k; i++)
+	v1[i] += v2[i];
+}
+
+/*---------------------------------------------------------------------------*\
+
+	FUNCTION....: norm()
+
+	AUTHOR......: David Rowe
+	DATE CREATED: 23/2/95
+
+	Divides each element in k dimensional vector v by n.
+
+\*---------------------------------------------------------------------------*/
+
+void norm(float v[], int k, long n)
+/*  float  v[];		ptr to start of vector		*/
+/*  int	   k;		dimension of vectors		*/
+/*  int	   n;		normalising factor		*/
+{
+    int	   i;
+
+    assert(n != 0);
+    for(i=0; i<k; i++)
+	v[i] /= n;
+}
+
+/*---------------------------------------------------------------------------*\
+
+	FUNCTION....: quantise()
+	AUTHOR......: David Rowe
+	DATE CREATED: 23/2/95
+
+	Quantises vec by choosing the nearest vector in codebook cb, and
+	returns the vector index.  The squared error of the quantised vector
+	is added to se.
+
+\*---------------------------------------------------------------------------*/
+
+long quantise(float  cb[],   // current VQ codebook
+              float  vec[],  // cb[][k];	current VQ codebook
+              int    k,      // dimension of vectors
+              int    m,      // size of codebook
+              int    st,     // start sample of vector to use in error calc
+              int    en,     // end sample of vector to use in error calc
+              float *beste,  // current squared error
+              float *se)     // accumulated squared error
+{
+   long	   besti;	/* best index so far		*/
+   long	   j;
+   int     i;
+   float   diff,e;
+
+   besti = 0;
+   *beste = 1E32;
+   for(j=0; j<m; j++) {
+	e = 0.0;
+	for(i=st; i<=en; i++) {
+	    diff = cb[j*k+i]-vec[i];
+	    e += diff*diff;
+	}
+	if (e < *beste) {
+	    *beste = e;
+	    besti = j;
+	}
+   }
+
+   *se += *beste;
+   
+   return(besti);
+}
+
-- 
cgit v1.2.3