From 1646066985a2f558254b29686cf22ce7babefb0c Mon Sep 17 00:00:00 2001 From: drowe67 Date: Thu, 31 Mar 2016 00:58:02 +0000 Subject: [PATCH] listening tests for new mask/circ/reampling work, sound spretty gd, working on VQ git-svn-id: https://svn.code.sf.net/p/freetel/code@2765 01035d8c-6547-0410-b346-abe4f91aad63 --- codec2-dev/octave/newamp.m | 41 +++++++ codec2-dev/octave/newamp_batch.m | 66 +++--------- codec2-dev/octave/newamp_fbf.m | 177 ++++--------------------------- 3 files changed, 79 insertions(+), 205 deletions(-) diff --git a/codec2-dev/octave/newamp.m b/codec2-dev/octave/newamp.m index 8b61229b..24da35ae 100644 --- a/codec2-dev/octave/newamp.m +++ b/codec2-dev/octave/newamp.m @@ -31,6 +31,47 @@ 1; melvq; + +function [maskdB_ maskdB_cyclic Dabs dk_] = decimate_in_freq(maskdB, cyclic=1, Nkeep=7) + + % Lets try to come up with a smoothed, cyclic model. Replace + % points from 3500 Hz to 4000Hz with a sequence that joins without + % a step to points at the 0Hz end of the spectrum. This will make + % it more cyclical and make the DFT happier, less high freq + % energy. Yes, happier is an extremely technical term. + + L = length(maskdB); + 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_cyclic = [maskdB(1:anchor) ppval(mask_pp, anchor+1:L)]; + + % Now DFT, truncating DFT coeffs to undersample + + if cyclic + D = fft(maskdB_cyclic); + else + D = fft(maskdB); + end + Dabs = abs(D); % this returned for plotting + D_ = D; % + 10*randn(1,L) + 10*j*randn(1,L); + D_(Nkeep+1:L-Nkeep) = 0; % truncate + d = ifft(D_); % back to spectrum at rate L + maskdB_ = real(d); + Dk_ = [0 D(2:Nkeep-1) real(D(Nkeep)) D(L-Nkeep+1:L)]; % build rate Nkeep vector for quantisation + dk_ = real(ifft(Dk_)); + + % OK now fix up last 500Hz, taper down 10dB at 4000Hz + + xpts = [ anchor-1 anchor L]; + ypts = [ maskdB_(anchor-1) maskdB_(anchor) (maskdB_(anchor)-10)]; + mask_pp = splinefit(xpts, ypts, 1); + maskdB_ = [maskdB_(1:anchor) ppval(mask_pp, anchor+1:L)]; + +endfunction + + % Create a "decimated mask" model using just a few samples of % critical band filter centre frequencies. For voiced speech, % we fit the amplitude of these samples to a straight line. diff --git a/codec2-dev/octave/newamp_batch.m b/codec2-dev/octave/newamp_batch.m index b649bafb..796a31d8 100644 --- a/codec2-dev/octave/newamp_batch.m +++ b/codec2-dev/octave/newamp_batch.m @@ -16,21 +16,18 @@ % process a whole file and write results -function [non_masked_f_log non_masked_amp_log] = newamp_batch(samname, optional_Am_out_name, optional_Aw_out_name) +function dk_log = newamp_batch(samname, optional_Am_out_name, optional_Aw_out_name) newamp; more off; max_amp = 80; - decimate_in_freq = 1; + dec_in_freq = 1; postfilter = 1; - decimate_in_time = 0; - synth_phase = 1; + dec_in_time = 0; + synth_phase = 0; freq_quant = 0; amp_quant = 0; - non_masked_f_log = []; - non_masked_m_log = []; - non_masked_amp_log = []; - best_min_mse_log = []; + dk_log = []; model_name = strcat(samname,"_model.txt"); model = load(model_name); @@ -59,6 +56,7 @@ function [non_masked_f_log non_masked_amp_log] = newamp_batch(samname, optional_ %pp_bw = gen_pp_bw; + sd_sum = 0; for f=1:frames printf("%d ", f); model_(f,2) = L = min([model(f,2) max_amp-1]); @@ -72,35 +70,23 @@ function [non_masked_f_log non_masked_amp_log] = newamp_batch(samname, optional_ mask_sample_freqs_kHz = (1:L)*Wo*4/pi; maskdB = mask_model(AmdB, Wo, L); - 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_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 + % save post filter positions for decode - non_masked_amp = masker_amps_dB ; - non_masked_amp_log = [non_masked_amp_log; non_masked_amp]; + maskdB_ = maskdB; + a_non_masked_m = find(AmdB > maskdB); + non_masked_m(f,1:length(a_non_masked_m)) = a_non_masked_m; - % Save this for decoder, so it knows where to apply post filter - % Basically the frequencies of the AbyS samples - - non_masked_m(f,:) = min(round(masker_freqs_kHz/(Wo*4/pi)),L); - non_masked_m_log = [non_masked_m_log; non_masked_m(f,:)]; - - non_masked_f_log = [non_masked_f_log; masker_freqs_kHz]; - maskdB_ = decmaskdB; - - else - % just approximate decimation in freq by using those mask samples we can hear - maskdB_ = maskdB; - a_non_masked_m = find(AmdB > maskdB); - non_masked_m(f,1:length(a_non_masked_m)) = a_non_masked_m; + if dec_in_freq + [maskdB_ tmp1 D dk_] = decimate_in_freq(maskdB, 1); + dk_log = [dk_log; dk_]; end + sd_sum += sum(maskdB - maskdB_); Am_ = zeros(1,max_amp); Am_ = 10 .^ (maskdB_(1:L-1)/20); model_(f,3:(L+1)) = Am_; end + printf("\nsd_sum: %5.2f\n", sd_sum/frames); % decoder loop ----------------------------------------------------- @@ -112,9 +98,8 @@ function [non_masked_f_log non_masked_amp_log] = newamp_batch(samname, optional_ maskdB_ = 20*log10(Am_); mask_sample_freqs_kHz = (1:L)*Wo*4/pi; - %maskdB_ = mask_model(maskdB_, Wo, L); - if decimate_in_time + if dec_in_time % decimate mask samples in time maskdB_ = decimate_frame_rate(model_, 4, f, frames, mask_sample_freqs_kHz); @@ -135,7 +120,6 @@ function [non_masked_f_log non_masked_amp_log] = newamp_batch(samname, optional_ m = max(find(non_masked_m(f,:) > 0)); a_non_masked_m = non_masked_m(f,1:m); - end % post filter - bump up samples by 6dB, reduce mask by same level to normalise gain @@ -171,22 +155,6 @@ function [non_masked_f_log non_masked_amp_log] = newamp_batch(samname, optional_ Aw(1:2:fft_enc*2) = cos(phase); Aw(2:2:fft_enc*2) = -sin(phase); - if 0 - % optional plotting to ensure amd and phase aligned on a few frames - figure(1) - clf; - subplot(211) - plot(mask_sample_freqs_kHz, maskdB_); - hold on; - plot(mask_sample_freqs_kHz, maskdB_pf,'g'); - hold off; - axis([0 4 0 70]) - - subplot(212) - plot((0:(fft_enc/2)-1)*8000/fft_enc, phase(1:fft_enc/2)) - axis([0 4000 -pi pi]) - k = kbhit(); - end fwrite(faw, Aw, "float32"); end end @@ -198,8 +166,6 @@ function [non_masked_f_log non_masked_amp_log] = newamp_batch(samname, optional_ printf("\n"); - printf("mse std: %5.2f\n", std(best_min_mse_log)); - endfunction diff --git a/codec2-dev/octave/newamp_fbf.m b/codec2-dev/octave/newamp_fbf.m index f7808143..17e6fbb4 100644 --- a/codec2-dev/octave/newamp_fbf.m +++ b/codec2-dev/octave/newamp_fbf.m @@ -14,35 +14,25 @@ % $ cd ~/codec2-dev/octave % octave:14> newamp_fbf("../build_linux/src/hts1a",50) + + function newamp_fbf(samname, f=10) - - more off; newamp; - phase_stuff = 0; - plot_not_masked = 0; + more off; plot_spectrum = 1; - freq_quant = 0; - amp_quant = 0; - decimate_in_freq = 0; + + % load up text files dumped from c2sim --------------------------------------- sn_name = strcat(samname,"_sn.txt"); Sn = load(sn_name); - sw_name = strcat(samname,"_sw.txt"); Sw = load(sw_name); - model_name = strcat(samname,"_model.txt"); model = load(model_name); [frames tmp] = size(model); - if phase_stuff - ak_name = strcat(samname,"_ak.txt"); - ak = load(ak_name); - end - - % pp_bw = gen_pp_bw; + % Keyboard loop -------------------------------------------------------------- - plot_all_masks = 0; k = ' '; do figure(1); @@ -59,67 +49,20 @@ function newamp_fbf(samname, f=10) Am = model(f,3:(L+2)); AmdB = 20*log10(Am); - % plotting - axis([1 4000 -20 80]); hold on; if plot_spectrum plot((1:L)*Wo*4000/pi, AmdB,";Am;r"); plot((1:L)*Wo*4000/pi, AmdB,"r+"); - %plot((0:255)*4000/256, Sw(f,:),";Sw;"); end [maskdB Am_freqs_kHz] = mask_model(AmdB, Wo, L); + [tmp1 tmp2 D] = decimate_in_freq(maskdB, 0); + [maskdB_ maskdB_cyclic D_cyclic dk_] = decimate_in_freq(maskdB, 1); - plot(Am_freqs_kHz*1000, maskdB, 'g'); - - % 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. - - 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)]; - - plot(Am_freqs_kHz*1000, maskdB_smooth, 'b'); - - 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 - - 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) - - 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; + plot(Am_freqs_kHz*1000, maskdB, ';mask;g'); + plot(Am_freqs_kHz*1000, maskdB_cyclic, ';mask cyclic;b'); + plot(Am_freqs_kHz*1000, maskdB_, ';mask trunated;c'); % 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 @@ -128,117 +71,41 @@ function newamp_fbf(samname, f=10) figure(3) clf - D = abs(fft(maskdB)); - en = floor(length(D)/2+1); + en = L/2+1; stem(D(2:en),'g') hold on; - D_smooth = abs(fft(maskdB_smooth)); - stem(D_smooth(2:en),'b') + stem(D_cyclic(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'); + plot(cumsum(D(2:en)/sum(D(2:en))),';cumsum;g'); hold on; - plot(cumsum(D_smooth(2:en)/sum(D_smooth(2:en))),'b'); + plot(cumsum(D_cyclic(2:en)/sum(D_cyclic(2:en))),';cumsum cyclic;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 - - %maskdB = decimate_frame_rate(maskdB, model, 4, f, frames, mask_sample_freqs_kHz); - %plot(mask_sample_freqs_kHz*1000, maskdB, 'k'); - - % optionally plot all masking curves - - if plot_all_masks - mask_sample_freqs_kHz = (1:L)*Wo*4/pi; - for m=1:L - maskdB = schroeder(m*Wo*4/pi, mask_sample_freqs_kHz) + AmdB(m); - plot(mask_sample_freqs_kHz*1000, maskdB, "k--"); - end - end - - hold off; - - if phase_stuff - - [phase Sdb s Aw] = determine_phase(model, f, ak(f,:)); - figure(3) - subplot(211) - plot(Sdb) - title('Mag (dB)'); - subplot(212) - plot(phase(1:256)) - hold on; - plot(angle(Aw(1:256))+0.5,'g') - hold off; - title('Phase (rads)'); - figure(4) - plot(s) - end + figure(5) + clf + stem(dk_) - % interactive menu - - printf("\rframe: %d menu: n-next b-back a-Am ", f); - if freq_quant - printf("F"); - else - printf("f"); - end - printf("-freq_quant "); - if amp_quant - printf("M"); - else - printf("m"); - end - printf("-amp_quant q-quit"); + % interactive menu ------------------------------------------ + printf("\rframe: %d menu: n-next b-back q-quit", f); fflush(stdout); k = kbhit(); + if (k == 'n') f = f + 1; endif if (k == 'b') f = f - 1; endif - if k == 'm' - if amp_quant == 0 - amp_quant = 3; - else - amp_quant = 0; - end - end - if k == 'f' - if freq_quant == 0 - freq_quant = 3; - else - freq_quant = 0; - end - end - if k == 'a' - if plot_spectrum == 0 - plot_spectrum = 1; - else - plot_spectrum = 0; - end - end until (k == 'q') printf("\n"); endfunction + -- 2.25.1