Tools: Tune: dcblock: Create set of cut-off frequencies

This patch modifies script example_dblock.m to create in
addition to existing default coefficients blob with parameter
R=0.98 a set of blobs those provide (-3 dB) cut-off frequencies
20, 30, 40, 50, 100, 200 Hz for 16 kHz and 48 kHz for various
usages. Human understandable parameters for blobs helps to select
the configuration for the needed dcblock usage.

The parameter value R for given frequencies is calculated
with iterative function dcblock_rval_calculate(). A closed
form equation might be possible to derive from the transfer
function. If such is found this function can be replaced with
quicker equation -- or with more advanced faster converging
iteration.

The topology blob export functions are also modified to add
comment line for exact build command.

Signed-off-by: Seppo Ingalsuo <[email protected]>

tools/tune/common/tplg2_write.m

@@ -1,4 +1,4 @@
-function tplg2_write(fn, blob8, component, comment)
+function tplg2_write(fn, blob8, component, comment, howto)
 
 % SPDX-License-Identifier: BSD-3-Clause
 %
@@ -9,6 +9,9 @@ function tplg2_write(fn, blob8, component, comment)
 if nargin < 4
 	comment = 'Exported Control Bytes';
 end
+if nargin < 5
+	howto = [];
+end
 
 %% Check that blob length is multiple of 32 bits
 n_blob = length(blob8);
@@ -23,6 +26,9 @@ function tplg2_write(fn, blob8, component, comment)
 fh = fopen(fn, 'w');
 nl = 8;
 fprintf(fh, '# %s %s\n', comment, date());
+if ~isempty(howto)
+	fprintf(fh, '# %s\n', howto);
+end
 fprintf(fh, 'Object.Base.data.\"%s\" {\n', component);
 fprintf(fh, '\tbytes \"\n');
 for i = 1:nl:n_blob

tools/tune/common/tplg_write.m

@@ -1,8 +1,11 @@
-function tplg_write(fn, blob8, name, comment)
+function tplg_write(fn, blob8, name, comment, howto)
 
 if nargin < 4
 	comment = 'Exported Control Bytes';
 end
+if nargin < 5
+	howto = [];
+end
 
 %% Pad blob length to multiple of four bytes
 n_orig = length(blob8);
@@ -14,6 +17,9 @@ function tplg_write(fn, blob8, name, comment)
 fh = fopen(fn, 'w');
 nl = 8;
 fprintf(fh, '# %s %s\n', comment, date());
+if ~isempty(howto)
+	fprintf(fh, '# %s\n', howto);
+end
 fprintf(fh, 'CONTROLBYTES_PRIV(%s_priv,\n', name);
 fprintf(fh, '`       bytes "');
 for i = 1:nl:n_new
@@ -33,5 +39,6 @@ function tplg_write(fn, blob8, name, comment)
 end
 fprintf(fh, ')\n');
 fclose(fh);
+fprintf('Blob size %d was written to file %s\n', n_new, fn);
 
 end

tools/tune/dcblock/example_dcblock.m

@@ -1,28 +1,61 @@
 function example_dcblock()
 
+% Default blob, about 150 Hz cut-off @ 48 kHz
+prm.fc = [];
+prm.fs = [];
+prm.R_coeffs = [0.98, 0.98, 0.98, 0.98, 0.98, 0.98, 0.98, 0.98];
+prm.id = "default";
+dcblock_blob_calculate(prm)
+
+% Generate a set of cut-off frequencies for 16 and 48 kHz rate
+prm.R_coeffs = [];
+for fs = [16e3 48e3]
+	for fc = [20 30 40 50 80 100 200]
+		prm.id = sprintf("%dhz_%dkhz", round(fc), round(fs / 1000));
+		prm.fs = fs;
+		prm.fc = fc;
+		dcblock_blob_calculate(prm)
+	end
+end
+
+end
+
+function dcblock_blob_calculate(prm)
+
 % Set the parameters here
-tplg1_fn = "../../topology/topology1/m4/dcblock_coef_default.m4"; % Control Bytes File
-tplg2_fn = "../../topology/topology2/include/components/dcblock/default.conf";
+tplg1_fn = sprintf("../../topology/topology1/m4/dcblock_coef_%s.m4", prm.id); % Control Bytes File
+tplg2_fn = sprintf("../../topology/topology2/include/components/dcblock/%s.conf", prm.id);
 % Use those files with sof-ctl to update the component's configuration
