first pass at constructing VQs on teh fly with individual amps
authordrowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Mon, 28 Aug 2017 23:41:53 +0000 (23:41 +0000)
committerdrowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Mon, 28 Aug 2017 23:41:53 +0000 (23:41 +0000)
git-svn-id: https://svn.code.sf.net/p/freetel/code@3355 01035d8c-6547-0410-b346-abe4f91aad63

codec2-dev/octave/newamp1_fbf.m
codec2-dev/octave/vq_construct_mg.m [new file with mode: 0644]

index 7c15b8936bcc29c4880693dcac5b8cbd0d603404..b4ee073f493f929bdc416cb1b9340bd452a89456 100644 (file)
@@ -24,41 +24,31 @@ function newamp1_fbf(samname, f=73, varargin)
   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"; 
   mask_en = 0;
-  nvec = 0;
+  nvq = 0; vq_start = [];
   quant_en = weight_en = 0;
 
-  % optional (split) VQ of rate K samples
-
-  ind = arg_exists(varargin, "vq");
-  if ind
-    nvec++;
-    vq_filename = varargin{ind+2};
-    x = load(vq_filename); vq = x.vq;
-    [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;
+  % specify one of more vqs, start index, vq file name, and search method
+
+  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;
+    avq_search =  varargin{ind+3};
+    if nvq < 2; vq_search = avq_search; else; vq_search = [vq_search; avq_search]; end;
+    printf("nvq %d vq_start: %d vq_filename: %s vq_search: %s\n", nvq, vq_start(nvq), avq_filename, avq_search);
+    anind = arg_exists(varargin(ind+1:length(varargin)), "vq");
+    if anind
+      ind += anind;
+    end
   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
 
+  ind = arg_exists(varargin, "construct");
+  vq_search = varargin{ind};
+  
   % optional exploration of phase
 
   ind = arg_exists(varargin, "phase");
@@ -119,133 +109,64 @@ function newamp1_fbf(samname, f=73, varargin)
     
     rate_K_vec_ = rate_K_vec_fit; 
  
-    if mask_en && nvec
+    if mask_en && nvq
       % 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);
       plot(rate_K_sample_freqs_kHz*1000, maskdB, ";mask dB;c+-");
     end
 
-    if nvec
-      if mask_en
-        target = rate_K_vec;
-        mask_thresh = 3;
-        ind = find (maskdB - target > mask_thresh);
-        target(ind) = maskdB(ind) - mask_thresh;
-        plot(rate_K_sample_freqs_kHz*1000, target, ";target;m+-");
-        target = target(vq_st:vq_en) - b;
-      else
-        target = rate_K_vec_fit(vq_st:vq_en);
-      end
-
-      weights = ones(1, vq_en - vq_st + 1);
-      if weight_en
-
-        % generate weighting.  if min is 20dB and max 40dB, weight at min
-        % is 1 and max 2, same for -10 to 10.  So gradient = 0.05, we only
-        % haveto calculate y intercept
-
-        gradient = 0.05; yint = 1 - gradient*min(target);
-        weights = gradient*target + yint;
-      end
+    if strcmp(vq_search, "construct")
+      target = rate_K_vec_fit;
+      [idx contrib errors b] = vq_construct_mg(target);
+      rate_K_vec_ = contrib;
     end
-    
-    if nvec
-      
-      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, "gain")
-        [idx contrib errors b] = vq_search_gain(vq, target, weights);
-      end
-
-      if strcmp(vq_search, "max")
-        [idx contrib errors b] = vq_search_max(vq, target);
-      end
+    if nvq
 
