From: drowe67 Date: Sun, 13 Aug 2017 04:00:15 +0000 (+0000) Subject: investigated and fixed a bug in Hilbert Transform based phase synthesis, hts1a quite... X-Git-Url: http://git.whiteaudio.com/gitweb/?a=commitdiff_plain;h=e5dc3d5febe46e6fab65ca6dcb673724a92bde54;p=freetel-svn-tracking.git investigated and fixed a bug in Hilbert Transform based phase synthesis, hts1a quite good with phase0, even after dec by 4 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 35a14ea3..4ff241fe 100644 --- 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 86473394..3a300f6a 100644 --- 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 index 00000000..8d7f96ce --- /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');