Fs = 8000; rate_K_sample_freqs_kHz = [0.1:0.1:4]; K = length(rate_K_sample_freqs_kHz);
quant_en = 0; vq_search = "gain";
mask_en = 0;
- nvec = 0;
+ nvq = 0; vq_start = [];
quant_en = weight_en = 0;
- % optional (split) VQ of rate K samples
-
- ind = arg_exists(varargin, "vq");
- if ind
- nvec++;
- vq_filename = varargin{ind+2};
- x = load(vq_filename); vq = x.vq;
- [vq_rows vq_cols] = size(vq); vq_st = varargin{ind+1}; vq_en = vq_st + vq_cols - 1;
- end
-
- % optional VQ of VQ gain coefficients
-
- ind = arg_exists(varargin, "vq_gain");
- if ind
- nvec++;
- vq_filename = varargin{ind+1};
- x = load(vq_filename); vq_gain = x.vq;
+ % specify one of more vqs, start index, vq file name, and search method
+
+ ind = anind = arg_exists(varargin, "vq");
+ while anind
+ nvq++;
+ vq_start = [vq_start varargin{ind+1}];
+ avq_filename = varargin{ind+2};
+ if nvq < 2; vq_filename = avq_filename; else; vq_filename = [vq_filename; avq_filename]; end;
+ avq_search = varargin{ind+3};
+ if nvq < 2; vq_search = avq_search; else; vq_search = [vq_search; avq_search]; end;
+ printf("nvq %d vq_start: %d vq_filename: %s vq_search: %s\n", nvq, vq_start(nvq), avq_filename, avq_search);
+ anind = arg_exists(varargin(ind+1:length(varargin)), "vq");
+ if anind
+ ind += anind;
+ end
end
- % different vq search algorithms
-
- ind = arg_exists(varargin, "vq_search");
- if ind
- vq_search = varargin{ind+1};
- end
-
fit_order = 0;
- ind = arg_exists(varargin, "noslope");
- if ind
- fit_order = 1;
- end
+ ind = arg_exists(varargin, "construct");
+ vq_search = varargin{ind};
+
% optional exploration of phase
ind = arg_exists(varargin, "phase");
rate_K_vec_ = rate_K_vec_fit;
- if mask_en && nvec
+ if mask_en && nvq
% experimental masking stuff that I can't seem to get to work
maskdB = determine_mask(rate_K_vec, rate_K_sample_freqs_kHz, rate_K_sample_freqs_kHz, bark_model=1);
plot(rate_K_sample_freqs_kHz*1000, maskdB, ";mask dB;c+-");
end
- if nvec
- if mask_en
- target = rate_K_vec;
- mask_thresh = 3;
- ind = find (maskdB - target > mask_thresh);
- target(ind) = maskdB(ind) - mask_thresh;
- plot(rate_K_sample_freqs_kHz*1000, target, ";target;m+-");
- target = target(vq_st:vq_en) - b;
- else
- target = rate_K_vec_fit(vq_st:vq_en);
- end
-
- weights = ones(1, vq_en - vq_st + 1);
- if weight_en
-
- % generate weighting. if min is 20dB and max 40dB, weight at min
- % is 1 and max 2, same for -10 to 10. So gradient = 0.05, we only
- % haveto calculate y intercept
-
- gradient = 0.05; yint = 1 - gradient*min(target);
- weights = gradient*target + yint;
- end
+ if strcmp(vq_search, "construct")
+ target = rate_K_vec_fit;
+ [idx contrib errors b] = vq_construct_mg(target);
+ rate_K_vec_ = contrib;
end
-
- if nvec
-
- if strcmp(vq_search, "mse")
- [idx contrib errors test_ g mg sl] = vq_search_mse(vq, target);
- rate_K_surface_fit_(f, vq_st:vq_en) = contrib;
- end
- if strcmp(vq_search, "gain")
- [idx contrib errors b] = vq_search_gain(vq, target, weights);
- end
-
- if strcmp(vq_search, "max")
- [idx contrib errors b] = vq_search_max(vq, target);
- end
+ if nvq
- if strcmp(vq_search, "sg")
- [idx contrib errors b] = vq_search_sg(vq, target);
- end
+ for i=1:nvq
+ avq_filename = char(cellstr(vq_filename)(i));
+ avq_search = char(cellstr(vq_search)(i));
+ x = load(avq_filename); vq = x.vq; [vq_rows vq_cols] = size(vq);
+ vq_st = vq_start(i); vq_en = vq_st + vq_cols - 1;
+ printf("\nsplit VQ: %d vq_filename: %s vq_search: %s vq_st: %d vq_en: %d nVec: %d\n", i, avq_filename, avq_search, vq_st, vq_en, vq_rows);
+ target = rate_K_vec_fit(vq_st:vq_en);
- if strcmp(vq_search, "mg")
- [idx contrib errors b] = vq_search_mg(vq, target);
- end
+ if strcmp(avq_search, "mse")
+ [idx contrib errors test_ g mg sl] = vq_search_mse(vq, target);
+ rate_K_surface_fit_(f, vq_st:vq_en) = contrib;
+ end
- if strcmp(vq_search, "slope")
- [idx contrib errors b_log] = vq_search_slope(vq, target, "closed_quant_slope");
- rate_K_surface_fit_(f, vq_st:vq_en) = contrib;
- printf(" mg: %3.2f sl: %3.2f g: %3.2f \n", b_log(1), b_log(2), b_log(3));
- if quant_en
- % set slope to 0
- contrib1 = contrib;
- contrib = b_log(1)*vq(idx,:) + b_log(3);
- rate_K_vec_(vq_en+1:K) -= b_log(2)*vq_cols;
+ if strcmp(avq_search, "gain")
+ [idx contrib errors b] = vq_search_gain(vq, target, weights);
end
- end
- if strcmp(vq_search, "para")
- printf("\n");
- [idx contrib errors b] = vq_search_para(vq, target);
+ if strcmp(avq_search, "max")
+ [idx contrib errors b] = vq_search_max(vq, target);
+ end
- k = 1:vq_cols; k2 = k.^2;
- para_target = k2*b(2) + k*b(3) + b(4);
- samples = [1 10 25];
- if quant_en
+ if strcmp(avq_search, "sg")
+ [idx contrib errors b] = vq_search_sg(vq, target);
+ end
-#{
- % search vq_gain for best match to gain coefficients
+ if strcmp(avq_search, "mg")
+ [idx contrib errors b] = vq_search_mg(vq, target);
+ end
- [nr nc] = size(vq_gain);
- d = g = zeros(nr,1);
- for r=1:nr
- g(r) = (sum(para_target) - sum(vq_gain(r,:)))/vq_cols;
- diff = para_target - (vq_gain(r,:) + g(r));
- d(r) = diff*diff';
+ if strcmp(avq_search, "slope")
+ [idx contrib errors b_log] = vq_search_slope(vq, target, "closed_quant_slope");
+ rate_K_surface_fit_(f, vq_st:vq_en) = contrib;
+ printf(" mg: %3.2f sl: %3.2f g: %3.2f \n", b_log(1), b_log(2), b_log(3));
+ if quant_en
+ % set slope to 0
+ contrib1 = contrib;
+ contrib = b_log(1)*vq(idx,:) + b_log(3);
+ rate_K_vec_(vq_en+1:K) -= b_log(2)*vq_cols;
end
- [dmin imin] = min(d);
-
-#}
- v = vq(idx,:);
-#{
- rng = 1:vq_cols-5;
- g = (sum(target(rng)) - sum(b(1)*v(rng)))/(vq_cols-5);
- v(vq_cols-5:vq_cols) += -10*(1:6);
- printf("g: %f\n", g);
- % recalc contrib
-#}
- para_target(1) = quantise([-20 -10 0 10], para_target(1));
- %para_target(10) = quantise([-6 +6], para_target(10));
- %para_target(10)
- b_ = polyfit([k(1) k(10) k(25)],
- [para_target(1) para_target(10) -10],
- 2);
-
- contrib1 = contrib;
- contrib = b(1)*v + b_(1)*k2 + b_(2)*k + b_(3);
- para = b_(1)*k2 + b_(2)*k + b_(3);
- %printf("imin: %d\n", imin);
end
- rate_K_surface_fit_(f, vq_st:vq_en) = contrib;
-
- end
-
- if strcmp(vq_search, "cubic")
- printf("\n");
- [idx contrib errors b] = vq_search_cubic(target);
- rate_K_surface_fit_(f, vq_st:vq_en) = contrib;
+ rate_K_vec_(vq_st:vq_en) = contrib;
end
-
- if strcmp(vq_search, "fourth")
- printf("\n");
- [idx contrib errors b] = vq_search_fourth(target);
- rate_K_surface_fit_(f, vq_st:vq_en) = contrib;
- % printf("g: %3.2f mg: %3.2f sl: %3.2f\n", g(idx), mg(idx), sl(idx));
- end
-
- rate_K_vec_(vq_st:vq_en) = contrib;
-
+
plot(rate_K_sample_freqs_kHz(vq_st:vq_en)*1000, contrib, 'm+-');
if strcmp(vq_search, "para")
plot(rate_K_sample_freqs_kHz(vq_st:vq_en)*1000, para_target, 'c+-');
sdL = std(abs(AmdB - AmdB_));
plot((1:L)*Wo*4000/pi, AmdB_,";AmdB bar;r+-");
- if nvec == 0
- l = sprintf(";errorx10 sd %3.2f dB;bk+-", sdL);
- plot((1:L)*Wo*4000/pi, 10*(AmdB - AmdB_), l);
+ if nvq == 0
+ l = sprintf(";error sd %3.2f dB;bk+-", sdL);
+ plot((1:L)*Wo*4000/pi, (AmdB - AmdB_), l);
end
hold off;
function ind = arg_exists(v, str)
ind = 0;
for i=1:length(v)
- if strcmp(v{i}, str)
+ if !ind && strcmp(v{i}, str)
ind = i;
end
end
--- /dev/null
+%----------------------------------------------------------------------
+% Construct a VQ entry using formant prototypes and gain/mag fitting
+
+#{
+ ptable - prototype table
+ - each row defines a formant prototype with 3 numbers: shape start end
+ - shape 1 - gaussian
+ 2 - symmetrical triangle
+ 3 - triangle with deep fall off at HF side
+ - start 0 - one after last prototype center position
+ n - start of search range in 100's of Hz, e.g. 20 -> 2000Hz
+ - end m - end of search range in 100's of Hz
+#}
+
+function [idx contrib errors b_log2] = vq_construct_mg(data)
+
+ ptable = [ 1 4 8 1 10;
+ 2 8 30 0 30;
+ 2 12 30 0 30;
+ 3 30 36 0 40];
+
+ protos = [-4 12 23 24 18 12 8 5 2;
+ 12 18 12 0 0 0 0 0 0;
+ 12 18 12 -12 -24 -36 0 0 0];
+ protoc = [4 2 2];
+ protol = [9 3 6];
+
+ nformants = rows(ptable);
+ [nRows nCols] = size(data);
+ K = nCols;
+
+ idx = zeros(nRows, nformants); b_log2 = [];
+
+ contrib = zeros(nRows, K);
+
+ #{
+ range of F1, specify range of each formant
+ specify shape
+
+ [X] table of shapes
+ [X] work out centre position
+ [X] another table to refer to each shape based on formant number
+ [X] const for number of formants
+ [ ] iterate over positions
+ [ ] gain/mag match at each position, MSE over entire vector
+ [ ] way to graphically peer inside what's going on
+ [ ] how to model flat in first 1kHz?
+ + cld perhaps exclude search entirely as an option
+ + or perhap sit will trun up low mags
+ #}
+
+ for r=1:nRows
+ % iterate across all formants
+
+ for f=1:nformants
+
+ % try formant centre at np points between st and en
+
+ shape = ptable(f,1); c_st = ptable(f,2); c_en = ptable(f,3);
+ if c_st == 0 % 0 code means start just after previous formant
+ c_st = idx(r,f-1) + 2;
+ end
+ np = c_en-c_st+1;
+ b_log = zeros(np, 2);
+ error = zeros(np, 1);
+
+ % we are matching a target segment between st and en
+
+ st = ptable(f,4); en = ptable(f,5);
+ if st == 0
+ st = c_st;
+ end
+ lt = en-st+1;
+ t = data(r,st:en);
+
+ printf("r: %d f: %d c_st: %d c_en: %d np: %d st: %d en: %d ------------------------\n",
+ r, f, c_st, c_en, np, st, en);
+
+ % test range of centres for formant
+
+ for c=c_st:c_en
+ v = zeros(1, K);
+ v_st = c - protoc(shape) + 1; v_en = c - protoc(shape) + protol(shape);
+ v(v_st:v_en) = protos(shape, 1:protol(shape));
+ v = v(st:en);
+ A = [v*v' sum(v); sum(v) lt];
+ d = [t*v' sum(t)]';
+ b = inv(A)*d;
+ b(1) = max(0.5,b(1));
+ diff = t - (b(1)*v + b(2));
+ p = c - c_st + 1;
+ b_log(p,:) = b;
+ error(p) = diff * diff';
+
+ %printf("r: %d f: %d c: %d e: %f b(1): %f b(2): %f \n", r, f, c, error(p), b(1), b(2));
+ end
+
+ % choose best for this formant
+
+ [mn p_min] = min(error);
+ c_min = p_min + c_st - 1
+ b = b_log(p_min,:);
+ v = zeros(1, K);
+ v_st = c_min - protoc(shape) + 1; v_en = c_min - protoc(shape) + protol(shape);
+ v(v_st:v_en) = protos(shape, 1:protol(shape));
+ contrib(r,st:en) = b(1)*v(st:en) + b(2);
+
+ % log index, mag, gain, might need a matric with one row per formant
+
+ idx(r,f) = c_min;
+ b_log2 = [b_log2; b];
+ end
+
+ errors(r) = sum((data(r,:) - contrib(r,:)) .^ 2);
+ end
+
+endfunction