about to clean up after a few weeks of amplitude and phase experiments
authordrowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Sat, 25 Aug 2012 04:35:58 +0000 (04:35 +0000)
committerdrowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Sat, 25 Aug 2012 04:35:58 +0000 (04:35 +0000)
git-svn-id: https://svn.code.sf.net/p/freetel/code@640 01035d8c-6547-0410-b346-abe4f91aad63

codec2-dev/src/ampexp.c
codec2-dev/src/c2sim.c
codec2-dev/src/phase.c
codec2-dev/src/phase.h
codec2-dev/src/quantise.c

index 8920cb36bf0d8f7a1fc3dd1c41cecceb6d4e61f8..4d3b70a854f469c9285f643a6960597f19c0ef5b 100644 (file)
@@ -162,14 +162,21 @@ struct AEXP *amp_experiment_create() {
     aexp->vq_var = 0.0;
     aexp->vq_var_n = 0;
 
-    aexp->vq1 = load("amp_1_20_1024.txt");
+    //aexp->vq1 = load("amp_1_80_1024a.txt");
+    //aexp->vq1 = load("../unittest/st1_10_1024.txt");
+    //aexp->vq1 = load("../unittest/amp41_80_1024.txt");
+    //aexp->vq1->offset = 40;
+    aexp->vq1 = load("../unittest/amp1_10_1024.txt");
     aexp->vq1->offset = 0;
-    aexp->vq2 = load("amp_21_40_1024.txt");
-    aexp->vq2->offset = 20;
-    aexp->vq3 = load("amp_41_60_1024.txt");
-    aexp->vq3->offset = 40;
-    aexp->vq4 = load("amp_61_80_128.txt");
-    aexp->vq4->offset = 60;
+    aexp->vq2 = load("../unittest/amp11_20_1024.txt");
+    aexp->vq2->offset = 10;
+
+    aexp->vq3 = load("../unittest/amp21_40_1024.txt");
+    aexp->vq3->offset = 20;
+    aexp->vq4 = load("../unittest/amp41_60_1024.txt");
+    aexp->vq4->offset = 40;
+    aexp->vq5 = load("../unittest/amp61_80_256.txt");
+    aexp->vq5->offset = 60;
 
     #ifdef CAND2_GS
     //aexp->vq1 = load("../unittest/t1_amp1_20_1024.txt");
@@ -327,17 +334,47 @@ static void print_sparse_pred_error(struct AEXP *aexp, MODEL *model, float mag_t
 }
 
 
-static void print_sparse_amp_error(struct AEXP *aexp, MODEL *model, float edB_thresh)
-{
-    int    m, index;
-    float  e, edB, enormdB, error, dWo;
-    float  sparse_pe[MAX_AMP];
+static float frame_energy(MODEL *model, float *enormdB) {
+    int   m;
+    float e, edB;
 
     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 */
+
+    #define VER_E0
+
+    #ifdef VER_E0
+    *enormdB = 10*log10(e/model->L); /* make high and low pitches have similar amps */
+    #endif
+
+    #ifdef VER_E1
+    e = 0.0;
+    for(m=1; m<=model->L; m++)
+       e += 10*log10(model->A[m]*model->A[m]);
+    *enormdB = e;
+    #endif
+
+    #ifdef VER_E2
+    e = 0.0;
+    for(m=1; m<=model->L; m++)
+       e += 10*log10(model->A[m]*model->A[m]);
+    *enormdB = e/model->L;
+    #endif
+    //printf("%f\n", enormdB);
+
+    return edB;
+}
+
+static void print_sparse_amp_error(struct AEXP *aexp, MODEL *model, float edB_thresh)
+{
+    int    m, index;
+    float  edB, enormdB, error, dWo, Am;
+    float  sparse_pe[MAX_AMP];
+
+    edB = frame_energy(model, &enormdB);
+    //printf("%f\n", enormdB);
     dWo = fabs((aexp->model_uq[2].Wo - aexp->model_uq[1].Wo)/aexp->model_uq[2].Wo);
     
     if ((edB > edB_thresh) && (dWo < 0.1)) {
@@ -358,7 +395,7 @@ static void print_sparse_amp_error(struct AEXP *aexp, MODEL *model, float edB_th
        
        for(m=0; m<MAX_AMP; m++)
            printf("%f ", sparse_pe[m]);
-           printf("\n");
+       printf("\n");
     }
 }
 
@@ -401,13 +438,14 @@ static int split_vq(float sparse_pe_out[], struct AEXP *aexp, struct codebook *v
     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);
+    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);
-           sparse_pe_out[j] = vq->cb[vq->k * vq_ind + i];
+           printf("%d ", j);
+           sparse_pe_in[j]  -= vq->cb[vq->k * vq_ind + i];
+           sparse_pe_out[j] += vq->cb[vq->k * vq_ind + i];
            non_zero++;
        }
     }
