investigated and fixed a bug in Hilbert Transform based phase synthesis, hts1a quite...

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

diff --git a/codec2-dev/octave/newamp1_batch.m b/codec2-dev/octave/newamp1_batch.m
index 35a14ea3758074246411a1b8a84509d56a57d4a2..4ff241fe8b20371ab9edffc027907b33fdbc988e 100644 (file)
--- a/codec2-dev/octave/newamp1_batch.m
+++ b/codec2-dev/octave/newamp1_batch.m
@@ -154,20 +154,20 @@ function [surface_no_mean surface] = newamp1_batch(input_prefix, varargin)
 
         % synthesis phase spectra from magnitiude spectra using minimum phase techniques
 
-        fft_enc = 256;
+        fft_enc = 512;
         phase = determine_phase(model_, f, fft_enc);
         assert(length(phase) == fft_enc);
 
-        % sample phase at centre of each harmonic, not 1st entry Hm[1] in octave Hm[0] in C
+        % sample phase at centre of each harmonic, not 1st entry Hm(1:2) in octave Hm[0] in C
         % is not used
 
         Hm = zeros(1, 2*max_amp);
         for m=1:L
           b = round(m*Wo*fft_enc/(2*pi));
-          Hm(2*m) = cos(phase(b));
-          Hm(2*m+1) = -sin(phase(b));
+          Hm(2*m+1) = cos(phase(b));
+          Hm(2*m+2) = sin(phase(b));
         end
-        fwrite(fhm, Hm, "float32");    
+        fwrite(fhm, Hm, "float32");
       end
     end
  
@@ -382,16 +382,9 @@ function [model_ rate_K_surface b_log] = experiment_const_freq(model, varargin)
       end
 
       if strcmp(vq_search, "slope")
-        [idx contrib errors b_log] = vq_search_slope(vq, rate_K_surface_no_mean(:,vq_st:vq_en), "closed_quant_slope");
+        [idx contrib errors b_log] = vq_search_slope(vq, rate_K_surface_no_mean(:,vq_st:vq_en),
+                                                     "closed_quant_slope", "open_quant_slope");
         b_log = [b_log energy' idx'];
-
-        for f=1:frames
-          target = rate_K_surface_no_mean(f,vq_st:vq_en);
-          b_log(f,2) = quantise([-1 -0.5 0.5 1], b_log(f,2));
-          b_log(f,3) = (sum(target) - sum(b_log(f,1)*vq(idx(f),:)+b_log(f,2)*(1:vq_cols)))/vq_cols;
-          contrib(f,:) = b_log(f,1)*vq(idx(f),:) + b_log(f,2)*(1:vq_cols) + b_log(f,3);
-          %rate_K_surface_no_mean_(f, vq_en+1:K) -= b_log(f,2)*vq_cols + b_log(f,3);
-        end
       end
 
       if strcmp(vq_search, "para")
@@ -459,7 +452,7 @@ function [model_ rate_K_surface b_log] = experiment_const_freq(model, varargin)
         hmg = hsl = zeros(1,frames);
         for f=1:frames
           hmg(f) = b_log(f, 1);
-          hsl(f) = b_log(f, 3);
+          hsl(f) = b_log(f, 2);
         end
         figure(fg++); clf; hist(hmg, 30); title('mag')
         figure(fg++); clf; hist(hsl, 30); title('slope')

diff --git a/codec2-dev/octave/newamp1_fbf.m b/codec2-dev/octave/newamp1_fbf.m
index 86473394acd5c0470771984788ad9a0071a17630..3a300f6a02ed0b11574ef3db54529b4803e43aa6 100644 (file)
--- a/codec2-dev/octave/newamp1_fbf.m
+++ b/codec2-dev/octave/newamp1_fbf.m
@@ -59,6 +59,14 @@ function newamp1_fbf(samname, f=73, varargin)
     fit_order = 1;
   end
 
+  % optional exploration of phase
+
+  ind = arg_exists(varargin, "phase");
+  phase_en = 0;
+  if ind
+    phase_en = 1;
+  end
+  
   if quant_en
     printf("quant_en: %d vq_filename: %s vq_st: %d vq_en: %d vq_search: %s\n", 
            quant_en, vq_filename, vq_st, vq_en, vq_search);
@@ -74,12 +82,18 @@ function newamp1_fbf(samname, f=73, varargin)
   model = load(model_name);
   [frames tmp] = size(model);
 
+  if phase_en
+    phase_name = strcat(samname,"_phase.txt");
+    phase = load(phase_name);
+  end
+  
   % Keyboard loop --------------------------------------------------------------
 
   k = ' ';
   do 
+    fg = 1;
     s = [ Sn(2*f-1,:) Sn(2*f,:) ];
-    figure(1); clf; plot(s); axis([1 length(s) -20000 20000]);
+    figure(fg++); clf; plot(s); axis([1 length(s) -20000 20000]);
 
     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;
@@ -94,7 +108,7 @@ function newamp1_fbf(samname, f=73, varargin)
 
     % plots ----------------------------------
   
-    figure(2); clf;
+    figure(fg++); clf;
     l = sprintf(";rate %d AmdB;g+-", L);
     plot((1:L)*Wo*4000/pi, AmdB, l);
     axis([1 4000 -20 80]);
@@ -249,14 +263,50 @@ function newamp1_fbf(samname, f=73, varargin)
     end
     hold off;
 
