first pass at quantising dct coeffs using kmeans designed quantisers, doesnt sound...

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

diff --git a/codec2-dev/octave/newamp.m b/codec2-dev/octave/newamp.m
index b6024723cef8a3adce1f636fffc069837f687699..1af8a4ace011d15ea18afd4b197de4faaea31b26 100644 (file)
--- a/codec2-dev/octave/newamp.m
+++ b/codec2-dev/octave/newamp.m
@@ -202,7 +202,7 @@ endfunction
 
 % generate a zig-zag linear to square mapping matrix
 
-function map = zigzag(nr,nc)
+function map = create_zigzag_map(nr,nc)
   map = zeros(nr, nc);
 
   state = 'zig';
@@ -244,6 +244,8 @@ function map = zigzag(nr,nc)
 endfunction
 
 
+% reshape matrix m as a vector v by reading out elements in zig-zag pattern
+
 function v = mtov_zigzag(m)
   [nr nc] = size(m);
   map = zigzag(nr,nc)
@@ -256,6 +258,68 @@ function v = mtov_zigzag(m)
 endfunction
 
 
+% extracts DCT information for rate K surface
+
+function unwrapped_dcts = dct_blocks(surf)
+  [frames K] = size(surf);
+
+  % break into 160ms blocks, 2D DCT, truncate, IDCT
+
+  Nt = 16;  % 16 blocks in time
+  Nblocks = floor(frames/Nt);
+  unwrapped_dcts = zeros(Nblocks,Nt*K);
+
+  for n=1:Nblocks
+    st = (n-1)*Nt+1; en = st + Nt - 1;
+    D = dct2(surf(st:en,:));
+    unwrapped_dcts(n,:) = reshape(D',1,Nt*K);
+  end
+endfunction
+
+
+% Determines a map for quantising 2D DCT coeffs in order of rms value
+% (ie std dev) of each coeff.  Those coeffs with the greatest change
+% need the most bits to quantise
+
+function [map rms_map mx unwrapped_dcts] = create_map_rms(rate_K_surface, nr, nc)
+  unwrapped_dcts = dct_blocks(rate_K_surface);
+  [mx mx_ind] = sort(std(unwrapped_dcts));
+  mx_ind = fliplr(mx_ind); mx = fliplr(mx);
+  map = rms_map = zeros(nr,nc);
+  for i=1:nr*nc
+    r = floor((mx_ind(i)-1)/nc) + 1;
+    c = mx_ind(i) - (r-1)*nc;
+    %printf("%d %d %d\n", i, r, c);
+    map(r,c) = i;
+    rms_map(r,c) = mx(i);
+  end
+endfunction
+
+
+% Design quantisers for each DCT coeff
+
+function [quantisers nbits] = design_quantisters(rate_K_surface, nr, nc, nbits_max)
+  [map rms_map mx unwrapped_dcts] = create_map_rms(rate_K_surface, nr, nc);
+  nq = nr*nc;
+  quantisers = zeros(nq, 2^nbits_max); nbits = zeros(nq,1);
+  for i=1:nq
+
+    % work out number of levels for this quantiser such that it is a
+    % power of 2 for integer number of bits
+
+    nlevels = (2^nbits_max);
+    nbits(i) = round(log2(nlevels));
+    nlevels = 2 .^ nbits(i);
+
+    if i <= 100
+      printf("%d %d\n", i, nlevels);
+      [idx, centers] = kmeans(unwrapped_dcts(:,i), nlevels);
+      quantisers(i,1:nlevels) = sort(centers);
+    end
+  end
+endfunction
+
+
 % Determine a phase spectra from a magnitude spectra
 % from http://www.dsprelated.com/showcode/20.php
 % Haven't _quite_ figured out how this works but have to start somewhere ....

diff --git a/codec2-dev/octave/newamp1_batch.m b/codec2-dev/octave/newamp1_batch.m
index 3dbbbc7cdfce2c1a08da0182a05ea3b6e54d0924..9901f7df86f94bb6cb6820178c60617e0c3a26c3 100644 (file)
--- a/codec2-dev/octave/newamp1_batch.m
+++ b/codec2-dev/octave/newamp1_batch.m
@@ -70,6 +70,8 @@ function surface = newamp1_batch(input_prefix, output_prefix)
 
   [model_ surface] = experiment_mel_freq(model, 0, 1, voicing);
 
+  %model_  = experiment_smoothed(model, 0);
+
   %model_ = experiment_dec_abys(model, 8, 1, 1, 1, voicing);
 
 #{
@@ -488,7 +490,6 @@ endfunction
 
 
 
-
 % ---------------------------------------------------------------------------------------
 % Various experiments tried during development
 
@@ -496,7 +497,11 @@ endfunction
 
 function [model_ rate_K_surface] = experiment_mel_freq(model, vq_en=0, plots=1, voicing)
   [frames nc] = size(model);
-  K = 20; Fs = 8000;
+  K = 20; Fs = 8000; correct_rate_K_en = 1;
+
+  quantisers = load("dct2quant.txt"); 
+  nbits = load("dct2quantnbits.txt");
+  map = load("dct2map.txt");
 
   for f=1:frames
     Wo = model(f,1);
@@ -505,15 +510,17 @@ function [model_ rate_K_surface] = experiment_mel_freq(model, vq_en=0, plots=1,
     AmdB = 20*log10(Am);
     Am_freqs_kHz = (1:L)*Wo*4/pi;
     [rate_K_vec rate_K_sample_freqs_kHz] = resample_const_rate_f_mel(model(f,:), K);
-    [tmp_ AmdB_] = resample_rate_L(model(f,:), rate_K_vec, rate_K_sample_freqs_kHz);
-    [rate_K_vec_corrected orig_error error nasty_error_log nasty_error_m_log] = correct_rate_K_vec(rate_K_vec, rate_K_sample_freqs_kHz, AmdB, AmdB_, K, Wo, L, Fs);
-    %rate_K_surface(f,:) = rate_K_vec;
-    rate_K_surface(f,:) = rate_K_vec_corrected;
+    if correct_rate_K_en
+      [tmp_ AmdB_] = resample_rate_L(model(f,:), rate_K_vec, rate_K_sample_freqs_kHz);
+      [rate_K_vec_corrected orig_error error nasty_error_log nasty_error_m_log] = correct_rate_K_vec(rate_K_vec, rate_K_sample_freqs_kHz, AmdB, AmdB_, K, Wo, L, Fs);
+      rate_K_surface(f,:) = rate_K_vec_corrected;
+    else
+      rate_K_surface(f,:) = rate_K_vec;
+    end
   end
 
   if plots
     figure(1); clf; mesh(rate_K_surface);
-    figure(2); clf; plot_dft_surface(rate_K_surface)
   end
 
   for f=1:frames
@@ -538,16 +545,50 @@ function [model_ rate_K_surface] = experiment_mel_freq(model, vq_en=0, plots=1,
     rate_K_surface = rate_K_surface_;
   end
 
+
+  % break into 160ms blocks, 2D DCT, truncate, IDCT
+
+  Nt = 16;  % 16 blocks in time
+  Nblocks = floor(frames/Nt);
+  unwrapped_dcts = zeros(Nblocks,Nt*K);
+
+  for n=1:Nblocks
+    st = (n-1)*Nt+1; en = st + Nt - 1;
+    D = dct2(rate_K_surface(st:en,:));
+
+    % move over surface and work out quantiser
+    % quantise, replace on map
+
+    E = zeros(Nt,K);
+    E = D;
+    for r=1:Nt
+      for c=1:K
+        quantiser_num = map(r,c);
+        nlevels = 2^nbits(quantiser_num);
+        if quantiser_num < 10
+          printf("r %d c %d quantiser_num %d nlevels %d\n", r,c,quantiser_num,nlevels);
+          levels = quantisers(quantiser_num,1:nlevels);
+          quant_out = quantise(levels, D(r,c));
+          printf("%f %f\n", D(r,c), quant_out );
+          E(r,c) = quant_out;
+        end
+      end
+    end
+
+    rate_K_surface(st:en,:) = idct2(E);
+  end
+
   model_ = resample_rate_L(model, rate_K_surface, rate_K_sample_freqs_kHz);
 
+  #{
   if plots
-    figure(3); clf; mesh(rate_K_surface_no_mean);
-    figure(4); clf; plot_dft_surface(rate_K_surface)
+    figure(2); clf; mesh(rate_K_surface_no_mean);
   end
 
   for f=1:frames
     rate_K_surface(f,:) -= mean(rate_K_surface(f,:));
   end
+  #}
 endfunction
 
 
@@ -938,6 +979,24 @@ function model_ =  experiment_filter_dec_filter(model)
 endfunction
 
 
+% smoothing using masking functions
+
+function model_ = experiment_smoothed(model, bark_model=1)
+  [frames nc] = size(model);
+
+  model_ = model;
+
+  for f=1:frames
+    Wo = model(f,1);
+    L = model(f,2);
+    Am = model(f,3:(L+2));
+    AmdB = 20*log10(Am);
+    [AmdB_ Am_freqs_kHz] = mask_model(AmdB, Wo, L, bark_model);
+    model_(f,3:(L+2)) = 10 .^ ((AmdB_)/20);
+  end
+
+endfunction
+
 #{
   
   My original idea was to used a 3-4 "resonators" to construct a