ham-ref/codec2
1
0
Fork 0
mirror of https://github.com/drowe67/codec2 synced 2026-08-14 19:31:34 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

rm misc directory

Browse source
This commit is contained in:
drowe67 2023-07-13 08:49:31 +09:30 • committed by David Rowe
parent 153781e91a
commit 3cd3a9f820
24 changed files with 9 additions and 3297 deletions
Download patch file Download diff file Expand all files Collapse all files

7
CMakeLists.txt
View file

@ -229,7 +229,6 @@ if(UNITTEST)
message(STATUS "Threads library flags: ${CMAKE_THREAD_LIBS_INIT}")
add_subdirectory(unittest)
add_subdirectory(misc)
endif(UNITTEST)
message(STATUS "Build type is: " ${CMAKE_BUILD_TYPE})
@ -361,7 +360,7 @@ endif()
# just the +ve freq exp(jwn) so output power is 0.5 input power
add_test(NAME test_quisk_filter
COMMAND sh -c "cd ${CMAKE_CURRENT_SOURCE_DIR}/octave;
${CMAKE_CURRENT_BINARY_DIR}/misc/mksine in.raw 1500 1;
${CMAKE_CURRENT_BINARY_DIR}/unittest/mksine in.raw 1500 1;
cat in.raw | ${CMAKE_CURRENT_BINARY_DIR}/unittest/tquisk_filter |
sox -t .s16 -r 8000 -c 1 - -t .s16 out.raw vol 2;
cd ${CMAKE_CURRENT_SOURCE_DIR}/octave;
@ -384,7 +383,7 @@ endif()
# check channel simulator measures correct Peak to Average Power Ratio (about 0dB) with a sine wave input signal
add_test(NAME test_ch_papr
COMMAND sh -c "cd ${CMAKE_CURRENT_BINARY_DIR};
./misc/mksine - 1000 10 | ./src/ch - /dev/null --ctest")
./unittest/mksine - 1000 10 | ./src/ch - /dev/null --ctest")
add_test(NAME test_codec2_700c_octave_port
COMMAND sh -c "
@ -1162,7 +1161,7 @@ endif(NOT APPLE)
add_test(NAME test_vq_mbest
COMMAND sh -c "./tvq_mbest; \
cat target.f32 | \
../misc/vq_mbest -k 4 -q vq1.f32,vq2.f32 --st 1 --en 2 --mbest 2 -v > /dev/null;"
./vq_mbest -k 4 -q vq1.f32,vq2.f32 --st 1 --en 2 --mbest 2 -v > /dev/null;"
WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/unittest
)
set_tests_properties(test_vq_mbest PROPERTIES PASS_REGULAR_EXPRESSION "MSE: 0.00")

47
misc/16_8_short.c
View file

@ -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;
}

35
misc/CMakeLists.txt
View file

@ -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)

59
misc/de.c
View file

@ -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;
}

117
misc/extract.c
View file

@ -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;
}

306
misc/ge_train.c
View file

@ -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;
}

183
misc/genlsp.c
View file

@ -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;
}

54
misc/mksine.c
View file

@ -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;
}

59
misc/pre.c
View file

@ -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;
}

45
misc/raw2h.c
View file

@ -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;
}

58
misc/speexnoisesup.c
View file

@ -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;
}

220
misc/tcodec2.c
View file

@ -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;
}

141
misc/tdec.c
View file

@ -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];
}
}

128
misc/timpulse.c
View file

@ -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;
}

151
misc/tinterp.c
View file

@ -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;
}

153
misc/tlininterp.c
View file

@ -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;
}

164
misc/tnlp.c
View file

@ -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;
}

53
misc/tprede.c
View file

@ -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;
}

214
misc/tquant.c
View file

@ -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;
}

109
misc/tsrc.c
View file

@ -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;
}

296
misc/vq_binary_switch.c
View file

@ -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;
}

302
misc/vq_mbest.c
View file

@ -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]);
}

399
misc/vqtrain.c
View file

