pager script to explore VQ entries

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

diff --git a/codec2-dev/octave/newamp1_batch.m b/codec2-dev/octave/newamp1_batch.m
index 0010c3776fd4d0868dfdc63303f6081bd94fced7..16a23c18d4a5bf2dae046c02018265bafeed61c2 100644 (file)
--- a/codec2-dev/octave/newamp1_batch.m
+++ b/codec2-dev/octave/newamp1_batch.m
@@ -40,16 +40,18 @@
 % In general, this function processes a bunch of amplitudes, we then
 % use c2sim to hear the results.  Bunch of different experiments below
 
-function surface = newamp1_batch(input_prefix, varargin)
+function [surface mean_f] = newamp1_batch(input_prefix, varargin)
   newamp;
   more off;
 
   max_amp = 160;
+  mean_f = [];
 
   % defaults
 
-  synth_phase = 1;
+  synth_phase = output = 1;
   output_prefix = input_prefix;
+  vq_type = "";
   vq_filename = "";
   mean_removal = 0;
   mode = "const";
@@ -72,11 +74,21 @@ function surface = newamp1_batch(input_prefix, varargin)
     if ind
       mean_removal = 1;
     end
-    ind = arg_exists(varargin, "vq");
+    ind = arg_exists(varargin, "vqh");
     if ind
+      vq_type = varargin{ind};
       vq_filename = varargin{ind+1};
+      printf("vq_type: %s vq_filename: %s\n", vq_type, vq_filename);
+    end
+    ind = arg_exists(varargin, "no_output");
+    if ind
+      output = 0;
+      synth_phase = 0;
     end
   end
+
+  printf("output_prefix: %s\nnomean: %d output: %d vq_type: %s\nvq_filename: %s\n", 
+         output_prefix, mean_removal, output, vq_type, vq_filename);
   
   model_name = strcat(input_prefix,"_model.txt");
   model = load(model_name);
@@ -99,7 +111,7 @@ function surface = 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, mean_removal, vq_filename);
+    [model_ surface mean_f] = experiment_const_freq(model, mean_removal, vq_type, vq_filename);
   end
   if strcmp(mode, 'piecewise')
     model_ = experiment_piecewise(model);
@@ -110,61 +122,69 @@ function surface = newamp1_batch(input_prefix, varargin)
 
   % ----------------------------------------------------
 
-  Am_out_name = sprintf("%s_am.out", output_prefix);
-  fam  = fopen(Am_out_name,"wb"); 
-
-  Wo_out_name = sprintf("%s_Wo.out", output_prefix);
-  fWo  = fopen(Wo_out_name,"wb"); 
+  if output
+    Am_out_name = sprintf("%s_am.out", output_prefix);
+    fam  = fopen(Am_out_name,"wb"); 
+ 
+    Wo_out_name = sprintf("%s_Wo.out", output_prefix);
+    fWo  = fopen(Wo_out_name,"wb"); 
   
-  if synth_phase
-    Hm_out_name = sprintf("%s_hm.out", output_prefix);
-    fhm = fopen(Hm_out_name,"wb"); 
-  end
+    if synth_phase
+     Hm_out_name = sprintf("%s_hm.out", output_prefix);
+     fhm = fopen(Hm_out_name,"wb"); 
+    end
 
-  for f=1:frames
-    %printf("%d ", f);   
-    Wo = model_(f,1); L = min([model_(f,2) max_amp-1]); Am = model_(f,3:(L+2));
-    assert(Wo*L < pi);
+    for f=1:frames
+      %printf("%d ", f);   
+      Wo = model_(f,1); L = min([model_(f,2) max_amp-1]); Am = model_(f,3:(L+2));
+      if Wo*L > pi
+        printf("Problem: %d  Wo*L > pi\n", f);   
+      end
 
-    Am_ = zeros(1,max_amp); Am_(2:L) = Am(1:L-1); fwrite(fam, Am_, "float32");
-    fwrite(fWo, Wo, "float32");
+      Am_ = zeros(1,max_amp); Am_(2:L) = Am(1:L-1); fwrite(fam, Am_, "float32");
+      fwrite(fWo, Wo, "float32");
 
-    if synth_phase
+      if synth_phase
 
-      % synthesis phase spectra from magnitiude spectra using minimum phase techniques
+        % synthesis phase spectra from magnitiude spectra using minimum phase techniques
 
-      fft_enc = 256;
-      phase = determine_phase(model_, f, fft_enc);
-      assert(length(phase) == fft_enc);
+        fft_enc = 256;
+        phase = determine_phase(model_, f, fft_enc);
+        assert(length(phase) == fft_enc);
 
