+
+#define N 160
+#define P 10
+
+#define LPC_FLOOR 0.0002 /* autocorrelation floor */
+#define LSP_DELTA1 0.2 /* grid spacing for LSP root searches */
+#define NDFT 256 /* DFT size for SD calculation */
+
+/* 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
+};
+
+/*---------------------------------------------------------------------------*\
+
+ find_aks_for_lsp()
+
+ This function takes a frame of samples, and determines the linear
+ prediction coefficients for that frame of samples. Modified version of
+ find_aks from lpc.c to include autocorrelation noise floor and lag window
+ to match Speex processing steps prior to LSP conversion.
+
+\*---------------------------------------------------------------------------*/
+
+void find_aks_for_lsp(
+ float Sn[], /* Nsam samples with order sample memory */
+ float a[], /* order+1 LPCs with first coeff 1.0 */
+ int Nsam, /* number of input speech samples */
+ int order, /* order of the LPC analysis */
+ float *E /* residual energy */
+)
+{
+ float Wn[N]; /* windowed frame of Nsam speech samples */
+ float R[P+1]; /* order+1 autocorrelation values of Sn[] */
+ int i;
+
+ hanning_window(Sn,Wn,Nsam);
+
+ autocorrelate(Wn,R,Nsam,order);
+ R[0] += LPC_FLOOR;
+ assert(order == 10); /* lag window only defined for order == 10 */
+ for(i=0; i<=order; i++)
+ R[i] *= lag_window[i];
+ levinson_durbin(R,a,order);
+
+ *E = 0.0;
+ for(i=0; i<=order; i++)
+ *E += a[i]*R[i];
+ if (*E < 0.0)
+ *E = 1E-12;
+}
+
+/*---------------------------------------------------------------------------*\
+
+ MAIN
+
+\*---------------------------------------------------------------------------*/
+
+int main(int argc, char *argv[])
+{
+ FILE *fin; /* input speech files */
+ short buf[N]; /* buffer of 16 bit speech samples */
+ float Sn[P+N]; /* input speech samples */
+ float E;
+ float ak[P+1]; /* LP coeffs */
+ float ak_[P+1]; /* quantised LP coeffs */
+ float lsp[P];
+ float lsp_[P]; /* quantised LSPs */
+ int roots; /* number of LSP roots found */
+ int frames; /* frames processed so far */
+ int i; /* loop variables */
+
+ SpeexBits bits;
+
+ float sd; /* SD for this frame */
+ float totsd; /* accumulated SD so far */
+ int gt2,gt4; /* number of frames > 2 and 4 dB SD */
+ int unstables; /* number of unstable LSP frames */
+
+ if (argc < 2) {
+ printf("usage: %s InputFile\n", argv[0]);
+ exit(0);
+ }
+
+ /* Open files */
+
+ if ((fin = fopen(argv[1],"rb")) == NULL) {
+ printf("Error opening input file: %s\n",argv[1]);
+ exit(0);
+ }
+
+ /* Initialise */
+
+ frames = 0;
+ for(i=0; i 2.0) gt2++;
+ if (sd > 4.0) gt4++;
+ totsd += sd;
+ }
+ else
+ unstables++;
+ }
+
+ fclose(fin);
+
+ printf("frames = %d Av sd = %3.2f dB", frames, totsd/frames);
+ printf(" >2 dB %3.2f%% >4 dB %3.2f%% unstables: %d\n",gt2*100.0/frames,
+ gt4*100.0/frames, unstables);
+
+ return 0;
+}
+
diff --git a/libs/libcodec2/unittest/t48_8 b/libs/libcodec2/unittest/t48_8
new file mode 100755
index 0000000000..9cf36606fc
Binary files /dev/null and b/libs/libcodec2/unittest/t48_8 differ
diff --git a/libs/libcodec2/unittest/t48_8.c b/libs/libcodec2/unittest/t48_8.c
new file mode 100644
index 0000000000..1eff979b20
--- /dev/null
+++ b/libs/libcodec2/unittest/t48_8.c
@@ -0,0 +1,112 @@
+/*
+ t48_8.c
+ David Rowe
+ May 10 2012
+
+ Unit test for 48 to 8 kHz sample rate conversion functions. I
+ evaluated output by plotting using Octave and looking for jaggies:
+
+ 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
+#include
+#include
+#include
+#include "fdmdv.h"
+
+#define N8 160 /* procssing buffer size at 8 kHz */
+#define N48 (N8*FDMDV_OS)
+#define MEM8 (FDMDV_OS_TAPS/FDMDV_OS)
+#define FRAMES 50
+#define TWO_PI 6.283185307
+#define FS 8000
+
+#define SINE
+
+int main() {
+ float in8k[MEM8 + N8];
+ float out48k[N48];
+ short out48k_short[N48];
+ FILE *f48;
+
+ float in48k[FDMDV_OS_TAPS + N48];
+ float out8k[N48];
+ short out8k_short[N8];
+ FILE *f8;
+
+ int i,f,t,t1;
+ float freq = 800.0;
+
+ f48 = fopen("out48.raw", "wb");
+ assert(f48 != NULL);
+ f8 = fopen("out8.raw", "wb");
+ assert(f8 != NULL);
+
+ /* clear filter memories */
+
+ for(i=0; i.
+*/
+
+#include
+#include
+#include
+#include
+#include
+
+#include "fdmdv_internal.h"
+#include "fdmdv.h"
+#include "octave.h"
+
+#define FRAMES 25
+#define CHANNEL_BUF_SIZE (10*M)
+
+int main(int argc, char *argv[])
+{
+ struct FDMDV *fdmdv;
+ int tx_bits[FDMDV_BITS_PER_FRAME];
+ COMP tx_symbols[NC+1];
+ COMP tx_baseband[NC+1][M];
+ COMP tx_fdm[M];
+ float channel[CHANNEL_BUF_SIZE];
+ int channel_count;
+ COMP rx_fdm[M+M/P];
+ float foff_coarse;
+ int nin, next_nin;
+ COMP rx_fdm_fcorr[M+M/P];
+ COMP rx_baseband[NC+1][M+M/P];
+ COMP rx_filt[NC+1][P+1];
+ float rx_timing;
+ float env[NT*P];
+ COMP rx_symbols[NC+1];
+ int rx_bits[FDMDV_BITS_PER_FRAME];
+ float foff_fine;
+ int sync_bit;
+ int fest_state;
+
+ int tx_bits_log[FDMDV_BITS_PER_FRAME*FRAMES];
+ COMP tx_symbols_log[(NC+1)*FRAMES];
+ COMP tx_baseband_log[(NC+1)][M*FRAMES];
+ COMP tx_fdm_log[M*FRAMES];
+ COMP pilot_baseband1_log[NPILOTBASEBAND*FRAMES];
+ COMP pilot_baseband2_log[NPILOTBASEBAND*FRAMES];
+ COMP pilot_lpf1_log[NPILOTLPF*FRAMES];
+ COMP pilot_lpf2_log[NPILOTLPF*FRAMES];
+ COMP S1_log[MPILOTFFT*FRAMES];
+ COMP S2_log[MPILOTFFT*FRAMES];
+ float foff_coarse_log[FRAMES];
+ float foff_log[FRAMES];
+ COMP rx_baseband_log[(NC+1)][(M+M/P)*FRAMES];
+ int rx_baseband_log_col_index;
+ COMP rx_filt_log[NC+1][(P+1)*FRAMES];
+ int rx_filt_log_col_index;
+ float env_log[NT*P*FRAMES];
+ float rx_timing_log[FRAMES];
+ COMP rx_symbols_log[NC+1][FRAMES];
+ float sig_est_log[NC+1][FRAMES];
+ float noise_est_log[NC+1][FRAMES];
+ int rx_bits_log[FDMDV_BITS_PER_FRAME*FRAMES];
+ float foff_fine_log[FRAMES];
+ int sync_bit_log[FRAMES];
+ int coarse_fine_log[FRAMES];
+ int nin_log[FRAMES];
+
+ FILE *fout;
+ int f,c,i,j;
+
+ fdmdv = fdmdv_create();
+ next_nin = M;
+ channel_count = 0;
+
+ rx_baseband_log_col_index = 0;
+ rx_filt_log_col_index = 0;
+
+ printf("sizeof FDMDV states: %d bytes\n", sizeof(struct FDMDV));
+
+ for(f=0; fprev_tx_symbols, tx_bits, &fdmdv->tx_pilot_bit);
+ memcpy(fdmdv->prev_tx_symbols, tx_symbols, sizeof(COMP)*(NC+1));
+ tx_filter(tx_baseband, tx_symbols, fdmdv->tx_filter_memory);
+ fdm_upconvert(tx_fdm, tx_baseband, fdmdv->phase_tx, fdmdv->freq);
+
+ /* --------------------------------------------------------*\
+ Channel
+ \*---------------------------------------------------------*/
+
+ nin = next_nin;
+ /*
+ if (f == 2)
+ nin = 120;
+ if (f == 3)
+ nin = 200;
+ if ((f !=2) && (f != 3))
+ nin = M;
+ */
+ /* add M tx samples to end of buffer */
+
+ assert((channel_count + M) < CHANNEL_BUF_SIZE);
+ for(i=0; icoarse_fine == COARSE)
+ fdmdv->foff = foff_coarse;
+ fdmdv_freq_shift(rx_fdm_fcorr, rx_fdm, fdmdv->foff, &fdmdv->foff_rect, &fdmdv->foff_phase_rect, nin);
+
+ /* baseband processing */
+
+ fdm_downconvert(rx_baseband, rx_fdm_fcorr, fdmdv->phase_rx, fdmdv->freq, nin);
+ rx_filter(rx_filt, rx_baseband, fdmdv->rx_filter_memory, nin);
+ rx_timing = rx_est_timing(rx_symbols, rx_filt, rx_baseband, fdmdv->rx_filter_mem_timing, env, fdmdv->rx_baseband_mem_timing, nin);
+ foff_fine = qpsk_to_bits(rx_bits, &sync_bit, fdmdv->phase_difference, fdmdv->prev_rx_symbols, rx_symbols);
+ snr_update(fdmdv->sig_est, fdmdv->noise_est, fdmdv->phase_difference);
+ memcpy(fdmdv->prev_rx_symbols, rx_symbols, sizeof(COMP)*(NC+1));
+
+ next_nin = M;
+
+ if (rx_timing > 2*M/P)
+ next_nin += M/P;
+
+ if (rx_timing < 0)
+ next_nin -= M/P;
+
+ fdmdv->coarse_fine = freq_state(sync_bit, &fdmdv->fest_state);
+ fdmdv->foff -= TRACK_COEFF*foff_fine;
+
+ /* --------------------------------------------------------*\
+ Log each vector
+ \*---------------------------------------------------------*/
+
+ memcpy(&tx_bits_log[FDMDV_BITS_PER_FRAME*f], tx_bits, sizeof(int)*FDMDV_BITS_PER_FRAME);
+ memcpy(&tx_symbols_log[(NC+1)*f], tx_symbols, sizeof(COMP)*(NC+1));
+ for(c=0; cpilot_baseband1, sizeof(COMP)*NPILOTBASEBAND);
+ memcpy(&pilot_baseband2_log[f*NPILOTBASEBAND], fdmdv->pilot_baseband2, sizeof(COMP)*NPILOTBASEBAND);
+ memcpy(&pilot_lpf1_log[f*NPILOTLPF], fdmdv->pilot_lpf1, sizeof(COMP)*NPILOTLPF);
+ memcpy(&pilot_lpf2_log[f*NPILOTLPF], fdmdv->pilot_lpf2, sizeof(COMP)*NPILOTLPF);
+ memcpy(&S1_log[f*MPILOTFFT], fdmdv->S1, sizeof(COMP)*MPILOTFFT);
+ memcpy(&S2_log[f*MPILOTFFT], fdmdv->S2, sizeof(COMP)*MPILOTFFT);
+ foff_coarse_log[f] = foff_coarse;
+ foff_log[f] = fdmdv->foff;
+
+ /* rx down conversion */
+
+ for(c=0; csig_est[c];
+ noise_est_log[c][f] = fdmdv->noise_est[c];
+ }
+ foff_fine_log[f] = foff_fine;
+ sync_bit_log[f] = sync_bit;
+
+ coarse_fine_log[f] = fdmdv->coarse_fine;
+ }
+
+
+ /*---------------------------------------------------------*\
+ Dump logs to Octave file for evaluation
+ by tfdmdv.m Octave script
+ \*---------------------------------------------------------*/
+
+ fout = fopen("tfdmdv_out.txt","wt");
+ assert(fout != NULL);
+ fprintf(fout, "# Created by tfdmdv.c\n");
+ octave_save_int(fout, "tx_bits_log_c", tx_bits_log, 1, FDMDV_BITS_PER_FRAME*FRAMES);
+ octave_save_complex(fout, "tx_symbols_log_c", tx_symbols_log, 1, (NC+1)*FRAMES, (NC+1)*FRAMES);
+ octave_save_complex(fout, "tx_baseband_log_c", (COMP*)tx_baseband_log, (NC+1), M*FRAMES, M*FRAMES);
+ octave_save_complex(fout, "tx_fdm_log_c", (COMP*)tx_fdm_log, 1, M*FRAMES, M*FRAMES);
+ octave_save_complex(fout, "pilot_lut_c", (COMP*)fdmdv->pilot_lut, 1, NPILOT_LUT, NPILOT_LUT);
+ octave_save_complex(fout, "pilot_baseband1_log_c", pilot_baseband1_log, 1, NPILOTBASEBAND*FRAMES, NPILOTBASEBAND*FRAMES);
+ octave_save_complex(fout, "pilot_baseband2_log_c", pilot_baseband2_log, 1, NPILOTBASEBAND*FRAMES, NPILOTBASEBAND*FRAMES);
+ octave_save_complex(fout, "pilot_lpf1_log_c", pilot_lpf1_log, 1, NPILOTLPF*FRAMES, NPILOTLPF*FRAMES);
+ octave_save_complex(fout, "pilot_lpf2_log_c", pilot_lpf2_log, 1, NPILOTLPF*FRAMES, NPILOTLPF*FRAMES);
+ octave_save_complex(fout, "S1_log_c", S1_log, 1, MPILOTFFT*FRAMES, MPILOTFFT*FRAMES);
+ octave_save_complex(fout, "S2_log_c", S2_log, 1, MPILOTFFT*FRAMES, MPILOTFFT*FRAMES);
+ octave_save_float(fout, "foff_log_c", foff_log, 1, FRAMES, FRAMES);
+ octave_save_float(fout, "foff_coarse_log_c", foff_coarse_log, 1, FRAMES, FRAMES);
+ octave_save_complex(fout, "rx_baseband_log_c", (COMP*)rx_baseband_log, (NC+1), rx_baseband_log_col_index, (M+M/P)*FRAMES);
+ octave_save_complex(fout, "rx_filt_log_c", (COMP*)rx_filt_log, (NC+1), rx_filt_log_col_index, (P+1)*FRAMES);
+ octave_save_float(fout, "env_log_c", env_log, 1, NT*P*FRAMES, NT*P*FRAMES);
+ octave_save_float(fout, "rx_timing_log_c", rx_timing_log, 1, FRAMES, FRAMES);
+ octave_save_complex(fout, "rx_symbols_log_c", (COMP*)rx_symbols_log, (NC+1), FRAMES, FRAMES);
+ octave_save_float(fout, "sig_est_log_c", (float*)sig_est_log, (NC+1), FRAMES, FRAMES);
+ octave_save_float(fout, "noise_est_log_c", (float*)noise_est_log, (NC+1), FRAMES, FRAMES);
+ octave_save_int(fout, "rx_bits_log_c", rx_bits_log, 1, FDMDV_BITS_PER_FRAME*FRAMES);
+ octave_save_float(fout, "foff_fine_log_c", foff_fine_log, 1, FRAMES, FRAMES);
+ octave_save_int(fout, "sync_bit_log_c", sync_bit_log, 1, FRAMES);
+ octave_save_int(fout, "coarse_fine_log_c", coarse_fine_log, 1, FRAMES);
+ octave_save_int(fout, "nin_log_c", nin_log, 1, FRAMES);
+ fclose(fout);
+
+ fdmdv_destroy(fdmdv);
+
+ return 0;
+}
+
diff --git a/libs/libcodec2/unittest/tfifo b/libs/libcodec2/unittest/tfifo
new file mode 100755
index 0000000000..bdbf6f9fb4
Binary files /dev/null and b/libs/libcodec2/unittest/tfifo differ
diff --git a/libs/libcodec2/unittest/tfifo.c b/libs/libcodec2/unittest/tfifo.c
new file mode 100644
index 0000000000..f52c4cf551
--- /dev/null
+++ b/libs/libcodec2/unittest/tfifo.c
@@ -0,0 +1,103 @@
+/*
+ tfifo.c
+ David Rowe
+ Nov 19 2012
+
+ Takes FIFOs, in particular thread safety.
+*/
+
+#include
+#include
+#include
+#include "fifo.h"
+
+#define FIFO_SZ 1024
+#define WRITE_SZ 10
+#define READ_SZ 8
+#define N_MAX 100
+#define LOOPS 1000000
+
+int run_thread = 1;
+struct FIFO *f;
+
+void writer(void);
+void *writer_thread(void *data);
+pthread_mutex_t mutex = PTHREAD_MUTEX_INITIALIZER;
+
+#define USE_THREADS
+//#define USE_MUTEX
+
+int main() {
+ pthread_t awriter_thread;
+ int i,j;
+ short read_buf[READ_SZ];
+ int n_out = 0;
+ int sucess;
+
+ f = fifo_create(FIFO_SZ);
+ #ifdef USE_THREADS
+ pthread_create(&awriter_thread, NULL, writer_thread, NULL);
+ #endif
+
+ for(i=0; i WRITE_SZ) {
+ for(i=0; i.
+*/
+
+#include
+#include
+#include
+#include
+#include
+#include "defines.h"
+#include "comp.h"
+#include "codec2.h"
+#include "quantise.h"
+#include "interp.h"
+#include "codec2_internal.h"
+
+float run_a_test(char raw_file_name[], int bit_to_corrupt)
+{
+ FILE *fin;
+ short buf[N];
+ struct CODEC2 *c2;
+ kiss_fft_cfg fft_fwd_cfg;
+ MODEL model;
+ float ak[LPC_ORD+1];
+ float lsps[LPC_ORD], e;
+ int lsp_indexes[LPC_ORD], found_bit;
+ float snr, snr_sum;
+ int frames, i, mask, index;
+
+ c2 = codec2_create(CODEC2_MODE_2400);
+ fft_fwd_cfg = kiss_fft_alloc(FFT_ENC, 0, NULL, NULL);
+
+ fin = fopen(raw_file_name, "rb");
+ assert(fin != NULL);
+
+ /* find bit we are corrupting */
+
+ found_bit = 0;
+ for(i=0; i lsp_bits(i))
+ bit_to_corrupt -= lsp_bits(i);
+ else {
+ index = i;
+ mask = (1 << bit_to_corrupt);
+ printf(" index: %d bit: %d mask: 0x%x ", index, bit_to_corrupt, mask);
+ found_bit = 1;
+ }
+ }
+ }
+ assert(found_bit == 1);
+
+ /* OK test a sample file, flipping bit */
+
+ snr_sum = 0.0;
+ frames = 0;
+ while(fread(buf, sizeof(short), N, fin) == N) {
+ analyse_one_frame(c2, &model, buf);
+ e = speech_to_uq_lsps(lsps, ak, c2->Sn, c2->w, LPC_ORD);
+ encode_lsps_scalar(lsp_indexes, lsps, LPC_ORD);
+
+ /* find and flip bit we are testing */
+
+ lsp_indexes[index] ^= mask;
+
+ /* decode LSPs and measure SNR */
+
+ decode_lsps_scalar(lsps, lsp_indexes, LPC_ORD);
+ check_lsp_order(lsps, LPC_ORD);
+ bw_expand_lsps(lsps, LPC_ORD);
+ lsp_to_lpc(lsps, ak, LPC_ORD);
+ aks_to_M2(fft_fwd_cfg, ak, LPC_ORD, &model, e, &snr, 0, 0, 1, 1, LPCPF_BETA, LPCPF_GAMMA);
+ snr_sum += snr;
+ frames++;
+ }
+
+ codec2_destroy(c2);
+
+ fclose(fin);
+
+ return snr_sum/frames;
+}
+
+int main(int argc, char *argv[]) {
+ int i;
+ int total_lsp_bits = 0;
+ float snr;
+
+ if (argc != 2) {
+ printf("usage: %s RawFile\n", argv[0]);
+ exit(1);
+ }
+
+ for(i=0; i
+#include
+#include
+#include
+#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
+#include
+#include
+#include
+#include
+
+#define N8 160 /* processing buffer size at 8 kHz */
+#define N48 ((int)N8*(48000/8000)) /* buf size assuming 48k max sample rate */
+
+int main(int argc, char *argv[]) {
+ FILE *f8k, *fout;
+ short in8k_short[N8];
+ float in8k[N8];
+ float out[N48];
+ short out_short[N48];
+ SRC_STATE *src;
+ SRC_DATA data;
+ int error;
+
+ if (argc != 4) {
+ printf("usage %s inputRawFile OutputRawFile OutputSamplerate\n", argv[0]);
+ exit(0);
+ }
+
+ f8k = fopen(argv[1], "rb");
+ assert(f8k != NULL);
+
+ fout = fopen(argv[2], "wb");
+ assert(fout != NULL);
+
+ src = src_new(SRC_SINC_FASTEST, 1, &error);
+ assert(src != NULL);
+
+ data.data_in = in8k;
+ data.data_out = out;
+ data.input_frames = N8;
+ data.output_frames = N48;
+ data.end_of_input = 0;
+ data.src_ratio = atof(argv[3])/8000;
+ printf("%f\n", data.src_ratio);
+
+ while(fread(in8k_short, sizeof(short), N8, f8k) == N8) {
+ src_short_to_float_array(in8k_short, in8k, N8);
+ src_process(src, &data);
+ printf("%d %d\n", (int)data.output_frames , (int)data.output_frames_gen);
+ assert(data.output_frames_gen <= N48);
+ src_float_to_short_array(out, out_short, data.output_frames_gen);
+ fwrite(out_short, sizeof(short), data.output_frames_gen, fout);
+ }
+
+ fclose(fout);
+ fclose(f8k);
+
+ return 0;
+}
diff --git a/libs/libcodec2/unittest/vq_train_jvm b/libs/libcodec2/unittest/vq_train_jvm
new file mode 100755
index 0000000000..81a28b4d91
Binary files /dev/null and b/libs/libcodec2/unittest/vq_train_jvm differ
diff --git a/libs/libcodec2/unittest/vq_train_jvm.c b/libs/libcodec2/unittest/vq_train_jvm.c
new file mode 100755
index 0000000000..9da091e085
--- /dev/null
+++ b/libs/libcodec2/unittest/vq_train_jvm.c
@@ -0,0 +1,486 @@
+/*---------------------------------------------------------------------------*\
+
+ FILE........: vq_train_jvm.c
+ AUTHOR......: Jean-Marc Valin
+ DATE CREATED: 21 Jan 2012
+
+ Multi-stage Vector Quantoser training program developed by Jean-Marc at
+ linux.conf.au 2012. Minor mods by David Rowe
+
+\*---------------------------------------------------------------------------*/
+
+/*
+ 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, 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 .
+*/
+
+
+#ifdef VALGRIND
+#include
+#endif
+
+#include
+#include
+#include
+#include
+
+#define MIN(a,b) ((a)<(b)?(a):(b))
+#define COEF 0.0f
+#define MAX_ENTRIES 16384
+
+void compute_weights(const float *x, float *w, int ndim)
+{
+ int i;
+ w[0] = MIN(x[0], x[1]-x[0]);
+ for (i=1;i.
+*/
+
+/*-----------------------------------------------------------------------*\
+
+ INCLUDES
+
+\*-----------------------------------------------------------------------*/
+
+#include
+#include
+#include
+#include
+#include
+#include
+
+typedef struct {
+ float real;
+ float imag;
+} COMP;
+
+/*-----------------------------------------------------------------------* \
+
+ DEFINES
+
+\*-----------------------------------------------------------------------*/
+
+#define DELTAQ 0.01 /* quiting distortion */
+#define MAX_STR 80 /* maximum string length */
+#define PI 3.141592654
+
+/*-----------------------------------------------------------------------*\
+
+ FUNCTION PROTOTYPES
+
+\*-----------------------------------------------------------------------*/
+
+void zero(COMP v[], int d);
+void acc(COMP v1[], COMP v2[], int d);
+void norm(COMP v[], int k);
+int quantise(COMP cb[], COMP vec[], int d, int e, float *se);
+void print_vec(COMP cb[], int d, int e);
+
+/*-----------------------------------------------------------------------* \
+
+ MAIN
+
+\*-----------------------------------------------------------------------*/
+
+int main(int argc, char *argv[]) {
+ int d,e; /* dimension and codebook size */
+ COMP *vec; /* current vector */
+ COMP *cb; /* vector codebook */
+ COMP *cent; /* centroids for each codebook entry */
+ int *n; /* number of vectors in this interval */
+ int J; /* number of vectors in training set */
+ int ind; /* index of current vector */
+ float se; /* total squared error for this iteration */
+ float var; /* variance */
+ float var_1; /* previous variance */
+ float delta; /* improvement in distortion */
+ FILE *ftrain; /* file containing training set */
+ FILE *fvq; /* file containing vector quantiser */
+ int ret;
+ int i,j, finished, iterations;
+ float b; /* equivalent number of bits */
+ float improvement;
+ float sd_vec, sd_element, sd_theory, bits_theory;
+ int var_n;
+
+ /* Interpret command line arguments */
+
+ if (argc != 5) {
+ printf("usage: %s TrainFile D(dimension) E(number of entries) VQFile\n", argv[0]);
+ exit(1);
+ }
+
+ /* 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 */
+
+ d = atoi(argv[2]);
+ e = atoi(argv[3]);
+ printf("\n");
+ printf("dimension D=%d number of entries E=%d\n", d, e);
+ vec = (COMP*)malloc(sizeof(COMP)*d);
+ cb = (COMP*)malloc(sizeof(COMP)*d*e);
+ cent = (COMP*)malloc(sizeof(COMP)*d*e);
+ n = (int*)malloc(sizeof(int)*e);
+ if (cb == NULL || cb == NULL || cent == NULL || vec == NULL) {
+ printf("Error in malloc.\n");
+ exit(1);
+ }
+
+ /* determine size of training set */
+
+ J = 0;
+ var_n = 0;
+ while(fread(vec, sizeof(COMP), d, ftrain) == (size_t)d) {
+ for(j=0; j 1) {
+ if (var > 0.0) {
+ delta = (var_1 - var)/var;
+ }
+ else
+ delta = 0;
+ if (delta < DELTAQ)
+ finished = 1;
+ }
+
+ if (!finished) {
+ /* determine new codebook from centroids */
+
+ for(i=0; i.
+*/
+
+/*-----------------------------------------------------------------------*\
+
+ INCLUDES
+
+\*-----------------------------------------------------------------------*/
+
+#include
+#include
+#include
+#include
+#include
+#include
+
+typedef struct {
+ float real;
+ float imag;
+} COMP;
+
+/*-----------------------------------------------------------------------* \
+
+ DEFINES
+
+\*-----------------------------------------------------------------------*/
+
+#define DELTAQ 0.01 /* quiting distortion */
+#define MAX_STR 80 /* maximum string length */
+
+/*-----------------------------------------------------------------------*\
+
+ FUNCTION PROTOTYPES
+
+\*-----------------------------------------------------------------------*/
+
+void zero(float v[], int d);
+void acc(float v1[], float v2[], int d);
+void norm(float v[], int k, int n[]);
+int quantise(float cb[], float vec[], int d, int e, float *se);
+void print_vec(float cb[], int d, int e);
+void split(float cb[], int d, int b);
+int gain_shape_quantise(float cb[], float vec[], int d, int e, float *se, float *best_gain);
+
+/*-----------------------------------------------------------------------* \
+
+ MAIN
+
+\*-----------------------------------------------------------------------*/
+
+int main(int argc, char *argv[]) {
+ int d,e; /* dimension and codebook size */
+ float *vec; /* current vector */
+ float *cb; /* vector codebook */
+ float *cent; /* centroids for each codebook entry */
+ int *n; /* number of vectors in this interval */
+ int J; /* number of vectors in training set */
+ int ind; /* index of current vector */
+ float se; /* total squared error for this iteration */
+ float var; /* variance */
+ float var_1; /* previous variance */
+ float delta; /* improvement in distortion */
+ FILE *ftrain; /* file containing training set */
+ FILE *fvq; /* file containing vector quantiser */
+ int ret;
+ int i,j, finished, iterations;
+ float sd;
+ int var_n, bits, b, levels;
+
+ /* Interpret command line arguments */
+
+ if (argc < 5) {
+ printf("usage: %s TrainFile D(dimension) B(number of bits) VQFile [error.txt file]\n", argv[0]);
+ exit(1);
+ }
+
+ /* 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 */
+
+ d = atoi(argv[2]);
+ bits = atoi(argv[3]);
+ e = 1< 1) {
+ if (var > 0.0) {
+ delta = (var_1 - var)/var;
+ }
+ else
+ delta = 0;
+ if (delta < DELTAQ)
+ finished = 1;
+ }
+
+ if (!finished) {
+ /* determine new codebook from centroids */
+
+ for(i=0; i