From: drowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Date: Fri, 10 Aug 2012 23:58:14 +0000 (+0000)
Subject: various experimental ampl vq schemes inclusing shape, gain shape, interpolation using VQ
X-Git-Url: http://git.whiteaudio.com/gitweb/?a=commitdiff_plain;h=6a3f94cef066d7f7d4fd9ebe7e97de731e7f9b2c;p=freetel-svn-tracking.git

various experimental ampl vq schemes inclusing shape, gain shape, interpolation using VQ

git-svn-id: https://svn.code.sf.net/p/freetel/code@610 01035d8c-6547-0410-b346-abe4f91aad63
---

diff --git a/codec2-dev/src/ampexp.c b/codec2-dev/src/ampexp.c
index dd159865..3cc7ab00 100644
--- a/codec2-dev/src/ampexp.c
+++ b/codec2-dev/src/ampexp.c
@@ -31,6 +31,7 @@
 #include <math.h>
 #include <stdio.h>
 #include <stdlib.h>
+#include <string.h>
 
 #include "ampexp.h"
 
@@ -56,8 +57,12 @@ struct AEXP {
     int              var_n;
     float            vq_var;
     int              vq_var_n;
-    struct codebook *vq;
-    int              offset;
+    struct codebook *vq1,*vq2,*vq3,*vq4,*vq5;
+
+    int              indexes[3];
+    MODEL            model[3];
+    float            mag[3];
+    MODEL            model_uq[3];
 };
 
 