+    if phase_en
+
+      % est phase using HT
+
+      Am_ = model(f,3:(L+2));
+      fft_enc = 512;
+      phase_est = determine_phase(model_, 1, fft_enc);
+      phase0 = zeros(1,L);
+      for m=1:L
+        b = round(m*Wo*fft_enc/(2*pi));
+        phase0(m) = phase_est(b);
+      end
+      
+      % plot amplitudes and phase for first 1kHz
+      
+      figure(fg++); clf;
+      subplot(211);
+      plot((1:L)*Wo*4000/pi, AmdB_(1:L),'+-');
+      subplot(212);
+      plot((1:L)*Wo*4000/pi, phase(f,1:L),'+-');
+      hold on;
+      plot((1:L)*Wo*4000/pi, phase0(1:L),'r+-');
+      hold off;
+      
+      % simple synthesis using sinusoidal parameters
+
+      figure(fg++); clf;
+      N = 320;
+      s = s_phase0 = zeros(1,N);
+      for m=1:L
+        s = s + Am_(m)*cos(m*Wo*(1:N) + phase(f,m));
+        s_phase0 = s_phase0 + Am_(m)*cos(m*Wo*(1:N) + phase0(m));
+      end
+      subplot(211); plot(s); subplot(212); plot(s_phase0,'g');
+    end
+    
     if quant_en
-      figure(4); clf;
+      figure(fg++); clf;
       plot(contrib1, 'b+-');
       hold on; plot(contrib,'r+'); hold off;
     end
     
     if weight_en
-      figure(3); clf;
+      figure(fg++); clf;
       subplot(211);
       plot((1:L)*Wo*4000/pi, AmdB,";AmdB;g+-");
       axis([1 4000 -20 80]);

diff --git a/codec2-dev/octave/tdetphase.m b/codec2-dev/octave/tdetphase.m
new file mode 100644 (file)
index 0000000..8d7f96c
--- /dev/null
+++ b/codec2-dev/octave/tdetphase.m
@@ -0,0 +1,84 @@
+% tdetphase.m
+% David Rowe August 2017
+%
+% Testing Hilbert Transform recover of phase from magnitude spectra
+
+newamp;
+Fs = 8000;
+
+w = 2*pi*500/Fs; gamma = 0.95
+ak = [1 -2*gamma*cos(w) gamma*gamma];
+Nfft = 512;
+
+% Test 1 - compare phase from freqz for 2nd order system (all pole filter)
+%        - uses internal test of determine_phase()
+
+h = freqz(1,ak,Nfft/2);
+
+% note dummy_model not used, as determine_phase() is used in test mode
+
+L = 20; Wo = pi/(L+1);
+dummy_model = [Wo L ones(1,L)];
+phase = determine_phase(dummy_model, 1, Nfft, ak);
+
+fg = 1;
+figure(fg++); clf;
+subplot(211); plot(20*log10(abs(h))); title('test 1');
+subplot(212); plot(angle(h)); hold on; plot(phase(1:Nfft/2),'g'); hold off;
+
+% Test 2 - feed in harmonic magnitudes
+
+F0 = 100; Wo = 2*pi*F0/Fs; L = floor(pi/Wo);
+Am = zeros(1,L);
+for m=1:L
+  b = round(m*Wo*Nfft/(2*pi));
+  Am(m) = abs(h(b));
+end
+AmdB = 20*log10(Am);
+model = [Wo L Am];
+[phase Gdbfk s] = determine_phase(model, 1, Nfft);
+
+fftx = (1:Nfft/2)*(Fs/Nfft);
+harmx = (1:L)*Wo*Fs/(2*pi);
+
+figure(fg++); clf;
+subplot(211); plot(fftx, Gdbfk(1:Nfft/2));
+subplot(212); plot(s(1:Nfft/2))
+
+figure(fg++); clf;
+subplot(211); plot(fftx, 20*log10(abs(h)));
+              hold on; plot(harmx, AmdB, 'g+'); plot(fftx, Gdbfk(1:Nfft/2), 'r'); hold off;
+subplot(212); plot(fftx, angle(h)); hold on; plot(fftx, phase(1:Nfft/2),'g'); hold off;
+
+% Test 3 - Use real harmonic amplitudes
+
+model = load("../build_linux/src/hts1a_model.txt");
+phase_orig = load("../build_linux/src/hts1a_phase.txt");
+
+f = 184;
+Wo = model(f,1); L = model(f,2); Am = model(f,3:L+2); AmdB = 20*log10(Am);
+[phase Gdbfk s] = determine_phase(model, f, Nfft);
+
+fftx = (1:Nfft/2)*(Fs/Nfft);
+harmx = (1:L)*Wo*Fs/(2*pi);
+
+figure(fg++); clf;
+subplot(211); plot(fftx, Gdbfk(1:Nfft/2));
+subplot(212); plot(s(1:Nfft/2))
+
+figure(fg++); clf;
+subplot(211); plot(harmx, AmdB, 'g+');
+              hold on; plot(fftx, Gdbfk(1:Nfft/2), 'r'); hold off;
+subplot(212); plot(fftx, phase(1:Nfft/2),'g');
+
+% synthesise using phases
+
+N = 320;
+s = s_phase = zeros(1,N);
+for m=1:L/4
+  s = s + Am(m)*cos(m*Wo*(1:N) + phase_orig(f,m));
+  b = round(m*Wo*Nfft/(2*pi));
+  s_phase = s_phase + Am(m)*cos(m*Wo*(1:N) + phase(b));
+end
+figure(fg++); clf;
+subplot(211); plot(s); subplot(212); plot(s_phase,'g');