From 8be3e55422ea970e4c917f434932efa4d7280fd5 Mon Sep 17 00:00:00 2001
From: drowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Date: Tue, 18 Jul 2017 03:33:17 +0000
Subject: [PATCH] batch and fbf working with slope based vq search

git-svn-id: https://svn.code.sf.net/p/freetel/code@3306 01035d8c-6547-0410-b346-abe4f91aad63
---
 codec2-dev/octave/newamp.m        |  20 ++++--
 codec2-dev/octave/newamp1_batch.m | 100 +++++++++++++++--------------
 codec2-dev/octave/newamp1_fbf.m   | 101 +++++++++++++-----------------
 3 files changed, 110 insertions(+), 111 deletions(-)

diff --git a/codec2-dev/octave/newamp.m b/codec2-dev/octave/newamp.m
index f31e81fd..805282b8 100644
--- a/codec2-dev/octave/newamp.m
+++ b/codec2-dev/octave/newamp.m
@@ -587,7 +587,7 @@ endfunction
 
 % Take a rate K surface and convert back to time varying rate L
 
-function [model_ AmdB_] = resample_rate_L(model, rate_K_surface, rate_K_sample_freqs_kHz, Fs=8000)
+function [model_ AmdB_] = resample_rate_L(model, rate_K_surface, rate_K_sample_freqs_kHz, Fs=8000, pad_end=0)
   max_amp = 160; K = columns(rate_K_surface);
   [frames col] = size(model);
 
