From: drowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Date: Thu, 5 Jan 2017 22:27:16 +0000 (+0000)
Subject: C versions of energy quant and post filter tested
X-Git-Url: http://git.whiteaudio.com/gitweb/?a=commitdiff_plain;h=8fb3014968850603ba4921dbfc46ba39e63290ce;p=freetel-svn-tracking.git

C versions of energy quant and post filter tested

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

diff --git a/codec2-dev/octave/newamp.m b/codec2-dev/octave/newamp.m
index e4f7a76d..a55af918 100644
--- a/codec2-dev/octave/newamp.m
+++ b/codec2-dev/octave/newamp.m
@@ -25,17 +25,18 @@ function y = interp_para(xp, yp, x)
   assert( (length(xp) >=3) && (length(yp) >= 3) );
 
   y = zeros(1,length(x));
+  k = 1;
   for i=1:length(x)
     xi = x(i);
 
     % k is index into xp of where we start 3 points used to form parabola
 
-    k = 1;
     while ((xp(k+1) < xi) && (k < (length(xp)-2))) 
       k++;
     end
     
     x1 = xp(k); y1 = yp(k); x2 = xp(k+1); y2 = yp(k+1); x3 = xp(k+2); y3 = yp(k+2);
+    %printf("k: %d i: %d xi: %f x1: %f y1: %f\n", k, i, xi, x1, y1);
 
     a = ((y3-y2)/(x3-x2)-(y2-y1)/(x2-x1))/(x3-x1);
     b = ((y3-y2)/(x3-x2)*(x2-x1)+(y2-y1)/(x2-x1)*(x3-x2))/(x3-x1);
@@ -178,6 +179,7 @@ function rate_K_sample_freqs_kHz = mel_sample_freqs_kHz(K)
   rate_K_sample_freqs_kHz = rate_K_sample_freqs_Hz/1000;
 endfunction
 
+
 function [rate_K_surface rate_K_sample_freqs_kHz] = resample_const_rate_f_mel(model, K) 
   rate_K_sample_freqs_kHz = mel_sample_freqs_kHz(K);
   rate_K_surface = resample_const_rate_f(model, rate_K_sample_freqs_kHz);