@@ -421,19 +459,19 @@ static void sparse_vq_pred_error(struct AEXP *aexp,
 )
 {
     int    m, index;
-    float  error, amp_dB, mag;
+    float  error, amp_dB, edB, enormdB;
     float  sparse_pe_in[MAX_AMP];
     float  sparse_pe_out[MAX_AMP];
     float  weights[MAX_AMP];
 
+    edB = frame_energy(model, &enormdB);
+
     for(m=0; m<MAX_AMP; m++) {
        sparse_pe_in[m] = 0.0;
        sparse_pe_out[m] = 0.0;
     }
 
-    mag = 0.0;
     for(m=1; m<=model->L; m++) {
-       mag += model->A[m]*model->A[m];
        assert(model->A[m] > 0.0);
        error = PRED_COEFF*20.0*log10(aexp->A_prev[m]) - 20.0*log10(model->A[m]);
 
@@ -442,7 +480,6 @@ static void sparse_vq_pred_error(struct AEXP *aexp,
        sparse_pe_in[index] = error;
        weights[index] = model->A[m];
     }
-    mag = 10.0*log10(mag/model->L);
 
     /* vector quantise */
         
@@ -464,7 +501,7 @@ static void sparse_vq_pred_error(struct AEXP *aexp,
     }
     #endif
 
-    if (mag > -100.0)
+    if (edB > -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);
@@ -484,19 +521,29 @@ static void sparse_vq_pred_error(struct AEXP *aexp,
 }
 
 
+static void split_error(struct AEXP *aexp, struct codebook *vq, float sparse_pe_in[], int ind)
+{
+    int i, j;
+
+    for(i=0, j=vq->offset; i<vq->k; i++,j++) {
+       if (sparse_pe_in[j] != 0.0) {
+           sparse_pe_in[j] -= vq->cb[vq->k * ind + i];
+       }
+    }
+}
+
+
 static void sparse_vq_amp(struct AEXP *aexp, MODEL *model)
 {
     int    m, index;
-    float  error, amp_dB, mag;
+    float  error, amp_dB, edB, enormdB;
     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;
+    edB = frame_energy(model, &enormdB);
+
+    aexp->mag[2] = enormdB;
    
     for(m=0; m<MAX_AMP; m++) {
        sparse_pe_in[m] = 0.0;
@@ -505,12 +552,12 @@ static void sparse_vq_amp(struct AEXP *aexp, MODEL *model)
 
     for(m=1; m<=model->L; m++) {
        assert(model->A[m] > 0.0);
-       error = 20.0*log10(model->A[m]) - mag;
+       error = 20.0*log10(model->A[m]) - enormdB;
 
        index = MAX_AMP*m*model->Wo/PI;
        assert(index < MAX_AMP);
        sparse_pe_in[index] = error;
-       weights[index] = model->A[m];
+       weights[index] = pow(model->A[m],0.8);
     }
 
     /* vector quantise */
@@ -518,23 +565,25 @@ static void sparse_vq_amp(struct AEXP *aexp, MODEL *model)
     for(m=0; m<MAX_AMP; m++) {
        sparse_pe_out[m] = sparse_pe_in[m];
     }
-
-    //#define SIM_VQ
-    #ifndef SIM_VQ
+    
+    for(m=0; m<80; m++)
+       sparse_pe_out[m] = 0;
+    
+    #define SPLIT
+    #ifdef SPLIT
     aexp->indexes[0][2] = split_vq(sparse_pe_out, aexp, aexp->vq1, weights, sparse_pe_in);
+    
     aexp->indexes[1][2] = split_vq(sparse_pe_out, aexp, aexp->vq2, weights, sparse_pe_in);
     aexp->indexes[2][2] = split_vq(sparse_pe_out, aexp, aexp->vq3, weights, sparse_pe_in);
     aexp->indexes[3][2] = 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;
-       }
-    }
+    aexp->indexes[4][2] = split_vq(sparse_pe_out, aexp, aexp->vq5, weights, sparse_pe_in);
+    #endif
+    //#define MULTISTAGE
+    #ifdef MULTISTAGE
+    aexp->indexes[0][2] = split_vq(sparse_pe_out, aexp, aexp->vq1, weights, sparse_pe_in);
+    aexp->indexes[1][2] = split_vq(sparse_pe_out, aexp, aexp->vq2, weights, sparse_pe_in);
+    aexp->indexes[2][2] = split_vq(sparse_pe_out, aexp, aexp->vq3, weights, sparse_pe_in);
+    //aexp->indexes[3][2] = split_vq(sparse_pe_out, aexp, aexp->vq4, weights, sparse_pe_in);
     #endif
 
     for(m=0; m<MAX_AMP; m++) {
@@ -547,7 +596,7 @@ static void sparse_vq_amp(struct AEXP *aexp, MODEL *model)
     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;
+       amp_dB = sparse_pe_out[index] + enormdB;
        model->A[m] = pow(10.0, amp_dB/20.0);
     }
     //exit(0);
@@ -758,6 +807,7 @@ static void vq_interp(struct AEXP *aexp, MODEL *model, int on)
        interp_split_vq(sparse_pe_out, aexp, aexp->vq2, sparse_pe_in, 1);
        interp_split_vq(sparse_pe_out, aexp, aexp->vq3, sparse_pe_in, 2);
        interp_split_vq(sparse_pe_out, aexp, aexp->vq4, sparse_pe_in, 3);
+       interp_split_vq(sparse_pe_out, aexp, aexp->vq5, sparse_pe_in, 4);
 
        for(m=1; m<=model->L; m++) {
            index = MAX_AMP*m*model->Wo/PI;
@@ -797,6 +847,172 @@ static void vq_interp(struct AEXP *aexp, MODEL *model, int on)
            aexp->indexes[j][i] = aexp->indexes[j][i+1];
        aexp->mag[i] = aexp->mag[i+1];
     }
+
+}
+
+
+/*
+  This functions tests theory that some bands can be combined together
+  due to less frequency resolution at higher frequencies.  This will
+  reduce the amount of information we need to encode.
+*/
+
+void smooth_samples(struct AEXP *aexp, MODEL *model, int mode)
+{
+    int    m, i, j, index, step, nav, v, en;
+    float  sparse_pe_in[MAX_AMP], av, amp_dB;
+    float  sparse_pe_out[MAX_AMP];
+    float  smoothed[MAX_AMP], smoothed_out[MAX_AMP];
+    float  weights[MAX_AMP];
+    float  edB, enormdB;
+
+    edB = frame_energy(model, &enormdB);
+    
+    for(m=0; m<MAX_AMP; m++) {
+       sparse_pe_in[m] = 0.0;
+       sparse_pe_out[m] = 0.0;
+    }
+
+    /* set up sparse array */
+
+    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_out[index] = sparse_pe_in[index] = 20.0*log10(model->A[m]) - enormdB;
+    }
+
+    /* now combine samples at high frequencies to reduce dimension */
+
+    step=4;
+    for(i=MAX_AMP/2,v=0; i<MAX_AMP; i+=step,v++) {
+
+       /* average over one band */
+
+       av = 0.0; nav = 0;
+       en = i+step;
+       if (en > (MAX_AMP-1))
+           en = MAX_AMP-1;
+       for(j=i; j<en; j++) {
+           if (sparse_pe_in[j] != 0.0) {
+               av += sparse_pe_in[j];
+               nav++;
+           }
+       }
+       if (nav) {
+           av /= nav;
+           smoothed[v] = av;
+           weights[v] = pow(10.0,av/20.0);
+           //weights[v] = 1.0;
+       }
+       else
+           smoothed[v] = 0.0;
+
+    }
+
+    if (mode == 1) {
+       for(i=0; i<v; i++)
+           printf("%5.2f ", smoothed[i]);
+       printf("\n");
+    }
+
+    if (mode == 2) {
+       for(i=0; i<v; i++)
+           smoothed_out[i] = 0;
+       split_vq(smoothed_out, aexp, aexp->vq1, weights, smoothed);
+       for(i=0; i<v; i++)
+           smoothed[i] = smoothed_out[i];
+    }
+       
+    /* set all samples to smoothed average */
+
+    step = 4;
+    for(i=MAX_AMP/2,v=0; i<MAX_AMP; i+=step,v++) {
+       en = i+step;
+       if (en > (MAX_AMP-1))
+           en = MAX_AMP-1;
+       for(j=i; j<en; j++)
+           sparse_pe_out[j] = smoothed[v];
+    }
+         
+    /* convert back to Am */
+    
+    for(m=1; m<=model->L; m++) {
+       index = MAX_AMP*m*model->Wo/PI;
+       assert(index < MAX_AMP);
+       amp_dB = sparse_pe_out[index] + enormdB;
+       //printf("%d %4.2f %4.2f\n", m, 10.0*log10(model->A[m]), amp_dB);
+       model->A[m] = pow(10.0, amp_dB/20.0);
+    }
+    
+}
+
+#define MAX_BINS 40
+static float bins[] = {
+    /*1000.0, 1200.0, 1400.0, 1600.0, 1800,*/
+    2000.0, 2400.0, 2800.0,
+    3000.0, 3400.0, 3600.0, 4000.0};
+
+void smooth_amp(struct AEXP *aexp, MODEL *model) {
+    int    m, i;
+    int    nbins;
+    int    b;
+    float  f;
+    float  av[MAX_BINS];
+    int    nav[MAX_BINS];
+
+    nbins = sizeof(bins)/sizeof(float);
+
+    /* clear all bins */
+
+    for(i=0; i<MAX_BINS; i++) {
+       av[i] = 0.0;
+       nav[i] = 0;
+    }
+
+    /* add amps into each bin */
+
+    for(m=1; m<=model->L; m++) {
+       f = m*model->Wo*FS/TWO_PI;
+       if (f > bins[0]) {
+
+           /* find bin  */
+
+           for(i=0; i<nbins; i++)
+               if ((f > bins[i]) && (f <= bins[i+1]))
+                   b = i;
+           assert(b < MAX_BINS);
+
+           av[b] += model->A[m]*model->A[m];
+           nav[b]++;
+       }
+           
+    }
+
+    /* use averages to est amps */
+
+    for(m=1; m<=model->L; m++) {
+       f = m*model->Wo*FS/TWO_PI;
+       if (f > bins[0]) {
+
+           /* find bin */
+
+           for(i=0; i<nbins; i++)
+               if ((f > bins[i]) && (f <= bins[i+1]))
+                   b = i;
+           assert(b < MAX_BINS);
+
+           /* add predicted phase error to this bin */
+
+           printf("L %d m %d f %4.f b %d\n", model->L, m, f, b);
+
+           printf(" %d: %4.3f -> ", m, 20*log10(model->A[m])); 
+           model->A[m] = sqrt(av[b]/nav[b]);
+           printf("%4.3f\n", 20*log10(model->A[m])); 
+       }
+    }
+    printf("\n");
 }
 
 /*---------------------------------------------------------------------------* \
@@ -809,16 +1025,18 @@ static void vq_interp(struct AEXP *aexp, MODEL *model, int on)
 
 void amp_experiment(struct AEXP *aexp, MODEL *model, char *arg) {
     int m,i;
-
+    
     memcpy(&aexp->model_uq[2], model, sizeof(MODEL));
 
-    //add_quant_noise(aexp, model, model->L/2, model->L, 3);
-    //sparse_vq_pred_error(aexp, model);
+    if (strcmp(arg, "qn") == 0) {
+       add_quant_noise(aexp, model, 1, model->L, 1);
+       update_snr_calc(aexp, &aexp->model_uq[2], model);
+   }
 
-    /* print training samplesthat can be > train.txt for training VQ */
+    /* print training samples that can be > train.txt for training VQ */
 
     if (strcmp(arg, "train") == 0) 