@@ -600,10 +600,20 @@ function [model_ AmdB_] = resample_rate_L(model, rate_K_surface, rate_K_sample_f
     
     % back down to rate L
 
-    AmdB_(f,1:L) = interp1([0 rate_K_sample_freqs_kHz Fs/2000], 
-                            [rate_K_surface(f,1) rate_K_surface(f,:) rate_K_surface(f,K)], 
-                            rate_L_sample_freqs_kHz, 
-                            "spline");
+    % dealing with end effects is an ongoing issue.....need a better solution
+
+    if pad_end
+      AmdB_(f,1:L) = interp1([0 rate_K_sample_freqs_kHz Fs/2000], 
+                             [rate_K_surface(f,1) rate_K_surface(f,:) rate_K_surface(f,K)], 
+                             rate_L_sample_freqs_kHz, 
+                             "spline");
+    else
+      AmdB_(f,1:L) = interp1([0 rate_K_sample_freqs_kHz], 
+                             [rate_K_surface(f,1) rate_K_surface(f,:)], 
+                             rate_L_sample_freqs_kHz, 
+                             "spline");
+    end
+
     %AmdB_(f,1:L) = interp_para(rate_K_sample_freqs_kHz, rate_K_surface(f,:), Fs/(2*1000), rate_L_sample_freqs_kHz);
     %printf("f: %d %f %f %f\n", f, rate_K_sample_freqs_kHz(1), rate_L_sample_freqs_kHz(1), AmdB_(1));
     model_(f,1) = Wo; model_(f,2) = L; model_(f,3:(L+2)) = 10 .^ (AmdB_(f, 1:L)/20);
diff --git a/codec2-dev/octave/newamp1_batch.m b/codec2-dev/octave/newamp1_batch.m
index 5eaddf19..cc568943 100644
--- a/codec2-dev/octave/newamp1_batch.m
+++ b/codec2-dev/octave/newamp1_batch.m
@@ -52,9 +52,10 @@ function [surface mean_f] = newamp1_batch(input_prefix, varargin)
   synth_phase = output = 1;
   output_prefix = input_prefix;
   vq_type = "";
-  vq_filename = "";
-  fit_order = 0;
+  vq_filename = ""; 
+  vq_search = "mse";
   mode = "const";
+  fit_order = 0;
 
   % parse variable argument list
 
@@ -70,19 +71,27 @@ function [surface mean_f] = newamp1_batch(input_prefix, varargin)
     if ind
       mode =  varargin{ind+1};
     end
-    ind = arg_exists(varargin, "nomean");
+    ind = arg_exists(varargin, "vq_search");
     if ind
-      fit_order = 0;
+      vq_search = varargin{ind+1};
     end
-    ind = arg_exists(varargin, "noslope");
+    ind = arg_exists(varargin, "vq");
     if ind
-      fit_order = 1;
+      vq_type = varargin{ind};
+      vq_filename = varargin{ind+1};
+      printf("vq_type: %s vq_filename: %s vq_search: %s\n", vq_type, vq_filename, vq_search);
     end
     ind = arg_exists(varargin, "vqh");
     if ind
       vq_type = varargin{ind};
       vq_filename = varargin{ind+1};
-      printf("vq_type: %s vq_filename: %s\n", vq_type, vq_filename);
+      printf("vq_type: %s vq_filename: %s vq_search: %s\n", vq_type, vq_filename, vq_search);
+    end
+    ind = arg_exists(varargin, "vql");
+    if ind
+      vq_type = varargin{ind};
+      vq_filename = varargin{ind+1};
+      printf("vq_type: %s vq_filename: %s vq_search: %s\n", vq_type, vq_filename, vq_search);
     end
     ind = arg_exists(varargin, "no_output");
     if ind
@@ -91,7 +100,10 @@ function [surface mean_f] = newamp1_batch(input_prefix, varargin)
     end
   end
 
-  printf("output_prefix: %s\nfit_order %d output: %d\n", output_prefix, fit_order, output);
+  printf("output: %d\n", output);
+  if (output)
+    printf("output_prefix: %s\n",  output_prefix);
+  end
   
   model_name = strcat(input_prefix,"_model.txt");
   model = load(model_name);
@@ -114,7 +126,7 @@ function [surface mean_f] = newamp1_batch(input_prefix, varargin)
     [model_ surface] = experiment_mel_freq(model, 0, 1, voicing);
   end
   if strcmp(mode, 'const')
-    [model_ surface] = experiment_const_freq(model, fit_order, vq_type, vq_filename);
+    [model_ surface] = experiment_const_freq(model, vq_type, vq_filename, vq_search);
   end
   if strcmp(mode, 'piecewise')
     model_ = experiment_piecewise(model);
@@ -133,8 +145,8 @@ function [surface mean_f] = newamp1_batch(input_prefix, varargin)
     fWo  = fopen(Wo_out_name,"wb"); 
   
     if synth_phase
-     Hm_out_name = sprintf("%s_hm.out", output_prefix);
-     fhm = fopen(Hm_out_name,"wb"); 
+      Hm_out_name = sprintf("%s_hm.out", output_prefix);
+      fhm = fopen(Hm_out_name,"wb"); 
     end
 
     for f=1:frames
@@ -186,14 +198,8 @@ function [surface mean_f] = newamp1_batch(input_prefix, varargin)
     end
   end
 
-  if fit_order == 0;
-    for f=1:frames
-      surface(f,:) -= mean(surface(f,:));
-    end
-  end
-
-  if fit_order == 1
-    surface = slope_and_mean_removal(surface);
+  for f=1:frames
+    surface(f,:) -= mean(surface(f,:));
   end
 
   printf("\n")
@@ -225,7 +231,7 @@ endfunction
 % Basic unquantised rate K linear sampling then back to rate L.  Used for generating 
 % training vectors and testing vector quntisers.
 
-function [model_ rate_K_surface] = experiment_const_freq(model, fit_order, vq_type, vq_filename)
+function [model_ rate_K_surface] = experiment_const_freq(model, vq_type, vq_filename, vq_search="mse")
   melvq;
   [frames nc] = size(model);
   Fs = 8000;
@@ -235,6 +241,14 @@ function [model_ rate_K_surface] = experiment_const_freq(model, fit_order, vq_ty
   rate_K_sample_freqs_kHz = [0.1:0.1:4];
   K = length(rate_K_sample_freqs_kHz);
 
+  % optional full band VQ
+
+  if strcmp(vq_type, "vq")
+    quant_en = 1;
+    x = load(vq_filename); vq = x.vq;
+    [vq_rows vq_cols] = size(vq); vq_st = 1; vq_en = vq_cols;
+  end
+ 
   if strcmp(vq_type, "vql")
     quant_en = 1;
     x = load(vq_filename); vq = x.vq;
@@ -256,42 +270,35 @@ function [model_ rate_K_surface] = experiment_const_freq(model, fit_order, vq_ty
   rate_K_surface = resample_const_rate_f(model, rate_K_sample_freqs_kHz, Fs);
 
   rate_K_surface_fit = zeros(frames,K);
-  b = slope = zeros(1,frames);
-
-  if fit_order == 0
-    for f=1:frames
-      b(f) = mean(rate_K_surface(f,:));
-      rate_K_surface_fit(f,:) = rate_K_surface(f,:) - b(f);
-    end
-  end
+  b = zeros(1,frames);
 
-  if fit_order == 1
-    for f=1:frames
-      [aslope ab] = linreg(1:K, rate_K_surface(f,:), K);
-      rate_K_surface_fit(f,:) = rate_K_surface(f,:) - aslope*(1:K) - ab;
-      slope(f) = aslope; b(f) = ab;
-    end
+  for f=1:frames
+    b(f) = mean(rate_K_surface(f,:));
+    rate_K_surface_fit(f,:) = rate_K_surface(f,:) - b(f);
   end
-
+ 
   % optional vector quantise
   
   if quant_en
  
     rate_K_surface_fit_ = rate_K_surface_fit;
-    res = zeros(frames,vq_cols); ind = [];
+    res = zeros(frames, vq_cols); ind = [];
 
     for f=1:frames
       target = rate_K_surface_fit(f, vq_st:vq_en);
 
-      [diff_weighted weights error g mn_ind] = search_vq_weighted(target, vq);
-      if (f>=73) && (f<=75)
-        printf("f: %d mn_ind: %d g: %3.2f sdK: %3.2f\n", f, mn_ind, g(mn_ind), error(mn_ind));
+      if strcmp(vq_search, "mse")
+        [idx contrib errors test_ g mg sl] = vq_search_mse(vq, target);
+        rate_K_surface_fit_(f, vq_st:vq_en) = contrib;
       end
-      rate_K_surface_fit_(f, vq_st:vq_en) = vq(mn_ind,:) + g(mn_ind);
 
-      res(f,vq_st:vq_en) = rate_K_surface_fit(f,vq_st:vq_en) - rate_K_surface_fit_(f,vq_st:vq_en);
-      %res(f,vq_st:vq_en) = diff_weighted(mn_ind,:);
-      ind = [ind mn_ind];
+      if strcmp(vq_search, "slope")
+        [idx contrib errors test_ g mg sl] = vq_search_slope(vq, target);
+        rate_K_surface_fit_(f, vq_st:vq_en) = contrib;
+      end
+      
+      res(f, vq_st:vq_en) = target - contrib;
+      ind = [ind idx];
     end
     
     figure(fg++); clf; mesh(res);
@@ -304,7 +311,7 @@ function [model_ rate_K_surface] = experiment_const_freq(model, fit_order, vq_ty
   end
 
   for f=1:frames
-    rate_K_surface_(f,:) = rate_K_surface_fit_(f,:) + slope(f)*(1:K) + b(f);
+    rate_K_surface_(f,:) = rate_K_surface_fit_(f,:) + b(f);
   end
 
   [model_ AmdB_] = resample_rate_L(model, rate_K_surface_, rate_K_sample_freqs_kHz, Fs);
@@ -313,16 +320,13 @@ function [model_ rate_K_surface] = experiment_const_freq(model, fit_order, vq_ty
   % in the rate K <-> L transition.  Can optionally plot distorted
   % frames
 
-  plot_sd_thresh = 10;
+  plot_sd_thresh = 5;
   sd = zeros(1,frames);
   for f=1:frames
     Wo = model(f,1);
     L = model(f,2);
     AmdB = 20*log10(model(f,3:(L+2)));
     sd(f) = std(AmdB(1:L) - AmdB_(f,1:L));
-    if (f>=73) && (f<=75)
-      printf("f: %d sdL: %3.2f\n", f, sd(f));
-    end
     if (sd(f) > plot_sd_thresh) && (fg < 10)
       printf("fg: %d f: %d\n", fg, f);
       figure(fg++); clf; plot((1:L)*Wo*4/pi, AmdB(1:L),'b+-'); hold on; plot((1:L)*Wo*4/pi, AmdB_(f,1:L),'r+-');
diff --git a/codec2-dev/octave/newamp1_fbf.m b/codec2-dev/octave/newamp1_fbf.m
index 3b9d29d8..82c98188 100644
--- a/codec2-dev/octave/newamp1_fbf.m
+++ b/codec2-dev/octave/newamp1_fbf.m
@@ -15,15 +15,25 @@
 %   octave:14> newamp1_fbf("../build_linux/src/hts1a",50)
 
 
-function newamp1_fbf(samname, f=10, varargin)
+function newamp1_fbf(samname, f=73, varargin)
   more off;
 
   newamp;
   melvq;
 
   Fs = 8000; rate_K_sample_freqs_kHz = [0.1:0.1:4]; K = length(rate_K_sample_freqs_kHz);
-  quant_en = 0;
+  quant_en = 0; vq_search = "mse";
 
+  % optional full band VQ
+
+  ind = arg_exists(varargin, "vq");
+  if ind
+    quant_en = 1;
+    vq_filename = varargin{ind+1};
+    x = load(vq_filename); vq = x.vq;
+    [vq_rows vq_cols] = size(vq); vq_st = 1; vq_en = vq_cols;
+  end
+  
   % optional split VQ low freq quantiser
 
   ind = arg_exists(varargin, "vql");
@@ -44,12 +54,24 @@ function newamp1_fbf(samname, f=10, varargin)
     [vq_rows vq_cols] = size(vq); vq_st = 11; vq_en = K;
   end
 
+  % different vq search algorithms
+
+  ind = arg_exists(varargin, "vq_search");
+  if ind
+    vq_search = varargin{ind+1};
+  end
+
   fit_order = 0;
   ind = arg_exists(varargin, "noslope");
   if ind
     fit_order = 1;
   end
 
+  if quant_en
+    printf("quant_en: %d vq_filename: %s vq_st: %d vq_en: %d vq_search: %s\n", 
+           quant_en, vq_filename, vq_st, vq_en, vq_search);
+  end
+
   % load up text files dumped from c2sim ---------------------------------------
 
   sn_name = strcat(samname,"_sn.txt");
@@ -60,14 +82,6 @@ function newamp1_fbf(samname, f=10, varargin)
   model = load(model_name);
   [frames tmp] = size(model);
 
-  voicing_name = strcat(samname,"_pitche.txt");
-  voicing = zeros(1,frames);
-  
-  if exist(voicing_name, "file") == 2
-    pitche = load(voicing_name);
-    voicing = pitche(:, 3);
-  end
-
   % Keyboard loop --------------------------------------------------------------
 
   k = ' ';
@@ -78,15 +92,13 @@ function newamp1_fbf(samname, f=10, varargin)
     Wo = model(f,1); L = model(f,2); Am = model(f,3:(L+2)); AmdB = 20*log10(Am);
     Am_freqs_kHz = (1:L)*Wo*4/pi;
 
+    % remove constant gain term
+
     rate_K_vec = resample_const_rate_f(model(f,:), rate_K_sample_freqs_kHz, Fs);
     if fit_order == 0
       slope = 0; b = mean(rate_K_vec); 
       rate_K_vec_fit = rate_K_vec - b;
     end
-    if fit_order == 1
-      [slope b] = linreg(1:K, rate_K_vec, K);
-      rate_K_vec_fit = rate_K_vec - slope*(1:K) - b;
-    end
 
     % plots ----------------------------------
   
@@ -99,67 +111,40 @@ function newamp1_fbf(samname, f=10, varargin)
     if quant_en
       target = rate_K_vec_fit(vq_st:vq_en);
       
-      [diff_weighted weights error g mn_ind] = search_vq_weighted(target, vq);
- 
-      rate_K_vec_fit_ = rate_K_vec_fit;
-      rate_K_vec_fit_(vq_st:vq_en) = vq(mn_ind,:) + g(mn_ind);
-      rate_K_vec_ = slope*(1:K) + rate_K_vec_fit_ + b;
+      if strcmp(vq_search, "mse")
+        [idx contrib errors test_ g mg sl] = vq_search_mse(vq, target);
+        rate_K_surface_fit_(f, vq_st:vq_en) = contrib;
+      end
+
+      if strcmp(vq_search, "slope")
+        [idx contrib errors test_ g mg sl] = vq_search_slope(vq, target);
+        rate_K_surface_fit_(f, vq_st:vq_en) = contrib;
+      end
+
+      rate_K_vec_ = rate_K_vec_fit; rate_K_vec_(vq_st:vq_en) = contrib;
+      rate_K_vec_ += b;
       [model_ AmdB_] = resample_rate_L(model(f,:), rate_K_vec_, rate_K_sample_freqs_kHz, Fs);
       AmdB_ = AmdB_(1:L);
 
       sdL = std(AmdB - AmdB_);
-      printf("f: %d mn_ind: %d g: %3.2f sdK: %3.2f sdL: %3.2f\n", 
-             f, mn_ind, g(mn_ind), error(mn_ind), sdL);
+      %printf("f: %d mn_ind: %d g: %3.2f sdK: %3.2f sdL: %3.2f\n", 
+      %       f, mn_ind, g(mn_ind), error(mn_ind), sdL);
 
-      plot(rate_K_sample_freqs_kHz(vq_st:vq_en)*1000, diff_weighted(mn_ind,:), ";diff;k+-");
+      plot(rate_K_sample_freqs_kHz(vq_st:vq_en)*1000, target - contrib, ";diff;k+-");
       plot((1:L)*Wo*4000/pi, AmdB_,";AmdB bar;r+-");
       hold off;
 
-      figure (4); clf; 
-      plot(rate_K_sample_freqs_kHz(vq_st:vq_en)*1000, weights, ";weights;k+-");  
-      axis([0 4000 0 max(weights)]);
-
-      figure (5); clf; plot(error);
-
-      % sort and plot top m matches
-
-      figure(3); clf;
-      m = 4;
-      [error_list index_list] = sort(error);
-
-      mse_list = error_list(1:m);
-      index_list = index_list(1:m);
-      for i=1:m
-        subplot(sqrt(m),sqrt(m),i);
-        indx = index_list(i);
-        y_offset = g(indx) + b;
-        if quant_en == 2
-           y_offset = 0;
-        end
-        l = sprintf("b+-;ind %d g %3.2f sd %3.2f;",indx, g(indx), error(indx));
-        plot(target + y_offset,l);
-        hold on;
-        if index_list(i) == mn_ind
-          plot(vq(indx,:) + y_offset,'r-+');
-        else
-          plot(vq(indx,:) + y_offset,'g+-');
-        end
-        if quant_en != 2
-          plot(diff_weighted(indx,:), "ko-");
-        end
-        hold off;
-      end
     end
 
     % interactive menu ------------------------------------------
 
-    printf("\rframe: %d  menu: n-next  b-back  q-quit  m-show_quant[%d] o-fit order[%d]", f, quant_en, fit_order);
+    printf("\rframe: %d  menu: n-next  b-back  q-quit  m-show_quant[%d]", f, quant_en);
     fflush(stdout);
     k = kbhit();
 
     if k == 'm'
       quant_en++;
-      if quant_en > 2
+      if quant_en == 2
         quant_en = 0;
       end
     endif
-- 
2.25.1