@@ -142,7 +147,7 @@ static struct codebook *load(const char * name)
 
 struct AEXP *amp_experiment_create() {
     struct AEXP *aexp;
-    int i;
+    int i,m;
 
     aexp = (struct AEXP *)malloc(sizeof(struct AEXP));
     assert (aexp != NULL);
@@ -156,9 +161,57 @@ struct AEXP *amp_experiment_create() {
     aexp->var_n = 0;
     aexp->vq_var = 0.0;
     aexp->vq_var_n = 0;
-    aexp->vq = load("../unittest/amp1_20_1024.txt");
-    //aexp->vq = load("hts1a_vq1_20_300.txt");
-    aexp->offset = 20;
+
+    aexp->vq1 = load("../unittest/amp1_20_1024.txt");
+    aexp->vq1->offset = 0;
+    #ifdef CAND2_GS
+    //aexp->vq1 = load("../unittest/t1_amp1_20_1024.txt");
+    //aexp->vq1 = load("../unittest/t2_amp1_20_1024.txt");
+    aexp->vq1 = load("../unittest/amp1_20_1024.txt");
+    aexp->vq1->offset = 0;
+    aexp->vq2 = load("../unittest/amp21_40_1024.txt");
+    aexp->vq2->offset = 20;
+    aexp->vq3 = load("../unittest/amp41_60_1024.txt");
+    aexp->vq3->offset = 40;
+    aexp->vq4 = load("../unittest/amp61_80_32.txt");
+    aexp->vq4->offset = 60;
+    #endif
+
+    //#define CAND2_GS
+    #ifdef CAND2_GS
+    aexp->vq1 = load("../unittest/amp1_20_1024.txt");
+    aexp->vq2 = load("../unittest/amp21_40_1024.txt");
+    aexp->vq3 = load("../unittest/amp41_80_1024.txt");
+    aexp->vq4 = load("../unittest/amp61_80_32.txt");
+    aexp->vq1->offset = 0;
+    aexp->vq2->offset = 20;
+    aexp->vq3->offset = 40;
+    aexp->vq4->offset = 60;
+    #endif
+
+    //#define CAND1
+    #ifdef CAND1
+    aexp->vq1 = load("../unittest/amp1_10_128.txt");
+    aexp->vq2 = load("../unittest/amp11_20_512.txt");
+    aexp->vq3 = load("../unittest/amp21_40_1024.txt");
+    aexp->vq4 = load("../unittest/amp41_60_1024.txt");
+    aexp->vq5 = load("../unittest/amp61_80_32.txt");
+    aexp->vq1->offset = 0;
+    aexp->vq2->offset = 10;
+    aexp->vq3->offset = 20;
+    aexp->vq4->offset = 40;
+    aexp->vq5->offset = 60;
+    #endif
+
+    for(i=0; i<3; i++) {
+	aexp->indexes[i] = 0;
+	aexp->mag[i] = 1.0;
+	aexp->model[i].Wo = TWO_PI*100.0/8000.0;
+	aexp->model[i].L = floor(PI/aexp->model[i].Wo);	
+	for(m=1; m<=MAX_AMP; m++)
+	    aexp->model[i].A[m] = 10.0;
+	aexp->model_uq[i] = aexp->model[i];
+    }
 
     return aexp;
 }
@@ -241,7 +294,7 @@ static void print_sparse_pred_error(struct AEXP *aexp, MODEL *model, float mag_t
     for(m=1; m<=model->L; m++)
 	mag += model->A[m]*model->A[m];
     mag = 10*log10(mag/model->L);
-    
+
     if (mag > mag_thresh) {
 	for(m=0; m<MAX_AMP; m++) {
 	    sparse_pe[m] = 0.0;
@@ -250,15 +303,51 @@ static void print_sparse_pred_error(struct AEXP *aexp, MODEL *model, float mag_t
 	for(m=1; m<=model->L; m++) {
 	    assert(model->A[m] > 0.0);
 	    error = PRED_COEFF*20.0*log10(aexp->A_prev[m]) - 20.0*log10(model->A[m]);
+	    //error = 20.0*log10(model->A[m]) - mag;
+
+	    index = MAX_AMP*m*model->Wo/PI;
+	    assert(index < MAX_AMP);
+	    sparse_pe[index] = error;
+	}
+
+	/* dump sparse amp vector */
+
+	for(m=0; m<MAX_AMP; m++)
+	    printf("%f ", sparse_pe[m]);
+	printf("\n");
+    }
+}
+
+
+static void print_sparse_amp_error(struct AEXP *aexp, MODEL *model, float edB_thresh)
+{
+    int    m, index;
+    float  e, edB, enormdB, error, level;
+    float  sparse_pe[MAX_AMP];
+
+    e = 0.0;
+    for(m=1; m<=model->L; m++)
+	e += model->A[m]*model->A[m];
+    edB = 10*log10(e);
+    enormdB = 10*log10(e/model->L); /* make high and low pitches have similar amps */
+
+    if (edB > edB_thresh) {
+	for(m=0; m<MAX_AMP; m++) {
+	    sparse_pe[m] = 0.0;
+	}
+
+	for(m=1; m<=model->L; m++) {
+	    assert(model->A[m] > 0.0);
+	    error = 20.0*log10(model->A[m]) - enormdB;
 
 	    index = MAX_AMP*m*model->Wo/PI;
 	    assert(index < MAX_AMP);
 	    sparse_pe[index] = error;
 	}
 
-	/* dump spare phase vector */
+	/* dump sparse amp vector */
 
-	for(m=0; m<MAX_AMP/4; m++)
+	for(m=0; m<MAX_AMP; m++)
 	    printf("%f ", sparse_pe[m]);
 	printf("\n");
     }
@@ -274,7 +363,7 @@ int vq_amp(float cb[], float vec[], float weights[], int d, int e, float *se)
    float   diff, metric, best_metric;
 
    besti = 0;
-   best_metric = 1E32;
+   best_metric = best_error = 1E32;
    for(j=0; j<e; j++) {
        metric = error = 0.0;
        for(i=0; i<d; i++) {
@@ -297,23 +386,24 @@ int vq_amp(float cb[], float vec[], float weights[], int d, int e, float *se)
 }
 
 
-static void split_vq(float sparse_pe_out[], struct AEXP *pexp, struct codebook *vq, float weights[], float sparse_pe_in[])
+static int split_vq(float sparse_pe_out[], struct AEXP *aexp, struct codebook *vq, float weights[], float sparse_pe_in[])
 {
     int i, j, non_zero, vq_ind;
     float se;
 
     vq_ind = vq_amp(vq->cb, &sparse_pe_in[vq->offset], &weights[vq->offset], vq->k, vq->m, &se);
     printf("\n offset %d k %d m %d vq_ind %d j: ", vq->offset, vq->k, vq->m, vq_ind);
-  
+    
     non_zero = 0;
     for(i=0, j=vq->offset; i<vq->k; i++,j++) {
 	if (sparse_pe_in[j] != 0.0) {
-	    //printf("%d ", j);
+	    printf("%d ", j);
 	    sparse_pe_out[j] = vq->cb[vq->k * vq_ind + i];
 	    non_zero++;
 	}
     }
-    pexp->vq_var_n += non_zero;
+    aexp->vq_var_n += non_zero;
+    return vq_ind;
 }
 
 
@@ -351,11 +441,13 @@ static void sparse_vq_pred_error(struct AEXP *aexp,
 	sparse_pe_out[m] = sparse_pe_in[m];
     }
 
-    split_vq(sparse_pe_out, aexp, aexp->vq, weights, sparse_pe_in);
-    #ifdef SIM_VQ
-    for(m=0; m<MAX_AMP; m++) {
+    //#define SIM_VQ
+    #ifndef SIM_VQ
+    split_vq(sparse_pe_out, aexp, aexp->vq1, weights, sparse_pe_in);
+    #else
+    for(m=aexp->vq->offset; m<aexp->vq->offset+aexp->vq->k; m++) {
 	if (sparse_pe_in[m] != 0.0) {
-	    float error = 4.0*(1.0 - 2.0*rand()/RAND_MAX);
+	    float error = 8*(1.0 - 2.0*rand()/RAND_MAX);
 	    aexp->vq_var += error*error;
 	    aexp->vq_var_n++;
 	    sparse_pe_out[m] = sparse_pe_in[m] + error;
@@ -363,7 +455,7 @@ static void sparse_vq_pred_error(struct AEXP *aexp,
     }
     #endif
 
-    if (mag > 40.0)
+    if (mag > -100.0)
 	for(m=0; m<MAX_AMP; m++) {
 	    if (sparse_pe_in[m] != 0.0)
 		aexp->vq_var += pow(sparse_pe_out[m] - sparse_pe_in[m], 2.0);
@@ -383,26 +475,272 @@ static void sparse_vq_pred_error(struct AEXP *aexp,
 }
 
 
-static void update_snr_calc(struct AEXP *aexp, MODEL *model, float before[])
+static void sparse_vq_amp(struct AEXP *aexp, MODEL *model)
+{
+    int    m, index;
+    float  error, amp_dB, mag;
+    float  sparse_pe_in[MAX_AMP];
+    float  sparse_pe_out[MAX_AMP];
+    float  weights[MAX_AMP];
+
+    mag = 0.0;
+    for(m=1; m<=model->L; m++)
+	mag += model->A[m]*model->A[m];
+    mag = 10*log10(mag/model->L);
+    aexp->mag[2] = mag;
+   
+    for(m=0; m<MAX_AMP; m++) {
+	sparse_pe_in[m] = 0.0;
+	sparse_pe_out[m] = 0.0;
+    }
+
+    for(m=1; m<=model->L; m++) {
+ 	assert(model->A[m] > 0.0);
+	error = 20.0*log10(model->A[m]) - mag;
+
+	index = MAX_AMP*m*model->Wo/PI;
+	assert(index < MAX_AMP);
+	sparse_pe_in[index] = error;
+	weights[index] = model->A[m];
+    }
+
+    /* vector quantise */
+        
+    for(m=0; m<MAX_AMP; m++) {
+	sparse_pe_out[m] = sparse_pe_in[m];
+    }
+
+    //#define SIM_VQ
+    #ifndef SIM_VQ
+    printf("hello %d\n", aexp->frames);
+    aexp->indexes[2] = split_vq(sparse_pe_out, aexp, aexp->vq1, weights, sparse_pe_in);
+    //split_vq(sparse_pe_out, aexp, aexp->vq2, weights, sparse_pe_in);
+    //split_vq(sparse_pe_out, aexp, aexp->vq3, weights, sparse_pe_in);
+    //split_vq(sparse_pe_out, aexp, aexp->vq4, weights, sparse_pe_in);
+    //split_vq(sparse_pe_out, aexp, aexp->vq5, weights, sparse_pe_in);
+    #else
+    for(m=aexp->vq->offset; m<aexp->vq->offset+aexp->vq->k; m++) {
+	if (sparse_pe_in[m] != 0.0) {
+	    float error = 0*(1.0 - 2.0*rand()/RAND_MAX);
+	    aexp->vq_var += error*error;
+	    aexp->vq_var_n++;
+	    sparse_pe_out[m] = sparse_pe_in[m] + error;
+	}
+    }
+    #endif
+
+    for(m=0; m<MAX_AMP; m++) {
+	if (sparse_pe_in[m] != 0.0)
+	    aexp->vq_var += pow(sparse_pe_out[m] - sparse_pe_in[m], 2.0);
+    }
+    
+    /* transform quantised amps back */
+
+    for(m=1; m<=model->L; m++) {
+	index = MAX_AMP*m*model->Wo/PI;
+	assert(index < MAX_AMP);
+	amp_dB = sparse_pe_out[index] + mag;
+	model->A[m] = pow(10.0, amp_dB/20.0);
+    }
+    //exit(0);
+}
+
+
+static void update_snr_calc(struct AEXP *aexp, MODEL *m1, MODEL *m2)
 {
     int m;
     float signal, noise, signal_dB;
 
-    signal = 0.0; noise = 0.0;
-    for(m=1; m<=model->L; m++) {	    
-	signal += before[m]*before[m];
-	noise  += pow(before[m] - model->A[m], 2.0);
+    assert(m1->L == m2->L);
+
+    signal = 0.0; noise = 1E-32;
+    for(m=1; m<=m1->L; m++) {	    
+	signal += m1->A[m]*m1->A[m];
+	noise  += pow(m1->A[m] - m2->A[m], 2.0);
 	//printf("%f %f\n", before[m], model->phi[m]);
     }
-    signal_dB = 10*log10(signal/model->L);
-    if (signal_dB > - 100.0) {
+    signal_dB = 10*log10(signal);
+    if (signal_dB > -100.0) {
 	aexp->snr += 10.0*log10(signal/noise);
 	aexp->snr_n++;
     }
 }
 
 
-/*---------------------------------------------------------------------------*\
+/* gain/shape vq search.  Returns index of best gain.  Gain is additive (as we use log quantisers) */
+
+int gain_shape_vq_amp(float cb[], float vec[], float weights[], int d, int e, float *se, float *best_gain)
+{
+   float   error;	/* current error		*/
+   int     besti;	/* best index so far		*/
+   float   best_error;	/* best error so far		*/
+   int	   i,j,m;
+   float   diff, metric, best_metric, gain, sumAm, sumCb;
+
+   besti = 0;
+   best_metric = best_error = 1E32;
+   for(j=0; j<e; j++) {
+
+       /* compute optimum gain */
+
+       sumAm = sumCb = 0.0;
+       m = 0;
+       for(i=0; i<d; i++) {
+	   if (vec[i] != 0.0) {
+	       m++;
+	       sumAm += vec[i];
+	       sumCb += cb[j*d+i];
+	   }
+       }
+       gain = (sumAm - sumCb)/m;
+       
+       /* compute error */
+
+       metric = error = 0.0;
+       for(i=0; i<d; i++) {
+	   if (vec[i] != 0.0) {
+	       diff = vec[i] - cb[j*d+i] - gain;
+	       error += diff*diff;
+	       metric += weights[i]*diff*diff;
+	   }
+       }
+       if (metric < best_metric) {
+	   best_error = error;
+	   best_metric = metric;
+	   *best_gain = gain;
+	   besti = j;
+       }
+   }
+
+   *se += best_error;
+
+   return(besti);
+}
+
+
+static void gain_shape_split_vq(float sparse_pe_out[], struct AEXP *aexp, struct codebook *vq, float weights[], float sparse_pe_in[], float *best_gain)
+{
+    int i, j, non_zero, vq_ind;
+    float se;
+
+    vq_ind = gain_shape_vq_amp(vq->cb, &sparse_pe_in[vq->offset], &weights[vq->offset], vq->k, vq->m, &se, best_gain);
+    printf("\n offset %d k %d m %d vq_ind %d gain: %4.2f j: ", vq->offset, vq->k, vq->m, vq_ind, *best_gain);
+  
+    non_zero = 0;
+    for(i=0, j=vq->offset; i<vq->k; i++,j++) {
+	if (sparse_pe_in[j] != 0.0) {
+	    printf("%d ", j);
+	    sparse_pe_out[j] = vq->cb[vq->k * vq_ind + i] + *best_gain;
+	    non_zero++;
+	}
+    }
+    aexp->vq_var_n += non_zero;
+}
+
+
+static void gain_shape_sparse_vq_amp(struct AEXP *aexp, MODEL *model)
+{
+    int    m, index;
+    float  error, amp_dB, best_gain;
+    float  sparse_pe_in[MAX_AMP];
+    float  sparse_pe_out[MAX_AMP];
+    float  weights[MAX_AMP];
+
+    for(m=0; m<MAX_AMP; m++) {
+	sparse_pe_in[m] = 0.0;
+	sparse_pe_out[m] = 0.0;
+    }
+
+    for(m=1; m<=model->L; m++) {
+ 	assert(model->A[m] > 0.0);
+
+	index = MAX_AMP*m*model->Wo/PI;
+	assert(index < MAX_AMP);
+	sparse_pe_in[index] = 20.0*log10(model->A[m]);
+	weights[index] = model->A[m];
+    }
+
+    /* vector quantise */
+        
+    for(m=0; m<MAX_AMP; m++) {
+	sparse_pe_out[m] = sparse_pe_in[m];
+    }
+
+    //gain_shape_split_vq(sparse_pe_out, aexp, aexp->vq1, weights, sparse_pe_in, &best_gain);
+    gain_shape_split_vq(sparse_pe_out, aexp, aexp->vq2, weights, sparse_pe_in, &best_gain);
+    //gain_shape_split_vq(sparse_pe_out, aexp, aexp->vq3, weights, sparse_pe_in, &best_gain);
+    //gain_shape_split_vq(sparse_pe_out, aexp, aexp->vq4, weights, sparse_pe_in, &best_gain);
+
+    for(m=0; m<MAX_AMP; m++) {
+	if (sparse_pe_in[m] != 0.0)
+	    aexp->vq_var += pow(sparse_pe_out[m] - sparse_pe_in[m], 2.0);
+    }
+    
+    /* transform quantised amps back */
+
+    for(m=1; m<=model->L; m++) {
+	index = MAX_AMP*m*model->Wo/PI;
+	assert(index < MAX_AMP);
+	amp_dB = sparse_pe_out[index];
+	model->A[m] = pow(10.0, amp_dB/20.0);
+    }
+    //exit(0);
+}
+
+
+static void vq_interp(struct AEXP *aexp, MODEL *model, int on)
+{
+    int              i, m, index;
+    float            amp_dB;
+    struct codebook *vq = aexp->vq1;
+ 
+    /* replace odd frames with interp */
+    /* once we get an even input frame we can interpolate and output odd */
+    /* use VQ to interpolate.  This assumes some correlation in adjacent VQ
+       samples */
+
+    memcpy(&aexp->model[2], model, sizeof(MODEL));
+
+    /* once we get an even input frame we have enough information to
+      replace prev odd frame with interpolated version */
+
+    if (on && ((aexp->frames % 2) == 0)) {
+
+	/* copy Wo, L, and phases */
+
+	memcpy(model, &aexp->model[1], sizeof(MODEL));
+	printf("mags: %4.2f %4.2f %4.2f Am: \n", aexp->mag[0], aexp->mag[1], aexp->mag[2]);
+
+	/* now replace Am by interpolation */
+
+	for(m=1; m<=model->L; m++) {
+	    index = MAX_AMP*m*model->Wo/PI;
+	    assert(index < MAX_AMP);
+
+	    if (index < vq->k) {
+		amp_dB  = 0.5*(aexp->mag[0] + vq->cb[vq->k * aexp->indexes[0] + index]);
+		amp_dB += 0.5*(aexp->mag[2] + vq->cb[vq->k * aexp->indexes[2] + index]);
+		printf("  %4.2f", 10.0*log10(model->A[m]));
+		//amp_dB = 10;
+		model->A[m] = pow(10.0, amp_dB/20.0);
+		printf("  %4.2f\n", 10.0*log10(model->A[m]));
+	    }
+	}
+	
+    }
+    else
+	memcpy(model, &aexp->model[1], sizeof(MODEL));
+
+    /* update memories */
+
+    for(i=0; i<2; i++) {
+	memcpy(&aexp->model[i], &aexp->model[i+1], sizeof(MODEL));
+	aexp->indexes[i] = aexp->indexes[i+1];
+	aexp->mag[i] = aexp->mag[i+1];
+    }
+}
+
+/*---------------------------------------------------------------------------*	\
 
   amp_experiment()
 
@@ -411,22 +749,31 @@ static void update_snr_calc(struct AEXP *aexp, MODEL *model, float before[])
 \*---------------------------------------------------------------------------*/
 
 void amp_experiment(struct AEXP *aexp, MODEL *model) {
-    int m;
-    float before[MAX_AMP];
+    int m,i;
 
-    for(m=1; m<=model->L; m++)
-	before[m] = model->A[m];
+    memcpy(&aexp->model_uq[2], model, sizeof(MODEL));
+
+    //for(m=1; m<=aexp->model[1].L; m++)
+    //	before[m] = aexp->model[1].A[m];
 
-    //print_sparse_pred_error(aexp, model, 40.0);
     //add_quant_noise(aexp, model, model->L/2, model->L, 3);
-    sparse_vq_pred_error(aexp, model);
+    //sparse_vq_pred_error(aexp, model);
+
+    //print_sparse_amp_error(aexp, model, 1);
+    sparse_vq_amp(aexp, model);
+    vq_interp(aexp, model, 1);
+
+    //gain_shape_sparse_vq_amp(aexp, model);
 
-    update_snr_calc(aexp, model, before);
+    printf("%d %d %f %f %f %f\n", model->L, aexp->model[0].L, model->Wo, aexp->model[0].Wo, model->A[1], aexp->model[0].A[1]);
+    update_snr_calc(aexp, &aexp->model_uq[1], model);
 
     /* update states */
 
     for(m=1; m<=model->L; m++)
 	aexp->A_prev[m] = model->A[m];	    
     aexp->frames++;
+    for(i=0; i<3; i++)
+	aexp->model_uq[i] = aexp->model_uq[i+1];
 }
 
diff --git a/codec2-dev/src/phase.c b/codec2-dev/src/phase.c
index b3d9e216..1442c152 100644
--- a/codec2-dev/src/phase.c
+++ b/codec2-dev/src/phase.c
@@ -960,7 +960,7 @@ static void sparse_vq_pred_error(struct PEXP     *pexp,
 	assert(index < MAX_AMP);
 	sparse_pe_in[index].real = cos(error);
 	sparse_pe_in[index].imag = sin(error);
-	//sparse_pe_out[index] = sparse_pe_in[index];
+	sparse_pe_out[index] = sparse_pe_in[index];
 	weights[index] = model->A[i];
 	printf("%d ", index);
     }
@@ -969,9 +969,9 @@ static void sparse_vq_pred_error(struct PEXP     *pexp,
         
     split_vq(sparse_pe_out, pexp, pexp->vq1, weights, sparse_pe_in);
     split_vq(sparse_pe_out, pexp, pexp->vq2, weights, sparse_pe_in);
-    split_vq(sparse_pe_out, pexp, pexp->vq3, weights, sparse_pe_in);
-    split_vq(sparse_pe_out, pexp, pexp->vq4, weights, sparse_pe_in);
-
+    //split_vq(sparse_pe_out, pexp, pexp->vq3, weights, sparse_pe_in);
+    //split_vq(sparse_pe_out, pexp, pexp->vq4, weights, sparse_pe_in);
+    
     /* transform quantised phases back */
 
     for(i=1; i<=model->L; i++) {
@@ -1130,12 +1130,9 @@ void phase_experiment(struct PEXP *pexp, MODEL *model) {
     //print_sparse_pred_error(pexp, model, 1, model->L, 40.0);
 
     sparse_vq_pred_error(pexp, model);
+ 
     //quant_phases(model, model->L/4+1, model->L, 3);
 
-    //sparse_vq_pred_error(pexp, model,  model->L/4+1,   model->L/2, pexp->vq2);
-    //sparse_vq_pred_error(pexp, model, model->L/2+1, 3*model->L/4, pexp->vq3);
-    //sparse_vq_pred_error(pexp, model, 3*model->L/4+21, model->L  , pexp->vq4);
-
     //predict_phases1(pexp, model, 1, model->L/4);
     //quant_phases(model, model->L/4+1, model->L, 3);