From: drowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Date: Fri, 28 Jul 2017 00:03:35 +0000 (+0000)
Subject: extended order to fit a parabola, sounds pretty good between 1000 and 4000 Hz
X-Git-Url: http://git.whiteaudio.com/gitweb/?a=commitdiff_plain;h=1f03e2bbe21f01458c38347eaa533c89196159df;p=freetel-svn-tracking.git

extended order to fit a parabola, sounds pretty good between 1000 and 4000 Hz

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

diff --git a/codec2-dev/octave/newamp1_batch.m b/codec2-dev/octave/newamp1_batch.m
index aa633f73..ef41d1e6 100644
--- a/codec2-dev/octave/newamp1_batch.m
+++ b/codec2-dev/octave/newamp1_batch.m
@@ -241,16 +241,18 @@ function [model_ rate_K_surface] = experiment_const_freq(model, varargin)
   mask_en = 0;
   vq_search = "gain";   % defaul to gain search method as it's our favourite atm
   nvq = 0;              % number of vector quantisers
-  vq_start = vq_filename = [];
+  vq_start = [];
   
   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 settle don "gain" for now
+  % 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
 
+  weight_en = arg_exists(varargin, "weight");
+
   ind = arg_exists(varargin, "vq_search");
   if ind
     vq_search = varargin{ind+1};
@@ -263,8 +265,12 @@ function [model_ rate_K_surface] = experiment_const_freq(model, varargin)
     nvq++;
     vq_start = [vq_start varargin{ind+1}];
     avq_filename =  varargin{ind+2};
-    vq_filename = [vq_filename; avq_filename];
-    printf("nvq %d vq_start: %d vq_filename: %s\n", nvq, vq_start(nvq), avq_filename);
+    if nvq < 2
+      vq_filename = avq_filename;
+    else
+      vq_filename = [vqfilname; avq_filename];
+    end
+    printf("nvq %d vq_start: %d vq_filename: %s weight_en: %d\n", nvq, vq_start(nvq), avq_filename, weight_en > 0);
     anind = arg_exists(varargin(ind+1:length(varargin)), "vq");
     if anind
       ind += anind;
@@ -323,20 +329,38 @@ function [model_ rate_K_surface] = experiment_const_freq(model, varargin)
       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);
 
+      weights = ones(frames, 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
+        % have to calculate y intercept
+
+        for f=1:frames
+          target = rate_K_surface_no_mean(f,vq_st:vq_en);
+          gradient = 0.05; yint = 1 - gradient*min(target);
+          weights(f,:) = gradient*target + yint;
+        end
+      end
+
       if strcmp(vq_search, "mse")
         [idx contrib errors test_ g mg sl] = vq_search_mse(vq, rate_K_surface_no_mean(:,vq_st:vq_en));
       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));
+        [idx contrib errors test_ g mg sl] = vq_search_gain(vq, rate_K_surface_no_mean(:,vq_st:vq_en), weights);
       end
 
       if strcmp(vq_search, "sg")
-        [idx contrib errors test_ g mg sl] = vq_search_gain(vq, rate_K_surface_no_mean(:,vq_st:vq_en));
+        [idx contrib errors test_ g mg sl] = vq_search_sg(vq, rate_K_surface_no_mean(:,vq_st:vq_en));
       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));
+        [idx contrib errors test_ g mg sl] = vq_search_slope(vq, rate_K_surface_no_mean(:,vq_st:vq_en), weights);
+      end
+
+      if strcmp(vq_search, "para")
+        [idx contrib errors test_ g mg sl] = vq_search_para(vq, rate_K_surface_no_mean(:,vq_st:vq_en), weights);
       end
 
       rate_K_surface_no_mean_(:, vq_st:vq_en) = contrib;
@@ -352,8 +376,8 @@ function [model_ rate_K_surface] = experiment_const_freq(model, varargin)
           hmg(f) = mg(f, idx(f));
           hsl(f) = sl(f, idx(f));
         end
-        figure(fg++); clf; hist(hmg, 30);
-        figure(fg++); clf; hist(hsl, 30);
+        figure(fg++); clf; hist(hmg, 30); title('mag')
+        figure(fg++); clf; hist(hsl, 30); title('slope')
       end
 
       sd_per_frame = std(res(:,vq_st:vq_en)');
diff --git a/codec2-dev/octave/newamp1_fbf.m b/codec2-dev/octave/newamp1_fbf.m
index 82c98188..0800a5aa 100644
--- a/codec2-dev/octave/newamp1_fbf.m
+++ b/codec2-dev/octave/newamp1_fbf.m
@@ -22,38 +22,21 @@ 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 = "mse";
+  quant_en = 0; vq_search = "gain";
+  mask_en = 0;
+  nvec = 0;
+  weight_en = 0;
 
-  % optional full band VQ
+  % optional VQ
 
   ind = arg_exists(varargin, "vq");
   if ind
-    quant_en = 1;
-    vq_filename = varargin{ind+1};
+    nvec++;
+    vq_filename = varargin{ind+2};
     x = load(vq_filename); vq = x.vq;
-    [vq_rows vq_cols] = size(vq); vq_st = 1; vq_en = vq_cols;
+    [vq_rows vq_cols] = size(vq); vq_st = varargin{ind+1}; vq_en = vq_st + vq_cols - 1;
   end
   
-  % optional split VQ low freq quantiser
-
-  ind = arg_exists(varargin, "vql");
-  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_st + vq_cols - 1;
-  end
-  
-  % optional split VQ high freq quantiser
-
-  ind = arg_exists(varargin, "vqh");
-  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 = 11; vq_en = K;
-  end
-
   % different vq search algorithms
 
   ind = arg_exists(varargin, "vq_search");
@@ -108,17 +91,56 @@ function newamp1_fbf(samname, f=73, varargin)
     hold on;
     plot(rate_K_sample_freqs_kHz*1000, rate_K_vec, ";rate K;b+-");
 
-    if quant_en
+    if mask_en
+      % 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 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 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 test_ g mg sl] = vq_search_gain(vq, target, weights);
+      end
+
       if strcmp(vq_search, "slope")
