Merge pull request #61 from singalsu/fix_matlab_iir_eq_proposal

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

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);
+        [eq.p_z, eq.p_p, eq.p_k] = eq_define_parametric_eq(eq.peq, eq.fs);
-        if length(eq.b_p) > 2*eq.iir_biquads_max+1
+        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.p_z, eq.p_p, eq.p_k] = tf2zp(1, 1);
-        eq.a_p = 1;
 end
-[eq.iir_eq_db, eq.iir_eq_ph, eq.iir_eq_gd] = eq_compute_response(eq.b_p, ...
+[eq.iir_eq_db, eq.iir_eq_ph, eq.iir_eq_gd] = eq_compute_response(eq.p_z, ...
-        eq.a_p, eq.f, eq.fs);
+        eq.p_p, eq.p_k, eq.f, eq.fs);
 
 
 %% FIR EQ computation

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 <seppo.ingalsuo@linux.intel.com>
 %
 
-h = freqz(b, a, f, fs);
+switch nargin
-m = 20*log10(abs(h));
+        case 3
-ph = 180/pi*angle(h);
+                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);
+                if length(b) == 1
-la = length(a);
+                        gd = zeros(1, length(f));
-if lb == 1 && la == 1
+                else
-	gd = zeros(1, length(f));
+			if exist('OCTAVE_VERSION', 'builtin')
-else
+				% grpdelay() has some issue so better to not show a plot
-	gd = 1/fs*grpdelay(b, a, f, fs);
+				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

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.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, 1, 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;
 

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.t_z = [];
-p.a_t = [];
+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.p_z = [];
-p.a_p = 0;
+p.p_p = [];
+p.p_k = [];
 p.iir_eq_db = [];
 p.iir_eq_ph = [];
 p.iir_eq_gd = [];

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_HP1, [z0, p0, k0] = butter(1, 2*f/fs, 'high');
-                        case PEQ_HP2, [b0, a0] = butter(2, 2*f/fs, 'high');
+                        case PEQ_HP2, [z0, p0, k0] = butter(2, 2*f/fs, 'high');
-                        case PEQ_HP4, [b0, a0] = butter(4, 2*f/fs, 'high');
+                        case PEQ_HP4, [z0, p0, k0] = butter(4, 2*f/fs, 'high');
-                        case PEQ_LP1, [b0, a0] = butter(1, 2*f/fs);
+                        case PEQ_LP1, [z0, p0, k0] = butter(1, 2*f/fs);
-                        case PEQ_LP2, [b0, a0] = butter(2, 2*f/fs);
+                        case PEQ_LP2, [z0, p0, k0] = butter(2, 2*f/fs);
-                        case PEQ_LP4, [b0, a0] = butter(4, 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

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
+%  z - zeros
-%  a - denominator coefficients
+%  p - poles
-%  bits - number of bits to quantize
+%  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')
+%  This a simple implementation of zp2sos() function. It is not used here due
-        [sos, gain] = tf2sos(eq_b, eq_a);
+%  to utilization of rather strong scaling and resulting low SNR with the
-else
+%  available word length in EQ in SOF. This poles and zeros allocation to
-        % TODO: tf2sos produces in Matlab EQs with zero output due to
+%  biquads is base only in ascending sort of angular frequency.
-        % very high gain variations within biquad. Need to investigate if
+sz = length(eq_z);
-        % this is incorrect usage and a bug here.
+sp = length(eq_p);
-        [sos, gain] = tf2sos(eq_b, eq_a, 'UP', Inf);
+sk = length(eq_k);
-        fprintf('Warning: Problems have been seen with some IIR designs.\n');
+nb = max(sz, sp)/2;
-        fprintf('Warning: Please check the result.\n');
+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

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

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 ~exist('OCTAVE_VERSION', 'builtin')
-if eq.enable_fir && eq.enable_iir
+	% Skip plot if running in Octave due to incorrect result
-        semilogx(eq.f, eq.tot_eq_gd * 1e3, ...
+	fh=figure(fn); fn = fn+1;
-                eq.f, eq.fir_eq_gd * 1e3, '--', ...
+	if eq.enable_fir && eq.enable_iir
-                eq.f, eq.iir_eq_gd * 1e3, '--');
+		semilogx(eq.f, eq.tot_eq_gd * 1e3, ...
-        legend('Combined','FIR','IIR');
+			 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

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];
+[z_pass, p_pass, k_pass] = tf2zp([1 0 0],[1 0 0]);
-a_pass = [1 0 0];
 
 %% Quantize and pack filter coefficients plus shifts etc.
-bq_pass = eq_iir_blob_quant(b_pass, a_pass);
+bq_pass = eq_iir_blob_quant(z_pass, p_pass, k_pass);
-bq_loud = eq_iir_blob_quant(eq_loud.b_p, eq_loud.a_p);
+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

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, ...