From 1ba4c2a24aa6367094c81f7e8832b9ca3f7ca25a Mon Sep 17 00:00:00 2001
From: drowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Date: Thu, 31 Mar 2016 05:07:40 +0000
Subject: [PATCH] VQ working for mask/circ/resampling idea, sounds OK on 14
 bit/sfr, using 10ms frames and orig phase

git-svn-id: https://svn.code.sf.net/p/freetel/code@2766 01035d8c-6547-0410-b346-abe4f91aad63
---
 codec2-dev/octave/newamp.m       | 33 ++++++++++++++++++++++++++++----
 codec2-dev/octave/newamp_batch.m | 14 ++++++++++----
 codec2-dev/octave/newamp_fbf.m   | 11 ++++++++++-
 3 files changed, 49 insertions(+), 9 deletions(-)

diff --git a/codec2-dev/octave/newamp.m b/codec2-dev/octave/newamp.m
index 24da35ae..f1dd6ae9 100644
--- a/codec2-dev/octave/newamp.m
+++ b/codec2-dev/octave/newamp.m
@@ -32,7 +32,7 @@
 melvq;
 
 
-function [maskdB_ maskdB_cyclic Dabs dk_] = decimate_in_freq(maskdB, cyclic=1, Nkeep=7)
+function [maskdB_ maskdB_cyclic Dabs dk_] = decimate_in_freq(maskdB, cyclic=1, k=7, vq)
 
     % Lets try to come up with a smoothed, cyclic model.  Replace
     % points from 3500 Hz to 4000Hz with a sequence that joins without
@@ -55,14 +55,39 @@ function [maskdB_ maskdB_cyclic Dabs dk_] = decimate_in_freq(maskdB, cyclic=1, N
       D = fft(maskdB);
     end
     Dabs = abs(D);                        % this returned for plotting
+
+if 0
     D_ = D; % + 10*randn(1,L) + 10*j*randn(1,L);
     D_(Nkeep+1:L-Nkeep) = 0;              % truncate
     d = ifft(D_);                         % back to spectrum at rate L
     maskdB_ = real(d);
-    Dk_ = [0 D(2:Nkeep-1) real(D(Nkeep)) D(L-Nkeep+1:L)];  % build rate Nkeep vector for quantisation
-    dk_ = real(ifft(Dk_));
+end
+
+    % truncate D to rate k, convert to 2k length real vector for quantisation and transmission
+
+    Dk = [0 D(2:k-1) real(D(k)) D(L-k+1:L)]; 
+    dk = real(ifft(Dk));
+    
+    % quantisation
 
-    % OK now fix up last 500Hz, taper down 10dB at 4000Hz
+    if nargin == 4
+       [res dk_ ] = mbest(vq, dk, 4);
+       std(dk_ - dk)
+    else
+       dk_ = dk;
+    end
+
+    % convert quantised dk back to rate L magnitude spectrum
+
+    Dk_ = fft(dk_);
+    D_ = zeros(1,L);
+    D_(1) = D(1);                         % lets assume energy comes through separately
+    D_(2:k-1) = Dk_(2:k-1);
+    D_(L-k+1:L) = Dk_(k+1:2*k);
+    d_ = ifft(D_);                        % back to spectrum at rate L
+    maskdB_ = real(d_);
+    
+    % Finally fix up last 500Hz, taper down 10dB at 4000Hz
 
     xpts = [ anchor-1 anchor L];
     ypts = [ maskdB_(anchor-1) maskdB_(anchor) (maskdB_(anchor)-10)];
diff --git a/codec2-dev/octave/newamp_batch.m b/codec2-dev/octave/newamp_batch.m
index 796a31d8..0659d66d 100644
--- a/codec2-dev/octave/newamp_batch.m
+++ b/codec2-dev/octave/newamp_batch.m
@@ -25,8 +25,7 @@ function dk_log = newamp_batch(samname, optional_Am_out_name, optional_Aw_out_na
   postfilter = 1;
   dec_in_time = 0;
   synth_phase = 0;
-  freq_quant = 0;
-  amp_quant = 0;
+  vq_en = 1;
   dk_log = [];
 
   model_name = strcat(samname,"_model.txt");
@@ -52,9 +51,12 @@ function dk_log = newamp_batch(samname, optional_Am_out_name, optional_Aw_out_na
     faw = fopen(Aw_out_name,"wb"); 
   end
 
+  if vq_en
+    load vq;
+  end
+
   % encoder loop ------------------------------------------------------
 
-  %pp_bw = gen_pp_bw;
 
   sd_sum = 0;
   for f=1:frames
@@ -77,7 +79,11 @@ function dk_log = newamp_batch(samname, optional_Am_out_name, optional_Aw_out_na
     non_masked_m(f,1:length(a_non_masked_m)) = a_non_masked_m;
 
     if dec_in_freq
-      [maskdB_ tmp1 D dk_] = decimate_in_freq(maskdB, 1);
+      if vq_en
+        [maskdB_ tmp1 D dk_] = decimate_in_freq(maskdB, 1, 7, vq);
+      else
+        [maskdB_ tmp1 D dk_] = decimate_in_freq(maskdB, 1);
+      end
       dk_log = [dk_log; dk_];
     end
     sd_sum += sum(maskdB - maskdB_);
diff --git a/codec2-dev/octave/newamp_fbf.m b/codec2-dev/octave/newamp_fbf.m
index 17e6fbb4..64f9ebf7 100644
--- a/codec2-dev/octave/newamp_fbf.m
+++ b/codec2-dev/octave/newamp_fbf.m
@@ -20,6 +20,7 @@ function newamp_fbf(samname, f=10)
   newamp;
   more off;
   plot_spectrum = 1;
+  vq_en = 1;
 
   % load up text files dumped from c2sim ---------------------------------------
 
@@ -31,6 +32,10 @@ function newamp_fbf(samname, f=10)
   model = load(model_name);
   [frames tmp] = size(model);
 
+  if vq_en
+    load vq;
+  end
+
   % Keyboard loop --------------------------------------------------------------
 
   k = ' ';
@@ -58,7 +63,11 @@ function newamp_fbf(samname, f=10)
 
     [maskdB Am_freqs_kHz] = mask_model(AmdB, Wo, L);
     [tmp1 tmp2 D] = decimate_in_freq(maskdB, 0);
-    [maskdB_ maskdB_cyclic D_cyclic dk_] = decimate_in_freq(maskdB, 1);
+    if vq_en
+      [maskdB_ maskdB_cyclic D_cyclic dk_] = decimate_in_freq(maskdB, 1, 7, vq);
+    else
+      [maskdB_ maskdB_cyclic D_cyclic dk_] = decimate_in_freq(maskdB, 1);
+    end
 
     plot(Am_freqs_kHz*1000, maskdB, ';mask;g');
     plot(Am_freqs_kHz*1000, maskdB_cyclic, ';mask cyclic;b');
-- 
2.25.1