rate_L_sample_freqs_kHz = (1:L)*Wo*Fs/(2000*pi);
- rate_K_surface(f,:) = interp1(rate_L_sample_freqs_kHz, AmdB, rate_K_sample_freqs_kHz, "spline", "extrap");
+ rate_K_surface(f,:) = interp1([0 rate_L_sample_freqs_kHz (Fs/2000)], [AmdB(1) AmdB AmdB(L)], rate_K_sample_freqs_kHz, "spline");
%rate_K_surface(f,:) = interp_para(rate_L_sample_freqs_kHz, AmdB, Fs/(2*1000), rate_K_sample_freqs_kHz);
%printf("%d\n", f);
endfunction
+function [diff_weighted weights error g slope min_ind] = search_vq_weighted(target, vq, weight_gain, fit_order)
+ [vq_rows vq_cols] = size(vq);
+
+ % find mse for each vector
+
+ error = g = slope = zeros(1, vq_rows);
+ diff = weights = diff_weighted = zeros(vq_rows, vq_cols);
+
+ for i=1:vq_rows
+
+ if fit_order == 0
+
+ % work out gain for best match
+
+ g(i) = sum(target - vq(i,:))/vq_cols;
+ slope(i) = 1;
+ end
+
+ if fit_order == 1
+ % work out linear + gradient fit
+ vi = vq(i,:);
+ A = [sum(vi) vq_cols; vi*vi' sum(vi)];
+ c = [sum(target) target*vi']';
+ b = inv(A)*c;
+ g(i) = b(2);
+ slope(i) = b(1);
+ end
+
+ % Find weighted difference. This allocated more importance
+ % (error) to samples with higher energy, and stops really low
+ % level harmonics from having any impact. Note addition in dB
+ % is multiplication in linear
+
+ diff(i,:) = target - slope(i)*vq(i,:) - g(i);
+ weights(i,:) = max(0.1, weight_gain .* (target + 20));
+
+ diff_weighted(i,:) = diff(i,:) .* weights(i,:);
+
+ % abs in dB is MSE in linear
+
+ error(i) = sum(abs(diff_weighted(i,:)));
+ end
+
+ [mn min_ind] = min(error)
+
+endfunction
+
+
% --------------------------------------------------------------------------------
% Experimental functions used for masking, piecewise models, not part of newamp1
% --------------------------------------------------------------------------------
% calculate and plot masking curve
- maskdB = zeros(1,length(mask_sample_freqs_kHz));
+ maskdB = -20*ones(1,length(mask_sample_freqs_kHz));
for m=1:length(masker_freqs_kHz)
maskdB = max(maskdB, schroeder(masker_freqs_kHz(m), mask_sample_freqs_kHz, bark_model) + masker_amps_dB(m));
%maskdB = max(maskdB, parabolic_resonator(masker_freqs_kHz(m), mask_sample_freqs_kHz) + masker_amps_dB(m));