about to refactor to support different experiments
authordrowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Mon, 12 Dec 2016 02:12:31 +0000 (02:12 +0000)
committerdrowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Mon, 12 Dec 2016 02:12:31 +0000 (02:12 +0000)
git-svn-id: https://svn.code.sf.net/p/freetel/code@2929 01035d8c-6547-0410-b346-abe4f91aad63

codec2-dev/octave/newamp1_batch.m

index 0fb04737d7acd68fa00cd0f6d693823c018acfac..603c892e1966ae90f006317d5fac19863b366e06 100644 (file)
 %   [ ] switches to turn on/off quantisation
 %   [ ] rename mask_sample_freqs, do we need "mask" any more
 
+#{
+    [ ] refactor
+#}
+
 function [fvec_log amps_log] = newamp1_batch(samname, optional_Am_out_name, optional_Aw_out_name)
   newamp;
   more off;
@@ -33,6 +37,58 @@ function [fvec_log amps_log] = newamp1_batch(samname, optional_Am_out_name, opti
   model = load(model_name);
   [frames nc] = size(model);
 
+  K = 40;
+  [rate_K_surface rate_K_sample_freqs_kHz] = resample_const_rate_f(model, K);
+
+  #{
+  So given surface, lets look at spectral content and see if we can reduce it while maintaining
+  speech quality.  First step is to dft across time and plot.
+  #}
+  
+  %Nf=20; Nf2 = floor(Nf/2)+1;
+  %b = fir1(Nf,0.25)
+  %b = zeros(1, Nf2);
+  Nf = 4; Nf2 = 5;
+  [b a]= cheby1(4, 1, 0.20);
+  dft_surface = zeros(frames,K);
+  rate_K_surface_filt = zeros(frames,K);
+  dft_surface_filt = zeros(frames,K);
+  for k=1:K
+    dft_surface(:,k) = fft(rate_K_surface(:,k).*hanning(frames));
+    rate_K_surface_filt(:,k) = filter(b, a, rate_K_surface(:,k));
+    dft_surface_filt(:,k) = fft(rate_K_surface_filt(:,k).*hanning(frames));
+  end
+  figure(1); clf;
+  mesh(abs(dft_surface))
+  figure(2); clf;
+  mesh(abs(dft_surface_filt))
+
+  Fs = 100; Ts = 1/Fs;
+  figure(3);
+  subplot(211);
+  h = freqz(b,a,Fs/2);
+  plot(1:Fs/2, 20*log10(abs(h)))
+  axis([1 Fs/2 -80 0])
+  ylabel('Gain (dB)')
+  grid;
+  subplot(212)
+  [g w] = grpdelay(b,a);
+  plot(w*Fs/pi, 1000*g*Ts)
+  %axis([1 Fs/2 0 0.5])
+  xlabel('Frequency (Hz)');
+  ylabel('Delay (ms)')
+  grid
+
+  figure(4); clf; stem(filter(b,a,[1 zeros(1,20)]))
+  
+  % adjust for time offset due to filtering
+
+  rate_K_surface_filt = [rate_K_surface_filt(Nf2:frames,:); zeros(Nf2, K)];
+
+  % back down to rate L
+
+  model_ = resample_rate_L(model, rate_K_surface_filt, rate_K_sample_freqs_kHz);
+
   if nargin == 2
     Am_out_name = optional_Am_out_name;
   else
@@ -41,13 +97,28 @@ function [fvec_log amps_log] = newamp1_batch(samname, optional_Am_out_name, opti
 
   fam  = fopen(Am_out_name,"wb"); 
 
+  Wo_out_name = sprintf("%s_Wo.out", samname);
+  fWo  = fopen(Wo_out_name,"wb"); 
+  
+  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));
+
+    Am_ = zeros(1,max_amp);
+    Am_(2:L) = Am(1:L-1);
+    fwrite(fam, Am_, "float32");
+    fwrite(fWo, Wo, "float32");
+  end
+#{
   % encoder loop ------------------------------------------------------
 
   fvec_log = []; amps_log = [];
 
-  % prime initial values so first pass iterpolator works.  TODO: improve this
-
   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));