-       print_sparse_amp_error(aexp, model, 50.0);
+       print_sparse_amp_error(aexp, model, 00.0);
 
     /* VQ of amplitudes, no interpolation (ie 10ms rate) */
 
@@ -838,12 +1056,31 @@ void amp_experiment(struct AEXP *aexp, MODEL *model, char *arg) {
 
     /* gain/shape VQ of amplitudes, 10ms rate (doesn't work that well) */
 
-   if (strcmp(arg, "gsvq") == 0) {
+    if (strcmp(arg, "gsvq") == 0) {
        gain_shape_sparse_vq_amp(aexp, model);
        vq_interp(aexp, model, 0);
        update_snr_calc(aexp, &aexp->model_uq[1], model);
     }
 
+    if (strcmp(arg, "smooth") == 0) {
+       smooth_samples(aexp, model, 0);
+       update_snr_calc(aexp, &aexp->model_uq[2], model);
+    }
+
+    if (strcmp(arg, "smoothtrain") == 0) {
+       smooth_samples(aexp, model, 1);
+       //update_snr_calc(aexp, &aexp->model_uq[2], model);
+    }
+
+    if (strcmp(arg, "smoothvq") == 0) {
+       smooth_samples(aexp, model, 2);
+       update_snr_calc(aexp, &aexp->model_uq[2], model);
+    }
+
+    if (strcmp(arg, "smoothamp") == 0) {
+       smooth_amp(aexp, model);
+       update_snr_calc(aexp, &aexp->model_uq[2], model);
+    }
 
     /* update states */
 
index 26531b11a56c8dc9782252a015ac2b23c1636e49..434f4f399a43a34c379d4b1f06440f95a96146aa 100644 (file)
@@ -79,6 +79,7 @@ int main(int argc, char *argv[])
     int   voiced1 = 0;
     char  out_file[MAX_STR];
     char  ampexp_arg[MAX_STR];
+    char  phaseexp_arg[MAX_STR];
     float snr;
     float sum_snr;
 
@@ -97,7 +98,7 @@ int main(int argc, char *argv[])
     int   postfilt;
     float bg_est;
 
-    int   hand_voicing = 0, phaseexp = 0, ampexp = 0;
+    int   hand_voicing = 0, phaseexp = 0, ampexp = 0, hi =0;
     FILE *fvoicing = 0;
 
     MODEL prev_model, interp_model;
@@ -136,12 +137,13 @@ int main(int argc, char *argv[])
         { "lspdt_mode", required_argument, NULL, 0 },
         { "lspjvm", no_argument, &lspjvm, 1 },
         { "phase0", no_argument, &phase0, 1 },
-        { "phaseexp", no_argument, &phaseexp, 1 },
+        { "phaseexp", required_argument, &phaseexp, 1 },
         { "ampexp", required_argument, &ampexp, 1 },
         { "postfilter", no_argument, &postfilt, 1 },
         { "hand_voicing", required_argument, &hand_voicing, 1 },
         { "dec", no_argument, &decimate, 1 },
         { "dt", no_argument, &dt, 1 },
+        { "hi", no_argument, &hi, 1 },
         { "dump_pitch_e", required_argument, &dump_pitch_e, 1 },
         { "sq_pitch_e", no_argument, &scalar_quant_Wo_e, 1 },
         { "vq_pitch_e", no_argument, &vector_quant_Wo_e, 1 },
@@ -236,6 +238,8 @@ int main(int argc, char *argv[])
                        optarg, strerror(errno));
                     exit(1);
                 }