-      % sample phase at centre of each harmonic, not 1st entry Hm[1] in octave Hm[0] in C
-      % is not used
+        % sample phase at centre of each harmonic, not 1st entry Hm[1] in octave Hm[0] in C
+        % is not used
 
-      Hm = zeros(1, 2*max_amp);
-      for m=1:L
-        b = round(m*Wo*fft_enc/(2*pi));
-        Hm(2*m) = cos(phase(b));
-        Hm(2*m+1) = -sin(phase(b));
+        Hm = zeros(1, 2*max_amp);
+        for m=1:L
+          b = round(m*Wo*fft_enc/(2*pi));
+          Hm(2*m) = cos(phase(b));
+          Hm(2*m+1) = -sin(phase(b));
+        end
+        fwrite(fhm, Hm, "float32");    
       end
-      fwrite(fhm, Hm, "float32");    
+      fclose(fhm);
     end
-  end
  
-  fclose(fam);
-  fclose(fWo);
-  if synth_phase
-    fclose(fhm);
-  end
+   fclose(fam);
+   fclose(fWo);
 
-  % save voicing file
+   % save voicing file
   
-  if exist("voicing_", "var")
-    v_out_name = sprintf("%s_v.txt", output_prefix);
-    fv  = fopen(v_out_name,"wt"); 
-    for f=1:length(voicing_)
-      fprintf(fv,"%d\n", voicing_(f));
+   if exist("voicing_", "var")
+     v_out_name = sprintf("%s_v.txt", output_prefix);
+     fv  = fopen(v_out_name,"wt"); 
+     for f=1:length(voicing_)
+       fprintf(fv,"%d\n", voicing_(f));
+     end
+     fclose(fv);
+    end
+  end
+
+  if mean_removal
+    for f=1:frames
+      surface(f,:) -= mean_f(f);
     end
-    fclose(fv);
   end
 
   printf("\n")
@@ -181,13 +201,32 @@ function ind = arg_exists(v, str)
     end
 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, mean_removal, vq_filename)
+function [model_ rate_K_surface mean_f] = experiment_const_freq(model, mean_removal, vq_type, vq_filename)
+  melvq;
+
   [frames nc] = size(model);
   Fs = 8000;
   fg = 1;
+  quant_en = 0;
+
+  rate_K_sample_freqs_kHz = [0.1:0.1:4];
+  K = length(rate_K_sample_freqs_kHz);
+
+  if strcmp(vq_type, "vql")
+    quant_en = 1;
+    x = load(vq_filename); vq = x.vq;
+    [vq_rows vq_cols] = size(vq); vq_st = 1; vq_en = vq_st + vq_cols - 1;
+  end
+
+  if strcmp(vq_type, "vqh")
+    quant_en = 1;
+    x = load(vq_filename); vq = x.vq;
+    [vq_rows vq_cols] = size(vq); vq_st = 11; vq_en = K;
+  end
 
   energy = zeros(1,frames);
   for f=1:frames
