some initial fbf tests on undersampling spectrum, cyclic extension etc

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

diff --git a/codec2-dev/octave/newamp.m b/codec2-dev/octave/newamp.m
index 3b96395d7e5d84625d0563999baf5e00c50104db..8b61229bdc4c734ef2039cca22df7838dd3df37c 100644 (file)
--- a/codec2-dev/octave/newamp.m
+++ b/codec2-dev/octave/newamp.m
@@ -157,6 +157,34 @@ function mse = search_mask(target, decmaskdB, mask_sample_freqs_kHz, AmdB, Wo, l
 end
 
 
+% simple resampling to a fixed length vector on mel scale
+% fixed freq grid resampling
+
+function [decmaskdB mel_sample_freqs_kHz mel_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)  
+
+  % set up mel sampling grid
+
+  Nmel = 20;
+  mel_st = freq2mel(Wo*4000/pi);
+  mel_en = freq2mel(L*Wo*4000/pi);
+  m_step = (mel_en-mel_st)/Nmel;
+  m = mel_st:m_step:mel_en;
+  mel_sample_freqs_kHz = mel2freq(m)/1000;
+
+  % resample on mel grid
+
+  mask_pp = splinefit(mask_sample_freqs_kHz, maskdB, L);
+  mel_masker_amps_dB = ppval(mask_pp, mel_sample_freqs_kHz);
+
+  % resample on Wo grid
+
+  mel_pp = splinefit(mel_sample_freqs_kHz, mel_masker_amps_dB, L);
+  decmaskdB = ppval(mel_pp, mask_sample_freqs_kHz);
+  size(decmaskdB)
+  min_error = 0; mse_log1=[]; best_min_mse=0; ind_log=[];
+end
+
+
 function mel = freq2mel(f)
   mel = 70*log10(1 + f/700);
 endfunction

diff --git a/codec2-dev/octave/newamp_batch.m b/codec2-dev/octave/newamp_batch.m
index 2dd4e8ca229e9af2404616e02207969e5e2273b8..b649bafb4f859bf98e2951d073aa35a7da44eefe 100644 (file)
--- a/codec2-dev/octave/newamp_batch.m
+++ b/codec2-dev/octave/newamp_batch.m
@@ -23,10 +23,10 @@ function [non_masked_f_log non_masked_amp_log] = newamp_batch(samname, optional_
   max_amp = 80;
   decimate_in_freq = 1;
   postfilter = 1;
-  decimate_in_time = 1;
+  decimate_in_time = 0;
   synth_phase = 1;
-  freq_quant = 3;
-  amp_quant = 3;
+  freq_quant = 0;
+  amp_quant = 0;
   non_masked_f_log = [];
   non_masked_m_log = [];
   non_masked_amp_log = [];
@@ -74,7 +74,7 @@ function [non_masked_f_log non_masked_amp_log] = newamp_batch(samname, optional_
 
     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_abys(maskdB, AmdB, Wo, L, mask_sample_freqs_kHz, freq_quant, amp_quant);
+      [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
 

diff --git a/codec2-dev/octave/newamp_fbf.m b/codec2-dev/octave/newamp_fbf.m
index 5634e05692a3582fa4b59060ccbefa3ba6522503..f78081431f908b96d71272d4967dc9aecd0cc51b 100644 (file)
--- a/codec2-dev/octave/newamp_fbf.m
+++ b/codec2-dev/octave/newamp_fbf.m
@@ -23,6 +23,7 @@ function newamp_fbf(samname, f=10)
   plot_spectrum = 1;
   freq_quant = 0;
   amp_quant = 0;
+  decimate_in_freq = 0;
 
   sn_name = strcat(samname,"_sn.txt");
   Sn = load(sn_name);
@@ -56,8 +57,6 @@ function newamp_fbf(samname, f=10)
     Wo = model(f,1);
     L = model(f,2);
     Am = model(f,3:(L+2));
-    %[h w] = freqz([1 -1],1,(1:L)*Wo); % pre-emphasise to reduce dynamic range
-    %Am = Am .* abs(h);
     AmdB = 20*log10(Am);
 
     % plotting
@@ -72,51 +71,92 @@ function newamp_fbf(samname, f=10)
 
     [maskdB Am_freqs_kHz] = mask_model(AmdB, Wo, L);
 
-    %plot(Am_freqs_kHz*1000, maskdB, 'g');
-    %plot(Am_freqs_kHz2*1000, maskdB2, 'r');
+    plot(Am_freqs_kHz*1000, maskdB, 'g');
 
-    % optionally plot synthesised spectrum (early simple model)
+    % 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.
 
-    if 0
-      AmdB_ = maskdB;
-      plot(Am_freqs_kHz*1000, AmdB_, 'g');
-      plot(Am_freqs_kHz*1000, AmdB_, 'g+');
-    end
+    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)];
 
-    % decimate in frequency
+    plot(Am_freqs_kHz*1000, maskdB_smooth, 'b');
 
-    mask_sample_freqs_kHz = (1:L)*Wo*4/pi;
-    [decmaskdB masker_freqs_kHz masker_amps_dB min_error mse_log1] = make_decmask_abys(maskdB, AmdB, Wo, L, mask_sample_freqs_kHz, freq_quant, amp_quant);
+    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
     
-    d = decmaskdB(2:L) - decmaskdB(1:L-1);
-    tp = [];
-    for m=1:L-1
-      if (d(m) > 0) && (d(m+1) < 0)
-        tp = [tp m+1];
+    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)
+      figure(2)
     
-    tonef_kHz = masker_freqs_kHz;
-    nlm = length(tonef_kHz);
-
-    plot(tonef_kHz*1000, zeros(1,nlm), ';AbyS Mask Freqs;bk+');
-    plot(mask_sample_freqs_kHz*1000, decmaskdB, ';AbyS Mask;m');
-    plot(tp*Wo*4000/pi, 5*ones(1,length(tp)), ';PF Freqs;r+');
-    legend('boxoff');
-    %plot(mask_sample_freqs_kHz*1000, min_error);
+      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;
 
-    if 0
+    % 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
+    % just a few coefficents to be non-zero
+
     figure(3)
     clf
-    plot((1:L)*Wo*4000/pi, mse_log1');
-    axis([0 4000 0 max(mse_log1(1,:))])
-    title('Error as a function of position for each AbyS stage');
-    end
+
+    D = abs(fft(maskdB));
+    en = floor(length(D)/2+1);
+    stem(D(2:en),'g')
+    hold on;
+    D_smooth = abs(fft(maskdB_smooth));
+    stem(D_smooth(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');
+    hold on;
+    plot(cumsum(D_smooth(2:en)/sum(D_smooth(2:en))),'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