progressing quantisation
authordrowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Mon, 29 Feb 2016 23:39:28 +0000 (23:39 +0000)
committerdrowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Mon, 29 Feb 2016 23:39:28 +0000 (23:39 +0000)
git-svn-id: https://svn.code.sf.net/p/freetel/code@2710 01035d8c-6547-0410-b346-abe4f91aad63

codec2-dev/octave/newamp.m
codec2-dev/octave/newamp_batch.m
codec2-dev/octave/newamp_fbf.m

index 6a22b86b02acb165e959a49ee5d1b71fe6e304e9..2fb053cda99b274a5f14fdc421b72c2f825b5e39 100644 (file)
@@ -135,7 +135,7 @@ endfunction
 % Ahh, takes me back to when I was a slip of a speech coder, playing
 % with my first CELP codec!
 
-function [decmaskdB dec_samples min_error mse_log1 mse_log2] = make_decmask_abys(maskdB, AmdB, Wo, L, mask_sample_freqs_kHz)
+function [decmaskdB masker_freqs_kHz min_error mse_log1 mse_log2] = make_decmask_abys(maskdB, AmdB, Wo, L, mask_sample_freqs_kHz, freq_quant, amp_quant)
 
     Nsamples = 4;
 
@@ -189,32 +189,40 @@ function [decmaskdB dec_samples min_error mse_log1 mse_log2] = make_decmask_abys
 
       decmaskdB = min_candidate;
       dec_samples = [dec_samples; AmdB(min_mse_m) min_mse_m];
+      masker_freqs_kHz = dec_samples(:,2)*Wo*4/pi;
       %printf("sample: %d min_mse_m: %d\n", sample, min_mse_m);
     end
 
-    % Differential Freq Quantisers - sounds acceptable
+    bits = [];
 
-    if 0
-      % sort into increasing freq order
+    % sort into increasing freq order
 
-      masker_amps_dB = dec_samples(:,1);
-      masker_freqs_kHz = dec_samples(:,2)*Wo*4/pi;
-      [fsrt fsrt_ind] = sort(masker_freqs_kHz);
-            
+    masker_amps_dB = dec_samples(:,1);
+    masker_freqs_kHz = dec_samples(:,2)*Wo*4/pi;
+    [fsrt fsrt_ind] = sort(masker_freqs_kHz);
+    masker_freqs_kHz = fsrt;
+    masker_amps_dB = masker_amps_dB(fsrt_ind);
+
+    % Differential Freq Quantisers - sounds acceptable
+
+    if freq_quant
+           
       % first freq quant to harmonic number m=1:8
 
       f0_kHz = Wo*4/pi;
-      masker_freqs_kHz(1) = quantise((1:8)*f0_kHz, fsrt(1));
-      
+      [masker_freqs_kHz(1) abits] = quantise((1:8)*f0_kHz, masker_freqs_kHz(1));
+      bits = [bits abits];
+     
       % then quantise differences
 
       for i=2:4
-        targ = fsrt(i) - masker_freqs_kHz(i-1);
-        masker_freqs_kHz(i) = masker_freqs_kHz(i-1) + quantise(0.2:0.2:1.6, targ);
+        targ = masker_freqs_kHz(i) - masker_freqs_kHz(i-1);
+        [q_freq abits] = quantise(0.2:0.2:1.6, targ);
+        bits = [bits abits];
+        masker_freqs_kHz(i) = masker_freqs_kHz(i-1) + q_freq;
       end
 
-      asrt = masker_amps_dB(fsrt_ind);
-      decmaskdB = determine_mask(asrt,  masker_freqs_kHz, mask_sample_freqs_kHz);
+       decmaskdB = determine_mask(masker_amps_dB,  masker_freqs_kHz, mask_sample_freqs_kHz);
     end
 
     if 0
@@ -279,34 +287,83 @@ function [decmaskdB dec_samples min_error mse_log1 mse_log2] = make_decmask_abys
       decmaskdB = determine_mask(masker_amps_dB,  masker_freqs_kHz, mask_sample_freqs_kHz);
     end
 
