refactoring
authordrowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Wed, 21 Dec 2016 07:58:37 +0000 (07:58 +0000)
committerdrowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Wed, 21 Dec 2016 07:58:37 +0000 (07:58 +0000)
git-svn-id: https://svn.code.sf.net/p/freetel/code@2940 01035d8c-6547-0410-b346-abe4f91aad63

codec2-dev/octave/newamp1_batch.m

index b189ddd30f134736eeae9dcc62ee6bdbf55949d0..63376553111da154770c9d0fd7f94ff3bd1414a9 100644 (file)
@@ -3,13 +3,20 @@
 % Copyright David Rowe 2016
 % This program is distributed under the terms of the GNU General Public License 
 % Version 2
-%
-% Octave script to batch process model parameters using the new
-% amplitude model.  Used for generating samples we can listen to.
-%
 
 #{
 
+  Octave script to batch process model parameters using the new
+  amplitude model, version 1.  Outputs anotehr set of model parameters
+  that can be fed to c2sim for listening tests.  The companion
+  newamp1_fbf.m script is used to visualise the processing frame by frame
+  c2sim -> dump files -> newamp1_batch.m -> output model params -> c2sim -> play
+  The newamp1_xxx scripts have evolved to (i) resample {Am} using a
+  mel frequency axis, (ii) 2 stage VQ the mean removed vector.  Seems to work
+  OK at 700 bit/s.
+
   Usage:
 
     ~/codec2-dev/build_linux/src$ ./c2sim ../../raw/hts1a.raw --dump hts1a
@@ -17,7 +24,9 @@
     octave:14> newamp1_batch("../build_linux/src/hts1a")
     ~/codec2-dev/build_linux/src$ ./c2sim ../../raw/hts1a.raw --amread hts1a_am.out -o - | play -t raw -r 8000 -s -2 -
  
-  Or with a little more processing, first dump energy and voicing, the import Wo, voicing, phase spectra:
+  Or with a little more processing, first dump energy and voicing, the
+  import Wo, voicing, phase spectra which simulates all the decoder
+  DSP.  We switch on lpc 10 just to dump voicing.
 
     $ ./c2sim ../../raw/vk5qi.raw --phase0 --postfilter --dump vk5qi --lpc 10 --dump_pitch_e vk5qi_pitche.txt
     octave:14> newamp1_batch("../build_linux/src/vk5qi", "../build_linux/src/vk5qi_am_Wo.out");
 #}
 
 
-% process a whole file and write results
-% TODO: 
-%   [ ] refactor decimate-in-time to avoid extra to/from model conversions
-%   [ ] switches to turn on/off quantisation
-%   [ ] rename mask_sample_freqs, do we need "mask" any more
-
-#{
-    [ ] refactor
-#}
-
 % In general, this function processes a bunch of amplitudes, we then
-% use c2sim to hear the results
+% use c2sim to hear the results.  Bunch of different experiments below
 
 function surface = newamp1_batch(samname, optional_Am_out_name, optional_Aw_out_name)
   newamp;
   more off;
 
   max_amp = 80;
-  postfilter = 0;
+  postfilter = 0;   % optional postfiler that runs on Am, not used atm
   synth_phase = 1;
 
   model_name = strcat(samname,"_model.txt");
