From: drowe67 Date: Fri, 21 Aug 2015 02:17:50 +0000 (+0000) Subject: wider bands in F1, improved handling of noise and hts2a, refactored into library... X-Git-Url: http://git.whiteaudio.com/gitweb/?a=commitdiff_plain;h=9cdcfada744a80a58dab1a4405acd0d93aecc1af;p=freetel-svn-tracking.git wider bands in F1, improved handling of noise and hts2a, refactored into library plus frame by frame and batch scripts git-svn-id: https://svn.code.sf.net/p/freetel/code@2278 01035d8c-6547-0410-b346-abe4f91aad63 --- diff --git a/codec2-dev/octave/newamp.m b/codec2-dev/octave/newamp.m index 2b36f743..a63cb9ea 100644 --- a/codec2-dev/octave/newamp.m +++ b/codec2-dev/octave/newamp.m @@ -1,11 +1,15 @@ +% newamp.m +% % Copyright David Rowe 2015 % This program is distributed under the terms of the GNU General Public License % Version 2 % -% Octave script to explore new ideas in amplitude modelling. +% Library of Octave functions to explore new ideas in amplitude +% (spectral envelope) modelling. See newamp_fby (frame by frame +% analysis) and newamp_batch (batch processing for listening tests) % % we don't care about -% + spectral tilt, in can vary on input and our quantiser shouldnt care. +% + spectral tilt, in can vary on input and our model shouldnt care. % We can vary it on output and the listener won't care % + absolute energy of entire signal % + harmonics beneath the masking curve @@ -15,100 +19,19 @@ % % TODO: % [ ] waterfall sounds -% [ ] tweak CB bandwidths at LF to be wider +% [X] tweak CB bandwidths at LF to be wider % [ ] way to output at various processing steps, like PF initial mask, pre PF % [ ] BPF at all? % [ ] phase model? Fit LPC, just swing phase at peaks? Try no phase tweaks 1; -% -% frame by frame analysis - -function newamp_frame(samname, f) - - 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); - - plot_all_masks = 0; - k = ' '; - do - figure(1); - clf; - s = [ Sn(2*f-1,:) Sn(2*f,:) ]; - size(s); - plot(s); - axis([1 length(s) -20000 20000]); - - figure(2); - Wo = model(f,1); - L = model(f,2); - Am = model(f,3:(L+2)); - AmdB = 20*log10(Am); - - % plotting - - plot((1:L)*Wo*4000/pi, AmdB,";Am;r"); - axis([1 4000 0 80]); - hold on; - plot((0:255)*4000/256, Sw(f,:),";Sw;"); - - [maskdB Am_freqs_kHz] = mask_model(AmdB, Wo, L); - plot(Am_freqs_kHz*1000, maskdB, 'g'); - - % Analysis by synthesis --------------------------------------- - - mask_sample_freqs_kHz = (1:L)*Wo*4/pi; - [newmaskdB local_maxima] = make_newmask(maskdB, AmdB, Wo, L, mask_sample_freqs_kHz); - - plot(local_maxima(:,2)*Wo*4000/pi, 70*ones(1,length(local_maxima)), 'r+'); - plot(mask_sample_freqs_kHz*1000, newmaskdB, 'm'); - - % 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; - - % interactive menu - - printf("\rframe: %d menu: n-next b-back p-png q-quit m-all", f); - fflush(stdout); - k = kbhit(); - if (k == 'n') - f = f + 1; - endif - if (k == 'b') - f = f - 1; - endif - if k == 'm' - if plot_all_masks == 0 - plot_all_masks = 1; - else - plot_all_masks = 0; - end - end - until (k == 'q') - printf("\n"); - -endfunction +% model mask using just a few samples of critical band filter centre frequencies +function [newmaskdB local_maxima_sort] = make_newmask(maskdB, AmdB, Wo, L, mask_sample_freqs_kHz) -% model mask using just a few samples of critical band filters + % find local maxima and sort in descending order of magnitude -function [newmaskdB local_maxima_sort] = make_newmask(maskdB, AmdB, Wo, L, mask_sample_freqs_kHz) local_maxima = []; if maskdB(1) > maskdB(2) local_maxima = [local_maxima; AmdB(1) 1]; @@ -125,10 +48,11 @@ function [newmaskdB local_maxima_sort] = make_newmask(maskdB, AmdB, Wo, L, mask_ end local_maxima_sort = flipud(sortrows(local_maxima,1)); - % construct new mask from sorted local maxima + % construct new mask from a small subset of sorted local maxima, + % this is the highly compressed model of the amplitude envelope + % that we send to the decoder masker_amps_dB = local_maxima_sort(1:min(4,nlm),1); - % masker_amps_dB = mean(masker_amps_dB)*ones(1,min(4,nlm)); masker_freqs_kHz = local_maxima_sort(1:min(4,nlm),2)*Wo*4/pi; newmaskdB = determine_mask(masker_amps_dB, masker_freqs_kHz, mask_sample_freqs_kHz); endfunction @@ -164,44 +88,6 @@ function [maskdB_pf Am_freqs_kHz peaks_kHz] = mask_model(AmdB, Wo, L) endfunction -% process a whole file and write results - -function newamp_batch(samname) - model_name = strcat(samname,"_model.txt"); - model = load(model_name); - [frames nc] = size(model); - max_amp = 80; - - Am_out_name = sprintf("%s_am.out", samname); - fam = fopen(Am_out_name,"wb"); - - for f=1:frames - - L = min([model(f,2) max_amp-1]); - Wo = model(f,1); - Am = model(f,3:(L+2)); - AmdB = 20*log10(Am); - - maskdB = mask_model(AmdB, Wo, L); - mask_sample_freqs_kHz = (1:L)*Wo*4/pi; - [newmaskdB local_maxima] = make_newmask(maskdB, AmdB, Wo, L, mask_sample_freqs_kHz); - - [nlm tmp] = size(local_maxima); - non_masked_m = local_maxima(1:min(4,nlm),2); - maskdB_pf = newmaskdB; - maskdB_pf(non_masked_m) = maskdB_pf(non_masked_m) + 6; - - Am_ = zeros(1,max_amp); - Am_(2:L) = 10 .^ (maskdB_pf(1:L-1)/20); - - % save to file - fwrite(fam, Am_, "float32"); - end - - fclose(fam); - -endfunction - % % Masking functions from http://www.perceptualentropy.com/coder.html#C % Thanks Jon Boley! @@ -209,13 +95,32 @@ endfunction % Calculate the Schroeder masking spectrum for a given frequency and SPL -function maskdB = schroeder(freq_tone_kHz, mask_sample_freqs_kHz) +function maskdB = schroeder(freq_tone_kHz, mask_sample_freqs_kHz, modified_bark_en=1) f_kHz = mask_sample_freqs_kHz; - A = 3.64*(f_kHz).^(-0.8) - 6.5*exp(-0.6*(f_kHz - 3.3).^2) + (10^(-3))*(f_kHz).^4; f_Hz = f_kHz*1000; % Schroeder Spreading Function - dz = bark(freq_tone_kHz*1000)-bark(f_Hz); + + if modified_bark_en == 1 + + % Modification by DR: Piecewise linear model that makes bands + % beneath 1.5kHz wider to match the width of F1 and + % "fill in" the spectra better for UV sounds. + + if freq_tone_kHz <= 0.5 + y = 0.5; + end + if (freq_tone_kHz > 0.5) && (freq_tone_kHz < 1.5) + y = 0.5*freq_tone_kHz + 0.25; + end + if freq_tone_kHz >= 1.5 + y = 1; + end + dz = y*(bark(freq_tone_kHz*1000) - bark(f_Hz)); + else + dz = bark(freq_tone_kHz*1000)-bark(f_Hz); + end + maskdB = 15.81 + 7.5*(dz+0.474) - 17.5*sqrt(1 + (dz+0.474).^2); endfunction @@ -228,31 +133,33 @@ function b=bark(f) endfunction -% return a sampling grid of frequencies in Hz given B equally spaced -% points on the bark scale - -function f_Hz = bark_freq_samples(B) - Fs2 = 4000; - bark_lut = bark(1:Fs2); - bark_grid_step = 1/B; - bark_grid = (bark_grid_step:bark_grid_step:1)*bark(Fs2); - f_Hz = interp1(bark_lut, 1:Fs2, bark_grid); -endfunction - - -% plot some masking curves +% plot some masking curves, used for working on masking filter changes function plot_masking Fs = 8000; - mask_sample_freqs_kHz = 0.1:0.1:(Fs/1000)/2; - maskdB = schroeder(1, mask_sample_freqs_kHz) + figure(1) + mask_sample_freqs_kHz = 0.1:0.1:(Fs/1000)/2; + maskdB = schroeder(0.5, mask_sample_freqs_kHz, 0); plot(mask_sample_freqs_kHz, maskdB); + hold on; + maskdB = schroeder(0.5, mask_sample_freqs_kHz, 1); + plot(mask_sample_freqs_kHz, maskdB,'g'); + + for f=1:0.5:3 + maskdB = schroeder(f, mask_sample_freqs_kHz, 0); + plot(mask_sample_freqs_kHz, maskdB); + maskdB = schroeder(f, mask_sample_freqs_kHz, 1); + plot(mask_sample_freqs_kHz, maskdB,'g'); + end + hold off; + axis([0.1 4 -30 0]) + + figure(2) + plot(mask_sample_freqs_kHz, bark(mask_sample_freqs_kHz*1000)) + hold on; + plot(mask_sample_freqs_kHz, modified_bark(mask_sample_freqs_kHz*1000),'g') + hold off; + grid endfunction -% Choose one of these to run - -more off; -%newamp_frame("../build_linux/src/hts1a", 121); -newamp_batch("../build_linux/src/hts1a"); -%plot_masking diff --git a/codec2-dev/octave/newamp_batch.m b/codec2-dev/octave/newamp_batch.m new file mode 100644 index 00000000..c23e3a2d --- /dev/null +++ b/codec2-dev/octave/newamp_batch.m @@ -0,0 +1,55 @@ +% newamp_batch.m +% +% Copyright David Rowe 2015 +% This program is distributed under the terms of the GNU General Public License +% Version 2 +% +% Octave script to batch process model parameters using the new +% amplitude model. Used for generating samples we can listen to. +% +% Usage: +% ~/codec2-dev/build_linux/src$ ./c2sim ../../raw/hts1a.raw --dump hts1a +% $ cd ~/codec2-dev/octave +% octave:14> newamp_batch("../build_linux/src/hts1a") +% ~/codec2-dev/build_linux/src$ ./c2sim ../../raw/hts1a.raw --amread hts1a_am.out -o - | play -t raw -r 8000 -s -2 - + +% process a whole file and write results + +function newamp_batch(samname) + newamp; + + model_name = strcat(samname,"_model.txt"); + model = load(model_name); + [frames nc] = size(model); + max_amp = 80; + + Am_out_name = sprintf("%s_am.out", samname); + fam = fopen(Am_out_name,"wb"); + + for f=1:frames + + L = min([model(f,2) max_amp-1]); + Wo = model(f,1); + Am = model(f,3:(L+2)); + AmdB = 20*log10(Am); + + maskdB = mask_model(AmdB, Wo, L); + mask_sample_freqs_kHz = (1:L)*Wo*4/pi; + [newmaskdB local_maxima] = make_newmask(maskdB, AmdB, Wo, L, mask_sample_freqs_kHz); + + [nlm tmp] = size(local_maxima); + non_masked_m = local_maxima(1:min(4,nlm),2); + maskdB_pf = newmaskdB; + maskdB_pf(non_masked_m) = maskdB_pf(non_masked_m) + 6; + + Am_ = zeros(1,max_amp); + Am_(2:L) = 10 .^ (maskdB_pf(1:L-1)/20); + + % save to file + fwrite(fam, Am_, "float32"); + end + + fclose(fam); + +endfunction + diff --git a/codec2-dev/octave/newamp_fbf.m b/codec2-dev/octave/newamp_fbf.m new file mode 100644 index 00000000..2b7ab5a5 --- /dev/null +++ b/codec2-dev/octave/newamp_fbf.m @@ -0,0 +1,96 @@ +% newamp_fbf.m +% +% Copyright David Rowe 2015 +% This program is distributed under the terms of the GNU General Public License +% Version 2 +% +% Interactive Octave script to explore frame by frame operation of new amplitude +% modelling model. +% +% Usage: +% ~/codec2-dev/build_linux/src$ ./c2sim ../../raw/hts1a.raw --dump hts1a +% $ cd ~/codec2-dev/octave +% octave:14> newamp_fbf("../build_linux/src/hts1a",50) + +function newamp_fbf(samname, f) + + newamp; + + 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); + + plot_all_masks = 0; + k = ' '; + do + figure(1); + clf; + s = [ Sn(2*f-1,:) Sn(2*f,:) ]; + size(s); + plot(s); + axis([1 length(s) -20000 20000]); + + figure(2); + Wo = model(f,1); + L = model(f,2); + Am = model(f,3:(L+2)); + AmdB = 20*log10(Am); + + % plotting + + plot((1:L)*Wo*4000/pi, AmdB,";Am;r"); + axis([1 4000 0 80]); + hold on; + plot((0:255)*4000/256, Sw(f,:),";Sw;"); + + [maskdB Am_freqs_kHz] = mask_model(AmdB, Wo, L); + plot(Am_freqs_kHz*1000, maskdB, 'g'); + + % Analysis by synthesis --------------------------------------- + + mask_sample_freqs_kHz = (1:L)*Wo*4/pi; + [newmaskdB local_maxima] = make_newmask(maskdB, AmdB, Wo, L, mask_sample_freqs_kHz); + + plot(local_maxima(:,2)*Wo*4000/pi, 70*ones(1,length(local_maxima)), 'r+'); + plot(mask_sample_freqs_kHz*1000, newmaskdB, 'm'); + + % 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; + + % interactive menu + + printf("\rframe: %d menu: n-next b-back p-png q-quit m-all", f); + fflush(stdout); + k = kbhit(); + if (k == 'n') + f = f + 1; + endif + if (k == 'b') + f = f - 1; + endif + if k == 'm' + if plot_all_masks == 0 + plot_all_masks = 1; + else + plot_all_masks = 0; + end + end + until (k == 'q') + printf("\n"); + +endfunction +