Attempt at correcting alias errors in rate K resampling, sig improvement for hst1a
authordrowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Tue, 6 Jun 2017 04:04:11 +0000 (04:04 +0000)
committerdrowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Tue, 6 Jun 2017 04:04:11 +0000 (04:04 +0000)
git-svn-id: https://svn.code.sf.net/p/freetel/code@3157 01035d8c-6547-0410-b346-abe4f91aad63

codec2-dev/octave/newamp.m
codec2-dev/octave/newamp1_batch.m
codec2-dev/octave/newamp1_fbf.m

index 010893425731489a11f703f6317d9462138a9e50..a1c64bc26c63c3b7d94855c4cf4b9f28aaa21f2b 100644 (file)
@@ -46,6 +46,43 @@ function y = interp_para(xp, yp, x)
 endfunction
 
 
+% choose largest sample in band, idea is we care more about finding
+% peaks, can handle some error in frequency. x are non linear
+% (arbitrary) sampling points in kHz
+
+function y = interp_largest(f0_Hz, AmdB, x_kHz)
+  L = length(AmdB);
+  x = x_kHz*1000;
+  y = zeros(1,length(x));
+  bw = x(2) - x(1);
+  k = 1;
+
+  for i=1:length(x)
+
+    % determine limits of this band
+
+    if i>1
+      bw = x(i) - x(i-1);
+    end
+    band_low = x(i) - bw/2; band_high = x(i) + bw/2;
+
+    % map band limits to harmonics    
+    if x(i) < f0_Hz
+      m_low = m_high = 1;
+    else
+      m_low = round(band_low/f0_Hz); m_high = round(band_high/f0_Hz)-1;
+      m_low = max(1, m_low); m_high = min(L, m_high); m_high = max(m_low, m_high);
+    end
+
+    printf("L: %d f0: %f i: %d band_low: %f band_high: %f m_low: %d m_high: %d\n",L, f0_Hz, i, band_low, band_high, m_low, m_high);
+    % find max in band
+
+    y(i) = max(AmdB(m_low:m_high));
+  end
+
+endfunction
+
 % simple linear interpolator
 
 function y = interp_linear(xp, yp, x)
@@ -172,7 +209,7 @@ endfunction
 
 
 function rate_K_sample_freqs_kHz = mel_sample_freqs_kHz(K)
-  mel_start = ftomel(200); mel_end = ftomel(3700); 
+  mel_start = ftomel(100); mel_end = ftomel(0.95*4000); 
   step = (mel_end-mel_start)/(K-1);
   mel = mel_start:step:mel_end;
   rate_K_sample_freqs_Hz = 700*((10 .^ (mel/2595)) - 1);
