Contents

% Example for multi-point equalization routine
%
% Sebastian J. Schlecht, Tuesday, 21. April 2020

clear; clc; close all;
rng(1);

Create some impulse responses from all pole models

n = 2^13;
numPoles = 50;
numOfIR = 10;

% create poles
poleAngles = rand(numPoles/2 , 1) * pi ;
poleMag = rand(numPoles/2 , 1) * 0.09 + 0.9;
polePos = poleMag .* exp(1i * poleAngles);
polePos = [polePos; conj(polePos)];

% simulated denominator
a = poly(polePos);

IR = [];
for it = 1:numOfIR
    % create randomised nominator
    zeroLengthFactor = 4; % scale length of nominator
    b = randn(numPoles * zeroLengthFactor, 1) ; b = b / norm(b);
    % compute simulated impulse response
    [h,t] = impz(b,a,n);
    IR = [IR, h];
end

Multi-point Equalization

[num, den] = MultiPointEQ(IR,numPoles * 2);

Plot results

figure(1); hold on; grid on;

% approximation frequency response
[h,w]  = freqz(num, den);
h = mag2db(abs(h));
p1 = plot(w, h, 'r','LineWidth',2);

% inverse frequency response
[h,w]  = freqz(den, 1);
h = mag2db(abs(h));
p2 = plot(w, h, 'g');

% last IR frequency response
for it = 1:numOfIR
    [h,w]  = freqz(IR(:,it), 1);
    h = mag2db(abs(h));
    p3 = plot(w, h);
end
uistack(p1,'top');

% labels
axis tight
xlabel('Frequency [rad]')
ylabel('Magnitude [dB]')
title('Magnitude Response of Multi-Point Equalisation')
legend([p1,p2,p3(1)],'Multi-point Approximation', 'Multi-point Inverse','Simulated IRs')
hold off;

Published with MATLAB® R2019b

Contents

function [num, den] = MultiPointEQ(IR,varargin)

%MultiPointEQ Estimates frequency response with allpole filter.
%   [num, den] = MultiPointEQ(IR,varargin) reads multiple impulse responses IR,
%   and approximates the magnitude response by an allpole model. The method
%   is taken from
%       Multiple-point equalization of room transfer functions by using common acoustical poles
%       by Haneda, Y. and Makino, S. and Kaneda, Y.
%       in IEEE Transactions on Speech and Audio Processing.
%
%   INPUT:
%       IR                          impulse response [ len, numCh ]
%       numPoles (optional)         filter order of the approximation
%
%   OUTPUT
%       num                         numerator of filter (always 1)
%       den                         denominator of filter [ 1 , numPoles ]
%
%
%
%
% Author: Dr.-Ing. Sebastian Jiro Schlecht,
% Aalto University, Finland
% email address: sebastian.schlecht@aalto.fi
% Website: sebastianjiroschlecht.com
% 21. April 2020; Last revision: 21. April 2020

input parser

p = inputParser;
defaultNumPoles = 100;

addRequired(p,'IR',@isnumeric);
addOptional(p,'numPoles',defaultNumPoles,@isnumeric);

parse(p,IR,varargin{:});

IR = p.Results.IR;
P = p.Results.numPoles;

[n,m] = size(IR);
if m > n % rotate IR if not correct
    error('Impulse Responses have to be column vectors')
end

Algortihm

W = [];
v = [];

for itm = 1:m
    H_i = [];
    h_i = IR(:,itm);
    for it = 1:P
        H_i = [H_i, [zeros(it - 1, 1); h_i; zeros(P - it, 1)]];
    end
    W = [W; H_i];
    v = [v; [h_i(2:end); zeros(P - 0, 1)]];
end

a = (W' * W) \ (W' * v);

Result

den = [1; -a];
num = 1;

Published with MATLAB® R2019b