## Called Functions

```% 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
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;
```
MultiPointEQ

## 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;

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;
```