From 2d41672da414f9892107241e07cf7ac40dde10e4 Mon Sep 17 00:00:00 2001
From: drowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Date: Fri, 11 Aug 2017 02:00:41 +0000
Subject: [PATCH] reas results with 2 bit slope quantiser 1-3k

git-svn-id: https://svn.code.sf.net/p/freetel/code@3351 01035d8c-6547-0410-b346-abe4f91aad63
---
 codec2-dev/octave/newamp1_batch.m   | 255 ++++++++++++++++++++++++++--
 codec2-dev/octave/newamp1_fbf.m     | 104 ++++++++++--
 codec2-dev/octave/vq_search_slope.m |  14 +-
 3 files changed, 340 insertions(+), 33 deletions(-)

diff --git a/codec2-dev/octave/newamp1_batch.m b/codec2-dev/octave/newamp1_batch.m
index 33db7099..35a14ea3 100644
--- a/codec2-dev/octave/newamp1_batch.m
+++ b/codec2-dev/octave/newamp1_batch.m
@@ -56,7 +56,8 @@ function [surface_no_mean surface] = newamp1_batch(input_prefix, varargin)
   vq_search = "mse";
   mode = "const";
   fit_order = 0;
-
+  mean_remove = 1;
+  
   % parse variable argument list
 
   if (length (varargin) > 0)
@@ -114,6 +115,10 @@ function [surface_no_mean surface] = newamp1_batch(input_prefix, varargin)
       end
     end
   end
+  if strcmp(mode, 'pred')
+    [model_ surface b_log] = experiment_const_freq_pred(model, varargin{:});
+    mean_remove = 0;
+  end
   if strcmp(mode, 'piecewise')
     model_ = experiment_piecewise(model);
   end
@@ -184,10 +189,14 @@ function [surface_no_mean surface] = newamp1_batch(input_prefix, varargin)
     end
   end
 
-  for f=1:frames
-    surface_no_mean(f,:) = surface(f,:) - mean(surface(f,:));
+  if mean_remove
+    for f=1:frames
+      surface_no_mean(f,:) = surface(f,:) - mean(surface(f,:));
+    end
+  else
+    surface_no_mean = surface;
   end
-
+  
   printf("\n")
 
 endfunction
@@ -223,7 +232,7 @@ function [model_ rate_K_surface b_log] = experiment_const_freq(model, varargin)
   Fs = 8000;
   fg = 1;
   mask_en = 0;
-  vq_search = "gain";   % defaul to gain search method as it's our favourite atm
+  vq_search = "gain";   % default to gain search method as it's our favourite atm
   nvq = 0;              % number of vector quantisers
   vq_start = [];
   quant_en = 0;
@@ -360,7 +369,7 @@ function [model_ rate_K_surface b_log] = experiment_const_freq(model, varargin)
       end
 
       if strcmp(vq_search, "gain")
-        [idx contrib errors test_ g mg sl] = vq_search_gain(vq, rate_K_surface_no_mean(:,vq_st:vq_en), weights);
+        [idx contrib errors b_log] = vq_search_gain(vq, rate_K_surface_no_mean(:,vq_st:vq_en), weights);
       end
 
       if strcmp(vq_search, "mg")
@@ -373,7 +382,16 @@ function [model_ rate_K_surface b_log] = experiment_const_freq(model, varargin)
       end
 
       if strcmp(vq_search, "slope")
-        [idx contrib errors test_ g mg sl] = vq_search_slope(vq, rate_K_surface_no_mean(:,vq_st:vq_en), weights);
+        [idx contrib errors b_log] = vq_search_slope(vq, rate_K_surface_no_mean(:,vq_st:vq_en), "closed_quant_slope");
+        b_log = [b_log energy' idx'];
+
+        for f=1:frames
+          target = rate_K_surface_no_mean(f,vq_st:vq_en);
+          b_log(f,2) = quantise([-1 -0.5 0.5 1], b_log(f,2));
+          b_log(f,3) = (sum(target) - sum(b_log(f,1)*vq(idx(f),:)+b_log(f,2)*(1:vq_cols)))/vq_cols;
+          contrib(f,:) = b_log(f,1)*vq(idx(f),:) + b_log(f,2)*(1:vq_cols) + b_log(f,3);
+          %rate_K_surface_no_mean_(f, vq_en+1:K) -= b_log(f,2)*vq_cols + b_log(f,3);
+        end
       end
 
       if strcmp(vq_search, "para")