@ -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(&cent[i*k], k);
n[i] = 0;
}
/* quantise training set */
se = 0.0; noutliers[0] = noutliers[1] = noutliers[2] = 0;
rewind(ftrain);
for(i=0; i<J; i++) {
ret = fread(vec, sizeof(float), k, ftrain);
assert(ret == k);
ind = quantise(cb, vec, k, m, st, en, &e, &se);
n[ind]++;
acc(&cent[ind*k], vec, k);
if (sqrt(e/(en-st+1)) > 1.0) noutliers[0]++;
if (sqrt(e/(en-st+1)) > 2.0) noutliers[1]++;
if (sqrt(e/(en-st+1)) > 3.0) noutliers[2]++;
}
var = se/(J*(en-st+1));
delta = (var_1-var)/var;
int n_min = J;
int n_max = 0;
for(i=0; i<m; i++) {
if (n[i] < n_min) n_min = n[i];
if (n[i] > n_max) n_max = n[i];
}
printf("\r It: %2ld, var: %5f sd: %f outliers > 1/2/3 dB = %3.2f/%3.2f/%3.2f Delta = %5.4f %d %d\n", j, var, sqrt(var),
(float)noutliers[0]/J, (float)noutliers[1]/J, (float)noutliers[2]/J, delta, n_min, n_max);
j++;
/* determine new codebook from centroids */
if (delta > deltaq_stop)
for(i=0; i<m; i++) {
if (n[i] != 0) {
norm(&cent[i*k], k, n[i]);
memcpy(&cb[i*k], &cent[i*k], k*sizeof(float));
}
}
} while (delta > deltaq_stop);
/* save VQ to disk */
fvq = fopen(argv[dx+3],"wt");
if (fvq == NULL) {
printf("Error opening VQ file: %s\n",argv[dx+3]);
exit(1);
}
fwrite(cb, sizeof(float), m*k, fvq);
/* optionally output residual error for next stage VQ */
if (fres != NULL) {
float res[k];
rewind(ftrain);
for(j=0; j<J; j++) {
ret = fread(vec, sizeof(float), k, ftrain);
ind = quantise(cb, vec, k, m, st, en, &e, &se);
for(i=0; i<k; i++) {
res[i] = vec[i] - cb[k*ind+i];
}
fwrite(res, sizeof(float), k, fres);
}
fclose(fres);
}
fclose(fvq);
fclose(ftrain);
free(vec);
free(n);
return 0;
}
/*-----------------------------------------------------------------------*\
FUNCTIONS
\*-----------------------------------------------------------------------*/
/*-----------------------------------------------------------------------*\
FUNCTION....: zero()
AUTHOR......: David Rowe
DATE CREATED: 23/2/95
Zeros a vector of length k.
\*-----------------------------------------------------------------------*/
void zero(float v[], int k)
/* float v[]; ptr to start of vector */
/* int k; lngth of vector */
{
int i;
for(i=0; i<k; i++)
v[i] = 0.0;
}
/*---------------------------------------------------------------------------*\
FUNCTION....: acc()
AUTHOR......: David Rowe
DATE CREATED: 23/2/95
Adds k dimensional vectors v1 to v2 and stores the result back in v1.
\*---------------------------------------------------------------------------*/
void acc(float v1[], float v2[], int k)
/* float v1[]; ptr to start of vector to accumulate */
/* float v2[]; ptr to start of vector to add */
/* int k; dimension of vectors */
{
int i;
for(i=0; i<k; i++)
v1[i] += v2[i];
}
/*---------------------------------------------------------------------------*\
FUNCTION....: norm()
AUTHOR......: David Rowe
DATE CREATED: 23/2/95
Divides each element in k dimensional vector v by n.
\*---------------------------------------------------------------------------*/
void norm(float v[], int k, long n)
/* float v[]; ptr to start of vector */
/* int k; dimension of vectors */
/* int n; normalising factor */
{
int i;
assert(n != 0);
for(i=0; i<k; i++)
v[i] /= n;
}
/*---------------------------------------------------------------------------*\
FUNCTION....: quantise()
AUTHOR......: David Rowe
DATE CREATED: 23/2/95
Quantises vec by choosing the nearest vector in codebook cb, and
returns the vector index. The squared error of the quantised vector
is added to se.
\*---------------------------------------------------------------------------*/
long quantise(float cb[], // current VQ codebook
float vec[], // cb[][k]; current VQ codebook
int k, // dimension of vectors
int m, // size of codebook
int st, // start sample of vector to use in error calc
int en, // end sample of vector to use in error calc
float *beste, // current squared error
float *se) // accumulated squared error
{
long besti; /* best index so far */
long j;
int i;
float diff,e;
besti = 0;
*beste = 1E32;
for(j=0; j<m; j++) {
e = 0.0;
for(i=st; i<=en; i++) {
diff = cb[j*k+i]-vec[i];
e += diff*diff;
}
if (e < *beste) {
*beste = e;
besti = j;
}
}
*se += *beste;
return(besti);
}

6
unittest/CMakeLists.txt
View file

@ -141,3 +141,9 @@ target_compile_options(golay23 PUBLIC -DGOLAY23_UNITTEST)
add_executable(golay23_runtime_tables ../src/golay23.c)
target_compile_options(golay23_runtime_tables PUBLIC -DGOLAY23_UNITTEST -DRUN_TIME_TABLES)
add_executable(mksine mksine.c)
target_link_libraries(mksine m)
add_executable(vq_mbest vq_mbest.c)
target_link_libraries(vq_mbest codec2)