+           } else if(strcmp(long_options[option_index].name, "phaseexp") == 0) {
+               strcpy(phaseexp_arg, optarg);
            } else if(strcmp(long_options[option_index].name, "ampexp") == 0) {
                strcpy(ampexp_arg, optarg);
            } else if(strcmp(long_options[option_index].name, "rate") == 0) {
@@ -364,19 +368,26 @@ int main(int argc, char *argv[])
        dump_Sn(Sn); dump_Sw(Sw); dump_model(&model);
         #endif
 
+       if (ampexp)
+           amp_experiment(aexp, &model, ampexp_arg);
+
        if (phaseexp) {
             #ifdef DUMP
            dump_phase(&model.phi[0], model.L);
             #endif
-           phase_experiment(pexp, &model);
+           phase_experiment(pexp, &model, phaseexp_arg);
             #ifdef DUMP
            dump_phase_(&model.phi[0], model.L);
             #endif
        }
 
-       if (ampexp)
-           amp_experiment(aexp, &model, ampexp_arg);
-           
+       if (hi) {
+           int m;
+           for(m=1; m<model.L/2; m++)
+               model.A[m] = 0.0;
+           for(m=3*model.L/4; m<=model.L; m++)
+               model.A[m] = 0.0;
+       }
 
        /*------------------------------------------------------------*\
 
index 1442c1528f949ad85c903d9ee64ef44f7e4bfe1a..f16fb852acde2e4e1bc2d4ffafde4920f57fd3bf 100644 (file)
@@ -217,7 +217,7 @@ void phase_synth_zero_order(
      I found that using just this frame's Wo improved quality for UV
      sounds compared to interpolating two frames Wo like this:
      
-     ex_phase[0] += (*prev_Wo+mode->Wo)*N/2;
+     ex_phase[0] += (*prev_Wo+model->Wo)*N/2;
   */
   
   ex_phase[0] += (model->Wo)*N;
@@ -225,19 +225,22 @@ void phase_synth_zero_order(
   r = TWO_PI/GLOTTAL_FFT_SIZE;
 
   for(m=1; m<=model->L; m++) {
-
+      
     /* generate excitation */
-
+           
     if (model->voiced) {
-       /* I think adding a little jitter helps improve low pitch
-          males like hts1a. This moves the onset of each harmonic
-          over at +/- 0.25 of a sample.
-       */
-        jitter = 0.25*(1.0 - 2.0*rand()/RAND_MAX);
+       float offs;
         b = floor(m*model->Wo/r + 0.5);
        if (b > ((GLOTTAL_FFT_SIZE/2)-1)) {
                b = (GLOTTAL_FFT_SIZE/2)-1;
        }
+
+       /* I think adding a little jitter helps improve low pitch
+          males like hts1a. This moves the onset of each harmonic
+          over +/- 0.25 of a sample.
+       */
+       jitter = 0.25*(1.0 - 2.0*rand()/RAND_MAX);
+
        Ex[m].real = cos(ex_phase[0]*m - jitter*model->Wo*m + glottal[b]);
        Ex[m].imag = sin(ex_phase[0]*m - jitter*model->Wo*m + glottal[b]);
     }
@@ -345,15 +348,18 @@ static struct codebook *load(const char * name)
 /* states for phase experiments */
 
 struct PEXP {
+    float            phi1;
     float            phi_prev[MAX_AMP];
     float            Wo_prev;
     int              frames;
     float            snr;
     float            var;
     int              var_n;
-    struct codebook *vq1,*vq2,*vq3,*vq4; 
+    struct codebook *vq1,*vq2,*vq3,*vq4,*vq5
     float            vq_var;
     int              vq_var_n;
+    MODEL            prev_model;
+    int              state;
 };
 
 
@@ -372,6 +378,7 @@ struct PEXP * phase_experiment_create() {
     pexp = (struct PEXP *)malloc(sizeof(struct PEXP));
     assert (pexp != NULL);
 
+    pexp->phi1 = 0;
     for(i=0; i<MAX_AMP; i++)
        pexp->phi_prev[i] = 0.0;
     pexp->Wo_prev = 0.0;
@@ -380,20 +387,30 @@ struct PEXP * phase_experiment_create() {
     pexp->var = 0.0;
     pexp->var_n = 0;
     
+    /* smoothed 10th order for 1st 1 khz */
+    //pexp->vq1 = load("../unittest/ph1_10_1024.txt");
+    //pexp->vq1->offset = 0;
+
     /* load experimental phase VQ */
 
-    pexp->vq1 = load("../unittest/testn1_20_1024.txt");
-    pexp->vq2 = load("../unittest/testn21_40_1024.txt");
-    pexp->vq3 = load("../unittest/testn41_60_32.txt");
-    pexp->vq4 = load("../unittest/testn61_80_32.txt");
+    //pexp->vq1 = load("../unittest/testn1_20_1024.txt");
+    pexp->vq1 = load("../unittest/test.txt");
+    //pexp->vq2 = load("../unittest/testn21_40_1024.txt");
+    pexp->vq2 = load("../unittest/test11_20_1024.txt");
+    pexp->vq3 = load("../unittest/test21_30_1024.txt");
+    pexp->vq4 = load("../unittest/test31_40_1024.txt");
+    pexp->vq5 = load("../unittest/test41_60_1024.txt");
     pexp->vq1->offset = 0;
-    pexp->vq2->offset = 20;
-    pexp->vq3->offset = 40;
-    pexp->vq4->offset = 60;
+    pexp->vq2->offset = 10;
+    pexp->vq3->offset = 20;
+    pexp->vq4->offset = 30;
+    pexp->vq5->offset = 40;
 
     pexp->vq_var = 0.0;
     pexp->vq_var_n = 0;
 
+    pexp->state = 0;
+
     return pexp;
 }
 
@@ -440,6 +457,7 @@ static void bubbleSort(struct AMPINDEX numbers[], int array_size)
   {
     for (j = 1; j <= i; j++)
     {
+       //printf("i %d j %d %f %f \n", i, j, numbers[j-1].amp, numbers[j].amp);
       if (numbers[j-1].amp < numbers[j].amp)
       {
         temp = numbers[j-1];
@@ -477,17 +495,83 @@ static void predict_phases(struct PEXP *pexp, MODEL *model, int start, int end)
     int i;
 
     for(i=start; i<=end; i++) {
-       model->phi[i] = pexp->phi_prev[i] + N*i*(model->Wo + pexp->Wo_prev)/2.0;
+       model->phi[i] = pexp->phi_prev[i] + N*i*model->Wo;
+    }
+   
+}
+static float refine_Wo(struct PEXP     *pexp, 
+                      MODEL           *model, 
+                      int              start, 
+                      int              end);
+
+/* Fancy state based phase prediction.  Actually works OK on most utterances,
+   but could use some tuning.  Breaks down a bit on mmt1. */
+
+static void predict_phases_state(struct PEXP *pexp, MODEL *model, int start, int end) {
+    int i, next_state;
+    float best_Wo, dWo;
+
+    //best_Wo = refine_Wo(pexp, model, start, end);
+    //best_Wo = (model->Wo + pexp->Wo_prev)/2.0;
+    best_Wo = model->Wo;
+
+    dWo = fabs(model->Wo - pexp->Wo_prev)/model->Wo;
+    next_state = pexp->state;
+    switch(pexp->state) {
+    case 0:
+       if (dWo < 0.1) {
+           /* UV -> V transition, so start with phases in lock.  They will
+              drift a bit over voiced track which is kinda what we want, so
+              we don't get clicky speech.
+           */
+           next_state = 1;
+           for(i=start; i<=end; i++)
+               pexp->phi_prev[i] = i*pexp->phi1;
+       }
+       
+       break;
+    case 1:
+       if (dWo > 0.1)
+           next_state = 0;
+       break;
     }
+    pexp->state = next_state;
+
+    if (pexp->state == 0)
+       for(i=start; i<=end; i++) {
+           model->phi[i] = PI*(1.0 - 2.0*rand()/RAND_MAX);
+       }
+    else
+       for(i=start; i<=end; i++) {
+           model->phi[i] = pexp->phi_prev[i] + N*i*best_Wo;
+       }
+    printf("state %d\n", pexp->state);
 }
 
 static void struct_phases(struct PEXP *pexp, MODEL *model, int start, int end) {
     int i;
 
-    model->phi[1] = pexp->phi_prev[1] + N*(model->Wo + pexp->Wo_prev)/2.0;
-
     for(i=start; i<=end; i++)
-       model->phi[i] = model->phi[1]*i;
+       model->phi[i] = pexp->phi1*i;
+   
+}
+
+
+static void predict_phases2(struct PEXP *pexp, MODEL *model, int start, int end) {
+    int i;
+    float pred, str, diff;
+
+    for(i=start; i<=end; i++) {
+       pred = pexp->phi_prev[i] + N*i*model->Wo;
+       str = pexp->phi1*i;
+       diff = str - pred;
+       diff = atan2(sin(diff), cos(diff));
+       if (diff > 0)
+           pred += PI/16;
+       else
+           pred -= PI/16;
+       model->phi[i] = pred;
+    }
    
 }
 
@@ -576,17 +660,6 @@ static void check_phase_quant(MODEL *model, float tol)
 }
 
 
-static void repeat_phases(MODEL *model, int period)
-{
-    int m,i;
-
-    for(m=period+1,i=1; m<=model->L; m++,i++) {
-       model->phi[m] = model->phi[m-1] + model->phi[i+1] - model->phi[i];
-    }
-
-}
-
-
 static float est_phi1(MODEL *model, int start, int end)
 {
     int m;
@@ -660,29 +733,48 @@ static void top_amp(struct PEXP *pexp, MODEL *model, int start, int end, int n_h
     int top, i, j;
     struct AMPINDEX sorted[MAX_AMP];
 
-    /* sort into acending order of amplitude */
+    /* sort into ascending order of amplitude */
 
+    printf("\n");
     for(i=start,j=0; i<end; i++,j++) {
        sorted[j].amp = model->A[i];
        sorted[j].index = i;
+       printf("%f ", model->A[i]);
     }
-    bubbleSort(&sorted[1], end-start);
+    bubbleSort(sorted, end-start);
+
+    printf("\n");
+    for(j=0; j<n_harm; j++)
+       printf("%d %f\n", j, sorted[j].amp);
 
     /* keep phase of top n_harm, predict others */
 
     for(i=start; i<end; i++) {         
        top = 0;
-       for(j=0; j<n_harm; j++)
+       for(j=0; j<n_harm; j++) {
            if (model->A[i] == sorted[j].amp) {
                top = 1;
                assert(i == sorted[j].index);
            }
-               
+       }
+       
+       #define ALTTOP
+       #ifdef ALTTOP
+       model->phi[i] = 0.0; /* make sure */
+       if (top) {
+           model->phi[i] = i*pexp->phi1;
+           removed++;
+       }
+       else {
+           model->phi[i] = PI*(1.0 - 2.0*(float)rand()/RAND_MAX); // note: try rand for higher harms
+           removed++;
+       }
+       #else
        if (!top) {
            model->phi[i] = 0.0; /* make sure */
            if (pred)  {
-               model->phi[i] = pexp->phi_prev[i] + i*N*(model->Wo + pexp->Wo_prev)/2.0;
-               //model->phi[i] = i*model->phi[1];
+               //model->phi[i] = pexp->phi_prev[i] + i*N*(model->Wo + pexp->Wo_prev)/2.0;
+               model->phi[i] = i*model->phi[1];
            }
            else
                model->phi[i] = PI*(1.0 - 2.0*(float)rand()/RAND_MAX); // note: try rand for higher harms
@@ -693,12 +785,13 @@ static void top_amp(struct PEXP *pexp, MODEL *model, int start, int end, int n_h
            quant_phase(&model->phi[i], -PI, PI, 3);        
            not_removed++;
        }
+       #endif
     }
     printf("dim: %d rem %d not_rem %d\n", end-start, removed, not_removed);
            
 }
 
-
+    
 static void limit_prediction_error(struct PEXP *pexp, MODEL *model, int start, int end, float limit) 
 {
     int   i;
@@ -822,6 +915,16 @@ static COMP cconj(COMP a)
     return res;
 }
 
+static COMP cadd(COMP a, COMP b)
+{
+    COMP res;
+
+    res.real = a.real + b.real;
+    res.imag = a.imag + b.imag;
+
+    return res;
+}
+
 static COMP cmult(COMP a, COMP b)
 {
     COMP res;
@@ -854,7 +957,7 @@ static int vq_phase(COMP cb[], COMP vec[], float weights[], int d, int e, float
                diffp = atan2(diffr.imag, diffr.real);
                error  += diffp*diffp;
                metric += weights[i]*weights[i]*diffp*diffp;
-               //metric diffp*diffp;
+               //metric += weights[i]*diffp*diffp;
                //metric = log10(weights[i]*fabs(diffp));
                //printf("diffp %f metric %f\n", diffp, metric);
                //if (metric < log10(PI/(8.0*sqrt(3.0))))
@@ -913,14 +1016,14 @@ static void split_vq(COMP sparse_pe_out[], struct PEXP *pexp, struct codebook *v
 {
     int i, j, non_zero, vq_ind;
     
-    printf("\n offset %d k %d m %d  j: ", vq->offset, vq->k, vq->m);
+    //printf("\n offset %d k %d m %d  j: ", vq->offset, vq->k, vq->m);
     vq_ind = vq_phase(vq->cb, &sparse_pe_in[vq->offset], &weights[vq->offset], vq->k, vq->m, &pexp->vq_var);
   
     non_zero = 0;
     for(i=0, j=vq->offset; i<vq->k; i++,j++) {
        //printf("%f ", atan2(sparse_pe[i].imag, sparse_pe[i].real));
        if ((sparse_pe_in[j].real != 0.0) && (sparse_pe_in[j].imag != 0.0)) {
-           printf("%d ", j);
+           //printf("%d ", j);
            sparse_pe_out[j] = vq->cb[vq->k * vq_ind + i];
            non_zero++;
        }
@@ -951,7 +1054,7 @@ static void sparse_vq_pred_error(struct PEXP     *pexp,
        sparse_pe_out[i].imag = 0.0;
     }
 
-    printf("\n");
+    //printf("\n");
     for(i=1; i<=model->L; i++) {
        pred = pexp->phi_prev[i] + N*i*best_Wo;
        error = pred - model->phi[i];
@@ -962,15 +1065,16 @@ static void sparse_vq_pred_error(struct PEXP     *pexp,
        sparse_pe_in[index].imag = sin(error);
        sparse_pe_out[index] = sparse_pe_in[index];
        weights[index] = model->A[i];
-       printf("%d ", index);
+       //printf("%d ", index);
     }
     
     /* vector quantise */
         
     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);
+    split_vq(sparse_pe_out, pexp, pexp->vq5, weights, sparse_pe_in);
     
     /* transform quantised phases back */
 
@@ -1110,6 +1214,370 @@ static void predict_phases1(struct PEXP *pexp, MODEL *model, int start, int end)
 }
 
 
+/*
+  This functions tests theory that some bands can be combined together
+  due to less frequency resolution at higher frequencies.  This will
+  reduce the amount of information we need to encode.
+*/
+
+void smooth_phase(struct PEXP *pexp, MODEL *model, int mode)
+{
+    int    m, i, j, index, step, v, en, nav, st;
+    COMP   sparse_pe_in[MAX_AMP], av;
+    COMP   sparse_pe_out[MAX_AMP];
+    COMP   smoothed[MAX_AMP];
+    float  best_Wo, pred, err;
+    float  weights[MAX_AMP];
+    float  avw, smoothed_weights[MAX_AMP];
+    COMP   smoothed_in[MAX_AMP], smoothed_out[MAX_AMP];
+
+    best_Wo = refine_Wo(pexp, model, 1, model->L);
+
+    for(m=0; m<MAX_AMP; m++) {
+       sparse_pe_in[m].real = sparse_pe_in[m].imag = 0.0;
+       sparse_pe_out[m].real = sparse_pe_out[m].imag = 0.0;
+    }
+
+    /* set up sparse array */
+
+    for(m=1; m<=model->L; m++) {
+       pred = pexp->phi_prev[m] + N*m*best_Wo;
+       err = model->phi[m] - pred;
+       err = atan2(sin(err),cos(err));
+
+       index = MAX_AMP*m*model->Wo/PI;
+       assert(index < MAX_AMP);
+       sparse_pe_in[index].real = model->A[m]*cos(err);
+       sparse_pe_in[index].imag = model->A[m]*sin(err);
+       sparse_pe_out[index] = sparse_pe_in[index];
+       weights[index] = model->A[m];
+    }
+
+    /* now combine samples at high frequencies to reduce dimension */
+
+    step = 2;
+    st = 0;
+    for(i=st,v=0; i<MAX_AMP; i+=step,v++) {
+
+       /* average over one band */
+
+       av.real = 0.0; av.imag = 0.0; avw = 0.0; nav = 0;
+       en = i+step;
+       if (en > (MAX_AMP-1))
+           en = MAX_AMP-1;
+       for(j=i; j<en; j++) {
+           if ((sparse_pe_in[j].real != 0.0) &&(sparse_pe_in[j].imag != 0.0) ) {
+               av = cadd(av, sparse_pe_in[j]);
+               avw += weights[j];
+               nav++;
+           }
+       }
+       if (nav) {
+           smoothed[v] = av;
+           smoothed_weights[v] = avw/nav;
+       }
+       else
+           smoothed[v].real = smoothed[v].imag = 0.0;
+    }
+
+    if (mode == 2) {
+       for(i=0; i<MAX_AMP; i++) {
+           smoothed_in[i] = smoothed[i];
+           smoothed_out[i] = smoothed_in[i];
+       }
+       split_vq(smoothed_out, pexp, pexp->vq1, smoothed_weights, smoothed_in);
+       for(i=0; i<MAX_AMP; i++)
+           smoothed[i] = smoothed_out[i];
+    }
+
+    /* set all samples to smoothed average */
+
+    for(i=st,v=0; i<MAX_AMP; i+=step,v++) {
+       en = i+step;
+       if (en > (MAX_AMP-1))
+           en = MAX_AMP-1;
+       for(j=i; j<en; j++)
+           sparse_pe_out[j] = smoothed[v];
+       if (mode == 1)
+           printf("%f ", atan2(smoothed[v].imag, smoothed[v].real));
+    }
+    if (mode == 1)
+       printf("\n");
+
+    /* convert back to Am */
+    
+    for(m=1; m<=model->L; m++) {
+       index = MAX_AMP*m*model->Wo/PI;
+       assert(index < MAX_AMP);
+       pred = pexp->phi_prev[m] + N*m*best_Wo;
+       err = atan2(sparse_pe_out[index].imag, sparse_pe_out[index].real);
+       model->phi[m] = pred + err;
+    }
+    
+}
+
+/*
+  Another version of a functions that tests the theory that some bands
+  can be combined together due to less frequency resolution at higher
+  frequencies.  This will reduce the amount of information we need to
+  encode.
+*/
+
+void smooth_phase2(struct PEXP *pexp, MODEL *model) {
+    float m;
+    float step;
+    int   a,b,h,i;
+    float best_Wo, pred, err, s,c, phi1_;
+
+    best_Wo = refine_Wo(pexp, model, 1, model->L);
+
+    step = (float)model->L/30;
+    printf("\nL: %d step: %3.2f am,bm: ", model->L, step);
+    for(m=(float)model->L/4; m<=model->L; m+=step) {
+       a = floor(m);
+       b = floor(m+step);
+       if (b > model->L) b = model->L;
+       h = b-a;
+
+       printf("%d,%d,(%d)  ", a, b, h);
+       c = s = 0.0;
+       if (h>1) {
+           for(i=a; i<b; i++) {
+               pred = pexp->phi_prev[i] + N*i*best_Wo;
+               err = model->phi[i] - pred;
+               c += cos(err); s += sin(err);
+           }
+           phi1_ = atan2(s,c);
+           for(i=a; i<b; i++) {
+               pred = pexp->phi_prev[i] + N*i*best_Wo;
+               printf("%d: %4.3f -> ", i, model->phi[i]); 
+               model->phi[i] = pred + phi1_;
+               model->phi[i] = atan2(sin(model->phi[i]),cos(model->phi[i]));
+               printf("%4.3f  ", model->phi[i]); 
+           }
+       }
+    }
+}
+
+
+#define MAX_BINS 40
+//static float bins[] = {2600.0, 2800.0, 3000.0, 3200.0, 3400.0, 3600.0, 3800.0, 4000.0};
+static float bins[] = {/*
+
+                       1000.0, 1100.0, 1200.0, 1300.0, 1400.0, 
+                      1500.0, 1600.0, 1700.0, 1800.0, 1900.0,*/
+
+    2000.0, 2400.0, 2800.0,
+    3000.0, 3400.0, 3600.0, 4000.0};
+
+void smooth_phase3(struct PEXP *pexp, MODEL *model) {
+    int    m, i;
+    int   nbins;
+    int   b;
+    float f, best_Wo, pred, err;
+    COMP  av[MAX_BINS];
+
+    nbins = sizeof(bins)/sizeof(float);
+    best_Wo = refine_Wo(pexp, model, 1, model->L);
+
+    /* clear all bins */
+
+    for(i=0; i<MAX_BINS; i++) {
+       av[i].real = 0.0;
+       av[i].imag = 0.0;
+    }
+
+    /* add phases into each bin */
+
+    for(m=1; m<=model->L; m++) {
+       f = m*model->Wo*FS/TWO_PI;
+       if (f > bins[0]) {
+
+           /* find bin  */
+
+           for(i=0; i<nbins; i++)
+               if ((f > bins[i]) && (f <= bins[i+1]))
+                   b = i;
+           assert(b < MAX_BINS);
+
+           /* est predicted phase from average */
+
+           pred = pexp->phi_prev[m] + N*m*best_Wo; 
+           err = model->phi[m] - pred;
+           av[b].real += cos(err); av[b].imag += sin(err);
+       }
+           
+    }
+
+    /* use averages to est phases */
+
+    for(m=1; m<=model->L; m++) {
+       f = m*model->Wo*FS/TWO_PI;
+       if (f > bins[0]) {
+
+           /* find bin */
+
+           for(i=0; i<nbins; i++)
+               if ((f > bins[i]) && (f <= bins[i+1]))
+                   b = i;
+           assert(b < MAX_BINS);
+
+           /* add predicted phase error to this bin */
+
+           printf("L %d m %d f %4.f b %d\n", model->L, m, f, b);
+
+           pred = pexp->phi_prev[m] + N*m*best_Wo;
+           err = atan2(av[b].imag, av[b].real);
+           printf(" %d: %4.3f -> ", m, model->phi[m]); 
+           model->phi[m] = pred + err;
+           model->phi[m] = atan2(sin(model->phi[m]),cos(model->phi[m]));
+           printf("%4.3f\n", model->phi[m]); 
+       }
+    }
+    printf("\n");
+}
+
+
+/* 
+   Try to code the phase of the largest amplitude in each band.  Randomise the
+   phase of the other harmonics. The theory is that only the largest harmonic
+   will be audible.
+*/
+
+void cb_phase1(struct PEXP *pexp, MODEL *model) {
+    int   m, i;
+    int   nbins;
+    int   b;
+    float f, best_Wo;
+    float max_val[MAX_BINS];
+    int   max_ind[MAX_BINS];
+
+    nbins = sizeof(bins)/sizeof(float);
+    best_Wo = refine_Wo(pexp, model, 1, model->L);
+
+    for(i=0; i<nbins; i++)
+       max_val[i] = 0.0;
+
+    /* determine max amplitude for each bin */
+
+    for(m=1; m<=model->L; m++) {
+       f = m*model->Wo*FS/TWO_PI;
+       if (f > bins[0]) {
+
+           /* find bin  */
+
+           for(i=0; i<nbins; i++)
+               if ((f > bins[i]) && (f <= bins[i+1]))
+                   b = i;
+           assert(b < MAX_BINS);
+
+           if (model->A[m] > max_val[b]) {
+               max_val[b] = model->A[m];
+               max_ind[b] = m;
+           }
+       }
+           
+    }
+
+    /* randomise phase of other harmonics */
+
+    for(m=1; m<=model->L; m++) {
+       f = m*model->Wo*FS/TWO_PI;
+       if (f > bins[0]) {
+
+           /* find bin */
+
+           for(i=0; i<nbins; i++)
+               if ((f > bins[i]) && (f <= bins[i+1]))
+                   b = i;
+           assert(b < MAX_BINS);
+
+           if (m != max_ind[b])
+               model->phi[m] = pexp->phi_prev[m] + N*m*best_Wo;
+       }
+    }
+}
+
+
+/* 
+   Theory is only the phase of the envelope of signal matters within a
+   Critical Band. So we estimate the position of an impulse that
+   approximates the envelope of the signal.
+*/
+
+void cb_phase2(struct PEXP *pexp, MODEL *model) {
+    int   st, m, i, a, b, step;
+    float diff,w,c,s,phi1_;
+    float A[MAX_AMP];
+
+    for(m=1; m<=model->L; m++) {
+       A[m] = model->A[m];
+       model->A[m] = 0;
+    }
+
+    st = 2*model->L/4;
+    step = 3;
+    model->phi[1] = pexp->phi_prev[1] + (pexp->Wo_prev+model->Wo)*N/2.0;
+
+    printf("L=%d ", model->L);
+    for(m=st; m<st+step*2; m+=step) {
+       a = m; b=a+step;
+       if (b > model->L)
+           b = model->L;
+       
+       c = s = 0;
+       for(i=a; i<b-1; i++) {
+           printf("diff %d,%d ", i, i+1);
+           diff = model->phi[i+1] - model->phi[i];
+           //w = (model->A[i+1] + model->A[i])/2; 
+           w = 1.0;
+           c += w*cos(diff); s += w*sin(diff);
+       }
+       phi1_ = atan2(s,c);
+       printf("replacing: ");
+       for(i=a; i<b; i++) {
+           //model->phi[i] = i*phi1_;
+           //model->phi[i] = i*model->phi[1];
+           //model->phi[i] = m*(pexp->Wo_prev+model->Wo)*N/2.0;
+           model->A[m] = A[m];
+           printf("%d ", i);
+       }
+       printf(" . ");
+    }
+    printf("\n");
+}
+
+
+static void smooth_phase4(MODEL *model) {
+    int    m;
+    float  phi_m, phi_m_1;
+
+    if (model->L > 25) {
+       printf("\nL %d\n", model->L);
+       for(m=model->L/2; m<=model->L; m+=2) {
+           if ((m+1) <= model->L) {
+               phi_m   = (model->phi[m] - model->phi[m+1])/2.0;
+               phi_m_1 = (model->phi[m+1] - model->phi[m])/2.0;
+               model->phi[m] = phi_m;
+               model->phi[m+1] = phi_m_1;
+               printf("%d %4.3f %4.3f  ", m, phi_m, phi_m_1);
+           }
+       }
+    }
+
+}
+
+/* try repeating last frame, just advance phases to account for time shift */
+
+static void repeat_phases(struct PEXP *pexp, MODEL *model) {
+    int m;
+
+    *model = pexp->prev_model;
+    for(m=1; m<=model->L; m++)
+       model->phi[m] += N*m*model->Wo;
+
+}
+
 /*---------------------------------------------------------------------------*\
 
   phase_experiment()
@@ -1118,29 +1586,93 @@ static void predict_phases1(struct PEXP *pexp, MODEL *model, int start, int end)
 
 \*---------------------------------------------------------------------------*/
 
-void phase_experiment(struct PEXP *pexp, MODEL *model) {
+void phase_experiment(struct PEXP *pexp, MODEL *model, char *arg) {
     int              m;
-    float            before[MAX_AMP];
+    float            before[MAX_AMP], best_Wo;
 
     assert(pexp != NULL);
     memcpy(before, &model->phi[0], sizeof(float)*MAX_AMP);
 
-    //quant_phases(model, 1, model->L, 3);
-    //update_variance_calc(pexp, model, before);
-    //print_sparse_pred_error(pexp, model, 1, model->L, 40.0);
+    if (strcmp(arg,"q3") == 0) { 
+       quant_phases(model, 1, model->L, 3);
+       update_snr_calc(pexp, model, before);
+       update_variance_calc(pexp, model, before);
+    }
 
-    sparse_vq_pred_error(pexp, model);
-    //quant_phases(model, model->L/4+1, model->L, 3);
+    if (strcmp(arg,"dec2") == 0) {
+       if ((pexp->frames % 2) != 0) {
+           predict_phases(pexp, model, 1, model->L);   
+           update_snr_calc(pexp, model, before);
+           update_variance_calc(pexp, model, before);
+       }
+    }
 
-    //predict_phases1(pexp, model, 1, model->L/4);
-    //quant_phases(model, model->L/4+1, model->L, 3);
+    if (strcmp(arg,"repeat") == 0) {
+       if ((pexp->frames % 2) != 0) {
+           repeat_phases(pexp, model); 
+           update_snr_calc(pexp, model, before);
+           update_variance_calc(pexp, model, before);
+       }
+    }
 
-    update_snr_calc(pexp, model, before);
-    update_variance_calc(pexp, model, before);
-    
-    //print_pred_error_sparse_wo_correction1(pexp, model,1, model->L, -100.0);
-    //print_sparse_pred_error(pexp, model, 1, model->L, 40.0);
+    if (strcmp(arg,"vq") == 0) {
+       sparse_vq_pred_error(pexp, model);
+       update_snr_calc(pexp, model, before);
+       update_variance_calc(pexp, model, before);
+    }
+
+    if (strcmp(arg,"pred") == 0) 
+       predict_phases_state(pexp, model, 1, model->L);
+
+    if (strcmp(arg,"pred1k") == 0) 
+       predict_phases(pexp, model, 1, model->L/4);
+
+    if (strcmp(arg,"smooth") == 0) {
+       smooth_phase(pexp, model,0);
+       update_snr_calc(pexp, model, before);
+    }
+    if (strcmp(arg,"smoothtrain") == 0) 
+       smooth_phase(pexp, model,1);
+    if (strcmp(arg,"smoothvq") == 0) {
+       smooth_phase(pexp, model,2);
+       update_snr_calc(pexp, model, before);
+    }
+
+    if (strcmp(arg,"smooth2") == 0) 
+       smooth_phase2(pexp, model);
+    if (strcmp(arg,"smooth3") == 0) 
+       smooth_phase3(pexp, model);
+    if (strcmp(arg,"smooth4") == 0) 
+       smooth_phase4(model);
+    if (strcmp(arg,"vqsmooth3") == 0)  {
+       sparse_vq_pred_error(pexp, model);
+       smooth_phase3(pexp, model);
+    }
+
+    if (strcmp(arg,"cb1") == 0) {
+       cb_phase1(pexp, model);
+       update_snr_calc(pexp, model, before);
+    }
+
+    if (strcmp(arg,"top") == 0) {
+       //top_amp(pexp, model, 1, model->L/4, 4, 1);
+       //top_amp(pexp, model, model->L/4, model->L/3, 4, 1);
+       //top_amp(pexp, model, model->L/3+1, model->L/2, 4, 1);
+       //top_amp(pexp, model, model->L/2, model->L, 6, 1);
+        //rand_phases(model, model->L/2, 3*model->L/4);
+       //struct_phases(pexp, model, model->L/2, 3*model->L/4);
+       //update_snr_calc(pexp, model, before);
+    }
+
+    if (strcmp(arg,"pred23") == 0) {
+       predict_phases2(pexp, model, model->L/2, model->L);
+       update_snr_calc(pexp, model, before);
+    }
+
+    if (strcmp(arg,"struct23") == 0) {
+       struct_phases(pexp, model, model->L/2, 3*model->L/4 );
+       update_snr_calc(pexp, model, before);
+    }
 
     /* normalise phases */
 
@@ -1149,13 +1681,28 @@ void phase_experiment(struct PEXP *pexp, MODEL *model) {
 
     /* update states */
 
+    //best_Wo = refine_Wo(pexp, model,  model->L/2, model->L);
+    pexp->phi1 += N*model->Wo;
+    
     for(m=1; m<=model->L; m++)
        pexp->phi_prev[m] = model->phi[m];          
     pexp->Wo_prev = model->Wo;
     pexp->frames++;
+    pexp->prev_model = *model;
 }
 
 #ifdef OLD_STUFF
+    //quant_phases(model, 1, model->L, 3);
+    //update_variance_calc(pexp, model, before);
+    //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);
+
+    //predict_phases1(pexp, model, 1, model->L/4);
+    //quant_phases(model, model->L/4+1, model->L, 3);
+
     //quant_phases(model, 1, model->L/8, 3);
 
     //update_snr_calc(pexp, model, before);
index d849f6b372438b5721f70b2739b51d8569fac7e6..8619c3fc9abfd3ca04f7b0b87e6dcf324a743224 100644 (file)
@@ -37,6 +37,6 @@ struct PEXP;
 
 struct PEXP * phase_experiment_create();
 void phase_experiment_destroy(struct PEXP *pexp);
-void phase_experiment(struct PEXP *pexp, MODEL *model);
+void phase_experiment(struct PEXP *pexp, MODEL *model, char *arg);
 
 #endif
index d8232925a14333a5057af95cfeb34e25b519d24b..5f24dd7337825857ddc43c8b8316eaf0c0d14ba2 100644 (file)
@@ -720,6 +720,14 @@ void aks_to_M2(
 
     signal += pow(model->A[m],2.0);
     noise  += pow(model->A[m] - Am,2.0);
+
+    /* this code improves perf of LPC model, in particular with phase0 */
+
+    if (Am > model->A[m])
+       Am *= 0.7;
+    if (Am < model->A[m])
+       Am *= 1.4;
+    
     model->A[m] = Am;
   }
   *snr = 10.0*log10(signal/noise);