@@ -430,14 +448,18 @@ function [model_ rate_K_surface b_log] = experiment_const_freq(model, varargin)
       res(:, vq_st:vq_en) = rate_K_surface_no_mean(:, vq_st:vq_en) - contrib;
       ind(:,i) = idx;
     
-      % histograms of higher order gain/shape params if we are in slope mode
+      % histograms of higher order gain/shape params
 
+      if strcmp(vq_search, "gain")
+        figure(fg++); clf; hist(b_log, 30); title('gain')
+      end
+      
       if strcmp(vq_search, "slope")
 
         hmg = hsl = zeros(1,frames);
         for f=1:frames
-          hmg(f) = mg(f, idx(f));
-          hsl(f) = sl(f, idx(f));
+          hmg(f) = b_log(f, 1);
+          hsl(f) = b_log(f, 3);
         end
         figure(fg++); clf; hist(hmg, 30); title('mag')
         figure(fg++); clf; hist(hsl, 30); title('slope')
@@ -463,6 +485,11 @@ function [model_ rate_K_surface b_log] = experiment_const_freq(model, varargin)
   % optional decimation
   
   if decimate
+    %lf = [1.00 0.75 0.50 0.25];
+    %lf = [1.00 1.00 0.50 0.00];
+    %lf = [1.00 1.00 0.00 0.00];
+    lf = [1.00 0.80 0.20 0.00];
+    
     for f=1:frames
 
       % determine frames that bracket the one we need to interp
@@ -475,7 +502,7 @@ function [model_ rate_K_surface b_log] = experiment_const_freq(model, varargin)
 
       % determine fraction of each frame to use
 
-      left_fraction  = 1 - mod((f-1),decimate)/decimate;
+      left_fraction  = lf(mod(f-1,decimate)+1);
       right_fraction = 1 - left_fraction;
 
       rate_K_surface_(f,:) = left_fraction*rate_K_surface_(left_f,:) + right_fraction*rate_K_surface_(right_f,:);
@@ -534,6 +561,212 @@ function [model_ rate_K_surface b_log] = experiment_const_freq(model, varargin)
    
 endfunction
 