@@ -208,7 +210,7 @@ function [rate_K_surface rate_K_sample_freqs_kHz] = resample_const_rate_f(model,
 
     AmdB_peak = max(AmdB);
     AmdB(find(AmdB < (AmdB_peak-50))) = AmdB_peak-50;
-
+    
     rate_L_sample_freqs_kHz = (1:L)*Wo*4/pi;
     
     %rate_K_surface(f,:) = interp1(rate_L_sample_freqs_kHz, AmdB, rate_K_sample_freqs_kHz, "spline", "extrap");
@@ -271,8 +273,44 @@ function vec = post_filter(vec, sample_freq_kHz, pf_gain = 1.5, voicing)
 endfunction
 
 
+% construct energy quantiser table, and save to text file to include in C
+
+function energy_q = create_energy_q
+    energy_q = 10 + 40/16*(0:15);
+endfunction
+
+function save_energy_q(fn)
+  energy_q = create_energy_q;
+  f = fopen(fn, "wt");
+  fprintf(f, "1 %d\n", length(energy_q));
+  for n=1:length(energy_q)
+    fprintf(f, "%f\n", energy_q(n));
+  end
+  fclose(f);
+endfunction
+
+
+% save's VQ in format that can be compiled by Codec 2 build system
+
+function save_vq(vqset, filenameprefix)
+  [Nvec order stages] = size(vqset);
+  for s=1:stages
+    fn = sprintf("%s_%d.txt", filenameprefix, s);
+    f = fopen(fn, "wt");
+    fprintf(f, "%d %d\n", order, Nvec);
+    for n=1:Nvec
+      for k=1:order
+        fprintf(f, "% 8.4f ", vqset(n,k,s));
+      end
+      fprintf(f, "\n");
+    end
+    fclose(f);
+  end
+endfunction
+
+
 % --------------------------------------------------------------------------------
-% Experimental functions used for masking, piecewise models
+% Experimental functions used for masking, piecewise models, not part of newamp1
 % --------------------------------------------------------------------------------
 
 
diff --git a/codec2-dev/octave/tnewamp1.m b/codec2-dev/octave/tnewamp1.m
index 59eb4304..af41622a 100644
--- a/codec2-dev/octave/tnewamp1.m
+++ b/codec2-dev/octave/tnewamp1.m
@@ -61,14 +61,36 @@ function tnewamp1(input_prefix)
   K = 20;
   [frames tmp] = size(rate_K_surface_c);
   [rate_K_surface sample_freqs_kHz] = resample_const_rate_f_mel(model(1:frames,:), K);
-rate_K_surface
+
+  melvq;
+  load train_120_vq; m=5;
+       
+  for f=1:frames
+    mean_f(f) = mean(rate_K_surface(f,:));
+    rate_K_surface_no_mean(f,:) = rate_K_surface(f,:) - mean_f(f);
+  end
+
+  [res rate_K_surface_no_mean_ ind] = mbest(train_120_vq, rate_K_surface_no_mean, m);
+
+  for f=1:frames
+    rate_K_surface_no_mean_(f,:) = post_filter(rate_K_surface_no_mean_(f,:), sample_freqs_kHz, 1.5);
+  end
+    
+  rate_K_surface_ = zeros(frames, K);
+  energy_q = create_energy_q;
+  for f=1:frames   
+    [mean_f_ indx] = quantise(energy_q, mean_f(f));
+    indexes(f,3) = indx - 1;
+    rate_K_surface_(f,:) = rate_K_surface_no_mean_(f,:) + mean_f_;
+  end
+
   figure(1);
-  mesh(rate_K_surface);
+  mesh(rate_K_surface_);
   figure(2);
-  mesh(rate_K_surface_c);
+  mesh(rate_K_surface__c);
   figure(3);
-  mesh(rate_K_surface - rate_K_surface_c);
-
+  mesh(rate_K_surface_ - rate_K_surface__c);
+  axis([1 K 1 frames -1 1])
   #{
 
   for f=1:frames
diff --git a/codec2-dev/src/CMakeLists.txt b/codec2-dev/src/CMakeLists.txt
index 5abc3e4f..fee648d9 100644
--- a/codec2-dev/src/CMakeLists.txt
+++ b/codec2-dev/src/CMakeLists.txt
@@ -55,12 +55,16 @@ set(CODEBOOKSLSPMELVQ
     ${D}/lspmelvq3.txt
 )
 
+set(CODEBOOKSGE ${D}/gecb.txt)
+
 set(CODEBOOKSNEWAMP1
     ${D}/train_120_1.txt
     ${D}/train_120_2.txt
 )
 
-set(CODEBOOKSGE ${D}/gecb.txt)
+set(CODEBOOKSNEWAMP1_ENERGY
+    ${D}/newamp1_energy_q.txt
+)
 
 # when crosscompiling we need a native executable
 if(CMAKE_CROSSCOMPILING)
@@ -138,6 +142,13 @@ add_custom_command(
     DEPENDS generate_codebook ${CODEBOOKSNEWAMP1}
 )
 
+# codebooknewamp1_energy.c
+add_custom_command(
+    OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/codebooknewamp1_energy.c
+    COMMAND generate_codebook newamp1_energy_cb ${CODEBOOKSNEWAMP1_ENERGY} > ${CMAKE_CURRENT_BINARY_DIR}/codebooknewamp1_energy.c
+    DEPENDS generate_codebook ${CODEBOOKSNEWAMP1_ENERGY}
+)
+
 #
 # codec2 library sources
 #
@@ -169,8 +180,9 @@ set(CODEC2_SRCS
     codebookjvm.c
     codebookmel.c
     codebooklspmelvq.c
-    codebooknewamp1.c
     codebookge.c
+    codebooknewamp1.c
+    codebooknewamp1_energy.c
     golay23.c
     freedv_api.c
     freedv_vhf_framing.c
diff --git a/codec2-dev/src/codebook/newamp1_energy_q.txt b/codec2-dev/src/codebook/newamp1_energy_q.txt
new file mode 100644
index 00000000..d2653f34
--- /dev/null
+++ b/codec2-dev/src/codebook/newamp1_energy_q.txt
@@ -0,0 +1,17 @@
+1 16
+10.000000
+12.500000
+15.000000
+17.500000
+20.000000
+22.500000
+25.000000
+27.500000
+30.000000
+32.500000
+35.000000
+37.500000
+40.000000
+42.500000
+45.000000
+47.500000
diff --git a/codec2-dev/src/defines.h b/codec2-dev/src/defines.h
index c1ebcdd9..ed9a0bbf 100644
--- a/codec2-dev/src/defines.h
+++ b/codec2-dev/src/defines.h
@@ -95,5 +95,6 @@ extern const struct lsp_codebook mel_cb[];
 extern const struct lsp_codebook ge_cb[];
 extern const struct lsp_codebook lspmelvq_cb[];
 extern const struct lsp_codebook newamp1vq_cb[];
+extern const struct lsp_codebook newamp1_energy_cb[];
 
 #endif
diff --git a/codec2-dev/src/quantise.c b/codec2-dev/src/quantise.c
index cc02bf00..05de19a5 100644
--- a/codec2-dev/src/quantise.c
+++ b/codec2-dev/src/quantise.c
@@ -2197,7 +2197,7 @@ void interp_para(float y[], float xp[], float yp[], int np, float x[], int n)
     
         x1 = xp[k]; y1 = yp[k]; x2 = xp[k+1]; y2 = yp[k+1]; x3 = xp[k+2]; y3 = yp[k+2];
 
-        printf("k: %d np: %d i: %d xi: %f x1: %f y1: %f\n", k, np, i, xi, x1, y1);
+        //printf("k: %d np: %d i: %d xi: %f x1: %f y1: %f\n", k, np, i, xi, x1, y1);
 
         a = ((y3-y2)/(x3-x2)-(y2-y1)/(x2-x1))/(x3-x1);
         b = ((y3-y2)/(x3-x2)*(x2-x1)+(y2-y1)/(x2-x1)*(x3-x2))/(x3-x1);
@@ -2262,7 +2262,7 @@ void resample_const_rate_f(MODEL *model, float rate_K_vec[], float rate_K_sample
             AmdB_peak = AmdB[m];
         }
         rate_L_sample_freqs_kHz[m] = m*model->Wo*4.0/M_PI;
-        printf("m: %d AmdB: %f AmdB_peak: %f  sf: %f\n", m, AmdB[m], AmdB_peak, rate_L_sample_freqs_kHz[m]);
+        //printf("m: %d AmdB: %f AmdB_peak: %f  sf: %f\n", m, AmdB[m], AmdB_peak, rate_L_sample_freqs_kHz[m]);
     }
     
     /* clip between peak and peak -50dB, to reduce dynamic range */
@@ -2277,4 +2277,110 @@ void resample_const_rate_f(MODEL *model, float rate_K_vec[], float rate_K_sample
 }
 
 
+/*---------------------------------------------------------------------------*\
+
+  FUNCTION....: rate_K_mbest_encode
+  AUTHOR......: David Rowe
+  DATE CREATED: Jan 2017
+
+  Two stage rate K newamp1 VQ quantiser using mbest search.
+
+\*---------------------------------------------------------------------------*/
+
+float rate_K_mbest_encode(int *indexes, float *x, float *xq, int ndim, int mbest_entries)
+{
+  int i, j, n1, n2;
+  const float *codebook1 = newamp1vq_cb[0].cb;
+  const float *codebook2 = newamp1vq_cb[1].cb;
+  struct MBEST *mbest_stage1, *mbest_stage2;
+  float target[ndim];
+  float w[ndim];
+  int   index[MBEST_STAGES];
+  float mse, tmp;
+
+  for(i=0; i<ndim; i++)
+      w[i] = 1.0;
+
+  mbest_stage1 = mbest_create(mbest_entries);
+  mbest_stage2 = mbest_create(mbest_entries);
+  for(i=0; i<MBEST_STAGES; i++)
+      index[i] = 0;
+
+  /* Stage 1 */
+
+  mbest_search(codebook1, x, w, ndim, newamp1vq_cb[0].m, mbest_stage1, index);
+  MBEST_PRINT("Stage 1:", mbest_stage1);
+
+  /* Stage 2 */
+
+  for (j=0; j<mbest_entries; j++) {
+      index[1] = n1 = mbest_stage1->list[j].index[0];
+      for(i=0; i<ndim; i++)
+	  target[i] = x[i] - codebook1[ndim*n1+i];
+      mbest_search(codebook2, target, w, ndim, newamp1vq_cb[1].m, mbest_stage2, index);
+  }
+  MBEST_PRINT("Stage 2:", mbest_stage2);
+
+  n1 = mbest_stage2->list[0].index[1];
+  n2 = mbest_stage2->list[0].index[0];
+  mse = 0.0;
+  for (i=0;i<ndim;i++) {
+      tmp = codebook1[ndim*n1+i] + codebook2[ndim*n2+i];
+      mse += (x[i]-tmp)*(x[i]-tmp);
+      xq[i] = tmp;
+  }
+
+  mbest_destroy(mbest_stage1);
+  mbest_destroy(mbest_stage2);
+
+  indexes[0] = n1; indexes[1] = n2;;
+
+  return mse;
+}
+
+
+/*---------------------------------------------------------------------------*\
 
+  FUNCTION....: post_filter
+  AUTHOR......: David Rowe
+  DATE CREATED: Jan 2017
+
+  Post Filter, has a big impact on speech quality after VQ.  When used
+  on a mean removed rate K vector, it raises formants, and supresses
+  anti-formants.  As it manipulates amplitudes, we normalise energy to
+  prevent clipping or large level variations.  pf_gain of 1.2 to 1.5
+  (dB) seems to work OK.  Good area for further investigations and
+  improvements in speech quality.
+
+\*---------------------------------------------------------------------------*/
+
+void post_filter_newamp1(float vec[], float sample_freq_kHz[], int K, float pf_gain)
+{
+    int k;
+
+    /*
+      vec is rate K vector describing spectrum of current frame lets
+      pre-emp before applying PF. 20dB/dec over 300Hz.  Postfilter
+      affects energy of frame so we measure energy before and after
+      and normalise.  Plenty of room for experiment here as well.
+    */
+    
+    float pre[K];
+    float e_before = 0.0;
+    float e_after = 0.0;
+    for(k=0; k<K; k++) {
+        pre[k] = 20.0*log10f(sample_freq_kHz[k]/0.3);
+        vec[k] += pre[k];
+        e_before += powf(10.0, 2.0*vec[k]/20.0);
+        vec[k] *= pf_gain;
+        e_after += powf(10.0, 2.0*vec[k]/20.0);        
+    }
+
+    float gain = e_after/e_before;
+    float gaindB = 10*log10f(gain);
+  
+    for(k=0; k<K; k++) {
+        vec[k] -= gaindB;
+        vec[k] -= pre[k];
+    }
+}
diff --git a/codec2-dev/src/quantise.h b/codec2-dev/src/quantise.h
index e6df842c..765b6efc 100644
--- a/codec2-dev/src/quantise.h
+++ b/codec2-dev/src/quantise.h
@@ -139,5 +139,7 @@ void interp_para(float y[], float xp[], float yp[], int np, float x[], int n);
 float ftomel(float fHz);
 void mel_sample_freqs_kHz(float rate_K_sample_freqs_kHz[], int K);
 void resample_const_rate_f(MODEL *model, float rate_K_vec[], float rate_K_sample_freqs_kHz[], int K);
+float rate_K_mbest_encode(int *indexes, float *x, float *xq, int ndim, int mbest_entries);
+void post_filter_newamp1(float vec[], float sample_freq_kHz[], int K, float pf_gain);
 
 #endif
diff --git a/codec2-dev/unittest/tnewamp1.c b/codec2-dev/unittest/tnewamp1.c
index 586ba9a0..eda9a18c 100644
--- a/codec2-dev/unittest/tnewamp1.c
+++ b/codec2-dev/unittest/tnewamp1.c
@@ -48,7 +48,7 @@ int main(int argc, char *argv[]) {
     MODEL model;
     void *nlp_states;
     float pitch, prev_uq_Wo;
-    int   i,m,f;
+    int   i,m,f,k;
     
     if (argc != 2) {
         printf("usage: ./tnewamp1 RawFile\n");
@@ -65,14 +65,20 @@ int main(int argc, char *argv[]) {
 
     int K = 20;
     float rate_K_sample_freqs_kHz[K];
-    float model_octave[FRAMES][MAX_AMP+2]; // model params in matrix format, useful for C <-> Octave  
-    float rate_K_surface[FRAMES][K];       // rate K vecs for each frame, form a surface that makes pretty graphs
+    float model_octave[FRAMES][MAX_AMP+2];    // model params in matrix format, useful for C <-> Octave  
+    float rate_K_surface[FRAMES][K];          // rate K vecs for each frame, form a surface that makes pretty graphs
+    float rate_K_surface_no_mean[FRAMES][K];  // mean removed surface  
+    float rate_K_surface_no_mean_[FRAMES][K]; // quantised mean removed surface  
+    float mean[FRAMES];
+    float mean_[FRAMES];
+    float rate_K_surface_[FRAMES][K];         // quantised rate K vecs for each frame
 
     for(f=0; f<FRAMES; f++)
         for(m=0; m<MAX_AMP+2; m++)
             model_octave[f][m] = 0.0;
 
     mel_sample_freqs_kHz(rate_K_sample_freqs_kHz, K);
+
     //for(int k=0; k<K; k++)
     //    printf("k: %d sf: %f\n", k, rate_K_sample_freqs_kHz[k]);
 
@@ -103,12 +109,33 @@ int main(int argc, char *argv[]) {
 	two_stage_pitch_refinement(&model, Sw);
 	estimate_amplitudes(&model, Sw, W, 1);
 
-        /* Resample at rate K ----------------------------------------*/
+        /* newamp1 processing ----------------------------------------*/
 
         resample_const_rate_f(&model, &rate_K_surface[f][0], rate_K_sample_freqs_kHz, K);
-        for(int k=0; k<K; k++)
-            printf("k: %d sf: %f sv: %f\n", k, rate_K_sample_freqs_kHz[k], rate_K_surface[f][k]);
-        printf("\n");
+        float sum = 0.0;
+        for(k=0; k<K; k++)
+            sum += rate_K_surface[f][k];
+        mean[f] = sum/K;
+        for(k=0; k<K; k++)
+            rate_K_surface_no_mean[f][k] = rate_K_surface[f][k] - mean[f];
+
+        int vq_indexes[2];
+        rate_K_mbest_encode(vq_indexes, &rate_K_surface_no_mean[f][0], &rate_K_surface_no_mean_[f][0], K, 5);
+
+        post_filter_newamp1(&rate_K_surface_no_mean_[f][0], rate_K_sample_freqs_kHz, K, 1.5);
+
+        int energy_index;
+        float w[1] = {1.0};
+        float se;
+        energy_index = quantise(newamp1_energy_cb[0].cb, &mean[f], w, newamp1_energy_cb[0].k, newamp1_energy_cb[0].m, &se);
+        mean_[f] = newamp1_energy_cb[0].cb[energy_index];
+
+        for(k=0; k<K; k++)
+            rate_K_surface_[f][k] = rate_K_surface_no_mean_[f][k] + mean_[f];
+
+        //for(int k=0; k<K; k++)
+        //    printf("k: %d sf: %f sv: %f\n", k, rate_K_sample_freqs_kHz[k], rate_K_surface[f][k]);
+        //printf("\n");
 
         /* log vectors */
  
@@ -127,6 +154,11 @@ int main(int argc, char *argv[]) {
     assert(fout != NULL);
     fprintf(fout, "# Created by tnewamp1.c\n");
     octave_save_float(fout, "rate_K_surface_c", (float*)rate_K_surface, FRAMES, K, K);
+    octave_save_float(fout, "mean_c", (float*)mean, FRAMES, 1, 1);
+    octave_save_float(fout, "rate_K_surface_no_mean_c", (float*)rate_K_surface_no_mean, FRAMES, K, K);
+    octave_save_float(fout, "rate_K_surface_no_mean__c", (float*)rate_K_surface_no_mean_, FRAMES, K, K);
+    octave_save_float(fout, "mean__c", (float*)mean_, FRAMES, 1, 1);
+    octave_save_float(fout, "rate_K_surface__c", (float*)rate_K_surface_, FRAMES, K, K);
     octave_save_float(fout, "model_c", (float*)model_octave, FRAMES, MAX_AMP+2, MAX_AMP+2);
     fclose(fout);