diff options
Diffstat (limited to 'misc')
| -rw-r--r-- | misc/16_8_short.c | 47 | ||||
| -rw-r--r-- | misc/CMakeLists.txt | 35 | ||||
| -rw-r--r-- | misc/de.c | 59 | ||||
| -rw-r--r-- | misc/extract.c | 117 | ||||
| -rw-r--r-- | misc/ge_train.c | 306 | ||||
| -rw-r--r-- | misc/genlsp.c | 183 | ||||
| -rw-r--r-- | misc/mksine.c | 54 | ||||
| -rw-r--r-- | misc/pre.c | 59 | ||||
| -rw-r--r-- | misc/raw2h.c | 45 | ||||
| -rw-r--r-- | misc/speexnoisesup.c | 58 | ||||
| -rw-r--r-- | misc/tcodec2.c | 220 | ||||
| -rw-r--r-- | misc/tdec.c | 141 | ||||
| -rw-r--r-- | misc/timpulse.c | 128 | ||||
| -rw-r--r-- | misc/tinterp.c | 151 | ||||
| -rw-r--r-- | misc/tlininterp.c | 153 | ||||
| -rw-r--r-- | misc/tnlp.c | 164 | ||||
| -rw-r--r-- | misc/tprede.c | 53 | ||||
| -rw-r--r-- | misc/tquant.c | 214 | ||||
| -rw-r--r-- | misc/tsrc.c | 109 | ||||
| -rw-r--r-- | misc/vq_binary_switch.c | 296 | ||||
| -rw-r--r-- | misc/vq_mbest.c | 302 | ||||
| -rw-r--r-- | misc/vqtrain.c | 399 |
22 files changed, 0 insertions, 3293 deletions
diff --git a/misc/16_8_short.c b/misc/16_8_short.c deleted file mode 100644 index a23722c..0000000 --- a/misc/16_8_short.c +++ /dev/null @@ -1,47 +0,0 @@ -/* - 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 deleted file mode 100644 index 55728b4..0000000 --- a/misc/CMakeLists.txt +++ /dev/null @@ -1,35 +0,0 @@ -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 deleted file mode 100644 index 9726743..0000000 --- a/misc/de.c +++ /dev/null @@ -1,59 +0,0 @@ -/* - 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 deleted file mode 100644 index b4f38db..0000000 --- a/misc/extract.c +++ /dev/null @@ -1,117 +0,0 @@ -/* - 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 deleted file mode 100644 index db786fc..0000000 --- a/misc/ge_train.c +++ /dev/null @@ -1,306 +0,0 @@ -/* - 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 deleted file mode 100644 index 536651f..0000000 --- a/misc/genlsp.c +++ /dev/null @@ -1,183 +0,0 @@ -/*--------------------------------------------------------------------------*\ - - 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 deleted file mode 100644 index fb5f822..0000000 --- a/misc/mksine.c +++ /dev/null @@ -1,54 +0,0 @@ -/* - 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 deleted file mode 100644 index 8016c3f..0000000 --- a/misc/pre.c +++ /dev/null @@ -1,59 +0,0 @@ -/* - 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 deleted file mode 100644 index cb17f7c..0000000 --- a/misc/raw2h.c +++ /dev/null @@ -1,45 +0,0 @@ -/* - 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 deleted file mode 100644 index 7328235..0000000 --- a/misc/speexnoisesup.c +++ /dev/null @@ -1,58 +0,0 @@ -/*----------------------------------------------------------------------------*\ - - 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 deleted file mode 100644 index ecb81d7..0000000 --- a/misc/tcodec2.c +++ /dev/null @@ -1,220 +0,0 @@ -/*---------------------------------------------------------------------------*\ - - 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 deleted file mode 100644 index 501edc6..0000000 --- a/misc/tdec.c +++ /dev/null @@ -1,141 +0,0 @@ -/* - 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 deleted file mode 100644 index 8227c0d..0000000 --- a/misc/timpulse.c +++ /dev/null @@ -1,128 +0,0 @@ -/* - 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 deleted file mode 100644 index 34f2eff..0000000 --- a/misc/tinterp.c +++ /dev/null @@ -1,151 +0,0 @@ -/*---------------------------------------------------------------------------*\ - - 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 deleted file mode 100644 index 7db6791..0000000 --- a/misc/tlininterp.c +++ /dev/null @@ -1,153 +0,0 @@ -/* - 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 deleted file mode 100644 index 35e2ea4..0000000 --- a/misc/tnlp.c +++ /dev/null @@ -1,164 +0,0 @@ -/*---------------------------------------------------------------------------*\ - - 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 deleted file mode 100644 index 4d3d09c..0000000 --- a/misc/tprede.c +++ /dev/null @@ -1,53 +0,0 @@ -/* - 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 deleted file mode 100644 index 2492123..0000000 --- a/misc/tquant.c +++ /dev/null @@ -1,214 +0,0 @@ -/*---------------------------------------------------------------------------*\ - - 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 deleted file mode 100644 index 6791b51..0000000 --- a/misc/tsrc.c +++ /dev/null @@ -1,109 +0,0 @@ -/* - 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 deleted file mode 100644 index 4270c3d..0000000 --- a/misc/vq_binary_switch.c +++ /dev/null @@ -1,296 +0,0 @@ -/* - 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 deleted file mode 100644 index a247f61..0000000 --- a/misc/vq_mbest.c +++ /dev/null @@ -1,302 +0,0 @@ -/* - 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 deleted file mode 100644 index 59e88f1..0000000 --- a/misc/vqtrain.c +++ /dev/null @@ -1,399 +0,0 @@ -/*--------------------------------------------------------------------------*\ - - 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(¢[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(¢[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(¢[i*k], k, n[i]); - memcpy(&cb[i*k], ¢[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); -} - |