+function model_ = experiment_dct(model)
+
+  [frames tmp] = size(model); max_amp = 160;
+
+  for f=1:frames
+    printf("%d ", f);   
+    Wo = model(f,1);
+    L = min([model(f,2) max_amp-1]);
+    Am = model(f,3:(L+2));
+    AmdB = 20*log10(Am);
+
+    % fit model
+
+    D = dct(AmdB);
+    E = zeros(1,L);
+    E(1:min(20,L)) = D(1:min(20,L));
+    AmdB_ = idct(D);
+
+    model_(f,1) = Wo; model_(f,2) = L; model_(f,3:(L+2)) = 10 .^ (AmdB_(1:L)/20);
+  end
+
+endfunction
+
+
+% Predictive VQ using rate K resampled vectors, and rate 2 decimation
+
+function [model_ rate_K_surface_pred_ b_log] = experiment_const_freq_pred(model, varargin)
+  melvq;
+  [frames nc] = size(model);
+  Fs = 8000;
+  fg = 1;
+  nvq = 0;              % number of vector quantisers
+  vq_start = [];
+  quant_en = 0;
+  b_log = [];
+  decimate = 2;
+  frames2 = floor(frames/decimate);
+  
+  rate_K_sample_freqs_kHz = [0.1:0.1:4];
+  K = length(rate_K_sample_freqs_kHz);
+  
+  % parse command line options
+
+  % set vq search algorithm, e.g. mse, gain, mag, slope.  We've settled on "gain" for now
+  % as slope required extra bits to quantise higher order parameters that offset advantages
+
+  ind = arg_exists(varargin, "vq_search");
+  if ind
+    vq_search = varargin{ind+1};
+  end
+
+  % specify one of more VQs
+
+  ind = anind = arg_exists(varargin, "vq");
+  while anind
+    nvq++;
+    vq_start = [vq_start varargin{ind+1}];
+    avq_filename =  varargin{ind+2};
+    if nvq < 2
+      vq_filename = avq_filename;
+    else
+      vq_filename = [vq_filename; avq_filename];
+    end
+    printf("nvq %d vq_start: %d vq_filename: %s\n", nvq, vq_start(nvq), avq_filename);
+    anind = arg_exists(varargin(ind+1:length(varargin)), "vq");
+    if anind
+      ind += anind;
+    end
+  end
+  
+  ind = arg_exists(varargin, "vq_gain");
+  if ind
+    avq_filename = varargin{ind+1};
+    x = load(avq_filename); vq_gain = x.vq;
+    quant_en = 1
+  end
+
+  ind = arg_exists(varargin, "decimate");
+  if ind
+    decimate = varargin{ind+1};
+  end
+   
+  % OK start processing ............................................
+
+  energy = zeros(1,frames/decimate);
+  for f=1:decimate:frames
+    L = model(f,2);
+    energy(f) = 10*log10(sum( model(f,3:(L+2)) .^ 2 ));
+  end
+
+  rate_K_surface = resample_const_rate_f(model, rate_K_sample_freqs_kHz, Fs);
+
+  % decimate and predict
+  
+  rate_K_surface_pred = zeros(frames2, K);
+  beta = 0.9;
+  g = 1;
+  for f=3:decimate:frames
+    rate_K_surface_pred(g,:) = rate_K_surface(f,:) - beta*rate_K_surface(f-2,:);
+    g++;
+  end
+  
+  rate_K_surface_ = zeros(frames,K);
+  rate_K_surface_(1,:) = rate_K_surface(1,:);
+  g = 1;
+
+  % optional vector quantise
+  
+  if nvq
+ 
+    % note we init with target (ideal) to fill in values not covered by this VQ
+
+    rate_K_surface_pred_ = rate_K_surface_pred;
+    res = zeros(frames2, K); ind = zeros(frames2, nvq);
+
+    % quantise using split VQs
+
+    for i=1:nvq
+      avq_filename = char(cellstr(vq_filename)(i));
+      x = load(avq_filename); vq = x.vq; [vq_rows vq_cols] = size(vq);
+      vq_st = vq_start(i); vq_en = vq_st + vq_cols - 1;
+      printf("split VQ: %d vq_filename: %s vq_st: %d vq_en: %d nVec: %d\n", i, avq_filename, vq_st, vq_en, vq_rows);
+
+      for f=3:decimate:frames
+        rate_K_vec_pred = rate_K_surface(f,:) - beta*rate_K_surface_(f-2,:);
+        
+        if strcmp(vq_search, "mg")
+          [idx contrib errors b_log] = vq_search_mg(vq, rate_K_vec_pred(vq_st:vq_en));
+        end
+        
+        rate_K_vec_pred_ = rate_K_vec_pred;
+        rate_K_vec_pred_(vq_st:vq_en) = contrib; 
+        rate_K_surface_(f,:) = beta*rate_K_surface_(f-2,:) + rate_K_vec_pred_;
+
+        res(:, vq_st:vq_en) = rate_K_surface_pred(:, vq_st:vq_en) - contrib;
+        ind(:,i) = idx;
+      end
+      
+      % histograms of higher order gain/shape params if we are in slope mode
+
+      sd_per_frame = std(res(:,vq_st:vq_en)');
+      t=sprintf("VQ %d", i); 
+      figure(fg++); subplot(211); plot(energy); title(t); subplot(212); plot(sd_per_frame); 
+      figure(fg++); subplot(211); hist(sd_per_frame); title(t); subplot(212); hist(ind(:,i),100);
+      printf("VQ rms SD: %3.2f\n", mean(sd_per_frame));
+    end
+ 
+    figure(fg++); clf; mesh(res);
+
+  else
+    rate_K_surface_pred_ = rate_K_surface_pred;
+    for f=3:decimate:frames
+      rate_K_surface_(f,:) = beta*rate_K_surface_(f-2,:) + rate_K_surface_pred_(g,:);
+      g++;
+    end
+  end
+
+  % re-interpolate back to 10ms rate
+  
+  if decimate
+    for f=1:frames
+
+      % determine frames that bracket the one we need to interp
+
+      left_f = decimate*floor((f-1)/decimate)+1; 
+      right_f = left_f + decimate;
+      if right_f > frames
+        right_f = left_f;
+      end
+
+      % determine fraction of each frame to use
+
+      left_fraction  = 1 - mod((f-1),decimate)/decimate;
+      right_fraction = 1 - left_fraction;
+
+      rate_K_surface_(f,:) = left_fraction*rate_K_surface_(left_f,:) + right_fraction*rate_K_surface_(right_f,:);
+      %printf("f: %d left_f: %d right_f: %d left_fraction: %3.2f right_fraction: %3.2f \n", f, left_f, right_f, left_fraction, right_fraction)
+    end
+  end
+  
+  [model_ AmdB_] = resample_rate_L(model, rate_K_surface_, rate_K_sample_freqs_kHz, Fs);
+
+  % Measure distortion between AmdB and AmdB_, ths includes distortion
+  % in the rate K <-> L transition.  Can optionally plot distorted
+  % frames
+
+  plot_sd_thresh = 5;
+  sd = zeros(1,frames2); g = 1;
+  for f=1:decimate:frames
+    Wo = model(f,1);
+    L = model(f,2);
+    AmdB = 20*log10(model(f,3:(L+2)));
+    sd(g) = std(AmdB(1:L) - AmdB_(f,1:L));
+    if (sd(g) > 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+-');
+      plot(rate_K_sample_freqs_kHz, rate_K_surface_(f,:), 'c+-'); hold off;
+    end
+    g++;
+  end
+  printf("rate K resampling SD: %3.2f\n", mean(sd));
+  figure(fg++); clf; subplot(211); plot(energy); subplot(212); plot(sd); title('sdL');
+  figure(fg++); clf; hist(sd);
+   
+endfunction
+
 
 function model_ = experiment_dct(model)
 
diff --git a/codec2-dev/octave/newamp1_fbf.m b/codec2-dev/octave/newamp1_fbf.m
index 224d3519..86473394 100644
--- a/codec2-dev/octave/newamp1_fbf.m
+++ b/codec2-dev/octave/newamp1_fbf.m
@@ -22,12 +22,12 @@ function newamp1_fbf(samname, f=73, varargin)
   melvq;
 
   Fs = 8000; rate_K_sample_freqs_kHz = [0.1:0.1:4]; K = length(rate_K_sample_freqs_kHz);
-  quant_en = 0; vq_search = "gain";
+  quant_en = 0; vq_search = "gain"; 
   mask_en = 0;
   nvec = 0;
-  weight_en = quant_en = 0;
+  quant_en = weight_en = 0;
 
-  % optional VQ
+  % optional (split) VQ of rate K samples
 
   ind = arg_exists(varargin, "vq");
   if ind
@@ -37,6 +37,15 @@ function newamp1_fbf(samname, f=73, varargin)
     [vq_rows vq_cols] = size(vq); vq_st = varargin{ind+1}; vq_en = vq_st + vq_cols - 1;
   end
   
+  % optional VQ of VQ gain coefficients
+
+  ind = arg_exists(varargin, "vq_gain");
+  if ind
+    nvec++;
+    vq_filename = varargin{ind+1};
+    x = load(vq_filename); vq_gain = x.vq;
+  end
+  
   % different vq search algorithms
 
   ind = arg_exists(varargin, "vq_search");
@@ -95,7 +104,7 @@ function newamp1_fbf(samname, f=73, varargin)
     % default to the ideal
     
     rate_K_vec_ = rate_K_vec_fit; 
-
+ 
     if mask_en && nvec
       % experimental masking stuff that I can't seem to get to work
       maskdB = determine_mask(rate_K_vec, rate_K_sample_freqs_kHz, rate_K_sample_freqs_kHz, bark_model=1);
@@ -134,24 +143,81 @@ function newamp1_fbf(samname, f=73, varargin)
       end
 
       if strcmp(vq_search, "gain")
-        [idx contrib errors test_ g mg sl] = vq_search_gain(vq, target, weights);
+        [idx contrib errors b] = vq_search_gain(vq, target, weights);
+      end
+
+      if strcmp(vq_search, "sg")
+        [idx contrib errors b] = vq_search_sg(vq, target);
       end
 
       if strcmp(vq_search, "slope")
-        [idx contrib errors test_ g mg sl] = vq_search_slope1(vq, target, weights);
+        [idx contrib errors b_log] = vq_search_slope(vq, target, "closed_quant_slope");
         rate_K_surface_fit_(f, vq_st:vq_en) = contrib;
-        % printf("g: %3.2f mg: %3.2f sl: %3.2f\n", g(idx), mg(idx), sl(idx));
+        printf(" mg: %3.2f sl: %3.2f g: %3.2f \n", b_log(1), b_log(2), b_log(3));
+        if quant_en
+          % set slope to 0
+          contrib1 = contrib;
+          contrib = b_log(1)*vq(idx,:) + b_log(3);
+          rate_K_vec_(vq_en+1:K) -= b_log(2)*vq_cols;
+        end
       end
 
       if strcmp(vq_search, "para")
         printf("\n");
         [idx contrib errors b] = vq_search_para(vq, target);
+
+        k = 1:vq_cols; k2 = k.^2;
+        para_target = k2*b(2) + k*b(3) + b(4);
+        samples = [1 10 25];
         if quant_en
-          % recalc gain using some index
-          g = (sum(target(1:vq_cols-10)) - sum(mg*vq(idx,1:vq_cols-10)))/length(target);
-          contrib = mg*vq(idx,:) + g;
-          contrib(vq_cols-4:vq_cols) += -3:-6:-30;
+
+#{
+          % search vq_gain for best match to gain coefficients
+
+          [nr nc] = size(vq_gain);
+          d = g = zeros(nr,1);
+          for r=1:nr
+            g(r) = (sum(para_target) - sum(vq_gain(r,:)))/vq_cols;
+            diff = para_target - (vq_gain(r,:) + g(r));
+            d(r) = diff*diff';
+          end
+          [dmin imin] = min(d);
+          
+#}
+          v = vq(idx,:);
+#{
+          rng = 1:vq_cols-5;
+          g = (sum(target(rng)) - sum(b(1)*v(rng)))/(vq_cols-5);
+          v(vq_cols-5:vq_cols) += -10*(1:6);
+          printf("g: %f\n", g);
+          % recalc contrib
+#}
+          para_target(1) = quantise([-20 -10 0 10], para_target(1));
+          %para_target(10) = quantise([-6 +6], para_target(10));
+          %para_target(10)
+          b_ = polyfit([k(1) k(10) k(25)],
+                       [para_target(1) para_target(10) -10],
+                       2);
+          
+          contrib1 = contrib;
+          contrib = b(1)*v + b_(1)*k2 + b_(2)*k + b_(3);
+          para = b_(1)*k2 + b_(2)*k + b_(3);
+          %printf("imin: %d\n", imin);
         end
+
+        rate_K_surface_fit_(f, vq_st:vq_en) = contrib;
+
+      end
+
+      if strcmp(vq_search, "cubic")
+        printf("\n");
+        [idx contrib errors b] = vq_search_cubic(target);
+        rate_K_surface_fit_(f, vq_st:vq_en) = contrib;
+      end
+
+      if strcmp(vq_search, "fourth")
+        printf("\n");
+        [idx contrib errors b] = vq_search_fourth(target);
         rate_K_surface_fit_(f, vq_st:vq_en) = contrib;
         % printf("g: %3.2f mg: %3.2f sl: %3.2f\n", g(idx), mg(idx), sl(idx));
       end
@@ -159,6 +225,12 @@ function newamp1_fbf(samname, f=73, varargin)
       rate_K_vec_(vq_st:vq_en) = contrib;
 
       plot(rate_K_sample_freqs_kHz(vq_st:vq_en)*1000, contrib, 'm+-');
+      if strcmp(vq_search, "para")
+        plot(rate_K_sample_freqs_kHz(vq_st:vq_en)*1000, para_target, 'c+-');
+        if quant_en
+          plot(rate_K_sample_freqs_kHz(vq_st:vq_en)*1000, para, 'r+-');
+        end
+      end
       l = sprintf(";diff vq sd = %3.2f;k+-", std(target - contrib));
       plot(rate_K_sample_freqs_kHz(vq_st:vq_en)*1000, target - contrib, l);
     end
@@ -177,6 +249,12 @@ function newamp1_fbf(samname, f=73, varargin)
     end
     hold off;
 
+    if quant_en
+      figure(4); clf;
+      plot(contrib1, 'b+-');
+      hold on; plot(contrib,'r+'); hold off;
+    end
+    
     if weight_en
       figure(3); clf;
       subplot(211);
@@ -198,9 +276,7 @@ function newamp1_fbf(samname, f=73, varargin)
 
     if k == 'w'
       quant_en++;
-      if quant_en == 2
-        quant_en = 0;
-      end
+      if quant_en == 2; quant_en = 0; end
     endif
     if k == 'n'
       f = f + 1;
diff --git a/codec2-dev/octave/vq_search_slope.m b/codec2-dev/octave/vq_search_slope.m
index 18aa99d0..f0c6cb82 100644
--- a/codec2-dev/octave/vq_search_slope.m
+++ b/codec2-dev/octave/vq_search_slope.m
@@ -1,7 +1,7 @@
 %----------------------------------------------------------------------
 % abs() search with a linear, ampl scaling, and slope term
 
-function [idx contrib errors b_log2] = vq_search_slope(vq, data)
+function [idx contrib errors b_log2] = vq_search_slope(vq, data, closed_quant_fn, open_quant_fn)
   [nVec nCols] = size(vq);
   nRows = rows(data);
   
@@ -23,20 +23,18 @@ function [idx contrib errors b_log2] = vq_search_slope(vq, data)
 
   for f=1:nRows
     target = data(f,:);
-    %target = 2*vq(1,:)+1;
 
     for i=1:nVec
       c = [sum(target) target*(1:nCols)' target*vq(i,:)' ]';
       b = inv(A(:,:,i))*c;
-
-      b(1) = max(0.5,b(1)); b(1) = min(1,b(1));
+      if nargin == 3;
+        b = feval(closed_quant_fn,b);
+      end;  
       b_log(i,:) = b; 
       
       diff(i,:) = target - (b(1)*vq(i,:) + b(2)*(1:nCols) + b(3));
-      %diff(i,nCols-5:nCols) *= 0.25;
 
       error(i) = diff(i,:) * diff(i,:)';
-      b_log(i,:) = b; 
 
       %printf("f: %d i: %d mg: %f g: %f sl: %f error: %f\n", f, i, b(1), b(2), b(3), error(i));
     end
@@ -47,8 +45,8 @@ function [idx contrib errors b_log2] = vq_search_slope(vq, data)
     b = b_log(min_ind,:);
     b_log2(f,:) = b;
     
-    printf("f: %d mg: %f sl: %f g: %f\n", f, b(1), b(2), b(3));
-    contrib(f,:) = test_(f,:) = 0.8*vq(min_ind,:) + b(2)*(1:nCols) + b(3);
+    %printf("f: %d i: %d mg: %f sl: %f g: %f\n", f, idx(f), b(1), b(2), b(3));
+    contrib(f,:) = b(1)*vq(min_ind,:) + b(2)*(1:nCols) + b(3);
   end
 
 endfunction
-- 
2.25.1