@@ -222,6 +259,64 @@ function [rate_K_surface rate_K_sample_freqs_kHz] = resample_const_rate_f(model,
 endfunction
 
 
+function [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)
+
+    % aliasing correction --------------------------------------
+
+    % The mel sample rate decreases as frequency increases. Look for
+    % any regions above 1000Hz where we have missed definition of a
+    % spectral peak (formant) due to aliasing.  Adjust the rate K
+    % sample levels to restore peaks.  Theory is that correct
+    % definition of a formant is less important than the frequency of
+    % the formant.  As long as we define a formant in that general
+    % frequency area it will sound OK.
+
+    Am_freqs_kHz = (1:L)*Wo*4/pi;
+
+    % Lets see where we have made an error
+
+    error = orig_error = AmdB - AmdB_;
+
+    Ncorrections = 3;      % maximum number of rate K samples to correct
+    error_thresh = 3;      % only worry about errors larger than thresh
+
+    start_m = floor(L*1000/(Fs/2));
+    error(1:start_m) = 0;  % first 1000Hz is densly sampled so ignore
+    nasty_error_m_log = []; nasty_error_log = [];
+
+
+    rate_K_vec_corrected = rate_K_vec;
+    for i=1:Ncorrections
+      [mx mx_m] = max(error);
+
+      if mx > error_thresh
+        nasty_error_log = [nasty_error_log mx];
+        nasty_error_m_log = [nasty_error_m_log mx_m];
+
+        % find closest rate K sample to nasty error
+
+        nasty_error_freq = mx_m*Wo*Fs/(2*pi*1000)
+        [tmp closest_k] = min(abs(rate_K_sample_freqs_kHz - nasty_error_freq));
+        rate_K_vec_corrected(closest_k) = AmdB(mx_m);
+
+        % zero out error in this region and look for another large error region
+
+        k = max(1, closest_k-1); 
+        rate_K_prev_sample_kHz = rate_K_sample_freqs_kHz(k)
+        k = min(K, closest_k+1); 
+        rate_K_next_sample_kHz = rate_K_sample_freqs_kHz(k)
+
+        [tmp st_m] = min(abs(Am_freqs_kHz - rate_K_prev_sample_kHz));
+        [tmp en_m] = min(abs(Am_freqs_kHz - rate_K_next_sample_kHz));
+        if closest_k == K
+         en_m = L;
+        end 
+        error(st_m:en_m) = 0;
+      end
+    end
+endfunction
+
+
 % Take a rate K surface and convert back to time varying rate L
 
 function [model_ AmdB_] = resample_rate_L(model, rate_K_surface, rate_K_sample_freqs_kHz)
@@ -1266,3 +1361,12 @@ function lmin = abys(AmdB_, AmdB, Wo, L, mask_sample_freqs_kHz)
 endfunction
 
 
+function rate_K_surface_no_slope = remove_slope(rate_K_surface)
+  [frames K] = size(rate_K_surface);
+  rate_K_surface_no_slope = zeros(frames,K);
+  for f=1:frames
+    [gradient intercept] = linreg(1:K, rate_K_surface(f,:), K);
+    printf("f: %d gradient: %f intercept: %f\n", f, gradient, intercept);
+    rate_K_surface_no_slope(f,:) = rate_K_surface(f,:) - (intercept + gradient*(1:K));
+  end
+endfunction
index 58793139a447ae263cc451b1e08364b8935628c6..3dbbbc7cdfce2c1a08da0182a05ea3b6e54d0924 100644 (file)
@@ -7,7 +7,7 @@
 #{
 
   Octave script to batch process model parameters using the new
-  amplitude model, version 1.  Outputs anotehr set of model parameters
+  amplitude model, version 1.  Outputs another 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
  
@@ -44,7 +44,7 @@ function surface = newamp1_batch(input_prefix, output_prefix)
   newamp;
   more off;
 
-  max_amp = 80;
+  max_amp = 160;
   postfilter = 0;   % optional postfiler that runs on Am, not used atm
   synth_phase = 1;
 
@@ -68,12 +68,16 @@ function surface = newamp1_batch(input_prefix, output_prefix)
   %model_ = experiment_filter(model);
   %model_ = experiment_filter_dec_filter(model);
 
-  %[model_ surface] = experiment_mel_freq(model, 1, 1, voicing);
+  [model_ surface] = experiment_mel_freq(model, 0, 1, voicing);
+
   %model_ = experiment_dec_abys(model, 8, 1, 1, 1, voicing);
 
+#{
   [model_ voicing_ indexes] = experiment_rate_K_dec(model, voicing); % encoder/decoder, lets toss away results except for indexes
-  %[model_ voicing_] = model_from_indexes(indexes);                   % decoder uses just indexes, outputs vecs for synthesis
   [model_ voicing_] = model_from_indexes_fbf(indexes);                   % decoder uses just indexes, outputs vecs for synthesis
+#}
+
+  %[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);
@@ -121,7 +125,7 @@ function surface = newamp1_batch(input_prefix, output_prefix)
 
       % synthesis phase spectra from magnitiude spectra using minimum phase techniques
 
-      fft_enc = 128;
+      fft_enc = 256;
       phase = determine_phase(model_, f, fft_enc);
       assert(length(phase) == fft_enc);
       %Aw = zeros(1, fft_enc*2); 
@@ -492,9 +496,21 @@ endfunction
 
 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);
-  
+  K = 20; Fs = 8000;
+
+  for f=1:frames
+    Wo = model(f,1);
+    L = model(f,2);
+    Am = model(f,3:(L+2));
+    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;
+  end
+
   if plots
     figure(1); clf; mesh(rate_K_surface);
     figure(2); clf; plot_dft_surface(rate_K_surface)