-blob3_fn = "../../ctl/ipc3/dcblock_coef.blob"; % Blob binary file
-alsa3_fn = "../../ctl/ipc3/dcblock_coef.txt"; % ALSA CSV format file
-blob4_fn = "../../ctl/ipc4/dcblock_coef.blob"; % Blob binary file
-alsa4_fn = "../../ctl/ipc4/dcblock_coef.txt"; % ALSA CSV format file
+blob3_fn = sprintf("../../ctl/ipc3/dcblock/coef_%s.blob", prm.id); % Blob binary file
+alsa3_fn = sprintf("../../ctl/ipc3/dcblock/coef_%s.txt", prm.id); % ALSA CSV format file
+blob4_fn = sprintf("../../ctl/ipc4/dcblock/coef_%s.blob", prm.id); % Blob binary file
+alsa4_fn = sprintf("../../ctl/ipc4/dcblock/coef_%s.txt", prm.id); % ALSA CSV format file
 
 endian = "little";
-R_coeffs = [0.98, 0.98, 0.98, 0.98, 0.98, 0.98, 0.98, 0.98];
+
+if isempty(prm.fc)
+	R_coeffs = prm.R_coeffs;
+else
+	channels = 8;
+	R = dcblock_rval_calculate(prm.fs, prm.fc);
+	R_coeffs = R * ones(1, channels);
+end
 
 addpath ./../common
 
 blob8 = dcblock_build_blob(R_coeffs, endian);
 blob8_ipc4 = dcblock_build_blob(R_coeffs, endian, 4);
 
 % Generate output files
-tplg_write(tplg1_fn, blob8, "DCBLOCK");
+tplg_write(tplg1_fn, blob8, "DCBLOCK", ...
+	   "Exported with script example_dcblock.m", ...
+	   "cd tools/tune/dcblock; octave example_dcblock.m");
 blob_write(blob3_fn, blob8);
 alsactl_write(alsa3_fn, blob8);
 
-tplg2_write(tplg2_fn, blob8_ipc4, "dcblock_config", "Exported with script example_dcblock.m");
+tplg2_write(tplg2_fn, blob8_ipc4, "dcblock_config", ...
+	    "Exported with script example_dcblock.m" , ...
+	    "cd tools/tune/dcblock; octave example_dcblock.m");
 blob_write(blob4_fn, blob8_ipc4);
 alsactl_write(alsa4_fn, blob8_ipc4);
 
@@ -36,3 +69,54 @@ function example_dcblock()
 rmpath ./../common
 
 end
+
+% Finds with iterative search parameter R for given cutoff frequency
+function R = dcblock_rval_calculate(fs, fc_target)
+
+if (fc_target / fs < 10 / 48e3 || fc_target / fs > 1000 / 48e3)
+	error("Illegal fc_target");
+end
+
+h_target = 1 / sqrt(2); % -3.01 dB
+R = 0.5;
+R_step = 0.005;
+sign = 1;
+err_prev = inf;
+w = 2 * pi * fc_target / fs;
+j = sqrt(-1);
+z = exp(j * w);
+
+% Iteration 0
+h = (1 - z^-1) / (1 - R * z^-1);
+a = abs(h);
+err_prev = (h_target - abs(h))^2;
+R = R + sign * R_step;
+
+% Do more iterations
+for n = 1 : 200
+	h = (1 - z^-1) / (1 - R * z^-1);
+	err = (h_target - abs(h))^2;
+	if (err > err_prev)
+		sign = -sign;
+		R_step = R_step / 2;
+	end
+	R = R + sign * R_step;
+	err_prev = err;
+end
+
+% Sane result?
+if R < eps || R > 1 - eps
+	error("Calculate of R iteration failed")
+end
+
+% Sane high-pass function?
+f = [1 fc_target fs/2];
+b = [ 1 -1 ]; a = [ 1 -R ];
+h = freqz(b, a, f, fs);
+h_db = 20*log10(abs(h));
+err_hfc = abs(20*log10(1/sqrt(2)) - h_db(2));
+if err_hfc > 0.01 || h_db(1) > -10 || h_db(3) < 0 || h_db(3) > 1
+	error("Failed high-pass response");
+end
+
+end