145 lines
4.9 KiB
Matlab
145 lines
4.9 KiB
Matlab
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(z, p, k)
|
|
%
|
|
% z - zeros
|
|
% p - poles
|
|
% k - gain
|
|
%
|
|
% iir_resp - vector to setup an IIR equalizer with number of sections, shifts,
|
|
% and quantized coefficients
|
|
%
|
|
|
|
%%
|
|
% Copyright (c) 2016, Intel Corporation
|
|
% All rights reserved.
|
|
%
|
|
% Redistribution and use in source and binary forms, with or without
|
|
% modification, are permitted provided that the following conditions are met:
|
|
% * Redistributions of source code must retain the above copyright
|
|
% notice, this list of conditions and the following disclaimer.
|
|
% * Redistributions in binary form must reproduce the above copyright
|
|
% notice, this list of conditions and the following disclaimer in the
|
|
% documentation and/or other materials provided with the distribution.
|
|
% * Neither the name of the Intel Corporation nor the
|
|
% names of its contributors may be used to endorse or promote products
|
|
% derived from this software without specific prior written permission.
|
|
%
|
|
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
|
|
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
|
|
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
|
|
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
|
|
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
|
|
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
|
|
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
|
|
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
|
|
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
|
|
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
|
|
% POSSIBILITY OF SUCH DAMAGE.
|
|
%
|
|
% Author: Seppo Ingalsuo <seppo.ingalsuo@linux.intel.com>
|
|
%
|
|
|
|
%% Settings
|
|
bits_iir = 32; % Q2.30
|
|
qf_iir = 30;
|
|
bits_gain = 16; % Q2.14
|
|
qf_gain = 14;
|
|
scale_max = -6; % dB, scale biquads peak gain to this
|
|
plot_pz = 0;
|
|
plot_fr = 0;
|
|
|
|
%% Convert IIR to 2nd order sections
|
|
% 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;
|
|
n_section = 7;
|
|
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);
|
|
tstr = sprintf('SOS %d poles and zeros', n);
|
|
title(tstr);
|
|
end
|
|
|
|
np = 1024;
|
|
[h, w] = freqz(b, a, np);
|
|
hm = max(abs(h));
|
|
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);
|
|
plot(w, 20*log10(abs(h))); grid on;
|
|
xlabel('Frequency (w)');
|
|
ylabel('Magnitude (dB)');
|
|
tstr = sprintf('SOS %d frequency response', n);
|
|
title(tstr);
|
|
end
|
|
|
|
%% Apply remaining gain at last section output
|
|
if n == nbr_sections
|
|
section_shift = -fix(log(gain)/log(2));
|
|
section_gain= gain/2^(-section_shift);
|
|
else
|
|
section_shift = 0;
|
|
section_gain = 1;
|
|
end
|
|
|
|
%% Note: Invert sign of a!
|
|
%% Note: a(1) is omitted, it's always 1
|
|
m = n_section_header+(n-1)*n_section+1;
|
|
iir_resp(m:m+1) = eq_coef_quant(-a(3:-1:2), bits_iir, qf_iir);
|
|
iir_resp(m+2:m+4) = eq_coef_quant( b(3:-1:1), bits_iir, qf_iir);
|
|
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
|