-      if strcmp(vq_search, "sg")
-        [idx contrib errors b] = vq_search_sg(vq, target);
-      end
+      for i=1:nvq
+        avq_filename = char(cellstr(vq_filename)(i));
+        avq_search = char(cellstr(vq_search)(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("\nsplit VQ: %d vq_filename: %s vq_search: %s vq_st: %d vq_en: %d nVec: %d\n", i, avq_filename, avq_search, vq_st, vq_en, vq_rows);
+        target = rate_K_vec_fit(vq_st:vq_en);
 
-      if strcmp(vq_search, "mg")
-        [idx contrib errors b] = vq_search_mg(vq, target);
-      end
+        if strcmp(avq_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 b_log] = vq_search_slope(vq, target, "closed_quant_slope");
-        rate_K_surface_fit_(f, vq_st:vq_en) = contrib;
-        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;
+        if strcmp(avq_search, "gain")
+          [idx contrib errors b] = vq_search_gain(vq, target, weights);
         end
-      end
 
-      if strcmp(vq_search, "para")
-        printf("\n");
-        [idx contrib errors b] = vq_search_para(vq, target);
+        if strcmp(avq_search, "max")
+          [idx contrib errors b] = vq_search_max(vq, target);
+        end
 
-        k = 1:vq_cols; k2 = k.^2;
-        para_target = k2*b(2) + k*b(3) + b(4);
-        samples = [1 10 25];
-        if quant_en
+        if strcmp(avq_search, "sg")
+          [idx contrib errors b] = vq_search_sg(vq, target);
+        end
 
-#{
-          % search vq_gain for best match to gain coefficients
+        if strcmp(avq_search, "mg")
+          [idx contrib errors b] = vq_search_mg(vq, target);
+        end
 
-          [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';
+        if strcmp(avq_search, "slope")
+          [idx contrib errors b_log] = vq_search_slope(vq, target, "closed_quant_slope");
+          rate_K_surface_fit_(f, vq_st:vq_en) = contrib;
+          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
-          [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;
+        rate_K_vec_(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
-
-      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+-');
@@ -265,9 +186,9 @@ function newamp1_fbf(samname, f=73, varargin)
     sdL = std(abs(AmdB - AmdB_));
 
     plot((1:L)*Wo*4000/pi, AmdB_,";AmdB bar;r+-");
-    if nvec == 0
-      l = sprintf(";errorx10 sd %3.2f dB;bk+-", sdL);
-      plot((1:L)*Wo*4000/pi, 10*(AmdB - AmdB_), l);
+    if nvq == 0
+      l = sprintf(";error sd %3.2f dB;bk+-", sdL);
+      plot((1:L)*Wo*4000/pi, (AmdB - AmdB_), l);
     end
     hold off;
 
@@ -357,7 +278,7 @@ endfunction
 function ind = arg_exists(v, str) 
    ind = 0;
    for i=1:length(v)
-      if strcmp(v{i}, str)
+      if !ind && strcmp(v{i}, str)
         ind = i;
       end     
     end
diff --git a/codec2-dev/octave/vq_construct_mg.m b/codec2-dev/octave/vq_construct_mg.m
new file mode 100644 (file)
index 0000000..843f812
--- /dev/null
@@ -0,0 +1,117 @@
+%----------------------------------------------------------------------
+% Construct a VQ entry using formant prototypes and gain/mag fitting
+
+#{
+  ptable - prototype table
+         - each row defines a formant prototype with 3 numbers: shape start end
+         - shape 1 - gaussian
+                 2 - symmetrical triangle
+                 3 - triangle with deep fall off at HF side
+         - start 0 - one after last prototype center position
+                 n - start of search range in 100's of Hz, e.g. 20 -> 2000Hz
+         - end   m - end of search range in 100's of Hz
+#}
+
+function [idx contrib errors b_log2] = vq_construct_mg(data)
+
+  ptable = [ 1  4  8 1 10;
+             2  8 30 0 30;
+             2 12 30 0 30;
+             3 30 36 0 40];
+  
+  protos = [-4  12 23  24  18  12 8 5 2;
+             12 18 12   0   0   0 0 0 0;
+             12 18 12 -12 -24 -36 0 0 0];
+  protoc = [4 2 2];
+  protol = [9 3 6];                      
+
+  nformants = rows(ptable);
+  [nRows nCols] = size(data);  
+  K = nCols;  
+  
+  idx = zeros(nRows, nformants); b_log2 = [];
+
+  contrib = zeros(nRows, K);
+
+  #{
+    range of F1, specify range of each formant
+    specify shape
+
+    [X] table of shapes
+    [X] work out centre position
+    [X] another table to refer to each shape based on formant number
+    [X] const for number of formants
+    [ ] iterate over positions
+    [ ] gain/mag match at each position, MSE over entire vector
+    [ ] way to graphically peer inside what's going on
+    [ ] how to model flat in first 1kHz?
+        + cld perhaps exclude search entirely as an option
+        + or perhap sit will trun up low mags
+  #}
+  
+  for r=1:nRows
+    % iterate across all formants
+    
+    for f=1:nformants
+
+      % try formant centre at np points between st and en
+      
+      shape = ptable(f,1); c_st = ptable(f,2); c_en = ptable(f,3); 
+      if c_st == 0  % 0 code means start just after previous formant
+        c_st = idx(r,f-1) + 2;
+      end
+      np = c_en-c_st+1;
+      b_log = zeros(np, 2);
+      error = zeros(np, 1);
+
+      % we are matching a target segment between st and en
+       
+      st = ptable(f,4); en = ptable(f,5);
+      if st == 0
+        st = c_st;
+      end
+      lt = en-st+1;
+      t = data(r,st:en);
+
+      printf("r: %d f: %d c_st: %d c_en: %d np: %d st: %d en: %d ------------------------\n",
+              r, f, c_st, c_en, np, st, en);
+
+      % test range of centres for formant
+      
+      for c=c_st:c_en
+        v = zeros(1, K);
+        v_st = c - protoc(shape) + 1; v_en = c - protoc(shape) + protol(shape);
+        v(v_st:v_en) = protos(shape, 1:protol(shape));
+        v = v(st:en);
+        A = [v*v' sum(v); sum(v)  lt];
+        d = [t*v' sum(t)]';
+        b = inv(A)*d;
+        b(1) = max(0.5,b(1));
+        diff = t - (b(1)*v + b(2));
+        p = c - c_st + 1;
+        b_log(p,:) = b; 
+        error(p) = diff * diff';
+
+        %printf("r: %d f: %d c: %d e: %f b(1): %f b(2): %f \n", r, f, c, error(p), b(1), b(2));
+      end
+
+      % choose best for this formant
+
+      [mn p_min] = min(error);
+      c_min = p_min + c_st - 1
+      b = b_log(p_min,:);
+      v = zeros(1, K);
+      v_st = c_min - protoc(shape) + 1; v_en = c_min - protoc(shape) + protol(shape);
+      v(v_st:v_en) = protos(shape, 1:protol(shape));
+      contrib(r,st:en) = b(1)*v(st:en) + b(2);
+      
+      % log index, mag, gain, might need a matric with one row per formant
+
+      idx(r,f) = c_min;
+      b_log2 = [b_log2; b];
+    end
+    
+    errors(r) = sum((data(r,:) - contrib(r,:)) .^ 2); 
+  end
+
+endfunction