From: drowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Date: Sun, 30 Jul 2017 08:13:22 +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=12e2e495dcc2ef3bda7eaf3b91de5c3dc6628ca3;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@3339 01035d8c-6547-0410-b346-abe4f91aad63
---

diff --git a/codec2-dev/octave/newamp1_batch.m b/codec2-dev/octave/newamp1_batch.m
index ef41d1e6..90935d75 100644
--- a/codec2-dev/octave/newamp1_batch.m
+++ b/codec2-dev/octave/newamp1_batch.m
@@ -71,30 +71,7 @@ function [surface mean_f] = newamp1_batch(input_prefix, varargin)
     if ind
       mode =  varargin{ind+1};
     end
-#{
-    ind = arg_exists(varargin, "vq_search");
-    if ind
-      vq_search = varargin{ind+1};
-    end
-    ind = arg_exists(varargin, "vq");
-    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, "vqh");
-    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, "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
       output = 0;
@@ -128,7 +105,12 @@ 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, varargin{:});
+    [model_ surface b_log] = experiment_const_freq(model, varargin{:});
+    ind = arg_exists(varargin, "vq_search");
+    if ind && strcmp(varargin{ind+1},"para")
+      fn = sprintf("%s_b_log.txt", output_prefix);
+      save(fn,"b_log");
+    end
   end
   if strcmp(mode, 'piecewise')
     model_ = experiment_piecewise(model);
@@ -233,7 +215,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, varargin)
+function [model_ rate_K_surface b_log] = experiment_const_freq(model, varargin)
   melvq;
   [frames nc] = size(model);
   Fs = 8000;
@@ -242,6 +224,8 @@ function [model_ rate_K_surface] = experiment_const_freq(model, varargin)
   vq_search = "gain";   % defaul to gain search method as it's our favourite atm
   nvq = 0;              % number of vector quantisers
   vq_start = [];
+  quant_en = 0;
+  b_log = [];
   
   rate_K_sample_freqs_kHz = [0.1:0.1:4];
   K = length(rate_K_sample_freqs_kHz);
@@ -268,7 +252,7 @@ function [model_ rate_K_surface] = experiment_const_freq(model, varargin)
     if nvq < 2
       vq_filename = avq_filename;
     else
-      vq_filename = [vqfilname; avq_filename];
+      vq_filename = [vq_filename; 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");
@@ -277,6 +261,8 @@ function [model_ rate_K_surface] = experiment_const_freq(model, varargin)
     end
   end
   
+  quant_en = arg_exists(varargin, "quant");
+
   % OK start processing .....
 
   energy = zeros(1,frames);
@@ -360,7 +346,19 @@ function [model_ rate_K_surface] = experiment_const_freq(model, varargin)
       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);
+        [idx contrib errors b_log] = vq_search_para(vq, rate_K_surface_no_mean(:,vq_st:vq_en));
+        b_log = [b_log energy' idx'];
+        
+        if quant_en
+          for f=1:frames
+            target = rate_K_surface_no_mean(f,vq_st:vq_en);
+            % recalc gain using some index
+            v = mg(f)*vq(idx(f),1:vq_cols-10);
+            g = (sum(target(1:vq_cols-10)) - sum(v))/length(target);
+            contrib(f,:) = mg(f)*vq(idx(f),:) + g;
+            contrib(f, vq_cols-4:vq_cols) += -3:-6:-30;
+          end
+        end
       end
 
       rate_K_surface_no_mean_(:, vq_st:vq_en) = contrib;
@@ -419,6 +417,7 @@ function [model_ rate_K_surface] = experiment_const_freq(model, varargin)
   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
 
 
diff --git a/codec2-dev/octave/newamp1_fbf.m b/codec2-dev/octave/newamp1_fbf.m
index 0800a5aa..224d3519 100644
--- a/codec2-dev/octave/newamp1_fbf.m
+++ b/codec2-dev/octave/newamp1_fbf.m
@@ -25,7 +25,7 @@ function newamp1_fbf(samname, f=73, varargin)
   quant_en = 0; vq_search = "gain";
   mask_en = 0;
   nvec = 0;
-  weight_en = 0;
+  weight_en = quant_en = 0;
 
   % optional VQ
 
@@ -79,46 +79,53 @@ function newamp1_fbf(samname, f=73, varargin)
 
     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;
+      slope = 0; meanf = mean(rate_K_vec); 
+      rate_K_vec_fit = rate_K_vec - meanf;
     end
 
     % plots ----------------------------------
   
