From: drowe67 Date: Tue, 18 Aug 2015 03:04:27 +0000 (+0000) Subject: first pass at new ampl quant mode using masking X-Git-Url: http://git.whiteaudio.com/gitweb/?a=commitdiff_plain;h=b55991bcecb5665fa6b7d0206b57b5dd4caea757;p=freetel-svn-tracking.git first pass at new ampl quant mode using masking git-svn-id: https://svn.code.sf.net/p/freetel/code@2274 01035d8c-6547-0410-b346-abe4f91aad63 --- diff --git a/codec2-dev/octave/melvq.m b/codec2-dev/octave/melvq.m index 4dbc8e03..b31d1cb6 100644 --- a/codec2-dev/octave/melvq.m +++ b/codec2-dev/octave/melvq.m @@ -124,7 +124,7 @@ function [res output_vecs ind] = mbest(vqset, input_vecs, m) endfunction -% Quantises a set of msl-lsps and saves back to disk so they can be read in by c2sim +% Quantises a set of mel-lsps and saves back to disk so they can be read in by c2sim % assumes we have a vq saved to disk called vq % % ~/codec2-dev/build_linux/src$ sox -r 8000 -s -2 ../../wav/vk5qi.wav -t raw -r 8000 -s -2 - sinc 300 sinc -2600 | ./c2sim - --lpc 6 --lpcpf --lspmel --dump vk5qi -o - | play -t raw -r 8000 -s -2 - vol 3 @@ -152,7 +152,7 @@ function ind = test_run(samplename) fclose(fmel_); end -ind = test_run("vk5qi"); +ind = test_run("hts1a"); %load "../build_linux/src/all_mel.txt" %vq = trainvq(all_mel, 64, 3); diff --git a/codec2-dev/octave/newamp.m b/codec2-dev/octave/newamp.m new file mode 100644 index 00000000..4461da5f --- /dev/null +++ b/codec2-dev/octave/newamp.m @@ -0,0 +1,167 @@ +% 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. +% +% we don't care about +% + spectral tilt, in can vary on input and our quantiser shouldnt care. +% We can vary it on output and the listener won't care +% + absolute energy of entire signal +% + harmonics beneath the masking curve +% we do care about: +% + clearly defined formant formation + +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); + + + Ew_on = 1; + 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 mask_sample_freqs_kHz] = determine_mask(AmdB, Wo, L); + non_masked_m = find(maskdB < AmdB); + plot(mask_sample_freqs_kHz*1000, maskdB, 'g'); + plot(non_masked_m*Wo*4000/pi, AmdB(non_masked_m),'b+'); + + hold off; + + % interactive menu + + printf("\rframe: %d menu: n-next b-back p-png q-quit", f); + fflush(stdout); + k = kbhit(); + if (k == 'n') + f = f + 1; + endif + if (k == 'b') + f = f - 1; + endif + + until (k == 'q') + printf("\n"); + +endfunction + + +function [maskdB mask_sample_freqs_kHz] = determine_mask(AmdB, Wo, L) + + % calculate and plot masking curve + + mask_sample_freqs_kHz = (1:L)*Wo*4/pi; + maskdB = zeros(1,L); + for m=1:L + fmasker_kHz = m*Wo*4/pi; + maskdB = max(maskdB, schroeder(fmasker_kHz, mask_sample_freqs_kHz) + AmdB(m)); + end +end + + +% 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); + + % zero out any harmonics beneath mask + + [maskdB mask_sample_freqs_kHz] = determine_mask(AmdB, Wo, L); + non_masked_m = find(maskdB < AmdB); + Am_ = zeros(1,max_amp); + Am_(2:L) = 10 .^ (maskdB(1:L-1)/20); + Am_(non_masked_m+1) *=2; + + % save to file + fwrite(fam, Am_, "float32"); + end + + fclose(fam); + +endfunction + +% +% Masking functions from http://www.perceptualentropy.com/coder.html#C +% Thanks Jon Boley! +% + +% Calculate the Schroeder masking spectrum for a given frequency and SPL + +function maskdB = schroeder(freq_tone_kHz, mask_sample_freqs_kHz) + 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); + maskdB = 15.81 + 7.5*(dz+0.474) - 17.5*sqrt(1 + (dz+0.474).^2); +endfunction + + +% Converts frequency to bark scale +% Frequency should be specified in Hertz + +function b=bark(f) + b = 13*atan(0.76*f/1000) + 3.5*atan((f/7500).^2); +endfunction + + +% plot some masking curves + +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) + plot(mask_sample_freqs_kHz, maskdB); +endfunction + +% Choose one of these to run + +more off; +newamp_frame("../build_linux/src/hts2a", 100); +%newamp_batch("../build_linux/src/kristoff"); +%plot_masking diff --git a/codec2-dev/octave/sd.m b/codec2-dev/octave/sd.m index d3005240..a5352928 100644 --- a/codec2-dev/octave/sd.m +++ b/codec2-dev/octave/sd.m @@ -13,9 +13,9 @@ function sd(raw_filename, dump_file_prefix, f) ak2_filename = sprintf("%s_ak_.txt", dump_file_prefix); ak1 = load(ak1_filename); ak2 = load(ak2_filename); - - [ak1_r, ak1_c] = size(ak1); - [ak2_r, ak2_c] = size(ak1); + + [ak1_r, ak1_c] = size(ak1) + [ak2_r, ak2_c] = size(ak2) frames = max([ak1_r ak2_r]); printf("%d frames\n", frames); sd = zeros(1,frames); @@ -84,7 +84,7 @@ function sd(raw_filename, dump_file_prefix, f) title('Average error across spectrum') mel_indexes_filename = sprintf("%s_mel_indexes.txt", dump_file_prefix); - if file_in_path(".", mel_indexes_filename) + if 0 %file_in_path(".", mel_indexes_filename) mel_indexes = load(mel_indexes_filename); figure(6) bins = [15, 7, 15, 7, 7, 7]; @@ -126,7 +126,7 @@ function sd(raw_filename, dump_file_prefix, f) end printf("\n"); - for l=1:6 + for l=1:10 plot([lsp1(fr,l)*4000/pi lsp1(fr,l)*4000/pi], [0 -10], 'r'); plot([lsp2(fr,l)*4000/pi lsp2(fr,l)*4000/pi], [-10 -20], 'b'); plot([mel(fr,l) mel(fr,l)], [0 10], 'g');