From 6bf6e35051ce75f93d2fc16ec084aa1e39f8a8f8 Mon Sep 17 00:00:00 2001 From: drowe67 Date: Wed, 28 Jun 2017 02:11:04 +0000 Subject: [PATCH] pager script to explore VQ entries git-svn-id: https://svn.code.sf.net/p/freetel/code@3261 01035d8c-6547-0410-b346-abe4f91aad63 --- codec2-dev/octave/newamp1_batch.m | 182 +++++++++++++++++++----------- codec2-dev/octave/newamp1_fbf.m | 68 ++++------- codec2-dev/octave/vq_pager.m | 47 ++++++++ 3 files changed, 185 insertions(+), 112 deletions(-) create mode 100644 codec2-dev/octave/vq_pager.m diff --git a/codec2-dev/octave/newamp1_batch.m b/codec2-dev/octave/newamp1_batch.m index 0010c377..16a23c18 100644 --- a/codec2-dev/octave/newamp1_batch.m +++ b/codec2-dev/octave/newamp1_batch.m @@ -40,16 +40,18 @@ % In general, this function processes a bunch of amplitudes, we then % use c2sim to hear the results. Bunch of different experiments below -function surface = newamp1_batch(input_prefix, varargin) +function [surface mean_f] = newamp1_batch(input_prefix, varargin) newamp; more off; max_amp = 160; + mean_f = []; % defaults - synth_phase = 1; + synth_phase = output = 1; output_prefix = input_prefix; + vq_type = ""; vq_filename = ""; mean_removal = 0; mode = "const"; @@ -72,11 +74,21 @@ function surface = newamp1_batch(input_prefix, varargin) if ind mean_removal = 1; end - ind = arg_exists(varargin, "vq"); + ind = arg_exists(varargin, "vqh"); if ind + vq_type = varargin{ind}; vq_filename = varargin{ind+1}; + printf("vq_type: %s vq_filename: %s\n", vq_type, vq_filename); + end + ind = arg_exists(varargin, "no_output"); + if ind + output = 0; + synth_phase = 0; end end + + printf("output_prefix: %s\nnomean: %d output: %d vq_type: %s\nvq_filename: %s\n", + output_prefix, mean_removal, output, vq_type, vq_filename); model_name = strcat(input_prefix,"_model.txt"); model = load(model_name); @@ -99,7 +111,7 @@ function surface = newamp1_batch(input_prefix, varargin) [model_ surface] = experiment_mel_freq(model, 0, 1, voicing); end if strcmp(mode, 'const') - [model_ surface] = experiment_const_freq(model, mean_removal, vq_filename); + [model_ surface mean_f] = experiment_const_freq(model, mean_removal, vq_type, vq_filename); end if strcmp(mode, 'piecewise') model_ = experiment_piecewise(model); @@ -110,61 +122,69 @@ function surface = newamp1_batch(input_prefix, varargin) % ---------------------------------------------------- - Am_out_name = sprintf("%s_am.out", output_prefix); - fam = fopen(Am_out_name,"wb"); - - Wo_out_name = sprintf("%s_Wo.out", output_prefix); - fWo = fopen(Wo_out_name,"wb"); + if output + Am_out_name = sprintf("%s_am.out", output_prefix); + fam = fopen(Am_out_name,"wb"); + + Wo_out_name = sprintf("%s_Wo.out", output_prefix); + fWo = fopen(Wo_out_name,"wb"); - if synth_phase - Hm_out_name = sprintf("%s_hm.out", output_prefix); - fhm = fopen(Hm_out_name,"wb"); - end + if synth_phase + Hm_out_name = sprintf("%s_hm.out", output_prefix); + fhm = fopen(Hm_out_name,"wb"); + end - 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)); - assert(Wo*L < pi); + 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)); + if Wo*L > pi + printf("Problem: %d Wo*L > pi\n", f); + end - Am_ = zeros(1,max_amp); Am_(2:L) = Am(1:L-1); fwrite(fam, Am_, "float32"); - fwrite(fWo, Wo, "float32"); + Am_ = zeros(1,max_amp); Am_(2:L) = Am(1:L-1); fwrite(fam, Am_, "float32"); + fwrite(fWo, Wo, "float32"); - if synth_phase + if synth_phase - % synthesis phase spectra from magnitiude spectra using minimum phase techniques + % synthesis phase spectra from magnitiude spectra using minimum phase techniques - fft_enc = 256; - phase = determine_phase(model_, f, fft_enc); - assert(length(phase) == fft_enc); + fft_enc = 256; + 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 - % is not used + % sample phase at centre of each harmonic, not 1st entry Hm[1] 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 = 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)); + end + fwrite(fhm, Hm, "float32"); end - fwrite(fhm, Hm, "float32"); + fclose(fhm); end - end - fclose(fam); - fclose(fWo); - if synth_phase - fclose(fhm); - end + fclose(fam); + fclose(fWo); - % save voicing file + % save voicing file - if exist("voicing_", "var") - v_out_name = sprintf("%s_v.txt", output_prefix); - fv = fopen(v_out_name,"wt"); - for f=1:length(voicing_) - fprintf(fv,"%d\n", voicing_(f)); + if exist("voicing_", "var") + v_out_name = sprintf("%s_v.txt", output_prefix); + fv = fopen(v_out_name,"wt"); + for f=1:length(voicing_) + fprintf(fv,"%d\n", voicing_(f)); + end + fclose(fv); + end + end + + if mean_removal + for f=1:frames + surface(f,:) -= mean_f(f); end - fclose(fv); end printf("\n") @@ -181,13 +201,32 @@ function ind = arg_exists(v, str) end endfunction + % Basic unquantised rate K linear sampling then back to rate L. Used for generating % training vectors and testing vector quntisers. -function [model_ rate_K_surface] = experiment_const_freq(model, mean_removal, vq_filename) +function [model_ rate_K_surface mean_f] = experiment_const_freq(model, mean_removal, vq_type, vq_filename) + melvq; + [frames nc] = size(model); Fs = 8000; fg = 1; + quant_en = 0; + + rate_K_sample_freqs_kHz = [0.1:0.1:4]; + K = length(rate_K_sample_freqs_kHz); + + if strcmp(vq_type, "vql") + quant_en = 1; + x = load(vq_filename); vq = x.vq; + [vq_rows vq_cols] = size(vq); vq_st = 1; vq_en = vq_st + vq_cols - 1; + end + + if strcmp(vq_type, "vqh") + quant_en = 1; + x = load(vq_filename); vq = x.vq; + [vq_rows vq_cols] = size(vq); vq_st = 11; vq_en = K; + end energy = zeros(1,frames); for f=1:frames @@ -195,38 +234,55 @@ function [model_ rate_K_surface] = experiment_const_freq(model, mean_removal, vq energy(f) = 10*log10(sum( model(f,3:(L+2)) .^ 2 )); end - melvq; - if length(vq_filename) - x = load(vq_filename); m=5; - vq = x.vq; - size(vq) - end - - rate_K_sample_freqs_kHz = [0.1:0.1:4]; - K = length(rate_K_sample_freqs_kHz); rate_K_surface = resample_const_rate_f(model, rate_K_sample_freqs_kHz, Fs); mean_f = zeros(1,frames); if mean_removal for f=1:frames mean_f(f) = mean(rate_K_surface(f,:)); - rate_K_surface(f,:) -= mean_f(f); + rate_K_surface_no_mean(f,:) = rate_K_surface(f,:) - mean_f(f); end end - if length(vq_filename) - [res rate_K_surface_ ind] = mbest(vq, rate_K_surface, m); - sd_per_frame = std(res'); + % optional vector quantise + + if quant_en + assert(mean_removal == 1); + + weight_gain = 0.1; % I like this vairable name as it is funny + + rate_K_surface_no_mean_ = zeros(frames, K); + res = zeros(frames,vq_cols); ind = []; + + for f=1:frames + target = rate_K_surface_no_mean(f, vq_st:vq_en); + + [diff_weighted weights error g mn_ind] = search_vq_weighted(target, vq, weight_gain); + + rate_K_surface_no_mean_(f,:) = rate_K_surface_no_mean(f,:); + rate_K_surface_no_mean_(f, vq_st:vq_en) = vq(mn_ind,:) + g(mn_ind); + + %res(f,vq_st:vq_en) = rate_K_surface_no_mean(f,vq_st:vq_en) - rate_K_surface_no_mean_(f,vq_st:vq_en); + res(f,vq_st:vq_en) = diff_weighted(mn_ind,:); + ind = [ind mn_ind]; + end + + figure(fg++); clf; mesh(res); + sd_per_frame = std(res(:,vq_st:vq_en)'); figure(fg++); subplot(211); plot(energy); subplot(212); plot(sd_per_frame); - figure(fg); hist(sd_per_frame); + figure(fg++); subplot(211); hist(sd_per_frame); subplot(212); hist(ind,100); printf("VQ rms SD: %3.2f\n", mean(sd_per_frame)); else - rate_K_surface_ = rate_K_surface; + if mean_removal + rate_K_surface_no_mean_ = rate_K_surface_no_mean; + else + rate_K_surface_ = rate_K_surface; + end end if mean_removal for f=1:frames - rate_K_surface_(f,:) += mean_f(f); + rate_K_surface_(f,:) = rate_K_surface_no_mean_(f,:) + mean_f(f); end end @@ -243,7 +299,7 @@ function [model_ rate_K_surface] = experiment_const_freq(model, mean_removal, vq L = model(f,2); AmdB = 20*log10(model(f,3:(L+2))); sd(f) = std(AmdB(1:L) - AmdB_(f,1:L)); - if sd(f) > plot_sd_thresh + if (sd(f) > plot_sd_thresh) && (fg < 10) printf("fg: %d f: %d\n", fg, f); figure(fg++); clf; plot((1:L)*Wo*4/pi, AmdB(1:L),'b+-'); hold on; plot((1:L)*Wo*4/pi, AmdB_(f,1:L),'r+-'); plot(rate_K_sample_freqs_kHz, rate_K_surface_(f,:), 'c+-'); hold off; diff --git a/codec2-dev/octave/newamp1_fbf.m b/codec2-dev/octave/newamp1_fbf.m index d9d6a128..65917c5b 100644 --- a/codec2-dev/octave/newamp1_fbf.m +++ b/codec2-dev/octave/newamp1_fbf.m @@ -83,63 +83,26 @@ function newamp1_fbf(samname, f=10, varargin) axis([1 4000 -20 80]); hold on; plot(rate_K_sample_freqs_kHz*1000, rate_K_vec, ";rate K;b+-"); - %if quant_en != 2 - % plot(rate_K_sample_freqs_kHz*1000, maskdB, ";mask;c+-"); - %end - if quant_en == 2 - plot(rate_K_sample_freqs_kHz*1000, maskdB-rate_K_vec, ";maskdB - rate K;m+-"); - end - #{ - masked_indx = find(maskdB >= rate_K_vec); - unmasked_indx = find(maskdB < rate_K_vec); - masked_indx = (vq_en-5:vq_en); - unmasked_indx = (vq_st:vq_en-6); - stem(rate_K_sample_freqs_kHz(unmasked_indx)*1000, 10*ones(1,length(unmasked_indx)), ";unmasked_index;ko-"); - #} if quant_en - error = g = zeros(1, vq_rows); - - % find mse for each vector - target = rate_K_vec_no_mean(vq_st:vq_en); - weight_gain = 0.1; - %masked_indx = find(maskdB(vq_en-5:vq_en) >= rate_K_vec(vq_en-5:vq_en)); - %masked_indx = find(masked_indx > (vq_cols-5)) - diff = weights = diff_weighted = zeros(vq_rows, vq_cols); - for i=1:vq_rows - - % work out gain for best match + weight_gain = 0.1; % I like this vairable name as it is funny - g(i) = sum(target - vq(i,:))/vq_cols; - - % Find weighted difference. This allocated more importance - % (error) to samples with higher energy, and stops really low - % level harmonics from having any impact. Note addition in dB - % is multiplication in linear - - diff(i,:) = target - vq(i,:) - g(i); - weights(i,:) = max(0.1,weight_gain .* (target + 20)); - - diff_weighted(i,:) = diff(i,:) .* weights(i,:); - - error(i) = sum(abs(diff_weighted(i,:))); - end - - [mn mn_ind] = min(error); + [diff_weighted weights error g mn_ind] = search_vq_weighted(target, vq, weight_gain); rate_K_vec_no_mean_ = rate_K_vec_no_mean; rate_K_vec_no_mean_(vq_st:vq_en) = vq(mn_ind,:) + g(mn_ind); rate_K_vec_ = rate_K_vec_no_mean_ + mean(rate_K_vec); [model_ AmdB_] = resample_rate_L(model(f,:), rate_K_vec_, rate_K_sample_freqs_kHz, Fs); AmdB_ = AmdB_(1:L); - if quant_en == 2 - plot(rate_K_sample_freqs_kHz(vq_st:vq_en)*1000, diff_weighted(mn_ind,:), ";diff;k+-"); - end + + plot(rate_K_sample_freqs_kHz(vq_st:vq_en)*1000, diff_weighted(mn_ind,:), ";diff;k+-"); plot((1:L)*Wo*4000/pi, AmdB_,";AmdB bar;r+-"); hold off; - figure (4); clf; plot(rate_K_sample_freqs_kHz(vq_st:vq_en)*1000, weights(mn_ind,:), ";weights;k+-"); + figure (4); clf; + plot(rate_K_sample_freqs_kHz(vq_st:vq_en)*1000, weights(mn_ind,:), ";weights;k+-"); + axis([0 4000 0 max(weights(mn_ind,:))]); % sort and plot top m matches @@ -152,15 +115,20 @@ function newamp1_fbf(samname, f=10, varargin) for i=1:m subplot(sqrt(m),sqrt(m),i); indx = index_list(i); - plot(target + mean(rate_K_vec),'b+-'); + y_offset = mean(rate_K_vec); + if quant_en == 2 + y_offset = 0; + end + plot(target + y_offset,'b+-'); hold on; if index_list(i) == mn_ind - plot(vq(indx,:) + g(indx) + mean(rate_K_vec),'r+-'); + plot(vq(indx,:) + g(indx) + y_offset,'r+-'); else - plot(vq(indx,:) + g(indx) + mean(rate_K_vec),'g+-'); + plot(vq(indx,:) + g(indx) + y_offset,'g+-'); + end + if quant_en != 2 + plot(diff_weighted(indx,:), "ko-"); end - plot(diff(indx,:), ";diff;k+-"); - plot(diff_weighted(indx,:), ";weighted diff;ko-"); hold off; legend("location", "southwest"); end @@ -198,3 +166,5 @@ function ind = arg_exists(v, str) end end endfunction + + diff --git a/codec2-dev/octave/vq_pager.m b/codec2-dev/octave/vq_pager.m new file mode 100644 index 00000000..f00de9e5 --- /dev/null +++ b/codec2-dev/octave/vq_pager.m @@ -0,0 +1,47 @@ +% vq_pager.m +% +% Interactive Octave script to inspect vectors, e.g. for newamp1.m vector quantisation +% +% Usage: +% +% octave:14> vq_pager(vq) + +function vq_pager(vq) + grid_sz = 4; + + [vq_rows vq_cols] = size(vq); + + r = 1; + + % Keyboard loop + + k = ' '; + do + figure(1); clf; + n_plot = min(grid_sz*grid_sz,vq_rows-r+1) + printf("r: %d n_plot: %d\n", r, n_plot); + for i = 1:n_plot; + subplot(grid_sz,grid_sz,i); + plot(vq(r+i-1,:)); + axis([1 vq_cols -20 20]); + end + + % interactive menu + + printf("\rr: %d menu: n-next b-back q-quit", r); + fflush(stdout); + k = kbhit(); + + if k == 'n' + r = r + grid_sz*grid_sz; + r = min(r, vq_rows); + endif + if k == 'b' + r = r - grid_sz*grid_sz; + r = max(r,1); + endif + until (k == 'q') + printf("\n"); + +endfunction + -- 2.25.1