listening tests for new mask/circ/reampling work, sound spretty gd, working on VQ
authordrowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Thu, 31 Mar 2016 00:58:02 +0000 (00:58 +0000)
committerdrowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Thu, 31 Mar 2016 00:58:02 +0000 (00:58 +0000)
git-svn-id: https://svn.code.sf.net/p/freetel/code@2765 01035d8c-6547-0410-b346-abe4f91aad63

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

index 8b61229bdc4c734ef2039cca22df7838dd3df37c..24da35aee540957437765123a6473d91dfc3ecf7 100644 (file)
 1;
 melvq;
 
+
+function [maskdB_ maskdB_cyclic Dabs dk_] = decimate_in_freq(maskdB, cyclic=1, Nkeep=7)
+
+    % Lets try to come up with a smoothed, cyclic model.  Replace
+    % points from 3500 Hz to 4000Hz with a sequence that joins without
+    % a step to points at the 0Hz end of the spectrum.  This will make
+    % it more cyclical and make the DFT happier, less high freq
+    % energy.  Yes, happier is an extremely technical term.
+
+    L = length(maskdB);
+    anchor = floor(7*L/8);
+    xpts = [ anchor-1 anchor L+1 L+2];
+    ypts = [ maskdB(anchor-1) maskdB(anchor) maskdB(1) maskdB(2)];
+    mask_pp = splinefit(xpts, ypts, 1);
+    maskdB_cyclic = [maskdB(1:anchor) ppval(mask_pp, anchor+1:L)];
+
+    % Now DFT, truncating DFT coeffs to undersample
+
+    if cyclic
+      D = fft(maskdB_cyclic);
+    else
+      D = fft(maskdB);
+    end
+    Dabs = abs(D);                        % this returned for plotting
+    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_));
+
+    % OK now fix up last 500Hz, taper down 10dB at 4000Hz
+
+    xpts = [ anchor-1 anchor L];
+    ypts = [ maskdB_(anchor-1) maskdB_(anchor) (maskdB_(anchor)-10)];
+    mask_pp = splinefit(xpts, ypts, 1);
+    maskdB_ = [maskdB_(1:anchor) ppval(mask_pp, anchor+1:L)];
+
+endfunction
+
+
 % Create a "decimated mask" model using just a few samples of
 % critical band filter centre frequencies.  For voiced speech,
 % we fit the amplitude of these samples to a straight line.
index b649bafb4f859bf98e2951d073aa35a7da44eefe..796a31d8cdba3735185c560c947636473006b560 100644 (file)
 
 % process a whole file and write results
 
-function [non_masked_f_log non_masked_amp_log] = newamp_batch(samname, optional_Am_out_name, optional_Aw_out_name)
+function dk_log = newamp_batch(samname, optional_Am_out_name, optional_Aw_out_name)
   newamp;
   more off;
 
   max_amp = 80;
-  decimate_in_freq = 1;
+  dec_in_freq = 1;
   postfilter = 1;
-  decimate_in_time = 0;
-  synth_phase = 1;
+  dec_in_time = 0;
+  synth_phase = 0;
   freq_quant = 0;
   amp_quant = 0;
-  non_masked_f_log = [];
-  non_masked_m_log = [];
-  non_masked_amp_log = [];
-  best_min_mse_log = [];
+  dk_log = [];
 
   model_name = strcat(samname,"_model.txt");
   model = load(model_name);
