--- /dev/null
+% dec_pager.m
+%
+% Interactive Octave script to explore variable decimation points for low rate codecs
+%
+% Usage:
+% Make sure codec2-dev is compiled with the -DDUMP option - see README for
+% instructions.
+% ~/codec2-dev/build_linux/src$ ./c2sim ../../raw/hts1a.raw --dump hts1a
+% $ cd ~/codec2-dev/octave
+% octave:14> dec_pager("../raw/hts1a.raw", "../build_linux/src/hts1a")
+
+function dec_pager(rawfilename, samname, f=1, varargin)
+ more off;
+
+ Fs = 8000; rate_K_sample_freqs_kHz = [0.1:0.1:4]; K = length(rate_K_sample_freqs_kHz);
+ Nsam = 80;
+
+ % Number of frames (time block) we consider for sample points. There are typically
+ % Nt/dec = 4 samples points per block
+
+ Nt = 16;
+
+ % read in optional vector that defines sampling points for
+ % decimation, each entry is the frame number of a sample e.g.
+ % [1 8 10 12 ....
+
+ decsamfile_en = 0;
+ ind = arg_exists(varargin, "decsamfile");
+ if ind
+ dec_filename = varargin{ind+1};
+ decvec = load(dec_filename); Ndv = length(decvec);
+ decsamfile_en = 1;
+ end
+
+ % load up raw samples and text files dumped from c2sim -----------------------
+
+ fs = fopen(rawfilename,"rb");
+ Sraw = fread(fs,Inf,"short");
+
+ 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);
+
+ % Keyboard loop --------------------------------------------------------------
+
+ k = ' ';
+ do
+ % compute rate K surface and mean energy
+
+ rate_K_surf = resample_const_rate_f(model(f:f+Nt-1,:), rate_K_sample_freqs_kHz, Fs);
+ rate_k_mean = mean(rate_K_surf,2);
+
+ offset = 2;
+ st = (f-1-offset)*Nsam + 1; en = st + Nt*Nsam - 1;
+ s = Sraw(st:en);
+ figure(1); clf;
+ subplot(211); plot(st:en, s); axis([st en -20000 20000]);
+ subplot(212); plot(f:(f+Nt-1), rate_k_mean); axis([f f+Nt-1 0 80]);
+ if decsamfile_en
+ hold on; stem(decvec, 60*ones(1,Ndv)); hold off;
+ end
+
+ % plots ----------------------------------
+
+ figure(2); clf;
+ mesh(rate_K_surf);
+ axis([1 K 1 Nt 0 80])
+
+ % interactive menu ------------------------------------------
+
+ printf("\rframe: %d menu: n-next b-back", f);
+ fflush(stdout);
+ k = kbhit();
+
+ if (k == 'n') && (f+2*Nt <= frames)
+ f = f + Nt;
+ endif
+ if k == 'b'
+ f = f - Nt;
+ endif
+ until (k == 'q')
+ 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
+
+
% In general, this function processes a bunch of amplitudes, we then
% use c2sim to hear the results. Bunch of different experiments below
-function [surface mean_f] = newamp1_batch(input_prefix, varargin)
+function [surface_no_mean surface] = newamp1_batch(input_prefix, varargin)
newamp;
more off;
end
for f=1:frames
- surface(f,:) -= mean(surface(f,:));
+ surface_no_mean(f,:) = surface(f,:) - mean(surface(f,:));
end
printf("\n")
vq_start = [];
quant_en = 0;
b_log = [];
+ decimate = 0; decimate_var = 0;
rate_K_sample_freqs_kHz = [0.1:0.1:4];
K = length(rate_K_sample_freqs_kHz);
if ind
avq_filename = varargin{ind+1};
x = load(avq_filename); vq_gain = x.vq;
- quant_en = 1;
+ quant_en = 1
end
- % OK start processing .....
+ quant_en = arg_exists(varargin, "quant");
+
+ ind = arg_exists(varargin, "decimate");
+ if ind
+ decimate = varargin{ind+1};
+ end
+
+ % read in optional vector that defines sampling points for
+ % decimation, each entry is the frame number of a sample e.g.
+ % [1 8 10 12 ....
+
+ ind = arg_exists(varargin, "decsamfile");
+ if ind
+ dec_filename = varargin{ind+1};
+ decvec = load(dec_filename); Ndv = length(decvec);
+ decimate_var = 1;
+ end
+
+ % OK start processing ............................................
energy = zeros(1,frames);
for f=1:frames
b_log = [b_log energy' idx'];
if quant_en
+ k = 1:vq_cols; k2 = k.^2;
+#{
+ [nr nc] = size(vq_gain);
for f=1:frames
v = vq(idx(f),:);
- k = 1:vq_cols; k2 = k.^2;
-#{
- target = rate_K_surface_no_mean(f,vq_st:vq_en);
+ b = b_log(f,:);
+
+ para_target = k2*b(2) + k*b(3) + b(4);
% search vq_gain for best match to gain coefficients
- [nr nc] = size(vq_gain);
d = g = zeros(nr,1);
for r=1:nr
- %con = vq_gain(r,1)*v + vq_gain(r,2)*k2 + vq_gain(r,3)*k;
- con = v + vq_gain(r,2)*k2 + vq_gain(r,3)*k;
- g(r) = (sum(target) - sum(con))/vq_cols;
- con += g(r);
- d(r) = (target-con)*(target-con)';
+ g(r) = (sum(para_target) - sum(vq_gain(r,:)))/vq_cols;
+ diff = para_target - (vq_gain(r,:) + g(r));
+ d(r) = diff*diff';
end
[dmin imin] = min(d);
- b_ = vq_gain(imin,:);
- printf("idx(%d): %d imin: %d b(1); %f b(2): %f b3: %f\n", f, idx(f), imin, b_(1), b_(2), b_(3));
-
+
% recalc contrib
- %contrib(f,:) = b_(1)*v + b_(2)*k2 + b_(3)*k + g(imin);
+ printf("f: %d imin: %d g: %f\n", f, imin, g(imin));
+ contrib(f,:) = b(1)*vq(idx(f),:) + vq_gain(imin,:) + g(imin);
#}
-
- contrib(f,:) = v + (-1/12.5)*b_log(f,3)*k2 + b_log(f,3)*k + b_log(f, 4);
- end
+ for f=1:frames
+ target = rate_K_surface_no_mean(f,vq_st:vq_en);
+ b = b_log(f,:);
+ para_target = k2*b(2) + k*b(3) + b(4);
+ para_target(1) = quantise([-20 -10 0 10], para_target(1));
+ para_target(10) = quantise([-3 +3], para_target(10));
+ b_ = polyfit([k(1) k(10) k(25)],
+ [para_target(1) para_target(10) -10],
+ 2);
+ v = vq(idx(f),:);
+ contrib(f,:) = b(1)*v + b_(1)*k2 + b_(2)*k + b_(3);
+ printf("f: %d pt1: %f pt10: %f b_(1): %f b_(2): %f b_(3): %f \n", f, para_target(1), para_target(10), b_(1), b_(2), b_(3));
+ end
end
end
b_log = [b_log energy' idx'];
end
- rate_K_surface_no_mean_(:, vq_st:vq_en) = contrib;
+ rate_K_surface_no_mean_(:, vq_st:vq_en) = contrib;
res(:, vq_st:vq_en) = rate_K_surface_no_mean(:, vq_st:vq_en) - contrib;
ind(:,i) = idx;
rate_K_surface_(f,:) = rate_K_surface_no_mean_(f,:) + meanf(f);
end
+ % optional decimation
+
+ if decimate
+ for f=1:frames
+
+ % determine frames that bracket the one we need to interp
+
+ left_f = decimate*floor((f-1)/decimate)+1;
+ right_f = left_f + decimate;
+ if right_f > frames
+ right_f = left_f;
+ end
+
+ % determine fraction of each frame to use
+
+ left_fraction = 1 - mod((f-1),decimate)/decimate;
+ right_fraction = 1 - left_fraction;
+
+ rate_K_surface_(f,:) = left_fraction*rate_K_surface_(left_f,:) + right_fraction*rate_K_surface_(right_f,:);
+ printf("f: %d left_f: %d right_f: %d left_fraction: %3.2f right_fraction: %3.2f \n", f, left_f, right_f, left_fraction, right_fraction)
+ end
+ end
+
+ if decimate_var
+ %decvec = 1:4:frames;
+ d = 1; Ndv = length(decvec); right_f = 1;
+ for f=1:frames
+
+ % determine frames that bracket the one we need to interp
+
+ if f == decvec(d)
+ left_f = right_f;
+ if d < Ndv
+ d++;
+ right_f = decvec(d);
+ decimate = decvec(d) - decvec(d-1);
+ end
+ end
+
+ % determine fraction of each frame to use
+
+ left_fraction = 1 - (f - left_f)/decimate;
+ right_fraction = 1 - left_fraction;
+
+ rate_K_surface_(f,:) = left_fraction*rate_K_surface_(left_f,:) + right_fraction*rate_K_surface_(right_f,:);
+ printf("f: %d dec: %d lf: %d rf: %d lfrac: %3.2f rfrac: %3.2f \n", f, decimate, left_f, right_f, left_fraction, right_fraction)
+ end
+ end
+
[model_ AmdB_] = resample_rate_L(model, rate_K_surface_, rate_K_sample_freqs_kHz, Fs);
% Measure distortion between AmdB and AmdB_, ths includes distortion