OrgACRN
/
soft
mirror of https://github.com/thesofproject/soft.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

EQ Tool: Fix compatibility with Matlab and improve accuracy of IIR filters 

This patch changes code to allocate filter poles and zeros to biquads
to allow seeing poor SNR due to strong scaling vs. available fixed point
filter word length. The issue was severe especially with Matlab version
of tf2sos() function. The replacement is straighforward z, p frequency
sort. The scaling of biquads was changed to max -6 dB instead of previous
-3 dB to be more safe against internal signal clipping.

The IIR design was changed entirely to be based on zeros and poles
(z, p, k) in intermediate storate instead of polynomial coefficients.
This improves the accuracy especially of low frequencies high-pass
designs.

The default max. IIR order is increased from 12 to 16.

Finally since the grpdelay() function in octave-signal package calculates
incorrect result the group delay computation and plotting is disabled
for now in Octave. It works in Matlab.

Signed-off-by: Seppo Ingalsuo <seppo.ingalsuo@linux.intel.com>
This commit is contained in:
Seppo Ingalsuo 2018-08-30 13:49:17 +03:00
parent 41861aac37
commit 206deaf457
10 changed files with 148 additions and 85 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

15
tune/eq/eq_compute.m
View File

@ -42,9 +42,9 @@ eq = preprocess_responses(eq);
%% Define target (e.g. speaker) response as parametric filter. This could also
% be numerical data interpolated to the grid.
if length(eq.parametric_target_response) > 0
[eq.b_t, eq.a_t] = eq_define_parametric_eq( ...
[eq.t_z, eq.t_p, eq.t_k] = eq_define_parametric_eq( ...
eq.parametric_target_response, eq.fs);
eq.t_db = eq_compute_response(eq.b_t, eq.a_t, eq.f, eq.fs);
eq.t_db = eq_compute_response(eq.t_z, eq.t_p, eq.t_k, eq.f, eq.fs);
end
if isempty(eq.t_db)
@ -64,16 +64,15 @@ eq.err_db = eq.t_db - eq.m_db;
%% Parametric IIR EQ definition
if eq.enable_iir
[eq.b_p, eq.a_p] = eq_define_parametric_eq(eq.peq, eq.fs);
if length(eq.b_p) > 2*eq.iir_biquads_max+1
[eq.p_z, eq.p_p, eq.p_k] = eq_define_parametric_eq(eq.peq, eq.fs);
if max(length(eq.p_z), length(eq.p_p)) > 2*eq.iir_biquads_max
error('Maximum number of IIR biquads is exceeded');
end
else
eq.b_p = 1;
eq.a_p = 1;
[eq.p_z, eq.p_p, eq.p_k] = tf2zp(1, 1);
end
[eq.iir_eq_db, eq.iir_eq_ph, eq.iir_eq_gd] = eq_compute_response(eq.b_p, ...
eq.a_p, eq.f, eq.fs);
[eq.iir_eq_db, eq.iir_eq_ph, eq.iir_eq_gd] = eq_compute_response(eq.p_z, ...
eq.p_p, eq.p_k, eq.f, eq.fs);
%% FIR EQ computation

47
tune/eq/eq_compute_response.m
View File

@ -1,4 +1,4 @@
function [m, ph, gd] = eq_compute_response(b, a, f, fs)
function [m, ph, gd] = eq_compute_response(z, p, k, f, fs)
%%
% Copyright (c) 2016, Intel Corporation
@ -30,16 +30,43 @@ function [m, ph, gd] = eq_compute_response(b, a, f, fs)
% Author: Seppo Ingalsuo <seppo.ingalsuo@linux.intel.com>
%
h = freqz(b, a, f, fs);
m = 20*log10(abs(h));
ph = 180/pi*angle(h);
switch nargin
case 3
b = z; % 1st arg
f = p; % 2nd arg
fs = k; % 3rd arg
h = freqz(b, 1, f, fs);
m = 20*log10(abs(h));
ph = 180/pi*angle(h);
lb = length(b);
la = length(a);
if lb == 1 && la == 1
gd = zeros(1, length(f));
else
gd = 1/fs*grpdelay(b, a, f, fs);
if length(b) == 1
gd = zeros(1, length(f));
else
if exist('OCTAVE_VERSION', 'builtin')
% grpdelay() has some issue so better to not show a plot
gd = NaN * zeros(1, length(f));
else
gd = 1/fs*grpdelay(b, 1, f, fs);
end
end
case 5
[b, a] = zp2tf(z, p, k);
h = freqz(b, a, f, fs);
m = 20*log10(abs(h));
ph = 180/pi*angle(h);
if length(z) == 0 && length(p) == 0
gd = zeros(1, length(f));
else
if exist('OCTAVE_VERSION', 'builtin')
% grpdelay() has some issue so better to not show a plot
gd = NaN * zeros(1, length(f));
else
gd = 1/fs*grpdelay(b, a, f, fs);
end
end
otherwise
error('Incorrect input parameters');
end
end

4
tune/eq/eq_compute_tot_response.m
View File

@ -31,8 +31,8 @@ function eq = eq_compute_tot_response(eq)
%
% FIR response and combined IIR+FIR EQ
[eq.iir_eq_db, eq.iir_eq_ph, eq.iir_eq_gd] = eq_compute_response(eq.b_p, eq.a_p, eq.f, eq.fs);
[eq.fir_eq_db, eq.fir_eq_ph, eq.fir_eq_gd] = eq_compute_response(eq.b_fir, 1, eq.f, eq.fs);
[eq.iir_eq_db, eq.iir_eq_ph, eq.iir_eq_gd] = eq_compute_response(eq.p_z, eq.p_p, eq.p_k, eq.f, eq.fs);
[eq.fir_eq_db, eq.fir_eq_ph, eq.fir_eq_gd] = eq_compute_response(eq.b_fir, eq.f, eq.fs);
eq.tot_eq_db = eq.iir_eq_db + eq.fir_eq_db;
eq.tot_eq_gd = eq.iir_eq_gd + eq.fir_eq_gd;

12
tune/eq/eq_defaults.m
View File

@ -58,7 +58,7 @@ p.fir_autoband = 1;
p.enable_fir = 0;
% IIR conf
p.iir_biquads_max = 6;
p.iir_biquads_max = 8;
p.enable_iir = 0;
% Initialize other fields those are computed later to allow use of struct
@ -73,15 +73,17 @@ p.num_responses = 0;
p.m_db_s = [];
p.gd_s_s = [];
p.logsmooth_noct = 0;
p.b_t = [];
p.a_t = [];
p.t_z = [];
p.t_p = [];
p.t_k = [];
p.t_db = [];
p.m_db_abs = 0;
p.m_db_offs = 0;
p.raw_m_noalign_db = [];
p.err_db = [];
p.b_p = 0;
p.a_p = 0;
p.p_z = [];
p.p_p = [];
p.p_k = [];
p.iir_eq_db = [];
p.iir_eq_ph = [];
p.iir_eq_gd = [];

32
tune/eq/eq_define_parametric_eq.m
View File

@ -1,4 +1,4 @@
function [b_t, a_t] = eq_define_parametric_eq(peq, fs)
function [z, p, k] = eq_define_parametric_eq(peq, fs)
%%
% Copyright (c) 2016, Intel Corporation
@ -36,20 +36,27 @@ PEQ_LS1 = 5; PEQ_LS2 = 6; PEQ_HS1 = 7; PEQ_HS2 = 8;
PEQ_PN2 = 9; PEQ_LP4 = 10; PEQ_HP4 = 11;
sp = size(peq);
b_t = 1; a_t = 1;
z = [];
p = [];
k = 1;
for i=1:sp(1)
type = peq(i,1);
f = peq(i,2);
g = peq(i,3);
Q = peq(i,4);
if f < fs/2
a0 = [];
b0 = [];
z0 = [];
p0 = [];
k0 = [];
switch peq(i,1)
case PEQ_HP1, [b0, a0] = butter(1, 2*f/fs, 'high');
case PEQ_HP2, [b0, a0] = butter(2, 2*f/fs, 'high');
case PEQ_HP4, [b0, a0] = butter(4, 2*f/fs, 'high');
case PEQ_LP1, [b0, a0] = butter(1, 2*f/fs);
case PEQ_LP2, [b0, a0] = butter(2, 2*f/fs);
case PEQ_LP4, [b0, a0] = butter(4, 2*f/fs);
case PEQ_HP1, [z0, p0, k0] = butter(1, 2*f/fs, 'high');
case PEQ_HP2, [z0, p0, k0] = butter(2, 2*f/fs, 'high');
case PEQ_HP4, [z0, p0, k0] = butter(4, 2*f/fs, 'high');
case PEQ_LP1, [z0, p0, k0] = butter(1, 2*f/fs);
case PEQ_LP2, [z0, p0, k0] = butter(2, 2*f/fs);
case PEQ_LP4, [z0, p0, k0] = butter(4, 2*f/fs);
case PEQ_LS1, [b0, a0] = low_shelf_1st(f, g, fs);
case PEQ_LS2, [b0, a0] = low_shelf_2nd(f, g, fs);
case PEQ_HS1, [b0, a0] = high_shelf_1st(f, g, fs);
@ -58,7 +65,14 @@ for i=1:sp(1)
otherwise
error('Unknown parametric EQ type');
end
b_t=conv(b_t, b0); a_t = conv(a_t, a0);
if length(a0) > 0
[z0, p0, k0] = tf2zp(b0, a0);
end
if length(k0) > 0
z = [z ; z0(:)];
p = [p ; p0(:)];
k = k * k0;
end
end
end
end

58
tune/eq/eq_iir_blob_quant.m
View File

@ -1,12 +1,12 @@
function iir_resp = eq_iir_blob_quant(eq_b, eq_a)
function iir_resp = eq_iir_blob_quant(eq_z, eq_p, eq_k)
%% Convert IIR coefficients to 2nd order sections and quantize
%
% iir_resp = eq_iir_blob_quant(b, a, bits)
% iir_resp = eq_iir_blob_quant(z, p, k)
%
% b - numerator coefficients
% a - denominator coefficients
% bits - number of bits to quantize
% z - zeros
% p - poles
% k - gain
%
% iir_resp - vector to setup an IIR equalizer with number of sections, shifts,
% and quantized coefficients
@ -47,21 +47,41 @@ bits_iir = 32; % Q2.30
qf_iir = 30;
bits_gain = 16; % Q2.14
qf_gain = 14;
scale_max = -3; % dB, scale biquad L-inf norm
scale_max = -6; % dB, scale biquads peak gain to this
plot_pz = 0;
plot_fr = 0;
%% Convert IIR to 2nd order sections
if exist('OCTAVE_VERSION', 'builtin')
[sos, gain] = tf2sos(eq_b, eq_a);
else
% TODO: tf2sos produces in Matlab EQs with zero output due to
% very high gain variations within biquad. Need to investigate if
% this is incorrect usage and a bug here.
[sos, gain] = tf2sos(eq_b, eq_a, 'UP', Inf);
fprintf('Warning: Problems have been seen with some IIR designs.\n');
fprintf('Warning: Please check the result.\n');
% This a simple implementation of zp2sos() function. It is not used here due
% to utilization of rather strong scaling and resulting low SNR with the
% available word length in EQ in SOF. This poles and zeros allocation to
% biquads is base only in ascending sort of angular frequency.
sz = length(eq_z);
sp = length(eq_p);
sk = length(eq_k);
nb = max(sz, sp)/2;
az = angle(eq_z);
ap = angle(eq_p);
[~, iz] = sort(abs(az));
[~, ip] = sort(abs(ap));
eq_z = eq_z(iz);
eq_p = eq_p(ip);
sos = zeros(nb, 6);
for i = 1:nb
j = 2*(i-1)+1;
if i == 1
[b, a] = zp2tf(eq_z(j:j+1), eq_p(j:j+1), eq_k);
else
[b, a] = zp2tf(eq_z(j:j+1), eq_p(j:j+1), 1);
end
sos(i,1:3) = b;
sos(i,4:6) = a;
end
gain = 1;
%% Convert 2nd order sections to SOF parameters format and scale the biquads
% with criteria below (Gain max -6 dB at any frequency). Then calculate
% scaling shifts and finally gain multiplier for output.
sz = size(sos);
nbr_sections = sz(1);
n_section_header = 2;
@ -70,10 +90,10 @@ iir_resp = int32(zeros(1,n_section_header+nbr_sections*n_section));
iir_resp(1) = nbr_sections;
iir_resp(2) = nbr_sections; % Note: All sections in series
scale_max_lin = 10^(scale_max/20);
for n=1:nbr_sections
b = sos(n,1:3);
a = sos(n,4:6);
if plot_pz
figure
zplane(b,a);
@ -84,14 +104,12 @@ for n=1:nbr_sections
np = 1024;
[h, w] = freqz(b, a, np);
hm = max(abs(h));
scale = 10^(scale_max/20)/hm;
scale = scale_max_lin/hm;
gain_remain = 1/scale;
gain = gain*gain_remain;
b = b * scale;
ma = max(abs(a));
mb = max(abs(b));
if plot_fr
figure
[h, w] = freqz(b, a, np);
@ -119,6 +137,8 @@ for n=1:nbr_sections
iir_resp(m+5) = section_shift;
iir_resp(m+6) = eq_coef_quant( section_gain, bits_gain, qf_gain);
%fprintf('sec=%d, shift=%d, gain=%f\n', n, section_shift, section_gain);
end
end

4
tune/eq/eq_norm.m
View File

@ -61,13 +61,13 @@ end
%% Adjust FIR and IIR gains if enabled
if eq.enable_fir && eq.enable_iir
eq.b_fir = eq.b_fir * g_fir * g_offs;
eq.b_p = eq.b_p * g_iir * g_offs;
eq.p_k = eq.p_k * g_iir * g_offs;
end
if eq.enable_fir && eq.enable_iir == 0
eq.b_fir = eq.b_fir * g_fir * g_offs;
end
if eq.enable_fir == 0 && eq.enable_iir
eq.b_p = eq.b_p * g_iir * g_offs;
eq.p_k = eq.p_k * g_iir * g_offs;
end
%% Re-compute response after adjusting gain

52
tune/eq/eq_plot.m
View File

@ -115,7 +115,7 @@ if length(eq.b_fir) > 1
end
%% IIR filter
if length(eq.b_p) > 1
if length(eq.p_z) > 1 || length(eq.p_p) > 1
% Response
fh=figure(fn); fn = fn+1;
semilogx(eq.f, eq.iir_eq_db);
@ -128,7 +128,7 @@ if length(eq.b_p) > 1
% Polar
fh=figure(fn); fn = fn+1;
zplane(eq.b_p, eq.a_p);
zplane(eq.p_z, eq.p_p);
grid on;
tstr = sprintf('IIR zeros and poles: %s', eq.name);
title(tstr);
@ -137,7 +137,8 @@ if length(eq.b_p) > 1
ti = 50e-3;
x = zeros(1, ti * round(eq.fs));
x(1) = 1;
y = filter(eq.b_p, eq.a_p, x);
sos = zp2sos(eq.p_z, eq.p_p, eq.p_k);
y = sosfilt(sos, x);
fh=figure(fn); fn = fn+1;
t = (0:(length(x)-1)) / eq.fs;
plot(t, y);
@ -149,28 +150,29 @@ if length(eq.b_p) > 1
end
%% Group delay
fh=figure(fn); fn = fn+1;
if eq.enable_fir && eq.enable_iir
semilogx(eq.f, eq.tot_eq_gd * 1e3, ...
eq.f, eq.fir_eq_gd * 1e3, '--', ...
eq.f, eq.iir_eq_gd * 1e3, '--');
legend('Combined','FIR','IIR');
if ~exist('OCTAVE_VERSION', 'builtin')
% Skip plot if running in Octave due to incorrect result
fh=figure(fn); fn = fn+1;
if eq.enable_fir && eq.enable_iir
semilogx(eq.f, eq.tot_eq_gd * 1e3, ...
eq.f, eq.fir_eq_gd * 1e3, '--', ...
eq.f, eq.iir_eq_gd * 1e3, '--');
legend('Combined','FIR','IIR');
end
if eq.enable_fir && eq.enable_iir == 0
semilogx(eq.f, eq.fir_eq_gd * 1e3);
legend('FIR');
end
if eq.enable_fir == 0 && eq.enable_iir
semilogx(eq.f, eq.iir_eq_gd * 1e3);
legend('IIR');
end
grid on;
xlabel('Frequency (Hz)');
ylabel('Group delay (ms)');
ax = axis; axis([eq.p_fmin eq.p_fmax ax(3:4)]);
tstr = sprintf('Filter group delay: %s', eq.name);
title(tstr);
end
if eq.enable_fir && eq.enable_iir == 0
semilogx(eq.f, eq.fir_eq_gd * 1e3);
legend('FIR');
end
if eq.enable_fir == 0 && eq.enable_iir
semilogx(eq.f, eq.iir_eq_gd * 1e3);
legend('IIR');
end
grid on;
xlabel('Frequency (Hz)');
ylabel('Group delay (ms)');
ax = axis; axis([eq.p_fmin eq.p_fmax ax(3:4)]);
tstr = sprintf('Filter group delay: %s', eq.name);
title(tstr);
end

7
tune/eq/example_iir_eq.m
View File

@ -41,12 +41,11 @@ fs = 48e3;
eq_loud = loudness_iir_eq(fs);
%% Define a passthru IIR EQ equalizer
b_pass = [1 0 0];
a_pass = [1 0 0];
[z_pass, p_pass, k_pass] = tf2zp([1 0 0],[1 0 0]);
%% Quantize and pack filter coefficients plus shifts etc.
bq_pass = eq_iir_blob_quant(b_pass, a_pass);
bq_loud = eq_iir_blob_quant(eq_loud.b_p, eq_loud.a_p);
bq_pass = eq_iir_blob_quant(z_pass, p_pass, k_pass);
bq_loud = eq_iir_blob_quant(eq_loud.p_z, eq_loud.p_p, eq_loud.p_k);
%% Build blob
platform_max_channels = 4; % Setup max 4 channels EQ

2
tune/eq/example_spk_eq.m
View File

@ -140,7 +140,7 @@ end
%% Export IIR part
if eq.enable_iir
bq_iir = eq_iir_blob_quant(eq.b_p, eq.a_p);
bq_iir = eq_iir_blob_quant(eq.p_z, eq.p_p, eq.p_k);
bm_iir = eq_iir_blob_merge(platform_max_channels, ...
num_responses, ...
assign_response, ...