From 206deaf457b4a9d7260ea4cd99220e7d4297099d Mon Sep 17 00:00:00 2001 From: Seppo Ingalsuo Date: Thu, 30 Aug 2018 13:49:17 +0300 Subject: [PATCH] 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 --- tune/eq/eq_compute.m | 15 ++++---- tune/eq/eq_compute_response.m | 47 +++++++++++++++++++------ tune/eq/eq_compute_tot_response.m | 4 +-- tune/eq/eq_defaults.m | 12 ++++--- tune/eq/eq_define_parametric_eq.m | 32 ++++++++++++----- tune/eq/eq_iir_blob_quant.m | 58 +++++++++++++++++++++---------- tune/eq/eq_norm.m | 4 +-- tune/eq/eq_plot.m | 52 ++++++++++++++------------- tune/eq/example_iir_eq.m | 7 ++-- tune/eq/example_spk_eq.m | 2 +- 10 files changed, 148 insertions(+), 85 deletions(-) diff --git a/tune/eq/eq_compute.m b/tune/eq/eq_compute.m index a5143b4..57a6293 100644 --- a/tune/eq/eq_compute.m +++ b/tune/eq/eq_compute.m @@ -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 diff --git a/tune/eq/eq_compute_response.m b/tune/eq/eq_compute_response.m index cf8a295..e8a4704 100644 --- a/tune/eq/eq_compute_response.m +++ b/tune/eq/eq_compute_response.m @@ -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 % -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 diff --git a/tune/eq/eq_compute_tot_response.m b/tune/eq/eq_compute_tot_response.m index 828b936..250c364 100644 --- a/tune/eq/eq_compute_tot_response.m +++ b/tune/eq/eq_compute_tot_response.m @@ -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; diff --git a/tune/eq/eq_defaults.m b/tune/eq/eq_defaults.m index cf16aea..880e4e3 100644 --- a/tune/eq/eq_defaults.m +++ b/tune/eq/eq_defaults.m @@ -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 = []; diff --git a/tune/eq/eq_define_parametric_eq.m b/tune/eq/eq_define_parametric_eq.m index 06136e2..2beecf2 100644 --- a/tune/eq/eq_define_parametric_eq.m +++ b/tune/eq/eq_define_parametric_eq.m @@ -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 diff --git a/tune/eq/eq_iir_blob_quant.m b/tune/eq/eq_iir_blob_quant.m index 823d256..0c19fb3 100644 --- a/tune/eq/eq_iir_blob_quant.m +++ b/tune/eq/eq_iir_blob_quant.m @@ -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 diff --git a/tune/eq/eq_norm.m b/tune/eq/eq_norm.m index fdb1cb5..47b40d9 100644 --- a/tune/eq/eq_norm.m +++ b/tune/eq/eq_norm.m @@ -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 diff --git a/tune/eq/eq_plot.m b/tune/eq/eq_plot.m index 2dc5743..95968c4 100644 --- a/tune/eq/eq_plot.m +++ b/tune/eq/eq_plot.m @@ -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 diff --git a/tune/eq/example_iir_eq.m b/tune/eq/example_iir_eq.m index 963d6b8..ddcc17e 100644 --- a/tune/eq/example_iir_eq.m +++ b/tune/eq/example_iir_eq.m @@ -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 diff --git a/tune/eq/example_spk_eq.m b/tune/eq/example_spk_eq.m index fd4cafc..7be22c3 100644 --- a/tune/eq/example_spk_eq.m +++ b/tune/eq/example_spk_eq.m @@ -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, ...