From: drowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Date: Thu, 22 Dec 2016 06:00:07 +0000 (+0000)
Subject: parabolic interp working ok on hts1a, hts2a, vk5qi
X-Git-Url: http://git.whiteaudio.com/gitweb/?a=commitdiff_plain;h=5e09443fd529eb7b8cfbbff47e8987ee53a0b043;p=freetel-svn-tracking.git

parabolic interp working ok on hts1a, hts2a, vk5qi

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

diff --git a/codec2-dev/octave/newamp.m b/codec2-dev/octave/newamp.m
index 5e152795..3609b9eb 100644
--- a/codec2-dev/octave/newamp.m
+++ b/codec2-dev/octave/newamp.m
@@ -14,9 +14,36 @@
 melvq; % mbest VQ functions
 
 % --------------------------------------------------------------------------------
-% Function used by rate K mel work
+% Functions used by rate K mel work
 % --------------------------------------------------------------------------------
 
+% General 2nd order parabolic interpolator.  Used splines orginally,
+% but this is much simpler and we don't need much accuracy.  Given two
+% vectors of points xp and yp, find interpolated values y at points x
+
+function y = interp_para(xp, yp, x)
+  assert( (length(xp) >=3) && (length(yp) >= 3) );
+
+  y = zeros(1,length(x));
+  for i=1:length(x)
+    xi = x(i);
+
+    % k is index into xp of where we start 3 points used to form parabola
+
+    k = 1;
+    while ((xp(k+1) < xi) && (k < (length(xp)-2))) 
+      k++;
+    end
+    
+    x1 = xp(k); y1 = yp(k); x2 = xp(k+1); y2 = yp(k+1); x3 = xp(k+2); y3 = yp(k+2);
+
+    a = ((y3-y2)/(x3-x2)-(y2-y1)/(x2-x1))/(x3-x1);
+    b = ((y3-y2)/(x3-x2)*(x2-x1)+(y2-y1)/(x2-x1)*(x3-x2))/(x3-x1);
+  
+    y(i) = a*(xi-x2)^2 + b*(xi-x2) + y2;
+  end
+endfunction
+
 
 % quantise input sample to nearest value in table, optionally return binary code
 
@@ -150,7 +177,6 @@ function [rate_K_surface rate_K_sample_freqs_kHz] = resample_const_rate_f(model,
   K = length(rate_K_sample_freqs_kHz);
   rate_K_surface = zeros(frames, K);
 
-
   for f=1:frames
     Wo = model(f,1);
     L = min([model(f,2) max_amp-1]);
@@ -166,8 +192,9 @@ function [rate_K_surface rate_K_sample_freqs_kHz] = resample_const_rate_f(model,
 
     rate_L_sample_freqs_kHz = (1:L)*Wo*4/pi;
     
-    rate_K_surface(f,:) = interp1(rate_L_sample_freqs_kHz, AmdB, rate_K_sample_freqs_kHz, "spline", "extrap");
-    
+    %rate_K_surface(f,:) = interp1(rate_L_sample_freqs_kHz, AmdB, rate_K_sample_freqs_kHz, "spline", "extrap");
+    rate_K_surface(f,:)  = interp_para(rate_L_sample_freqs_kHz, AmdB, rate_K_sample_freqs_kHz);    
+
     %printf("\r%d/%d", f, frames);
   end
   %printf("\n");
@@ -188,7 +215,8 @@ function [model_ AmdB_] = resample_rate_L(model, rate_K_surface, rate_K_sample_f
     
     % back down to rate L
 
-    AmdB_ = interp1(rate_K_sample_freqs_kHz, rate_K_surface(f,:), rate_L_sample_freqs_kHz, "spline", 0);
+    % AmdB_ = interp1(rate_K_sample_freqs_kHz, rate_K_surface(f,:), rate_L_sample_freqs_kHz, "spline", 0);
+    AmdB_ = interp_para([ 0 rate_K_sample_freqs_kHz 4], [0 rate_K_surface(f,:) 0], rate_L_sample_freqs_kHz);
 
     model_(f,1) = Wo; model_(f,2) = L; model_(f,3:(L+2)) = 10 .^ (AmdB_(1:L)/20);
    end
diff --git a/codec2-dev/octave/newamp1_batch.m b/codec2-dev/octave/newamp1_batch.m
index 63376553..ef537e01 100644
--- a/codec2-dev/octave/newamp1_batch.m
+++ b/codec2-dev/octave/newamp1_batch.m
@@ -19,6 +19,8 @@
 
   Usage:
 
+    build codec2 with -DDUMP - see codec2-dev/README, then:
+
     ~/codec2-dev/build_linux/src$ ./c2sim ../../raw/hts1a.raw --dump hts1a
     $ cd ~/codec2-dev/octave
     octave:14> newamp1_batch("../build_linux/src/hts1a")
@@ -155,7 +157,13 @@ endfunction
 % -----------------------------------------------------------------------------------------
 % Linear decimator/interpolator that operates at rate K, includes VQ, post filter, and Wo/E
 % quantisation.  Evolved from abys decimator below.  Simulates the entire encoder/decoder.
- 
+
+#{
+   todo: 
+     [ ] M=8 nomenclature, make it closer to existing 25ms C modes
+     [ ] frame by frame-ish operation, close to C implementations
+#}
+
 function [model_ voicing_ indexes] = experiment_rate_K_dec(model, voicing)
   max_amp = 80;
   [frames nc] = size(model);
@@ -171,6 +179,9 @@ function [model_ voicing_ indexes] = experiment_rate_K_dec(model, voicing)
   [surface sample_freqs_kHz] = resample_const_rate_f_mel(model, K);
   target_surface = surface;
 
+  figure(1);
+  mesh(surface);
+
   % VQ rate K surface.  TODO: If we are decimating by M/2=4 we really
   % only need to do this every 4th frame.
 
@@ -181,6 +192,8 @@ function [model_ voicing_ indexes] = experiment_rate_K_dec(model, voicing)
     mean_f(f) = mean(surface(f,:));
     surface_no_mean(f,:) = surface(f,:) - mean_f(f);
   end
+  figure(2);
+  mesh(surface_no_mean);
 
   [res surface_no_mean_ ind] = mbest(train_120_vq, surface_no_mean, m);
   indexes(:,1:2) = ind;
@@ -188,6 +201,8 @@ function [model_ voicing_ indexes] = experiment_rate_K_dec(model, voicing)
   for f=1:frames
     surface_no_mean_(f,:) = post_filter(surface_no_mean_(f,:), sample_freqs_kHz, 1.5);
   end
+  figure(3);
+  mesh(surface_no_mean_);
     
   surface_ = zeros(frames, K);
   energy_q = 10 + 40/16*(0:15);
@@ -198,6 +213,9 @@ function [model_ voicing_ indexes] = experiment_rate_K_dec(model, voicing)
     surface_(f,:) = surface_no_mean_(f,:) + mean_f_;
   end
 
+  figure();
+  mesh(surface_);
+
   % break into segments of M frames.  We have 3 samples in M frame
   % segment spaced M/2 apart and interpolate the rest.  This evolved
   % from AbyS scheme below but could be simplified to simple linear