aboutsummaryrefslogtreecommitdiff
path: root/misc
diff options
context:
space:
mode:
Diffstat (limited to 'misc')
-rw-r--r--misc/16_8_short.c47
-rw-r--r--misc/CMakeLists.txt35
-rw-r--r--misc/de.c59
-rw-r--r--misc/extract.c117
-rw-r--r--misc/ge_train.c306
-rw-r--r--misc/genlsp.c183
-rw-r--r--misc/mksine.c54
-rw-r--r--misc/pre.c59
-rw-r--r--misc/raw2h.c45
-rw-r--r--misc/speexnoisesup.c58
-rw-r--r--misc/tcodec2.c220
-rw-r--r--misc/tdec.c141
-rw-r--r--misc/timpulse.c128
-rw-r--r--misc/tinterp.c151
-rw-r--r--misc/tlininterp.c153
-rw-r--r--misc/tnlp.c164
-rw-r--r--misc/tprede.c53
-rw-r--r--misc/tquant.c214
-rw-r--r--misc/tsrc.c109
-rw-r--r--misc/vq_binary_switch.c296
-rw-r--r--misc/vq_mbest.c302
-rw-r--r--misc/vqtrain.c399
22 files changed, 3293 insertions, 0 deletions
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);
+}
+
generated by cgit v1.2.3 (git 2.46.0) at 2026-01-27 05:35:10 +0000