-    figure(2); clf; 
-    plot((1:L)*Wo*4000/pi, AmdB,";AmdB;g+-");
+    figure(2); clf;
+    l = sprintf(";rate %d AmdB;g+-", L);
+    plot((1:L)*Wo*4000/pi, AmdB, l);
     axis([1 4000 -20 80]);
     hold on;
     plot(rate_K_sample_freqs_kHz*1000, rate_K_vec, ";rate K;b+-");
 
-    if mask_en
+    % 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);
       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
+    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
+      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
+        % 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;
+        gradient = 0.05; yint = 1 - gradient*min(target);
+        weights = gradient*target + yint;
+      end
     end
-
+    
     if nvec
       
       if strcmp(vq_search, "mse")
@@ -138,24 +145,38 @@ function newamp1_fbf(samname, f=73, varargin)
 
       if strcmp(vq_search, "para")
         printf("\n");
-        [idx contrib errors test_ g mg sl] = vq_search_para(vq, target, weights);
+        [idx contrib errors b] = vq_search_para(vq, target);
+        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;
+        end
         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;
-      rate_K_vec_ += b;
-      [model_ AmdB_] = resample_rate_L(model(f,:), rate_K_vec_, rate_K_sample_freqs_kHz, Fs);
-      AmdB_ = AmdB_(1:L);
+      rate_K_vec_(vq_st:vq_en) = contrib;
 
-      sdL = std(AmdB - AmdB_);
       plot(rate_K_sample_freqs_kHz(vq_st:vq_en)*1000, contrib, 'm+-');
-      l = sprintf(";diff sd = %3.2f;k+-", sdL);
+      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);
-      plot((1:L)*Wo*4000/pi, AmdB_,";AmdB bar;r+-");
-      hold off;
     end
 
+    % And .... back to rate L
+    
+    rate_K_vec_ += meanf;
+    [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_));
+
+    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);
+    end
+    hold off;
+
     if weight_en
       figure(3); clf;
       subplot(211);
@@ -171,14 +192,14 @@ function newamp1_fbf(samname, f=73, varargin)
 
     % interactive menu ------------------------------------------
 
-    printf("\rframe: %d  menu: n-next  b-back  q-quit  w-weight[%d]", f, weight_en);
+    printf("\rframe: %d  menu: n-next  b-back  q-quit  w-quant[%d]", f, quant_en);
     fflush(stdout);
     k = kbhit();
 
     if k == 'w'
-      weight_en++;
-      if weight_en == 2
-        weight_en = 0;
+      quant_en++;
+      if quant_en == 2
+        quant_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
index 6455f519..f01bf775 100644
--- a/codec2-dev/octave/vq_search_para.m
+++ b/codec2-dev/octave/vq_search_para.m
@@ -1,11 +1,11 @@
 %----------------------------------------------------------------------
 % abs() search with a linear, ampl scaling, and slope term
 
-function [idx contrib errors test_ g mg sl] = vq_search_para(vq, data)
+function [idx contrib errors b_log2] = vq_search_para(vq, data)
   [nVec nCols] = size(vq);
   nRows = rows(data);
   
-  g = mg = sl = zeros(nRows, nVec);
+  g = mg = sl = zeros(nRows, 1);
   diff  = zeros(nVec, nCols);
   idx = errors = zeros(1, nRows);
   error = zeros(1, nVec);
@@ -19,6 +19,7 @@ function [idx contrib errors test_ g mg sl] = vq_search_para(vq, data)
   
   
   b_log = zeros(nVec, 4);
+  b_log2 = [];
   
   for i=1:nVec
     v = vq(i,:);
@@ -35,7 +36,9 @@ function [idx contrib errors test_ g mg sl] = vq_search_para(vq, data)
     for i=1:nVec
       v = vq(i,:);      
       c = [t*v' t*k2' t*k' sum(t)]';
-      b = inv(A(:,:,i))*c;    
+      b = inv(A(:,:,i))*c;
+      % b(1) = max(b(1),0.5);
+      % b(2) = max(b(2),-0.2); b(2) = min(b(2),0.1);
       diff(i,:) = t - (b(1)*v + b(2)*k2 + b(3)*k + b(4));
       b_log(i,:) = b; 
       error(i) = diff(i,:) * diff(i,:)';
@@ -50,7 +53,8 @@ function [idx contrib errors test_ g mg sl] = vq_search_para(vq, data)
     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);
+    contrib(f,:) = b(1)*v + b(2)*k2 + b(3)*k + b(4);
+    b_log2 = [b_log2; b];
   end
 
 endfunction