From ef95c917e22c3063196ac651996871b9d7708a6a Mon Sep 17 00:00:00 2001 From: drowe67 Date: Tue, 13 Jun 2017 21:15:03 +0000 Subject: [PATCH] first pass at quantising dct coeffs using kmeans designed quantisers, doesnt sound too gd git-svn-id: https://svn.code.sf.net/p/freetel/code@3188 01035d8c-6547-0410-b346-abe4f91aad63 --- codec2-dev/octave/newamp.m | 66 +++++++++++++++++++++++++- codec2-dev/octave/newamp1_batch.m | 77 +++++++++++++++++++++++++++---- 2 files changed, 133 insertions(+), 10 deletions(-) diff --git a/codec2-dev/octave/newamp.m b/codec2-dev/octave/newamp.m index b6024723..1af8a4ac 100644 --- a/codec2-dev/octave/newamp.m +++ b/codec2-dev/octave/newamp.m @@ -202,7 +202,7 @@ endfunction % generate a zig-zag linear to square mapping matrix -function map = zigzag(nr,nc) +function map = create_zigzag_map(nr,nc) map = zeros(nr, nc); state = 'zig'; @@ -244,6 +244,8 @@ function map = zigzag(nr,nc) endfunction +% reshape matrix m as a vector v by reading out elements in zig-zag pattern + function v = mtov_zigzag(m) [nr nc] = size(m); map = zigzag(nr,nc) @@ -256,6 +258,68 @@ function v = mtov_zigzag(m) endfunction +% extracts DCT information for rate K surface + +function unwrapped_dcts = dct_blocks(surf) + [frames K] = size(surf); + + % break into 160ms blocks, 2D DCT, truncate, IDCT + + Nt = 16; % 16 blocks in time + Nblocks = floor(frames/Nt); + unwrapped_dcts = zeros(Nblocks,Nt*K); + + for n=1:Nblocks + st = (n-1)*Nt+1; en = st + Nt - 1; + D = dct2(surf(st:en,:)); + unwrapped_dcts(n,:) = reshape(D',1,Nt*K); + end +endfunction + + +% Determines a map for quantising 2D DCT coeffs in order of rms value +% (ie std dev) of each coeff. Those coeffs with the greatest change +% need the most bits to quantise + +function [map rms_map mx unwrapped_dcts] = create_map_rms(rate_K_surface, nr, nc) + unwrapped_dcts = dct_blocks(rate_K_surface); + [mx mx_ind] = sort(std(unwrapped_dcts)); + mx_ind = fliplr(mx_ind); mx = fliplr(mx); + map = rms_map = zeros(nr,nc); + for i=1:nr*nc + r = floor((mx_ind(i)-1)/nc) + 1; + c = mx_ind(i) - (r-1)*nc; + %printf("%d %d %d\n", i, r, c); + map(r,c) = i; + rms_map(r,c) = mx(i); + end +endfunction + + +% Design quantisers for each DCT coeff + +function [quantisers nbits] = design_quantisters(rate_K_surface, nr, nc, nbits_max) + [map rms_map mx unwrapped_dcts] = create_map_rms(rate_K_surface, nr, nc); + nq = nr*nc; + quantisers = zeros(nq, 2^nbits_max); nbits = zeros(nq,1); + for i=1:nq + + % work out number of levels for this quantiser such that it is a + % power of 2 for integer number of bits + + nlevels = (2^nbits_max); + nbits(i) = round(log2(nlevels)); + nlevels = 2 .^ nbits(i); + + if i <= 100 + printf("%d %d\n", i, nlevels); + [idx, centers] = kmeans(unwrapped_dcts(:,i), nlevels); + quantisers(i,1:nlevels) = sort(centers); + end + end +endfunction + + % Determine a phase spectra from a magnitude spectra % from http://www.dsprelated.com/showcode/20.php % Haven't _quite_ figured out how this works but have to start somewhere .... diff --git a/codec2-dev/octave/newamp1_batch.m b/codec2-dev/octave/newamp1_batch.m index 3dbbbc7c..9901f7df 100644 --- a/codec2-dev/octave/newamp1_batch.m +++ b/codec2-dev/octave/newamp1_batch.m @@ -70,6 +70,8 @@ function surface = newamp1_batch(input_prefix, output_prefix) [model_ surface] = experiment_mel_freq(model, 0, 1, voicing); + %model_ = experiment_smoothed(model, 0); + %model_ = experiment_dec_abys(model, 8, 1, 1, 1, voicing); #{ @@ -488,7 +490,6 @@ endfunction - % --------------------------------------------------------------------------------------- % Various experiments tried during development @@ -496,7 +497,11 @@ endfunction function [model_ rate_K_surface] = experiment_mel_freq(model, vq_en=0, plots=1, voicing) [frames nc] = size(model); - K = 20; Fs = 8000; + K = 20; Fs = 8000; correct_rate_K_en = 1; + + quantisers = load("dct2quant.txt"); + nbits = load("dct2quantnbits.txt"); + map = load("dct2map.txt"); for f=1:frames Wo = model(f,1); @@ -505,15 +510,17 @@ function [model_ rate_K_surface] = experiment_mel_freq(model, vq_en=0, plots=1, AmdB = 20*log10(Am); Am_freqs_kHz = (1:L)*Wo*4/pi; [rate_K_vec rate_K_sample_freqs_kHz] = resample_const_rate_f_mel(model(f,:), K); - [tmp_ AmdB_] = resample_rate_L(model(f,:), rate_K_vec, rate_K_sample_freqs_kHz); - [rate_K_vec_corrected orig_error error nasty_error_log nasty_error_m_log] = correct_rate_K_vec(rate_K_vec, rate_K_sample_freqs_kHz, AmdB, AmdB_, K, Wo, L, Fs); - %rate_K_surface(f,:) = rate_K_vec; - rate_K_surface(f,:) = rate_K_vec_corrected; + if correct_rate_K_en + [tmp_ AmdB_] = resample_rate_L(model(f,:), rate_K_vec, rate_K_sample_freqs_kHz); + [rate_K_vec_corrected orig_error error nasty_error_log nasty_error_m_log] = correct_rate_K_vec(rate_K_vec, rate_K_sample_freqs_kHz, AmdB, AmdB_, K, Wo, L, Fs); + rate_K_surface(f,:) = rate_K_vec_corrected; + else + rate_K_surface(f,:) = rate_K_vec; + end end if plots figure(1); clf; mesh(rate_K_surface); - figure(2); clf; plot_dft_surface(rate_K_surface) end for f=1:frames @@ -538,16 +545,50 @@ function [model_ rate_K_surface] = experiment_mel_freq(model, vq_en=0, plots=1, rate_K_surface = rate_K_surface_; end + + % break into 160ms blocks, 2D DCT, truncate, IDCT + + Nt = 16; % 16 blocks in time + Nblocks = floor(frames/Nt); + unwrapped_dcts = zeros(Nblocks,Nt*K); + + for n=1:Nblocks + st = (n-1)*Nt+1; en = st + Nt - 1; + D = dct2(rate_K_surface(st:en,:)); + + % move over surface and work out quantiser + % quantise, replace on map + + E = zeros(Nt,K); + E = D; + for r=1:Nt + for c=1:K + quantiser_num = map(r,c); + nlevels = 2^nbits(quantiser_num); + if quantiser_num < 10 + printf("r %d c %d quantiser_num %d nlevels %d\n", r,c,quantiser_num,nlevels); + levels = quantisers(quantiser_num,1:nlevels); + quant_out = quantise(levels, D(r,c)); + printf("%f %f\n", D(r,c), quant_out ); + E(r,c) = quant_out; + end + end + end + + rate_K_surface(st:en,:) = idct2(E); + end + model_ = resample_rate_L(model, rate_K_surface, rate_K_sample_freqs_kHz); + #{ if plots - figure(3); clf; mesh(rate_K_surface_no_mean); - figure(4); clf; plot_dft_surface(rate_K_surface) + figure(2); clf; mesh(rate_K_surface_no_mean); end for f=1:frames rate_K_surface(f,:) -= mean(rate_K_surface(f,:)); end + #} endfunction @@ -938,6 +979,24 @@ function model_ = experiment_filter_dec_filter(model) endfunction +% smoothing using masking functions + +function model_ = experiment_smoothed(model, bark_model=1) + [frames nc] = size(model); + + model_ = model; + + for f=1:frames + Wo = model(f,1); + L = model(f,2); + Am = model(f,3:(L+2)); + AmdB = 20*log10(Am); + [AmdB_ Am_freqs_kHz] = mask_model(AmdB, Wo, L, bark_model); + model_(f,3:(L+2)) = 10 .^ ((AmdB_)/20); + end + +endfunction + #{ My original idea was to used a 3-4 "resonators" to construct a -- 2.25.1