-        [idx contrib errors test_ g mg sl] = vq_search_slope(vq, target);
+        [idx contrib errors test_ g mg sl] = vq_search_slope1(vq, target, weights);
         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
+
+      if strcmp(vq_search, "para")
+        printf("\n");
+        [idx contrib errors test_ g mg sl] = vq_search_para(vq, target, weights);
+        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_ = rate_K_vec_fit; rate_K_vec_(vq_st:vq_en) = contrib;
@@ -127,25 +149,36 @@ function newamp1_fbf(samname, f=73, varargin)
       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);
-
-      plot(rate_K_sample_freqs_kHz(vq_st:vq_en)*1000, target - contrib, ";diff;k+-");
+      plot(rate_K_sample_freqs_kHz(vq_st:vq_en)*1000, contrib, 'm+-');
+      l = sprintf(";diff sd = %3.2f;k+-", sdL);
+      plot(rate_K_sample_freqs_kHz(vq_st:vq_en)*1000, target - contrib, l);
       plot((1:L)*Wo*4000/pi, AmdB_,";AmdB bar;r+-");
       hold off;
+    end
 
+    if weight_en
+      figure(3); clf;
+      subplot(211);
+      plot((1:L)*Wo*4000/pi, AmdB,";AmdB;g+-");
+      axis([1 4000 -20 80]);
+      hold on;
+      plot(rate_K_sample_freqs_kHz*1000, rate_K_vec, ";rate K;b+-");
+      hold off;
+      subplot(212);
+      plot(rate_K_sample_freqs_kHz(vq_st:vq_en)*1000, weights);
+      axis([1 4000 0 8]);
     end
 
     % interactive menu ------------------------------------------
 
-    printf("\rframe: %d  menu: n-next  b-back  q-quit  m-show_quant[%d]", f, quant_en);
+    printf("\rframe: %d  menu: n-next  b-back  q-quit  w-weight[%d]", f, weight_en);
     fflush(stdout);
     k = kbhit();
 
-    if k == 'm'
-      quant_en++;
-      if quant_en == 2
-        quant_en = 0;
+    if k == 'w'
+      weight_en++;
+      if weight_en == 2
+        weight_en = 0;
       end
     endif
     if k == 'n'
diff --git a/codec2-dev/octave/vq_search_para.m b/codec2-dev/octave/vq_search_para.m
new file mode 100644
index 00000000..6455f519
--- /dev/null
+++ b/codec2-dev/octave/vq_search_para.m
@@ -0,0 +1,56 @@
+%----------------------------------------------------------------------
+% abs() search with a linear, ampl scaling, and slope term
+
+function [idx contrib errors test_ g mg sl] = vq_search_para(vq, data)
+  [nVec nCols] = size(vq);
+  nRows = rows(data);
+  
+  g = mg = sl = zeros(nRows, nVec);
+  diff  = zeros(nVec, nCols);
+  idx = errors = zeros(1, nRows);
+  error = zeros(1, nVec);
+  contrib = zeros(nRows, nCols);
+  test_ = zeros(nVec, nCols);
+
+  weights = ones(1,nCols);
+
+  A = zeros(4,4,nVec);
+  K = nCols; k=(1:K); k2 = k.^2; 
+  
+  
+  b_log = zeros(nVec, 4);
+  
+  for i=1:nVec
+    v = vq(i,:);
+    
+    A(:,:,i) = [v*v'    k2*v'    k*v'    sum(v); ...
+                k2*v'   k2*k2'   k*k2'   k*k';   ...                
+                k*v'    k*k2'    k*k'    sum(k); ...
+                sum(v)  k*k'     sum(k)  K       ];
+  end
+  
+  for f=1:nRows
+    t = data(f,:);
+
+    for i=1:nVec
+      v = vq(i,:);      
+      c = [t*v' t*k2' t*k' sum(t)]';
+      b = inv(A(:,:,i))*c;    
+      diff(i,:) = t - (b(1)*v + b(2)*k2 + b(3)*k + b(4));
+      b_log(i,:) = b; 
+      error(i) = diff(i,:) * diff(i,:)';
+
+      %printf("f: %d i: %d e: %f b(1): %f b(2): %f b(3): %f b(4): %f\n", f, i, error(i), b(1), b(2), b(3), b(4));
+    end
+    
+    [mn min_ind] = min(error);
+    errors(f) = mn; 
+    idx(f) = min_ind(1);
+    b = b_log(min_ind,:);
+    v = vq(min_ind,:);
+    
+    printf("f: %d i: %d b(1): %f b(2): %f b(3): %f b(4): %f\n", f, idx(f), b(1), b(2), b(3), b(4));
+    contrib(f,:) = test_(f,:) = b(1)*v + b(2)*k2 + b(3)*k + b(4);
+  end
+
+endfunction