@@ -1057,3 +1073,4 @@ function plot_dft_surface(surface)
   Fs = 100; dF = Fs/frames;
   mesh(1:K, (1:frames/2)*dF, 20*log10(abs(dft_surface(1:frames/2,:,:))));
 endfunction
+
index 2946d5de0c5776c3057db7d184e099a98a7afe3d..fe7bb2ba3cc19c7da361d4443a06f185a170f3cb 100644 (file)
 function newamp1_fbf(samname, f=10)
   newamp;
   more off;
-  quant_en = 0; pf_en = 0; plot_phase = 1;
+  quant_en = 0; pf_en = 0; plot_phase = 0;
   melvq;
-
-  K=20; load train_120_vq; m=5; 
+  Fs = 8000; K=20; load train_120_vq; m=5; 
 
   % load up text files dumped from c2sim ---------------------------------------
 
@@ -66,9 +66,14 @@ function newamp1_fbf(samname, f=10)
     AmdB = 20*log10(Am);
     Am_freqs_kHz = (1:L)*Wo*4/pi;
 
-    % plots for mel sampling
+    % resample at rate K
 
     [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);
+    [tmp_ AmdB_corrected] = resample_rate_L(model(f,:), rate_K_vec, rate_K_sample_freqs_kHz);
+
+    % quantisation simulation ---------------------------------
 
     mean_f = mean(rate_K_vec);
     rate_K_vec_no_mean = rate_K_vec - mean_f;
@@ -85,6 +90,7 @@ function newamp1_fbf(samname, f=10)
 
     rate_K_vec_ = rate_K_vec_no_mean_ + mean_f;
     [model_ AmdB_] = resample_rate_L(model(f,:), rate_K_vec_, rate_K_sample_freqs_kHz);
+    %[model_ AmdB_largest] = resample_rate_L(model(f,:), rate_K_largest, rate_K_sample_freqs_kHz);
 
     % plots ----------------------------------
 
@@ -94,17 +100,20 @@ function newamp1_fbf(samname, f=10)
 
     axis([1 4000 -20 80]);
     hold on;
+
     plot((1:L)*Wo*4000/pi, AmdB,";Am;b+-");
-    plot(rate_K_sample_freqs_kHz*1000, rate_K_vec, ';rate K mel;g+-');
-    if quant_en >= 1
-      plot((1:L)*Wo*4000/pi, AmdB_,";Am quant;k+-");
-    end
-    if quant_en == 2
-      plot(rate_K_sample_freqs_kHz*1000, rate_K_vec_, ';rate K mel quant;r+-');   
-    end
+    stem(rate_K_sample_freqs_kHz*1000, rate_K_vec, 'g');   
+    plot(rate_K_sample_freqs_kHz*1000, rate_K_vec, ';rate K;g:');   
+    plot(rate_K_sample_freqs_kHz*1000, rate_K_vec_corrected, ';rate K;g-', 'linewidth', 2);   
+
+    plot((1:L)*Wo*4000/pi, -10 + orig_error ,";AmdB - AmdB_ error;k+-");
+    plot((1:L)*Wo*4000/pi, -10 + error ,";AmdB - AmdB_ error;r-");
+    plot(nasty_error_m_log*Wo*4000/pi, AmdB(nasty_error_m_log), 'ro', 'markersize', 10);
+
 
     hold off;
 
+    #{
     figure(3);
     clf;
     title('Frequency Domain 2');
@@ -116,6 +125,7 @@ function newamp1_fbf(samname, f=10)
     plot(rate_K_sample_freqs_kHz*1000, rate_K_vec_no_mean_, ';rate K mel no mean quant;r+-');
     end
     hold off;
+    #}
 
     if plot_phase
       phase_512 = determine_phase(model_, 1, 512);
@@ -135,7 +145,7 @@ function newamp1_fbf(samname, f=10)
 
     if k == 'm'
       quant_en++;
-      if quant_en > 2
+      if quant_en > 1
         quant_en = 0;
       end
     endif