@@ -60,30 +131,116 @@ function [fvec_log amps_log] = newamp1_batch(samname, optional_Am_out_name, opti
     fvec_log = [fvec_log; fvec];
     amps_log = [amps_log; amps];
 
-    model_(f,1) = Wo; model_(f,2) = L; model_(f,3:(L+2)) = 10 .^ (AmdB_(1:L)/20);
+    model_(f,1) = Wo; model_(f,2) = L; model_(f,3:(L+2)) = 10 .^ (AmdB(1:L)/20);
   end
 
   for f=1:frames
+    if 0
     if (f > 1) && (e(f) > (e(f-1)+3))
         decimate = 2;
       else
         decimate = 4;
     end
+    end
     printf("%d ", decimate);   
-
+    
     AmdB_ = decimate_frame_rate(model_, decimate, f, frames);
     L = length(AmdB_);
 
     Am_ = zeros(1,max_amp);
     Am_(2:L) = 10 .^ (AmdB_(1:L-1)/20);  % C array doesnt use A[0]
     fwrite(fam, Am_, "float32");
-    Am_ = zeros(1,max_amp);
+    fwrite(fWo, Wo, "float32");
   end
+#}
 
   fclose(fam);
+  fclose(fWo);
   printf("\n")
 
-  figure(1); clf;
-  plot(e,'+-');
 endfunction
   
+
+#{
+  [X] resample A(mWo,t) to fixed rate, creating a surface
+  [X] resample surface at half frame increments
+      + this will test resampling in freq is OK
+  [X] listen to a few files and determine if any distortion
+  [X] experiment with filtering in time
+  [ ] decimate to 25Hz
+  [ ] resample to 100Hz
+  [ ] listen to a few files and compare to just filtering at 12Hz
+#}
+
+function [rate_K_surface rate_K_sample_freqs_kHz] = resample_const_rate_f(model, K=50)
+
+  % convert rate L=pi/Wo amplitude samples to fixed rate K
+
+  max_amp = 80;
+  [frames col] = size(model);
+  rate_K_sample_freqs_kHz = (1:K)*4/K;
+  rate_K_surface = zeros(frames, K);
+
+  for f=1:frames
+    Wo = model(f,1);
+    L = min([model(f,2) max_amp-1]);
+    Am = model(f,3:(L+2));
+    AmdB = 20*log10(Am);
+    rate_L_sample_freqs_kHz = (1:L)*Wo*4/pi;
+    
+    rate_K_surface(f,:) = interp1(rate_L_sample_freqs_kHz, AmdB, rate_K_sample_freqs_kHz, "spline", "extrap");
+  end
+endfunction
+
+
+function model_ = resample_rate_L(model, rate_K_surface, rate_K_sample_freqs_kHz, K=50)
+  max_amp = 80;
+  [frames col] = size(model);
+
+  model_ = zeros(frames, max_amp+3);
+  for f=1:frames-1
+    Wo = (model(f,1) + model(f+1,1))/2;
+    L = min(pi/Wo, max_amp-1);
+    rate_L_sample_freqs_kHz = (1:L)*Wo*4/pi;
+    
+    % back down to rate L
+
+    AmdB_ = interp1(rate_K_sample_freqs_kHz, rate_K_surface(f,:), rate_L_sample_freqs_kHz, "spline", "extrap");
+
+    % todo a way to save all model params ... Am, Wo, L, (v?)  Start with current facility then
+    % make it better
+
+    model_(f,1) = Wo; model_(f,2) = L; model_(f,3:(L+2)) = 10 .^ (AmdB_(1:L)/20);
+   end
+endfunction
+
+
+function model_ = resample_half_frame_offset(model, rate_K_surface, rate_K_sample_freqs_kHz, K=50)
+  max_amp = 80;
+  [frames col] = size(model);
+
+  % Check distortion by shifting half a frame and returning model for
+  % synthesis Hmm how to handle Wo?  Well lets assume it's OK and
+  % average it.  If it sounds OK then it must be OK as well and Am
+  % interpolation.  If not then will do some more thinking.
+
+  model_ = zeros(frames, max_amp+3);
+  for f=1:frames-1
+    Wo = (model(f,1) + model(f+1,1))/2;
+    L = min(pi/Wo, max_amp-1);
+    rate_L_sample_freqs_kHz = (1:L)*Wo*4/pi;
+    
+    % interpolate at half frame offset
+
+    half = 0.5*rate_K_surface(f,:) + 0.5*rate_K_surface(f+1,:);
+    
+    % back down to rate L
+
+    AmdB_ = interp1(rate_K_sample_freqs_kHz, half, rate_L_sample_freqs_kHz, "spline", "extrap");
+
+    % todo a way to save all model params ... Am, Wo, L, (v?)  Start with current facility then
+    % make it better
+
+    model_(f,1) = Wo; model_(f,2) = L; model_(f,3:(L+2)) = 10 .^ (AmdB_(1:L)/20);
+   end
+endfunction