@@ -195,38 +234,55 @@ function [model_ rate_K_surface] = experiment_const_freq(model, mean_removal, vq
     energy(f) = 10*log10(sum( model(f,3:(L+2)) .^ 2 ));
   end
 
-  melvq;
-  if length(vq_filename) 
-    x = load(vq_filename); m=5;
-    vq = x.vq;
-    size(vq)
-  end
-
-  rate_K_sample_freqs_kHz = [0.1:0.1:4];
-  K = length(rate_K_sample_freqs_kHz);
   rate_K_surface = resample_const_rate_f(model, rate_K_sample_freqs_kHz, Fs);
 
   mean_f = zeros(1,frames);
   if mean_removal
     for f=1:frames
       mean_f(f) = mean(rate_K_surface(f,:));
-      rate_K_surface(f,:) -= mean_f(f);
+      rate_K_surface_no_mean(f,:) = rate_K_surface(f,:) - mean_f(f);
     end
   end
 
-  if length(vq_filename) 
-    [res rate_K_surface_ ind] = mbest(vq, rate_K_surface, m);
-    sd_per_frame = std(res');
+  % optional vector quantise
+  
+  if quant_en
+    assert(mean_removal == 1);
+ 
+    weight_gain = 0.1; % I like this vairable name as it is funny
+
+    rate_K_surface_no_mean_ = zeros(frames, K);
+    res = zeros(frames,vq_cols); ind = [];
+
+    for f=1:frames
+      target = rate_K_surface_no_mean(f, vq_st:vq_en);
+
+      [diff_weighted weights error g mn_ind] = search_vq_weighted(target, vq, weight_gain);
+
+      rate_K_surface_no_mean_(f,:) = rate_K_surface_no_mean(f,:);
+      rate_K_surface_no_mean_(f, vq_st:vq_en) = vq(mn_ind,:) + g(mn_ind);
+
+      %res(f,vq_st:vq_en) = rate_K_surface_no_mean(f,vq_st:vq_en) - rate_K_surface_no_mean_(f,vq_st:vq_en);
+      res(f,vq_st:vq_en) = diff_weighted(mn_ind,:);
+      ind = [ind mn_ind];
+    end
+    
+    figure(fg++); clf; mesh(res);
+    sd_per_frame = std(res(:,vq_st:vq_en)');
     figure(fg++); subplot(211); plot(energy); subplot(212); plot(sd_per_frame);
-    figure(fg); hist(sd_per_frame);
+    figure(fg++); subplot(211); hist(sd_per_frame); subplot(212); hist(ind,100);
     printf("VQ rms SD: %3.2f\n", mean(sd_per_frame));
   else
-    rate_K_surface_ = rate_K_surface;
+    if mean_removal
+       rate_K_surface_no_mean_ = rate_K_surface_no_mean;
+    else
+      rate_K_surface_ = rate_K_surface;
+    end
   end
 
   if mean_removal
     for f=1:frames
-      rate_K_surface_(f,:) += mean_f(f);
+      rate_K_surface_(f,:) = rate_K_surface_no_mean_(f,:) + mean_f(f);
     end
   end
 
@@ -243,7 +299,7 @@ function [model_ rate_K_surface] = experiment_const_freq(model, mean_removal, vq
     L = model(f,2);
     AmdB = 20*log10(model(f,3:(L+2)));
     sd(f) = std(AmdB(1:L) - AmdB_(f,1:L));
-    if sd(f) > plot_sd_thresh
+    if (sd(f) > plot_sd_thresh) && (fg < 10)
       printf("fg: %d f: %d\n", fg, f);
       figure(fg++); clf; plot((1:L)*Wo*4/pi, AmdB(1:L),'b+-'); hold on; plot((1:L)*Wo*4/pi, AmdB_(f,1:L),'r+-');
       plot(rate_K_sample_freqs_kHz, rate_K_surface_(f,:), 'c+-'); hold off;

diff --git a/codec2-dev/octave/newamp1_fbf.m b/codec2-dev/octave/newamp1_fbf.m
index d9d6a12822d8bcad501bd7db7ebd6bb8f6be3994..65917c5be7eb53d5ed26421648b3082ce3cbcdc1 100644 (file)
--- a/codec2-dev/octave/newamp1_fbf.m
+++ b/codec2-dev/octave/newamp1_fbf.m
@@ -83,63 +83,26 @@ function newamp1_fbf(samname, f=10, varargin)
     axis([1 4000 -20 80]);
     hold on;
     plot(rate_K_sample_freqs_kHz*1000, rate_K_vec, ";rate K;b+-");
-    %if quant_en != 2
-    %  plot(rate_K_sample_freqs_kHz*1000, maskdB, ";mask;c+-");
-    %end
-    if quant_en == 2
-      plot(rate_K_sample_freqs_kHz*1000, maskdB-rate_K_vec, ";maskdB - rate K;m+-");
-    end
-    #{
-    masked_indx = find(maskdB >= rate_K_vec);
-    unmasked_indx = find(maskdB < rate_K_vec);
-    masked_indx = (vq_en-5:vq_en);
-    unmasked_indx = (vq_st:vq_en-6);
-    stem(rate_K_sample_freqs_kHz(unmasked_indx)*1000, 10*ones(1,length(unmasked_indx)), ";unmasked_index;ko-");
-    #}
 
     if quant_en
-      error = g = zeros(1, vq_rows);
-
-      % find mse for each vector
-
       target = rate_K_vec_no_mean(vq_st:vq_en);
-      weight_gain = 0.1;
-      %masked_indx = find(maskdB(vq_en-5:vq_en) >= rate_K_vec(vq_en-5:vq_en));
-      %masked_indx = find(masked_indx > (vq_cols-5))
-      diff =  weights = diff_weighted = zeros(vq_rows, vq_cols);
-      for i=1:vq_rows
-
-        % work out gain for best match
+      weight_gain = 0.1; % I like this vairable name as it is funny
 
-        g(i) = sum(target - vq(i,:))/vq_cols;
-
-        % Find weighted difference.  This allocated more importance
-        % (error) to samples with higher energy, and stops really low
-        % level harmonics from having any impact.  Note addition in dB
-        % is multiplication in linear
-
-        diff(i,:) = target - vq(i,:) - g(i);
-        weights(i,:) = max(0.1,weight_gain .* (target + 20));
-        
-        diff_weighted(i,:) = diff(i,:) .* weights(i,:);
-
-        error(i) = sum(abs(diff_weighted(i,:)));
-      end
-
-      [mn mn_ind] = min(error);
+      [diff_weighted weights error g mn_ind] = search_vq_weighted(target, vq, weight_gain);
 
       rate_K_vec_no_mean_ = rate_K_vec_no_mean;
       rate_K_vec_no_mean_(vq_st:vq_en) = vq(mn_ind,:) + g(mn_ind);
       rate_K_vec_ = rate_K_vec_no_mean_ + mean(rate_K_vec);
       [model_ AmdB_] = resample_rate_L(model(f,:), rate_K_vec_, rate_K_sample_freqs_kHz, Fs);
       AmdB_ = AmdB_(1:L);
-      if quant_en == 2
-        plot(rate_K_sample_freqs_kHz(vq_st:vq_en)*1000, diff_weighted(mn_ind,:), ";diff;k+-");
-      end
+
+      plot(rate_K_sample_freqs_kHz(vq_st:vq_en)*1000, diff_weighted(mn_ind,:), ";diff;k+-");
       plot((1:L)*Wo*4000/pi, AmdB_,";AmdB bar;r+-");
       hold off;
 
-      figure (4); clf; plot(rate_K_sample_freqs_kHz(vq_st:vq_en)*1000, weights(mn_ind,:), ";weights;k+-");
+      figure (4); clf; 
+      plot(rate_K_sample_freqs_kHz(vq_st:vq_en)*1000, weights(mn_ind,:), ";weights;k+-");  
+      axis([0 4000 0 max(weights(mn_ind,:))]);
 
       % sort and plot top m matches
 
@@ -152,15 +115,20 @@ function newamp1_fbf(samname, f=10, varargin)
       for i=1:m
         subplot(sqrt(m),sqrt(m),i);
         indx = index_list(i);
-        plot(target + mean(rate_K_vec),'b+-');
+        y_offset = mean(rate_K_vec);
+        if quant_en == 2
+           y_offset = 0;
+        end
+        plot(target + y_offset,'b+-');
         hold on;
         if index_list(i) == mn_ind
-          plot(vq(indx,:) + g(indx) + mean(rate_K_vec),'r+-');
+          plot(vq(indx,:) + g(indx) + y_offset,'r+-');
         else
-          plot(vq(indx,:) + g(indx) + mean(rate_K_vec),'g+-');
+          plot(vq(indx,:) + g(indx) + y_offset,'g+-');
+        end
+        if quant_en != 2
+          plot(diff_weighted(indx,:), "ko-");
         end
-        plot(diff(indx,:), ";diff;k+-");
-        plot(diff_weighted(indx,:), ";weighted diff;ko-");
         hold off;
         legend("location", "southwest");
       end
@@ -198,3 +166,5 @@ function ind = arg_exists(v, str)
       end     
     end
 endfunction
+
+

diff --git a/codec2-dev/octave/vq_pager.m b/codec2-dev/octave/vq_pager.m
new file mode 100644 (file)
index 0000000..f00de9e
--- /dev/null
+++ b/codec2-dev/octave/vq_pager.m
@@ -0,0 +1,47 @@
+% vq_pager.m
+%
+% Interactive Octave script to inspect vectors, e.g. for newamp1.m vector quantisation
+%
+% Usage:
+%
+%   octave:14> vq_pager(vq)
+
+function vq_pager(vq)
+  grid_sz = 4;
+
+  [vq_rows vq_cols] = size(vq);
+  
+  r = 1;
+
+  % Keyboard loop 
+
+  k = ' ';
+  do 
+    figure(1); clf;
+    n_plot = min(grid_sz*grid_sz,vq_rows-r+1)
+    printf("r: %d n_plot: %d\n", r, n_plot);
+    for i = 1:n_plot;
+      subplot(grid_sz,grid_sz,i);
+      plot(vq(r+i-1,:));
+      axis([1 vq_cols -20 20]);
+    end
+
+    % interactive menu 
+
+    printf("\rr: %d  menu: n-next  b-back  q-quit", r);
+    fflush(stdout);
+    k = kbhit();
+
+    if k == 'n'
+      r = r + grid_sz*grid_sz;
+      r = min(r, vq_rows);
+    endif
+    if k == 'b'
+      r = r - grid_sz*grid_sz;
+      r = max(r,1);
+    endif
+  until (k == 'q')
+  printf("\n");
+
+endfunction
+