-    % fit straight line to amplitudes (sounds bad, still a bug somewhere)
+    % Amplitude quantisation by fitting a straight line -------------------------
+
+    if amp_quant
+
+      % Fit straight line
+
+      f = masker_freqs_kHz*1000;
+      [gradient intercept] = linreg(f, masker_amps_dB, 4);
+      % use quantised gradient to take into account quantisation
+      % errors in rest of quantisation
 
-    if 1
-      f = dec_samples(:,2)*Wo*4000/pi;
-      [gradient intercept] = linreg(f, dec_samples(:,1), 4);
-      masker_amps_dB_lin = f*gradient + intercept;
-      masker_amps_dB_lin_delta = dec_samples(:,1) - masker_amps_dB_lin;
       gradient_ = quantise(-0.04:0.005:0.04, gradient);
-      printf("gradient: %f gradient_: %f\n", gradient, gradient_);
-      %masker_amps_dB = f*gradient_ + masker_amps_dB_lin_delta + intercept;
-      masker_amps_dB = f*gradient_ + intercept;
-      masker_freqs_kHz = f/1000;
+      
+      % determine deltas, or errors in straight line fit
+
+      masker_amps_dB_lin = f*gradient_ + intercept;
+      masker_amps_dB_lin_delta = masker_amps_dB - masker_amps_dB_lin;
+
+      % quantise the deltas
+
+      masker_amps_dB_lin_delta_ = zeros(4,1);
+      for i=1:4
+        masker_amps_dB_lin_delta_(i) = quantise(-15:5:20, masker_amps_dB_lin_delta(i));
+      end
+      %masker_amps_dB_lin_delta
+      %masker_amps_dB_lin_delta_
+
+      masker_amps_dB = f*gradient_ + masker_amps_dB_lin_delta_ + intercept;
+      %decmaskdB = determine_mask(masker_amps_dB,  masker_freqs_kHz, mask_sample_freqs_kHz);
       decmaskdB = determine_mask(masker_amps_dB,  masker_freqs_kHz, mask_sample_freqs_kHz);
     end
 endfunction
 
 
-% quantise input sample to nearest value in table 
+% quantise input sample to nearest value in table, optionally return bianry code
+
+function [quant_out bits] = quantise(levels, quant_in)
+
+  % find closest quantiser level
 
-function quant_out = quantise(levels, quant_in)
   best_se = 1E32;
   for i=1:length(levels)
     se = (levels(i) - quant_in)^2;
     if se < best_se
       quant_out = levels(i);
       best_se = se;
+      best_i = i;
     end
   end
+
+  % convert index to binary bits
+
+  numbits = ceil(log2(length(levels)));
+  bits = zeros(1, numbits);
+  for b=1:numbits
+    bits(b) = bitand(best_i-1,2^(numbits-b)) != 0;
+  end
+
+endfunction
+
+
+function masker_freqs_kHz = unquantise_freqs(bits, Wo)
+  for i=1:4
+    st = (i-1)*3+1; en=i*3; 
+    index(i) = bits(st:en) * [4 2 1]' + 1;
+  end
+  f0_kHz = Wo*4/pi;
+
+  masker_freqs_kHz(1) = index(1)*f0_kHz;
+     
+  % then unquantise differences
+
+  q_freqs = 0.2:0.2:1.6;
+
+  for i=2:4
+    masker_freqs_kHz(i) = masker_freqs_kHz(i-1) + q_freqs(index(i));
+  end
 endfunction
 
 
index aab1986598d3fd567f05b40897d62c104ebf6181..d77a3ca3efc312f4651d7dc3eef2e4690da9448b 100644 (file)
 
 % process a whole file and write results
 
-function newamp_batch(samname, optional_Am_out_name, optional_Aw_out_name)
+function [non_masked_f_log non_masked_amp_log] = newamp_batch(samname, optional_Am_out_name, optional_Aw_out_name)
   newamp;
   more off;
 
   max_amp = 80;
   decimate_in_freq = 1;
   postfilter = 1;