@@ -59,6 +56,7 @@ function [non_masked_f_log non_masked_amp_log] = newamp_batch(samname, optional_
 
   %pp_bw = gen_pp_bw;
 
+  sd_sum = 0;
   for f=1:frames
     printf("%d ", f);   
     model_(f,2) = L = min([model(f,2) max_amp-1]);
@@ -72,35 +70,23 @@ function [non_masked_f_log non_masked_amp_log] = newamp_batch(samname, optional_
     mask_sample_freqs_kHz = (1:L)*Wo*4/pi;
     maskdB = mask_model(AmdB, Wo, L);
 
-    if decimate_in_freq
-      % decimate mask samples in frequency
-      [decmaskdB masker_freqs_kHz masker_amps_dB min_error mse_log1 best_min_mse ind_log] = make_decmask_mel(maskdB, AmdB, Wo, L, mask_sample_freqs_kHz, freq_quant, amp_quant);
-      best_min_mse_log = [best_min_mse_log best_min_mse];
-      ind_log
+    % save post filter positions for decode
 
-      non_masked_amp = masker_amps_dB ;
-      non_masked_amp_log = [non_masked_amp_log; non_masked_amp];
+    maskdB_ = maskdB;
+    a_non_masked_m = find(AmdB > maskdB);
+    non_masked_m(f,1:length(a_non_masked_m)) = a_non_masked_m;
 
-      % Save this for decoder, so it knows where to apply post filter
-      % Basically the frequencies of the AbyS samples
-
-      non_masked_m(f,:) = min(round(masker_freqs_kHz/(Wo*4/pi)),L);
-      non_masked_m_log = [non_masked_m_log; non_masked_m(f,:)];
-
-      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
-      maskdB_ = maskdB;
-      a_non_masked_m = find(AmdB > maskdB);
-      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);
+      dk_log = [dk_log; dk_];
     end
+    sd_sum += sum(maskdB - maskdB_);
 
     Am_ = zeros(1,max_amp);
     Am_ = 10 .^ (maskdB_(1:L-1)/20); 
     model_(f,3:(L+1)) = Am_;
   end
+  printf("\nsd_sum: %5.2f\n", sd_sum/frames);
 
   % decoder loop -----------------------------------------------------
 
@@ -112,9 +98,8 @@ function [non_masked_f_log non_masked_amp_log] = newamp_batch(samname, optional_
 
     maskdB_ = 20*log10(Am_);
     mask_sample_freqs_kHz = (1:L)*Wo*4/pi;
-    %maskdB_ = mask_model(maskdB_, Wo, L);
 
-    if decimate_in_time
+    if dec_in_time
       % decimate mask samples in time
       maskdB_ = decimate_frame_rate(model_, 4, f, frames, mask_sample_freqs_kHz);
 
@@ -135,7 +120,6 @@ function [non_masked_f_log non_masked_amp_log] = newamp_batch(samname, optional_
 
       m = max(find(non_masked_m(f,:) > 0)); 
       a_non_masked_m = non_masked_m(f,1:m);
-
     end
 
     % post filter - bump up samples by 6dB, reduce mask by same level to normalise gain
@@ -171,22 +155,6 @@ function [non_masked_f_log non_masked_amp_log] = newamp_batch(samname, optional_
       Aw(1:2:fft_enc*2) = cos(phase);
       Aw(2:2:fft_enc*2) = -sin(phase);
 
-      if 0
-         % optional plotting to ensure amd and phase aligned on a few frames
-         figure(1)
-         clf;
-         subplot(211)
-         plot(mask_sample_freqs_kHz, maskdB_);
-         hold on;
-         plot(mask_sample_freqs_kHz, maskdB_pf,'g');
-         hold off;
-         axis([0 4 0 70])
-
-         subplot(212)
-         plot((0:(fft_enc/2)-1)*8000/fft_enc, phase(1:fft_enc/2))
-         axis([0 4000 -pi pi])
-         k = kbhit();         
-      end 
       fwrite(faw, Aw, "float32");    
     end
   end
@@ -198,8 +166,6 @@ function [non_masked_f_log non_masked_amp_log] = newamp_batch(samname, optional_
 
   printf("\n");
 
-  printf("mse std: %5.2f\n", std(best_min_mse_log));
-
 endfunction
 
 
index f78081431f908b96d71272d4967dc9aecd0cc51b..17e6fbb4d55f5399d6e74ad829799660b1afcbe1 100644 (file)
 %   $ cd ~/codec2-dev/octave
 %   octave:14> newamp_fbf("../build_linux/src/hts1a",50)
 
+
+
 function newamp_fbf(samname, f=10)
-  
-  more off;
   newamp;
-  phase_stuff = 0;
-  plot_not_masked = 0;
+  more off;
   plot_spectrum = 1;
-  freq_quant = 0;
-  amp_quant = 0;
-  decimate_in_freq = 0;
+
+  % load up text files dumped from c2sim ---------------------------------------
 
   sn_name = strcat(samname,"_sn.txt");
   Sn = load(sn_name);
-
   sw_name = strcat(samname,"_sw.txt");
   Sw = load(sw_name);
-
   model_name = strcat(samname,"_model.txt");
   model = load(model_name);
   [frames tmp] = size(model);
 
-  if phase_stuff
-    ak_name = strcat(samname,"_ak.txt");
-    ak = load(ak_name);
-  end
-
-  % pp_bw = gen_pp_bw;
+  % Keyboard loop --------------------------------------------------------------
 
-  plot_all_masks = 0;
   k = ' ';
   do 
     figure(1);
@@ -59,67 +49,20 @@ function newamp_fbf(samname, f=10)
     Am = model(f,3:(L+2));
     AmdB = 20*log10(Am);
 
-    % plotting
-
     axis([1 4000 -20 80]);
     hold on;
     if plot_spectrum
       plot((1:L)*Wo*4000/pi, AmdB,";Am;r");
       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);
+    [tmp1 tmp2 D] = decimate_in_freq(maskdB, 0);
+    [maskdB_ maskdB_cyclic D_cyclic dk_] = decimate_in_freq(maskdB, 1);
 
-    plot(Am_freqs_kHz*1000, maskdB, 'g');
-
-    % Lets try to come up with a smoothed model.  Replace points from 3500 Hz on with
-    % those than match low end of spectrum.  This will make it more cyclical and make
-    % DFT happier, less high freq energy.  Need a better way to describe that.
-
-    anchor = floor(7*L/8);
-    xpts = [ anchor-1 anchor L+1 L+2];
-    ypts = [ maskdB(anchor-1) maskdB(anchor) maskdB(1) maskdB(2)];
-    mask_pp = splinefit(xpts, ypts, 1);
-    maskdB_smooth = [maskdB(1:anchor) ppval(mask_pp, anchor+1:L)];
-
-    plot(Am_freqs_kHz*1000, maskdB_smooth, 'b');
-
-    if decimate_in_freq
-      % decimate in frequency
-
-      mask_sample_freqs_kHz = (1:L)*Wo*4/pi;
-      [decmaskdB masker_freqs_kHz masker_amps_dB min_error mse_log1] = make_decmask_mel(maskdB, AmdB, Wo, L, mask_sample_freqs_kHz, freq_quant, amp_quant);
-
-    % find turning points - prototype for finding PF freqs when we decimate in time
-    
-    if 0
-      d = decmaskdB(2:L) - decmaskdB(1:L-1);
-      size(d)
-      tp = [];
-      for m=1:L-1
-        if (d(m) > 0) && (d(m+1) < 0)
-          tp = [tp m+1];
-        end
-      end
-    end
-
-      figure(2)
-    
-      tonef_kHz = masker_freqs_kHz;
-      nlm = length(tonef_kHz);
-
-      plot(tonef_kHz*1000, zeros(1,nlm), ';AbyS Mask Freqs;bk+');
-      if decimate_in_freq
-        plot(mask_sample_freqs_kHz*1000, decmaskdB, ';AbyS Mask;m');
-      end
-      %plot(tp*Wo*4000/pi, 5*ones(1,length(tp)), ';PF Freqs;r+');
-      legend('boxoff');
-      %plot(mask_sample_freqs_kHz*1000, min_error);
-
-    end
-
-    hold off;
+    plot(Am_freqs_kHz*1000, maskdB, ';mask;g');
+    plot(Am_freqs_kHz*1000, maskdB_cyclic, ';mask cyclic;b');
+    plot(Am_freqs_kHz*1000, maskdB_, ';mask trunated;c');
 
     % lets get a feel for the "spectrum" of the smoothed spectral envelope
     % this will give us a feel for how hard it is to code, ideally we would like
@@ -128,117 +71,41 @@ function newamp_fbf(samname, f=10)
     figure(3)
     clf
 
-    D = abs(fft(maskdB));
-    en = floor(length(D)/2+1);
+    en = L/2+1;
     stem(D(2:en),'g')
     hold on;
-    D_smooth = abs(fft(maskdB_smooth));
-    stem(D_smooth(2:en),'b')
+    stem(D_cyclic(2:en),'b')
     hold off;
 
     % let plot the cumulative amount of energy in each DFT
 
     figure(4)
     clf
-    plot(cumsum(D(2:en)/sum(D(2:en))),'g');
+    plot(cumsum(D(2:en)/sum(D(2:en))),';cumsum;g');
     hold on;
-    plot(cumsum(D_smooth(2:en)/sum(D_smooth(2:en))),'b');
+    plot(cumsum(D_cyclic(2:en)/sum(D_cyclic(2:en))),';cumsum cyclic;b');
     hold off;
     axis([1 L 0 1])
 
-    % try truncating DFT coeffs and see what that looks like
-
-    D = fft(maskdB_smooth);
-    keep = 7;
-    D(keep+1:L-keep) = 0;
-    d = ifft(D);
-
-    figure(2)
-    hold on;
-    plot(Am_freqs_kHz*1000, real(d), 'c')
-    hold off;
-
-    % decimated in time
-
-    %maskdB = decimate_frame_rate(maskdB, model, 4, f, frames, mask_sample_freqs_kHz);
-    %plot(mask_sample_freqs_kHz*1000, maskdB, 'k');
-
-    % optionally plot all masking curves
-
-    if plot_all_masks
-      mask_sample_freqs_kHz = (1:L)*Wo*4/pi;
-      for m=1:L
-        maskdB = schroeder(m*Wo*4/pi, mask_sample_freqs_kHz) + AmdB(m);
-        plot(mask_sample_freqs_kHz*1000, maskdB, "k--");
-      end
-    end
-
-    hold off;
-
-    if phase_stuff
-
-      [phase Sdb s Aw] = determine_phase(model, f, ak(f,:));
-      figure(3)
-      subplot(211)
-      plot(Sdb)
-      title('Mag (dB)');
-      subplot(212)
-      plot(phase(1:256))
-      hold on;
-      plot(angle(Aw(1:256))+0.5,'g')
-      hold off;
-      title('Phase (rads)');
-      figure(4)
-      plot(s)
-    end
+    figure(5)
+    clf
+    stem(dk_)
 
-    % interactive menu
-
-    printf("\rframe: %d  menu: n-next  b-back a-Am ", f);
-    if freq_quant
-      printf("F");
-    else
-      printf("f");
-    end
-    printf("-freq_quant ");
-    if amp_quant
-      printf("M");
-    else
-      printf("m");
-    end
-    printf("-amp_quant q-quit");
+    % interactive menu ------------------------------------------
 
+    printf("\rframe: %d  menu: n-next  b-back q-quit", f);
     fflush(stdout);
     k = kbhit();
+
     if (k == 'n')
       f = f + 1;
     endif
     if (k == 'b')
       f = f - 1;
     endif
-    if k == 'm'
-      if amp_quant == 0
-         amp_quant = 3;
-      else
-         amp_quant = 0;
-      end
-    end
-    if k == 'f'
-      if freq_quant == 0
-        freq_quant = 3;
-      else
-        freq_quant = 0;
-      end
-    end
-    if k == 'a'
-      if plot_spectrum == 0
-         plot_spectrum = 1;
-      else
-         plot_spectrum = 0;
-      end
-    end
   until (k == 'q')
   printf("\n");
 
 endfunction
 
+