@@ -66,8 +65,9 @@ function surface = newamp1_batch(samname, optional_Am_out_name, optional_Aw_out_
 
   %[model_ surface] = experiment_mel_freq(model, 1, 1, voicing);
   %model_ = experiment_dec_abys(model, 8, 1, 1, 1, voicing);
-  [amodel_ avoicing_ indexes] = experiment_rate_K_dec(model, voicing);
-  [model_ voicing_] = model_from_indexes(indexes);
+
+  [amodel_ avoicing_ indexes] = experiment_rate_K_dec(model, voicing); % encoder, toss away results except for indexes
+  [model_ voicing_] = model_from_indexes(indexes);                     % decoder uses just indexes, outputs vecs for synthesis
 
   %model_ = experiment_dec_linear(model_);
   %model_ = experiment_energy_rate_linear(model, 1, 0);
@@ -152,193 +152,9 @@ function surface = newamp1_batch(samname, optional_Am_out_name, optional_Aw_out_
 endfunction
  
 
-% experiment to resample freq axis on mel scale, then optionally vq
-
-function [model_ rate_K_surface] = experiment_mel_freq(model, vq_en=0, plots=1, voicing)
-  [frames nc] = size(model);
-  K = 20; 
-  [rate_K_surface rate_K_sample_freqs_kHz] = resample_const_rate_f_mel(model, K);
-  
-  if plots
-    figure(1); clf; mesh(rate_K_surface);
-    figure(2); clf; plot_dft_surface(rate_K_surface)
-  end
-
-  for f=1:frames
-    mean_f(f) = mean(rate_K_surface(f,:));
-    rate_K_surface_no_mean(f,:) = rate_K_surface(f,:) - mean_f(f);
-  end
-
-  if vq_en
-    melvq;
-    load train_120_vq; m=5;
-       
-    [res rate_K_surface_ ind] = mbest(train_120_vq, rate_K_surface_no_mean, m);
-
-    for f=1:frames
-        rate_K_surface_(f,:) = post_filter(rate_K_surface_(f,:), 1.5, voicing(f));
-    end
-
-    for f=1:frames
-      rate_K_surface_(f,:) += mean_f(f);
-    end
-
-    rate_K_surface = rate_K_surface_;
-  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)
-  end
-
-  for f=1:frames
-    rate_K_surface(f,:) -= mean(rate_K_surface(f,:));
-  end
-endfunction
-
-
-% mel spaced sampling, differential in time VQ.  Curiously, couldn't
-% get very good results out of this, I suspect a bug
-
-function [model_ rate_K_surface_diff] = experiment_mel_diff_freq(model, vq_en=0)
-  [frames nc] = size(model);
-  K = 20; 
-  [rate_K_surface rate_K_sample_freqs_kHz] = resample_const_rate_f_mel(model, K);
-  
-  if vq_en
-    melvq;
-    load surface_diff_vq; m=5;
-  end
-
-  for f=1:frames
-    mean_f(f,:) = mean(rate_K_surface(f,:));
-    rate_K_surface_no_mean(f,:) = rate_K_surface(f,:) - mean_f(f,:);
-  end
-
-  rate_K_surface_no_mean_ = zeros(frames, K);
-  rate_K_surface_no_mean_diff = zeros(frames, K);
-  rate_K_surface_(1,:) = rate_K_surface_diff(1,:) = zeros(1, K);
-
-  for f=2:frames
-    rate_K_surface_diff(f,:) = rate_K_surface_no_mean(f,:) - 0.8*rate_K_surface_no_mean_(f-1,:);
-    if vq_en
-      [res arate_K_surface_diff_ ind] = mbest(surface_diff_vq, rate_K_surface_diff(f,:), m);
-      rate_K_surface_diff_(f,:) = arate_K_surface_diff_;
-    else
-      rate_K_surface_diff_(f,:) = rate_K_surface_diff(f,:);
-    end
-    rate_K_surface_no_mean_(f,:) = 0.8*rate_K_surface_no_mean_(f-1,:) + rate_K_surface_diff_(f,:);
-  end
-  
-  for f=1:frames
-    rate_K_surface_(f,:) = rate_K_surface_no_mean_(f,:) + mean_f(f,:);
-  end
-
-  model_ = resample_rate_L(model, rate_K_surface_, rate_K_sample_freqs_kHz);
-
-endfunction
-
-
-% try vq with open and closed loop mean removal, turns out they give
-% identical results lol
-
-function [model_ rate_K_surface] = experiment_closed_loop_mean(model)
-  [frames nc] = size(model);
-  K = 15; 
-  mel_start = ftomel(300); mel_end = ftomel(3000); 
-  step = (mel_end-mel_start)/(K-1);
-  mel = mel_start:step:mel_end;
-  rate_K_sample_freqs_Hz = 700*((10 .^ (mel/2595)) - 1);
-  rate_K_sample_freqs_kHz = rate_K_sample_freqs_Hz/1000;
-
-  rate_K_surface = resample_const_rate_f(model, rate_K_sample_freqs_kHz);
-  
-  load surface_vq; m=1;
-   
-  for f=1:frames
-    amean = mean(rate_K_surface(f,:));
-    rate_K_target_no_mean = rate_K_surface(f,:) - amean;
-    mse_open_loop(f) = search_vq2(surface_vq(:,:,1), rate_K_target_no_mean, m, 0);
-    mse_closed_loop(f) = search_vq2(surface_vq(:,:,1), rate_K_surface(f,:), m, 1);
-  end
-
-  printf("rms open loop..: %f\nrms closed loop: %f\n", sqrt(mean(mse_open_loop)), sqrt(mean(mse_closed_loop)));
-
-  % just return model_ as we have to so nothing breaks, it's not actually useful
-
-  model_ = resample_rate_L(model, rate_K_surface, rate_K_sample_freqs_kHz);
-    
-endfunction
-
-
-% Experiment with 10ms update rate for energy but 40ms update for spectrum,
-% using linear interpolation for spectrum.
-
-function model_c = experiment_energy_rate_linear(model, vq_en, plot_en)
-  max_amp = 80;
-  [frames nc] = size(model);
-
-  % 10ms mel freq modelling and VQ
-
-  model_ = experiment_mel_freq(model, vq_en, plot_en);
-
-  % Remove energy.  Hmmmm, this is done on Ams rather than surface but that's
-  % similar I guess
-
-  e = zeros(1,frames);
-  model_a = zeros(frames,max_amp+3);
-  for f=1:frames
-    L = min([model_(f,2) max_amp-1]);
-    Am_ = model_(f,3:(L+2));
-    AmdB_ = 20*log10(Am_);
-    mean_f(f) = mean(AmdB_);
-    AmdB_ -= mean_f(f);
-    model_a(f,1) = model_(f,1); model_a(f,2) = L; model_a(f,3:(L+2)) = 10 .^ (AmdB_(1:L)/20);
-  end
-
-  % linear interp after removing energy (mean AmdB)
-
-  model_b = experiment_dec_linear(model_a);
-
-  % add back in energy
-
-  model_c = zeros(frames,max_amp+3);
-  for f=1:frames
-    L = min([model_b(f,2) max_amp-1]);
-    Am_ = model_b(f,3:(L+2));
-    AmdB_ = 20*log10(Am_);
-    AmdB_ += mean_f(f);
-    model_c(f,1) = model_b(f,1); model_c(f,2) = L; model_c(f,3:(L+2)) = 10 .^ (AmdB_(1:L)/20);
-  end
-
-endfunction
-
-
-% conventional decimation in time without any filtering, then linear
-% interpolation.  Linear interpolation is a two-tap (weak) form of fir
-% filtering that may have problems with signals with high freq
-% components, for example after quantisation noise is added.  Need to
-% look into this some more.
-function model_ = experiment_dec_linear(model)
-  newamp;
-  max_amp = 80;
-
-  [frames nc] = size(model);
-  model_ = zeros(frames, max_amp+3);
-  decimate = 4;
-  for f=1:frames    
-    AmdB_ = decimate_frame_rate(model, decimate, f, frames);
-    L = length(AmdB_);
-    model_(f,1) = model(f,1); model_(f,2) = L; model_(f,3:(L+2)) = 10 .^ (AmdB_(1:L)/20);
-  end
-endfunction
-
-
+% -----------------------------------------------------------------------------------------
 % Linear decimator/interpolator that operates at rate K, includes VQ, post filter, and Wo/E
-% quantisation.  Evolved from abys decimator below.
+% quantisation.  Evolved from abys decimator below.  Simulates the entire encoder/decoder.
  
 function [model_ voicing_ indexes] = experiment_rate_K_dec(model, voicing)
   max_amp = 80;
@@ -475,7 +291,9 @@ function [model_ voicing_ indexes] = experiment_rate_K_dec(model, voicing)
 endfunction
 
 
-% Stand alone decoder that takes indexes and creates model_
+% ---------------------------------------------------------------------------------------
+% Stand alone decoder that takes indexes and creates model_, models
+% decoder and an important step in proving everything works
  
 function [model_ voicing_] = model_from_indexes(indexes)
   max_amp = 80;
@@ -578,6 +396,194 @@ function [model_ voicing_] = model_from_indexes(indexes)
 endfunction
 
 
+% ---------------------------------------------------------------------------------------
+% Various experiments tried during development
+
+% experiment to resample freq axis on mel scale, then optionally vq
+
+function [model_ rate_K_surface] = experiment_mel_freq(model, vq_en=0, plots=1, voicing)
+  [frames nc] = size(model);
+  K = 20; 
+  [rate_K_surface rate_K_sample_freqs_kHz] = resample_const_rate_f_mel(model, K);
+  
+  if plots
+    figure(1); clf; mesh(rate_K_surface);
+    figure(2); clf; plot_dft_surface(rate_K_surface)
+  end
+
+  for f=1:frames
+    mean_f(f) = mean(rate_K_surface(f,:));
+    rate_K_surface_no_mean(f,:) = rate_K_surface(f,:) - mean_f(f);
+  end
+
+  if vq_en
+    melvq;
+    load train_120_vq; m=5;
+       
+    [res rate_K_surface_ ind] = mbest(train_120_vq, rate_K_surface_no_mean, m);
+
+    for f=1:frames
+        rate_K_surface_(f,:) = post_filter(rate_K_surface_(f,:), 1.5, voicing(f));
+    end
+
+    for f=1:frames
+      rate_K_surface_(f,:) += mean_f(f);
+    end
+
+    rate_K_surface = rate_K_surface_;
+  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)
+  end
+
+  for f=1:frames
+    rate_K_surface(f,:) -= mean(rate_K_surface(f,:));
+  end
+endfunction
+
+
+% mel spaced sampling, differential in time VQ.  Curiously, couldn't
+% get very good results out of this, I suspect a bug
+
+function [model_ rate_K_surface_diff] = experiment_mel_diff_freq(model, vq_en=0)
+  [frames nc] = size(model);
+  K = 20; 
+  [rate_K_surface rate_K_sample_freqs_kHz] = resample_const_rate_f_mel(model, K);
+  
+  if vq_en
+    melvq;
+    load surface_diff_vq; m=5;
+  end
+
+  for f=1:frames
+    mean_f(f,:) = mean(rate_K_surface(f,:));
+    rate_K_surface_no_mean(f,:) = rate_K_surface(f,:) - mean_f(f,:);
+  end
+
+  rate_K_surface_no_mean_ = zeros(frames, K);
+  rate_K_surface_no_mean_diff = zeros(frames, K);
+  rate_K_surface_(1,:) = rate_K_surface_diff(1,:) = zeros(1, K);
+
+  for f=2:frames
+    rate_K_surface_diff(f,:) = rate_K_surface_no_mean(f,:) - 0.8*rate_K_surface_no_mean_(f-1,:);
+    if vq_en
+      [res arate_K_surface_diff_ ind] = mbest(surface_diff_vq, rate_K_surface_diff(f,:), m);
+      rate_K_surface_diff_(f,:) = arate_K_surface_diff_;
+    else
+      rate_K_surface_diff_(f,:) = rate_K_surface_diff(f,:);
+    end
+    rate_K_surface_no_mean_(f,:) = 0.8*rate_K_surface_no_mean_(f-1,:) + rate_K_surface_diff_(f,:);
+  end
+  
+  for f=1:frames
+    rate_K_surface_(f,:) = rate_K_surface_no_mean_(f,:) + mean_f(f,:);
+  end
+
+  model_ = resample_rate_L(model, rate_K_surface_, rate_K_sample_freqs_kHz);
+
+endfunction
+
+
+% try vq with open and closed loop mean removal, turns out they give
+% identical results lol
+
+function [model_ rate_K_surface] = experiment_closed_loop_mean(model)
+  [frames nc] = size(model);
+  K = 15; 
+  mel_start = ftomel(300); mel_end = ftomel(3000); 
+  step = (mel_end-mel_start)/(K-1);
+  mel = mel_start:step:mel_end;
+  rate_K_sample_freqs_Hz = 700*((10 .^ (mel/2595)) - 1);
+  rate_K_sample_freqs_kHz = rate_K_sample_freqs_Hz/1000;
+
+  rate_K_surface = resample_const_rate_f(model, rate_K_sample_freqs_kHz);
+  
+  load surface_vq; m=1;
+   
+  for f=1:frames
+    amean = mean(rate_K_surface(f,:));
+    rate_K_target_no_mean = rate_K_surface(f,:) - amean;
+    mse_open_loop(f) = search_vq2(surface_vq(:,:,1), rate_K_target_no_mean, m, 0);
+    mse_closed_loop(f) = search_vq2(surface_vq(:,:,1), rate_K_surface(f,:), m, 1);
+  end
+
+  printf("rms open loop..: %f\nrms closed loop: %f\n", sqrt(mean(mse_open_loop)), sqrt(mean(mse_closed_loop)));
+
+  % just return model_ as we have to so nothing breaks, it's not actually useful
+
+  model_ = resample_rate_L(model, rate_K_surface, rate_K_sample_freqs_kHz);
+    
+endfunction
+
+
+% Experiment with 10ms update rate for energy but 40ms update for spectrum,
+% using linear interpolation for spectrum.
+
+function model_c = experiment_energy_rate_linear(model, vq_en, plot_en)
+  max_amp = 80;
+  [frames nc] = size(model);
+
+  % 10ms mel freq modelling and VQ
+
+  model_ = experiment_mel_freq(model, vq_en, plot_en);
+
+  % Remove energy.  Hmmmm, this is done on Ams rather than surface but that's
+  % similar I guess
+
+  e = zeros(1,frames);
+  model_a = zeros(frames,max_amp+3);
+  for f=1:frames
+    L = min([model_(f,2) max_amp-1]);
+    Am_ = model_(f,3:(L+2));
+    AmdB_ = 20*log10(Am_);
+    mean_f(f) = mean(AmdB_);
+    AmdB_ -= mean_f(f);
+    model_a(f,1) = model_(f,1); model_a(f,2) = L; model_a(f,3:(L+2)) = 10 .^ (AmdB_(1:L)/20);
+  end
+
+  % linear interp after removing energy (mean AmdB)
+
+  model_b = experiment_dec_linear(model_a);
+
+  % add back in energy
+
+  model_c = zeros(frames,max_amp+3);
+  for f=1:frames
+    L = min([model_b(f,2) max_amp-1]);
+    Am_ = model_b(f,3:(L+2));
+    AmdB_ = 20*log10(Am_);
+    AmdB_ += mean_f(f);
+    model_c(f,1) = model_b(f,1); model_c(f,2) = L; model_c(f,3:(L+2)) = 10 .^ (AmdB_(1:L)/20);
+  end
+
+endfunction
+
+
+% conventional decimation in time without any filtering, then linear
+% interpolation.  Linear interpolation is a two-tap (weak) form of fir
+% filtering that may have problems with signals with high freq
+% components, for example after quantisation noise is added.  Need to
+% look into this some more.
+function model_ = experiment_dec_linear(model)
+  newamp;
+  max_amp = 80;
+
+  [frames nc] = size(model);
+  model_ = zeros(frames, max_amp+3);
+  decimate = 4;
+  for f=1:frames    
+    AmdB_ = decimate_frame_rate(model, decimate, f, frames);
+    L = length(AmdB_);
+    model_(f,1) = model(f,1); model_(f,2) = L; model_(f,3:(L+2)) = 10 .^ (AmdB_(1:L)/20);
+  end
+endfunction
+
+
 
 % Experimental AbyS decimator that chooses best frames to match
 % surface based on AbyS approach.  Can apply post filter at different
@@ -828,7 +834,14 @@ endfunction
 
 
 #{
-todo: get this working again
+  
+  My original idea was to used a 3-4 "resonators" to construct a
+  piecewise model of the spectrum.  Kind of got distracted by the
+  surface and mel sampling that ended up working OK.  This method was
+  working OK, soem issues with backgorund noise but rather easy to
+  quantise.
+
+  todo: get this working again
 
 function model_ = experiment_piecewise(model)
   % encoder loop ------------------------------------------------------
@@ -875,9 +888,6 @@ endfunction
 #}
 
 
-
-
-
 % early test, devised to test rate K<->L changes along frequency axis
 
 function model_ = resample_half_frame_offset(model, rate_K_surface, rate_K_sample_freqs_kHz)
@@ -941,6 +951,9 @@ function [mse_list index_list] = search_vq2(vq, target, m, closed_loop_dc = 0)
 
 endfunction
 
+
+% helper function to plot surfaces
+
 function plot_dft_surface(surface)
   [frames K] = size(surface);