-  decimate_in_time = 0;
+  decimate_in_time = 1;
   synth_phase = 0;
+  freq_quant = 1;
+  amp_quant = 1;
+  non_masked_f_log = [];
+  non_masked_m_log = [];
+  non_masked_amp_log = [];
 
   model_name = strcat(samname,"_model.txt");
   model = load(model_name);
@@ -39,10 +44,10 @@ function newamp_batch(samname, optional_Am_out_name, optional_Aw_out_name)
   end
   if nargin >= 3
     Aw_out_name = optional_Aw_out_name;
+    faw = fopen(Aw_out_name,"wb"); 
   end
    
   fam  = fopen(Am_out_name,"wb"); 
-  faw = fopen(Aw_out_name,"wb"); 
 
   for f=1:frames
     printf("%d ", f);
@@ -58,9 +63,14 @@ function newamp_batch(samname, optional_Am_out_name, optional_Aw_out_name)
     maskdB = mask_model(AmdB, Wo, L);
 
     if decimate_in_freq
-      % decimate in mask samples in frequency
-      [decmaskdB decmasksamples] = make_decmask_abys(maskdB, AmdB, Wo, L, mask_sample_freqs_kHz);
-      non_masked_m = decmasksamples(:,2);
+      % decimate mask samples in frequency
+      [decmaskdB masker_freqs_kHz] = make_decmask_abys(maskdB, AmdB, Wo, L, mask_sample_freqs_kHz, freq_quant, amp_quant);
+      non_masked_amp = decmaskdB;
+      non_masked_amp_log = [non_masked_amp_log; non_masked_amp'];
+      Wo*4/pi
+      non_masked_m = round(masker_freqs_kHz/(Wo*4/pi))
+      non_masked_m_log = [non_masked_m_log; non_masked_m'];
+      non_masked_f_log = [non_masked_f_log; masker_freqs_kHz];
       maskdB_ = decmaskdB;
     else
       % just approximate decimation in freq by using those mask samples we can hear
@@ -70,7 +80,7 @@ function newamp_batch(samname, optional_Am_out_name, optional_Aw_out_name)
 
     if decimate_in_time
       % decimate mask samples in time
-      maskdB_ = decimate_frame_rate(maskdB_, model, decimate, f, frames, mask_sample_freqs_kHz);
+      maskdB_ = decimate_frame_rate(maskdB_, model, 4, f, frames, mask_sample_freqs_kHz);
     end
 
     % post filter - bump up samples by 6dB, reduce mask by same level to normalise gain
@@ -112,10 +122,11 @@ function newamp_batch(samname, optional_Am_out_name, optional_Aw_out_name)
   end
 
   fclose(fam);
-  fclose(faw);
+  if synth_phase
+    fclose(faw);
+  end
 
   printf("\n");
-
 endfunction
 
 
index 1f6e69390bf234873f9df42686f51524c2950066..6635de44039f1892cf53d41dda71b9d852116d6d 100644 (file)
 %   $ cd ~/codec2-dev/octave
 %   octave:14> newamp_fbf("../build_linux/src/hts1a",50)
 
-function newamp_fbf(samname, f)
+function newamp_fbf(samname, f=10)
   
   more off;
   newamp;
   phase_stuff = 0;
   plot_not_masked = 0;
-  plot_spectrum = 0;
+  plot_spectrum = 1;
+  freq_quant = 0;
+  amp_quant = 0;
 
   sn_name = strcat(samname,"_sn.txt");
   Sn = load(sn_name);
@@ -62,19 +64,12 @@ function newamp_fbf(samname, f)
     hold on;
     if plot_spectrum
       plot((1:L)*Wo*4000/pi, AmdB,";Am;r");
-      plot((1:L)*Wo*4000/pi, AmdB,";Am;r+");
-      plot((0:255)*4000/256, Sw(f,:),";Sw;");
+      plot((1:L)*Wo*4000/pi, AmdB,"r+");
+      %plot((0:255)*4000/256, Sw(f,:),";Sw;");
     end
 
     [maskdB Am_freqs_kHz] = mask_model(AmdB, Wo, L);
-    plot(Am_freqs_kHz*1000, maskdB, 'g');
-
-    % optionally show harmonics that are not masked
-
-    not_masked_m = find(maskdB < AmdB);
-    if plot_not_masked
-      plot(not_masked_m*Wo*4000/pi, 70*ones(1,length(not_masked_m)), 'bk+');
-    end
+    %plot(Am_freqs_kHz*1000, maskdB, 'g');
 
     % optionally plot synthesised spectrum (early simple model)
 
@@ -88,28 +83,21 @@ function newamp_fbf(samname, f)
     % decimate in frequency
 
     mask_sample_freqs_kHz = (1:L)*Wo*4/pi;
-    [decmaskdB local_maxima min_error mse_log1 mse_log2] = make_decmask_abys(maskdB, AmdB, Wo, L, mask_sample_freqs_kHz);
+    [decmaskdB masker_freqs_kHz min_error mse_log1 mse_log2] = make_decmask_abys(maskdB, AmdB, Wo, L, mask_sample_freqs_kHz, freq_quant, amp_quant);
     
-    [nlm tmp] = size(local_maxima(:,2));
-    nlm = min(nlm,4);
-    tonef_kHz = local_maxima(1:nlm,2)*Wo*4/pi;
-    toneamp_dB = local_maxima(1:nlm,1);
+    tonef_kHz = masker_freqs_kHz;
+    nlm = length(tonef_kHz);
 
-    plot(tonef_kHz*1000, 70*ones(1,nlm), 'bk+');
-    plot(mask_sample_freqs_kHz*1000, decmaskdB, 'm');
+    plot(tonef_kHz*1000, 20*ones(1,nlm), ';AbyS Mask Freqs;bk+');
+    plot(mask_sample_freqs_kHz*1000, decmaskdB, ';AbyS Mask;m');
     plot(mask_sample_freqs_kHz*1000, min_error);
+    hold off;
 
     figure(3)
     clf
     plot((1:L)*Wo*4000/pi, mse_log1');
     axis([0 4000 0 max(mse_log1(1,:))])
-    title('Basis 1 MSE as a function of position for each stage');
-
-    % fit a line to amplitudes
-
-    %[m b] = linreg(tonef_kHz, toneamp_dB, nlm);
-    %plot(tonef_kHz*1000, tonef_kHz*m + b, "bk");
-    %plot(tonef_kHz*1000, 60 + toneamp_dB - (tonef_kHz*m + b), "r+");
+    title('Error as a function of position for each AbyS stage');
 
     % decimated in time
 
@@ -147,7 +135,7 @@ function newamp_fbf(samname, f)
 
     % interactive menu
 
-    printf("\rframe: %d  menu: n-next  b-back m-allmasks o-notmasked s-spectrum q-quit", f);
+    printf("\rframe: %d  menu: n-next  b-back a-Am f-freq_quant m-amp_quant q-quit", f);
     fflush(stdout);
     k = kbhit();
     if (k == 'n')
@@ -157,20 +145,20 @@ function newamp_fbf(samname, f)
       f = f - 1;
     endif
     if k == 'm'
-      if plot_all_masks == 0
-         plot_all_masks = 1;
+      if amp_quant == 0
+         amp_quant = 1;
       else
-         plot_all_masks = 0;
+         amp_quant = 0;
       end
     end
-    if k == 'o'
-      if plot_not_masked == 0
-         plot_not_masked = 1;
+    if k == 'f'
+      if freq_quant == 0
+        freq_quant = 1;
       else
-         plot_not_masked = 0;
+        freq_quant = 0;
       end
     end
-    if k == 's'
+    if k == 'a'
       if plot_spectrum == 0
          plot_spectrum = 1;
       else