From 120435acd51f007cf5cd5a8da9af92ff01b3ded7 Mon Sep 17 00:00:00 2001 From: drowe67 Date: Thu, 27 Aug 2009 09:52:09 +0000 Subject: [PATCH] progressing scalar and VQ quantisation git-svn-id: https://svn.code.sf.net/p/freetel/code@33 01035d8c-6547-0410-b346-abe4f91aad63 --- codec2/README.txt | 10 +- codec2/TODO.txt | 3 +- codec2/octave/plamp.m | 40 ++++- codec2/script/menu.sh | 3 + codec2/src/Makefile | 6 +- codec2/src/dump.c | 17 +++ codec2/src/dump.h | 1 + codec2/src/quantise.c | 89 +++++++++++- codec2/src/sinedec.c | 10 +- codec2/unittest/Makefile | 11 +- codec2/unittest/extract.c | 121 ++++++++++++++++ codec2/unittest/genlsp.c | 123 ++++++++++++++++ codec2/unittest/vqtrain.c | 297 ++++++++++++++++++++++++++++++++++++++ 13 files changed, 719 insertions(+), 12 deletions(-) create mode 100644 codec2/unittest/extract.c create mode 100644 codec2/unittest/genlsp.c create mode 100644 codec2/unittest/vqtrain.c diff --git a/codec2/README.txt b/codec2/README.txt index 47b74261..a62ee9a4 100644 --- a/codec2/README.txt +++ b/codec2/README.txt @@ -41,12 +41,18 @@ LPC Modelling $ ./sinedec ../raw/hts1a.raw hts1a.mdl --lpc 10 - hts1a_lpc10.raw -Disucss why LPC modelling works so well when Am recovered via RMS method +Discuss why LPC modelling works so well when Am recovered via RMS method (Section 5.1 of thesis). Equal area model of LPC spectra versus harmonic? Seems to work remarkably well, especially compared to sampling. SNRs up to 30dB on female frames. -Octave Scripts +m=1 harmonic problem for males when LPC modelled. The amplitude of this harmonic +comes up by as much as 30dB after LP modelling as (I think) LPC spectra must +have zero derivative at DC. This means it's poor at modelling very low freq +harmonics which unfortunately ear is very sensitive to. Consider automatic +lowering for 20dB of this harmonic or maybe an extra few bits to quantise error. + +Octave Scripts -------------- pl.m - plot a segment from a raw file diff --git a/codec2/TODO.txt b/codec2/TODO.txt index 93f4140d..f7c0cbbb 100644 --- a/codec2/TODO.txt +++ b/codec2/TODO.txt @@ -10,4 +10,5 @@ TODO for codec2 [X] GPL2 notice in each file [ ] Replace Numerical Recipes in C (NRC) four1.c and four1.h with Gnu Science Lib (GSL) SL FFT as NRC code has restrictive licencing -[ ] change speech files from raw 16 bit files to wave files \ No newline at end of file +[ ] A way to handle m=1 harmonic for males when LPC modelling +[ ] Is BW expansion and Rk noise floor required before LSP quant diff --git a/codec2/octave/plamp.m b/codec2/octave/plamp.m index daaed16b..b7b479d0 100644 --- a/codec2/octave/plamp.m +++ b/codec2/octave/plamp.m @@ -16,6 +16,10 @@ function plamp(samname, f) modelq = load(modelq_name); pw_name = strcat(samname,"_pw.txt"); Pw = load(pw_name); + lsp_name = strcat(samname,"_lsp.txt"); + if (file_in_path(".",lsp_name)) + lsp = load(lsp_name); + endif do figure(1); @@ -39,11 +43,17 @@ function plamp(samname, f) noise = (Am-Amq) * (Am-Amq)'; snr = 10*log10(signal/noise); Am_err_label = sprintf(";Am_err SNR %4.2f dB;",snr); - %plot((1:L)*Wo*4000/pi, 20*log10(Amq) - 20*log10(Am),";Am_err;" ); plot((1:L)*Wo*4000/pi, 20*log10(Amq) - 20*log10(Am), Am_err_label); + if (file_in_path(".",lsp_name)) + for l=1:10 + plot([lsp(f,l)*4000/pi lsp(f,l)*4000/pi], [60 80], 'r'); + endfor + endif hold off; - printf("\rframe: %d menu: n-next b-back q-quit ", f); + % interactive menu + + printf("\rframe: %d menu: n-next b-back p-png q-quit ", f); fflush(stdout); k = kbhit(); if (k == 'n') @@ -52,6 +62,32 @@ function plamp(samname, f) if (k == 'b') f = f - 1; endif + + % optional print to PNG + + if (k == 'p') + + pngname = sprintf("%s_%d",samname,f); + + % small image + + __gnuplot_set__ terminal png size 420,300 + ss = sprintf("__gnuplot_set__ output \"%s.png\"", pngname); + eval(ss) + replot; + + % larger image + + __gnuplot_set__ terminal png size 800,600 + ss = sprintf("__gnuplot_set__ output \"%s_large.png\"", pngname); + eval(ss) + replot; + + % for some reason I need this to stop large plot getting wiped + __gnuplot_set__ output "/dev/null" + + endif + until (k == 'q') printf("\n"); diff --git a/codec2/script/menu.sh b/codec2/script/menu.sh index ce9a89d9..78f0fd08 100755 --- a/codec2/script/menu.sh +++ b/codec2/script/menu.sh @@ -59,6 +59,9 @@ do stty cbreak # or stty raw readchar=`dd if=/dev/tty bs=1 count=1 2>/dev/null` stty -cbreak + if [ $readchar == 'q' ] ; then + readchar=0 + fi if [ $readchar -ne 0 ] ; then play -f s -r 8000 -s w ${file[$readchar]} $dsp > /dev/null fi diff --git a/codec2/src/Makefile b/codec2/src/Makefile index b889ad04..6acf8144 100644 --- a/codec2/src/Makefile +++ b/codec2/src/Makefile @@ -1,9 +1,11 @@ CC=gcc -CFLAGS=-g -Wall +CFLAGS=-g -Wall -I. -I../src -I../speex -Wall -g -DFLOATING_POINT -DVAR_ARRAYS SINENC_OBJ = sinenc.o globals.o initenc.o four1.o refine.o spec.o SINEDEC_OBJ = sinedec.o globals.o initenc.o initdec.o four1.o synth.o \ - quantise.o lpc.o dump.o refine.o + quantise.o lpc.o dump.o refine.o ../speex/lsp.o \ + ../speex/quant_lsp.o ../speex/bits.o ../speex/lsp_tables_nb.o \ + ../speex/high_lsp_tables.c all: sinenc sinedec diff --git a/codec2/src/dump.c b/codec2/src/dump.c index 5a824676..55698f40 100644 --- a/codec2/src/dump.c +++ b/codec2/src/dump.c @@ -36,6 +36,7 @@ static FILE *fsw; static FILE *fmodel; static FILE *fqmodel; static FILE *fpw; +static FILE *flsp; void dump_on(char prefix[]) { char s[MAX_STR]; @@ -61,13 +62,19 @@ void dump_on(char prefix[]) { sprintf(s,"%s_pw.txt", prefix); fpw = fopen(s, "wt"); assert(fpw != NULL); + + sprintf(s,"%s_lsp.txt", prefix); + flsp = fopen(s, "wt"); + assert(flsp != NULL); } void dump_off(){ fclose(fsn); + fclose(fsw); fclose(fmodel); fclose(fqmodel); fclose(fpw); + fclose(flsp); } void dump_Sn(float Sn[]) { @@ -133,4 +140,14 @@ void dump_Pw(COMP Pw[]) { fprintf(fpw,"\n"); } +void dump_lsp(float lsp[]) { + int i; + + if (!dumpon) return; + + for(i=0; i<10; i++) + fprintf(flsp,"%f\t",lsp[i]); + fprintf(flsp,"\n"); +} + diff --git a/codec2/src/dump.h b/codec2/src/dump.h index 8e4ccc61..61300c28 100644 --- a/codec2/src/dump.h +++ b/codec2/src/dump.h @@ -36,5 +36,6 @@ void dump_Sw(COMP Sw[]); void dump_model(MODEL *m); void dump_quantised_model(MODEL *m); void dump_Pw(COMP Pw[]); +void dump_lsp(float lsp[]); #endif diff --git a/codec2/src/quantise.c b/codec2/src/quantise.c index 2a7e380f..b4a5b49e 100644 --- a/codec2/src/quantise.c +++ b/codec2/src/quantise.c @@ -24,13 +24,88 @@ Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */ +#include #include "sine.h" #include "quantise.h" #include "lpc.h" #include "dump.h" +#include +#include +#include #define MAX_ORDER 20 +#define LPC_FLOOR 0.0002 /* autocorrelation floor */ +#define LSP_DELTA1 0.2 /* grid spacing for LSP root searches */ + +/* Speex lag window */ + +const float lag_window[11] = { + 1.00000, 0.99716, 0.98869, 0.97474, 0.95554, 0.93140, 0.90273, 0.86998, + 0.83367, 0.79434, 0.75258 +}; + +/*---------------------------------------------------------------------------*\ + + quantise_uniform + + Simulates uniform quantising of a float. + +\*---------------------------------------------------------------------------*/ + +void quantise_uniform(float *val, float min, float max, int bits) +{ + int levels = 1 << (bits-1); + float norm; + int index; + + /* hard limit to quantiser range */ + + printf("min: %f max: %f val: %f ", min, max, val[0]); + if (val[0] < min) val[0] = min; + if (val[0] > max) val[0] = max; + + norm = (*val - min)/(max-min); + printf("%f norm: %f ", val[0], norm); + index = fabs(levels*norm + 0.5); + + *val = min + index*(max-min)/levels; + + printf("index %d val_: %f\n", index, val[0]); +} + +/*---------------------------------------------------------------------------*\ + + lsp_quantise + + Differential lsp quantiser + +\*---------------------------------------------------------------------------*/ + +void lsp_quantise( + float lsp[], + float lsp_[], + int order +) +{ + int i; + float dlsp[MAX_ORDER]; + float dlsp_[MAX_ORDER]; + + dlsp[0] = lsp[0]; + for(i=1; i {Am} LPC decode */ return snr; diff --git a/codec2/src/sinedec.c b/codec2/src/sinedec.c index bd013436..aa4f937a 100644 --- a/codec2/src/sinedec.c +++ b/codec2/src/sinedec.c @@ -76,6 +76,7 @@ int main(int argc, char *argv[]) float sum_snr; int lpc_model, order; + int lsp, lsp_quantiser; int dump; if (argc < 3) { @@ -135,9 +136,12 @@ int main(int argc, char *argv[]) } dump = switch_present("--dump",argc,argv); - if (dump) { + if (dump) dump_on(argv[dump+1]); - } + + lsp = switch_present("--lsp",argc,argv); + if (lsp) + lsp_quantiser = atoi(argv[lsp+1]); /* Initialise ------------------------------------------------------------*/ @@ -168,7 +172,7 @@ int main(int argc, char *argv[]) /* optional LPC model amplitudes */ if (lpc_model) { - snr = lpc_model_amplitudes(Sn, &model, order, 0); + snr = lpc_model_amplitudes(Sn, &model, order, lsp_quantiser); sum_snr += snr; } diff --git a/codec2/unittest/Makefile b/codec2/unittest/Makefile index eb9d67e4..0dd6d98b 100644 --- a/codec2/unittest/Makefile +++ b/codec2/unittest/Makefile @@ -1,6 +1,6 @@ CFLAGS = -I. -I../src -I../speex -Wall -g -DFLOATING_POINT -DVAR_ARRAYS -all: genres lsptest +all: genres lsptest genlsp extract vqtrain genres: genres.o ../src/lpc.o gcc $(CFLAGS) -o genres genres.o ../src/lpc.o -lm @@ -12,5 +12,14 @@ LSP_TEST_OBJ = lsptest.o ../src/lpc.o ../speex/lsp.o sd.o ../src/four1.o \ lsptest: $(LSP_TEST_OBJ) gcc $(CFLAGS) -o lsptest $(LSP_TEST_OBJ) -lm +genlsp: genlsp.o ../src/lpc.o ../speex/lsp.o + gcc $(CFLAGS) -o genlsp genlsp.o ../src/lpc.o ../speex/lsp.o -lm + +extract: extract.o + gcc $(CFLAGS) -o extract extract.o -lm + +vqtrain: vqtrain.o + gcc $(CFLAGS) -o vqtrain vqtrain.o -lm + %.o : %.c $(CC) -c $(CFLAGS) $< -o $@ diff --git a/codec2/unittest/extract.c b/codec2/unittest/extract.c new file mode 100644 index 00000000..59c6267a --- /dev/null +++ b/codec2/unittest/extract.c @@ -0,0 +1,121 @@ +/*--------------------------------------------------------------------------*\ + + FILE........: extract.c + AUTHOR......: David Rowe + DATE CREATED: 23/2/95 + + This program extracts a float file of vectors from a text file + of vectors. The float files are easier to process quickly + during VQ training. A subset of the text file VQ may be + extracted to faciltate split VQ of scaler VQ design. + +\*--------------------------------------------------------------------------*/ + +/* + 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 General Public License version 2, 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 General Public License + along with this program; if not, write to the Free Software + Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. +*/ + +#define MAX_STR 256 /* maximum string length */ + +#include +#include +#include + +void scan_line(FILE *fp, float f[], int n); + +int main(int argc, char *argv[]) { + FILE *ftext; /* text file of vectors */ + FILE *ffloat; /* float file of vectors */ + int st,en; /* start and end values of vector to copy */ + float *buf; /* ptr to vector read from ftext */ + long lines; /* lines read so far */ + + if (argc != 5) { + printf("usage: extract TextFile FloatFile start end\n"); + exit(0); + } + + /* read command line arguments and open files */ + + ftext = fopen(argv[1],"rt"); + if (ftext == NULL) { + printf("Error opening text file: %s\n",argv[1]); + exit(1); + } + + ffloat = fopen(argv[2],"wb"); + if (ffloat == NULL) { + printf("Error opening float file: %s\n",argv[2]); + exit(1); + } + + st = atoi(argv[3]); + en = atoi(argv[4]); + + buf = (float*)malloc(en*sizeof(float)); + if (buf == NULL) { + printf("Error in malloc()\n"); + exit(1); + } + + lines = 0; + while(!feof(ftext)) { + scan_line(ftext, buf, en); + fwrite(&buf[st-1], sizeof(float), en-st+1, ffloat); + printf("\r%ld lines",lines++); + } + printf("\n"); + + /* clean up and exit */ + + free(buf); + fclose(ftext); + fclose(ffloat); + + return 0; +} + +/*---------------------------------------------------------------------------*\ + + FUNCTION....: scan_line() + + AUTHOR......: David Rowe + DATE CREATED: 20/9/96 + + This function reads a vector of floats from a line in a text file. + +\*---------------------------------------------------------------------------*/ + +void scan_line(FILE *fp, float f[], int n) +/* FILE *fp; file ptr to text file */ +/* float f[]; array of floats to return */ +/* int n; number of floats in line */ +{ + char s[MAX_STR]; + char *ps,*pe; + int i; + + fgets(s,MAX_STR,fp); + ps = pe = s; + for(i=0; i +#include +#include +#include "lpc.h" /* LPC analysis functions */ +#include "lsp.h" /* LSP encode/decode functions */ + +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 E; /* frame energy */ + long af; /* number frames with "active" speech */ + float Eres; /* LPC residual energy */ + int i; + int roots; + int unstables; + + /* Initialise ------------------------------------------------------*/ + + if (argc != 3) { + printf("usage: gentest RawFile LSPTextFile\n"); + exit(0); + } + + /* 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); + } + + for(i=0; i THRESH) { + af++; + printf("Active Frame: %ld unstables: %d\n",af, unstables); + + find_aks(Sn, ak, N, P, &Eres); + roots = lpc_to_lsp(&ak[1], P , lsp, 10, LSP_DELTA1, NULL); + if (roots == P) { + for(i=0; i +#include +#include +#include +#include + +/*-----------------------------------------------------------------------*\ + + DEFINES + +\*-----------------------------------------------------------------------*/ + +#define DELTAQ 0.01 /* quiting 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, 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 se; /* squared error for this iteration */ + float Dn,Dn_1; /* current and previous iterations distortion */ + float delta; /* improvement in distortion */ + FILE *ftrain; /* file containing training set */ + FILE *fvq; /* file containing vector quantiser */ + + /* Interpret command line arguments */ + + if (argc != 5) { + printf("usage: vqtrain TrainFile K M VQFile\n"); + exit(0); + } + + /* Open training file */ + + ftrain = fopen(argv[1],"rb"); + if (ftrain == NULL) { + printf("Error opening training database file: %s\n",argv[1]); + exit(1); + } + + /* determine k and m, and allocate arrays */ + + k = atol(argv[2]); + m = atol(argv[3]); + printf("dimension K=%ld number of entries M=%ld\n", 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 (cb == NULL || cb == NULL || cent == NULL || vec == NULL) { + printf("Error in malloc.\n"); + exit(1); + } + + /* determine size of training set */ + + J = 0; + while(fread(vec, sizeof(float), k, ftrain) == k) + J++; + printf("J=%ld entries in training set\n", J); + + /* set up initial codebook state from samples of training set */ + + rewind(ftrain); + fread(cb, sizeof(float), k*m, ftrain); + + /* main loop */ + + Dn = 1E32; + j = 1; + do { + Dn_1 = Dn; + + /* zero centroids */ + + for(i=0; i DELTAQ) + for(i=0; i DELTAQ); + + /* save codebook to disk */ + + fvq = fopen(argv[4],"wt"); + if (fvq == NULL) { + printf("Error opening VQ file: %s\n",argv[4]); + exit(1); + } + + for(j=0; j