From: drowe67 Date: Sat, 25 Aug 2012 04:35:58 +0000 (+0000) Subject: about to clean up after a few weeks of amplitude and phase experiments X-Git-Url: http://git.whiteaudio.com/gitweb/?a=commitdiff_plain;h=fe898b0466968f3bec36c39d1f7178f9a4c07f16;p=freetel-svn-tracking.git about to clean up after a few weeks of amplitude and phase experiments git-svn-id: https://svn.code.sf.net/p/freetel/code@640 01035d8c-6547-0410-b346-abe4f91aad63 --- diff --git a/codec2-dev/src/ampexp.c b/codec2-dev/src/ampexp.c index 8920cb36..4d3b70a8 100644 --- a/codec2-dev/src/ampexp.c +++ b/codec2-dev/src/ampexp.c @@ -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; mcb, &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; ik; 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; mL; 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; mvq_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; ik; 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; mL; 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; mindexes[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; mvq->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; mL; 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; mL; 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-1)) + en = MAX_AMP-1; + for(j=i; jvq1, weights, smoothed); + for(i=0; i (MAX_AMP-1)) + en = MAX_AMP-1; + for(j=i; jL; 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; iL; m++) { + f = m*model->Wo*FS/TWO_PI; + if (f > bins[0]) { + + /* find bin */ + + for(i=0; i 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 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 */ diff --git a/codec2-dev/src/c2sim.c b/codec2-dev/src/c2sim.c index 26531b11..434f4f39 100644 --- a/codec2-dev/src/c2sim.c +++ b/codec2-dev/src/c2sim.c @@ -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, &exp, 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; mWo)*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; iphi_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; iA[i]; sorted[j].index = i; + printf("%f ", model->A[i]); } - bubbleSort(&sorted[1], end-start); + bubbleSort(sorted, end-start); + + printf("\n"); + for(j=0; jA[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; ik; 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; mL; 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-1)) + en = MAX_AMP-1; + for(j=i; jvq1, smoothed_weights, smoothed_in); + for(i=0; i (MAX_AMP-1)) + en = MAX_AMP-1; + for(j=i; jL; 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; iphi_prev[i] + N*i*best_Wo; + err = model->phi[i] - pred; + c += cos(err); s += sin(err); + } + phi1_ = atan2(s,c); + for(i=a; iphi_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; iL; m++) { + f = m*model->Wo*FS/TWO_PI; + if (f > bins[0]) { + + /* find bin */ + + for(i=0; i 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 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; iL; m++) { + f = m*model->Wo*FS/TWO_PI; + if (f > bins[0]) { + + /* find bin */ + + for(i=0; i 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 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 model->L) + b = model->L; + + c = s = 0; + for(i=a; iphi[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; iphi[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); diff --git a/codec2-dev/src/phase.h b/codec2-dev/src/phase.h index d849f6b3..8619c3fc 100644 --- a/codec2-dev/src/phase.h +++ b/codec2-dev/src/phase.h @@ -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 diff --git a/codec2-dev/src/quantise.c b/codec2-dev/src/quantise.c index d8232925..5f24dd73 100644 --- a/codec2-dev/src/quantise.c +++ b/codec2-dev/src/quantise.c @@ -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);