first pass at indep amp and gain estimation for each formant

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

diff --git a/codec2-dev/octave/newamp1_fbf.m b/codec2-dev/octave/newamp1_fbf.m
index b4ee073f493f929bdc416cb1b9340bd452a89456..0858a5fa91d24ec6e10a2d826b990746d2f19c82 100644 (file)
--- a/codec2-dev/octave/newamp1_fbf.m
+++ b/codec2-dev/octave/newamp1_fbf.m
@@ -47,7 +47,13 @@ function newamp1_fbf(samname, f=73, varargin)
   fit_order = 0;
 
   ind = arg_exists(varargin, "construct");
-  vq_search = varargin{ind};
+  if ind
+    vq_search = varargin{ind};
+  end
+  ind = arg_exists(varargin, "construct_indep");
+  if ind
+    vq_search = varargin{ind};
+  end
   
   % optional exploration of phase
 
@@ -116,8 +122,12 @@ function newamp1_fbf(samname, f=73, varargin)
     end
 
     if strcmp(vq_search, "construct")
-      target = rate_K_vec_fit;
-      [idx contrib errors b] = vq_construct_mg(target);
+      [idx contrib errors b] = vq_construct_mg(rate_K_vec);
+      rate_K_vec_ = contrib;
+    end
+
+    if strcmp(vq_search, "construct_indep")
+      [idx contrib errors b] = vq_construct_indep_mg(rate_K_vec);
       rate_K_vec_ = contrib;
     end
 
@@ -180,7 +190,9 @@ function newamp1_fbf(samname, f=73, varargin)
 
     % And .... back to rate L
     
-    rate_K_vec_ += meanf;
+    if (strcmp(vq_search, "construct") == 0) && (strcmp(vq_search, "construct_indep") == 0)
+      rate_K_vec_ += meanf;
+    end
     [model_ AmdB_] = resample_rate_L(model(f,:), rate_K_vec_, rate_K_sample_freqs_kHz, Fs);
     AmdB_ = AmdB_(1:L);
     sdL = std(abs(AmdB - AmdB_));

diff --git a/codec2-dev/octave/vq_construct_indep_mg.m b/codec2-dev/octave/vq_construct_indep_mg.m
new file mode 100644 (file)
index 0000000..037f9cc
--- /dev/null
+++ b/codec2-dev/octave/vq_construct_indep_mg.m
@@ -0,0 +1,118 @@
+% vq_construct_indep_mag.m
+% David Rowe
+% Sep 2017
+%
+%----------------------------------------------------------------------
+% Construct a VQ entry using formant prototypes and gain/mag fitting
+% indep amplitude and gain for each formant
+
+#{
+  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_indep_mg(data)
+
+#{
+  ptable = [ 1  4  8;
+             2  8 30;
+             2 12 30;
+             3 30 36];
+#}
+
+  ptable = [ 1  4  8 ];
+  
+  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);
+  
+  for r=1:nRows
+    % iterate across all formants
+    
+    % set up target
+      
+    t = data(r,:);
+
+    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);
+
+      printf("r: %d f: %d c_st: %d c_en: %d np: %d ------------------------\n", r, f, c_st, c_en, np);
+
+      % test range of centres for formant
+      
+      for c=c_st:c_en
+
+        % construct current prototype vector and gain window vector
+
+        p_st = c - protoc(shape) + 1; p_en = c - protoc(shape) + protol(shape);
+        p = zeros(1,K);
+        p(p_st:p_en) = protos(shape, 1:protol(shape));
+        g = zeros(1,K);
+        g(p_st:p_en) = ones(1, protol(shape));
+        
+        A = [p*p' sum(p); sum(p) protol(shape) ];
+        d = [t*p' t*g']';
+        b = inv(A)*d;
+
+        v = b(1)*p + b(2)*g;
+        diff = t - v;
+        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,:);
+      printf("f: %d cmin: %d b(1): %f b(2): %f -------------\n", f, c_min, b(1), b(2));
+      
+      % sum this formats contribution
+      
+      v = zeros(1, K);
+      v_st = c_min - protoc(shape) + 1; v_en = c_min - protoc(shape) + protol(shape);
+      v(v_st:v_en) = b(1)*protos(shape, 1:protol(shape)) + b(2);
+      contrib(r,:) += v;
+
+      % update residual target for next formant search
+      
+      t = data(r,:) - contrib(r,:);
+      
+      % log index, mag, gain
+
+      idx(r,f) = c_min;
+      b_log2 = [b_log2; b];
+    end
+    
+    errors(r) = sum(t .^ 2); 
+  end
+
+endfunction