From c8203bf2c7712df79eddf93d9a3e7e24bedc4061 Mon Sep 17 00:00:00 2001 From: drowe67 Date: Thu, 13 Jul 2017 04:37:37 +0000 Subject: [PATCH] first pass. Works, but has some high freq artefacts, probably the interpolator losing it git-svn-id: https://svn.code.sf.net/p/freetel/code@3299 01035d8c-6547-0410-b346-abe4f91aad63 --- codec2-dev/octave/c2wideband_batch.m | 302 +++++++++++++++++++++++++++ 1 file changed, 302 insertions(+) create mode 100644 codec2-dev/octave/c2wideband_batch.m diff --git a/codec2-dev/octave/c2wideband_batch.m b/codec2-dev/octave/c2wideband_batch.m new file mode 100644 index 00000000..c6031519 --- /dev/null +++ b/codec2-dev/octave/c2wideband_batch.m @@ -0,0 +1,302 @@ +% c2wideband_batch.m +% +% Copyright David Rowe 2017 +% This program is distributed under the terms of the GNU General Public License +% Version 2 + +#{ + + Octave script to batch process model parameters for wideband Codec 2. + + Outputs a set of model parameters that can be fed to c2sim for + listening tests. The companion c2wideband_fbf.m script is used to + visualise the processing frame by frame + + c2sim -> dump files -> c2wideband_batch.m -> output model params -> c2sim -> play + + Usage: + + Build codec2 with -DDUMP (see codec2-dev/README), then generate dump files: + + ~/codec2-dev/build_linux/src$ ./c2sim ~/Desktop/c2_hd/speech_orig_16k.wav --Fs 16000 --dump speech + + Start Octave and generate the map file, this only needs to be done once: + + $ cd ~/codec2-dev/octave + $ octave + octave:1> c2wideband_batch("../build_linux/src/speech", "mode", "generate map") + + Then to run bathc simulation and generate output speech: + + octave:1> c2wideband_batch("../build_linux/src/speech"); + + ~/codec2-dev/build_linux/src$ ./c2sim ~/Desktop/c2_hd/speech_orig_16k.wav --Fs 16000 --phase0 --postfilter --amread speech_am.out --hmread speech_hm.out -o | play -t raw -r 16000 -s -2 - +#} + +#{ + TODO: + [ ] how to generate dct coeff readout map for wideband + + maybe sep function + + save map as disk file for sucessive runs + + document process, special mode +#} + +function [surface mean_f] = c2wideband_batch(input_prefix, varargin) + newamp; + more off; + + max_amp = 160; + mean_f = []; + + % defaults + + synth_phase = output = 1; + output_prefix = input_prefix; + fit_order = 0; + mode = "dct2"; + + % parse variable argument list + + if (length (varargin) > 0) + + ind = arg_exists(varargin, "mode"); + if ind + mode = varargin{ind+1}; + end + + % check for the "output_prefix" option + + ind = arg_exists(varargin, "output_prefix"); + if ind + output_prefix = varargin{ind+1}; + end + + ind = arg_exists(varargin, "no_output"); + if ind + output = 0; + synth_phase = 0; + end + end + + if output && !strcmp(mode,"generate map") + printf("output_prefix: %s\n", output_prefix); + end + + model_name = strcat(input_prefix,"_model.txt"); + model = load(model_name); + [frames nc] = size(model); + + % Choose experiment to run test here ----------------------- + + if strcmp(mode, "generate map") + generate_map(model, K=30, "c2wideband_map"); + output = 0; + end + if strcmp(mode, "dct2") + [model_ surface] = experiment_rate_K_dct2(model, 1); + end + + % ---------------------------------------------------- + + if output + Am_out_name = sprintf("%s_am.out", output_prefix); + fam = fopen(Am_out_name,"wb"); + + 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)); + 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"); + + if synth_phase + + % 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); + + % 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)); + end + fwrite(fhm, Hm, "float32"); + end + end + + fclose(fam); + if synth_phase + fclose(fhm); + end + end % if output ..... + + printf("\n") + +endfunction + + +function ind = arg_exists(v, str) + ind = 0; + for i=1:length(v) + if strcmp(v{i}, str) + ind = i; + end + end +endfunction + + +% Create "map" from a training database. The map tells us which order +% to read out and quantise DCT coeffs. Could be approximated by a +% zig-zag pattern. +% +% Example 3x3 map in a zig-zag pattern, so (1,1) is first coeff, (1,2) +% 2nd .... +% +% 1 2 6 +% 2 5 7 +% 4 8 9 +% +% TODO: [ ] Come up with a better name than +% [ ] Script to convert this to a C header file. + +function generate_map(model, K, map_filename) + newamp; + + Fs = 16000; dec = 2; Nt = 16; + + [frames nc] = size(model); + surface = resample_const_rate_f_mel(model, K, Fs); + [nr nc] = size(surface); + asurf = surface(1:dec:nr,:); + [map rms_map mx mx_ind unwrapped_dcts] = create_map_rms(asurf, Nt, K); + + %printf("non zero coeffs: %d\n", sum(sum(map == 1))); + figure(2); clf; + mesh(map); + + printf("generating map file: %s\n", map_filename); + save("-ascii", map_filename, "map"); +endfunction + + +% --------------------------------------------------------------------------------------- +% rate K mel-resampling, high end correction, and DCT experiment workhorse + +function [model_ rate_K_surface] = experiment_rate_K_dct2(model, plots=1) + newamp; + [frames nc] = size(model); + K = 30; Fs = 16000; correct_rate_K_en = 1; + + % Resample variable rate L vectors to fixed length rate K. We have + % left high end correction out for now, this is less of an issue + % with a higher K + + [rate_K_surface rate_K_sample_freqs_kHz] = resample_const_rate_f_mel(model, K, Fs); + + % break into 160ms blocks, 2D DCT, truncate, IDCT + + Tf = 0.01; % frame period in seconds + Nt = 16; % number of 10ms frames blocks in time + dec = 2; % decimation factor + dist_dB = 2; % use enough coefficients to get this distortion ond DCT coeffs + + Nblocks = floor(frames/(Nt*dec)); + printf("frames: %d Nblocks: %d\n", frames, Nblocks); + + % map that defines order we read out and quantise DCT coeffs + + map = load("c2wideband_map"); + + % init a bunch of output variables + + rate_K_surface_ = zeros(frames, K); + sumnz = zeros(1,Nblocks); + dct2_sd = zeros(1,Nblocks); + + % create arrays to reverse map quantiser_num to r,c Luts + + rmap = cmap = zeros(1,Nt*K); + for r=1:Nt + for c=1:K + quantiser_num = map(r,c); + rmap(quantiser_num) = r; + cmap(quantiser_num) = c; + end + end + + for n=1:Nblocks + st = (n-1)*dec*Nt+1; en = st + dec*Nt - 1; + %printf("st: %d en: %d\n", st, en); + D = dct2(rate_K_surface(st:dec:en,:)); + + % So D is the 2D block of DCT coeffs at the encoder. We want to + % create a quantised version at the "decoder" E. This loop copies + % DCTs coeffs from D to E, until we get beneath a distortion + % threshold. + + % This is essentially variable rate quantisation, but gives us + % some idea of the final bit rate. In practice we will also need + % to constrain the total number of bits (ie bit rate), and + % quantise each coefficient. + + % Turns out than mean SD (across many blocks/frames) is about the + % same in the DCT domain as the rate K domain. So we can just + % measure MSE between D and E to estimate mean SD of the rate K + % vectors after quantisation. + + E = mapped = zeros(Nt,K); + + qn = 0; + adct2_sd = mean(std(D-E)); + while adct2_sd > dist_dB + qn++; + E(rmap(qn), cmap(qn)) = 1*round(D(rmap(qn), cmap(qn))/1); + adct2_sd = mean(std(D-E)); + %printf("qn %d %f\n", qn, adct2_sd); + end + sumnz(n) = qn; + + % note neat trick to interpolate to 10ms frames despite dec to 20ms, this means + % we don't need a separate decode side interpolator. + + rate_K_surface_(st:en,:) = idct2([sqrt(dec)*E; zeros(Nt*(dec-1), K)]); + + dct2_sd(n) = mean(std(D-E)); + end + + if plots + figure(4); clf; plot(sumnz); hold on; + plot([1 length(sumnz)],[mean(sumnz) mean(sumnz)]); hold off; title('Non Zero'); + figure(5); clf; plot(dct2_sd); title('DCT SD'); + end + printf("average dct spectral distortion: %3.2f dB\n", mean(dct2_sd)); + printf("mean number of coeffs/DCT: %3.2f/%d\n", mean(sumnz), Nt*K); + printf("coeffs/second: %3.2f\n", mean(sumnz)/(Nt*Tf*dec)); + printf("bits/s: %3.2f\n", 2.9*mean(sumnz)/(Nt*Tf*dec)); + + % prevent /0 errors at end of run + + rate_K_surface_(dec*Nt*Nblocks+1:frames,:) = rate_K_surface(dec*Nt*Nblocks+1:frames,:); + model_ = resample_rate_L(model, rate_K_surface_, rate_K_sample_freqs_kHz, Fs); + + dist = std((rate_K_surface_(1:dec:frames,:) - rate_K_surface(1:dec:frames,:))'); + + if plots + figure(1); clf; plot(dist); title('Rate K SD'); + printf("Rate K spectral distortion mean: %3.2f dB var: %3.2f\n", mean(dist), var(dist)); + end +endfunction + -- 2.25.1