EQ: Add tool for acoustical frequency response measurement
This patch adds into SOFT/tune/eq a script mls_freq_resp.m that can be used as low budget or quick developer tool to measure device frequency response. The configuration files mls_play_config.txt and mls_rec_config.txt configure the audio devices for local or remote ssh aplay/arecord execution. The MLS measurement with a short time window, e.g. less than 10 ms can partially achieve the same as anechoic measurement. However this is not intended to replace usage any professional audio measurement. If using this use with care! Note also that as such this script measures the response of a speaker - microphone system. One response need to be known and under calibration to get reliable results for other. Even ADC and DAC responses may impact. Signed-off-by: Seppo Ingalsuo <seppo.ingalsuo@linux.intel.com>
This commit is contained in:
Seppo Ingalsuo
parent
e37aa3f422
commit
1cbba6c5e8
|
|
@ -0,0 +1,98 @@
|
|||
function [f, m_db, b] = mls_calc_resp(csv_fn, ref_wfn, meas_wfn, t_tot, np, f_lo, f_hi)
|
||||
|
||||
%% Calculate frequency response from MLS recordings
|
||||
%
|
||||
% [f, m, b] = calc_freq_resp(csv_fn, ref_fn, meas_fn, t_tot, np, f_lo, f_hi)
|
||||
%
|
||||
% Inputs parameters:
|
||||
% csv_fn - File name for CSV format response output
|
||||
% ref_fn - Reference MLS wave file
|
||||
% meas_fn - Captured MSL wave file
|
||||
% t_tot - Time window of MLS analysis in s, e.g. 5 - 10 ms
|
||||
% np - Number of frequency points
|
||||
% f_lo - Lower frequency limit for analysis
|
||||
% f_hi - Upper frequency limit for analysis
|
||||
|
||||
%%
|
||||
% Copyright (c) 2018, 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>
|
||||
%
|
||||
|
||||
%% Read audio files
|
||||
[xr, fs] = audioread(ref_wfn);
|
||||
[xm, fsm] = audioread(meas_wfn);
|
||||
sm = size(xm);
|
||||
nch = sm(2);
|
||||
|
||||
if (fs ~= fsm)
|
||||
error('Reference and measured sample rates do not match.');
|
||||
end
|
||||
|
||||
ind_imp = zeros(1, nch);
|
||||
for ch=1:nch
|
||||
[acor(:,ch),lag] = xcorr(xm(:,ch), xr, 'coeff');
|
||||
[~,I] = max(abs(acor(:,ch)));
|
||||
timeDiff = lag(I);
|
||||
% Crosscorrelation maximum points to max of impulse response.
|
||||
% Use it beginning of impulse response
|
||||
ind_imp(ch) = find(lag == timeDiff);
|
||||
end
|
||||
% Keep the channels time aligned by using average of max as reference
|
||||
ind_imp_avg = round(sum(ind_imp)/nch);
|
||||
|
||||
b = [];
|
||||
nwin_half = round(t_tot*fs);
|
||||
nwin_tot = 2*nwin_half;
|
||||
win = hanning(nwin_tot);
|
||||
for ch = 1:nch
|
||||
% All channels are windowed from the same location
|
||||
p1 = ind_imp_avg-nwin_half;
|
||||
p2 = p1 + nwin_tot-1;
|
||||
acor_windowed = acor(p1:p2,ch) .* win;
|
||||
b(:,ch) = acor_windowed;
|
||||
end
|
||||
|
||||
|
||||
%% Compute frequency response
|
||||
f = logspace(log10(f_lo),log10(f_hi),np);
|
||||
for ch=1:nch
|
||||
h(:,ch) = freqz(b(:,ch), 1, f, fs);
|
||||
end
|
||||
m_db = 20*log10(abs(h));
|
||||
|
||||
fh = fopen(csv_fn, 'w');
|
||||
for i = 1:np
|
||||
fprintf(fh, '%10.3f', f(i));
|
||||
for j = 1:nch
|
||||
fprintf(fh, ', %8.3f', m_db(i, j));
|
||||
end
|
||||
fprintf(fh, '\n');
|
||||
end
|
||||
fclose(fh);
|
||||
|
||||
end
|
||||
|
|
@ -0,0 +1,306 @@
|
|||
function [f, m_db] = mls_freq_resp(id)
|
||||
%% Measure frequency response with MLS test signal
|
||||
%
|
||||
% [f, m] = mls_freq_resp(id)
|
||||
%
|
||||
% Input parameters
|
||||
% id - A string identifier for test case. An id 'selftest' is for special
|
||||
% usage. It calculates response of filtered MLS signal and computes
|
||||
% the measurement vs. known. Deviation is reported as error. It can
|
||||
% be useful if the internal MLS measurement parameters are adjusted.
|
||||
%
|
||||
% Output parameters
|
||||
% f - Frequency vector in Hz
|
||||
% m - Measured magnitude responses in dB
|
||||
%
|
||||
% Configuration (edit these):
|
||||
% mls_play_config.txt
|
||||
% mls_rec_config.txt
|
||||
%
|
||||
% The script will return also a text CSV format file with name mls-<id>.txt.
|
||||
%
|
||||
|
||||
%%
|
||||
% Copyright (c) 2018, 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
|
||||
np = 1024; % Number of frequency points to use
|
||||
f_lo = 100; % Lower frequency limit for analysis
|
||||
f_hi = 20e3; % Upper frequency limit for analysis
|
||||
t_tot = 10e-3; % MLS analysis window length in s
|
||||
t_mls_s = 1.0; % MLS test signal length in s
|
||||
a_mls_db = -10; % MLS test signal amplitude in dB
|
||||
fs = 48e3; % Sample rate in Hz
|
||||
bits = 16; % Audio format to use (bits)
|
||||
fmt = 'S16_LE'; % Audio format to use (ALSA)
|
||||
dir = '/tmp'; % Directory for temporary files
|
||||
capture_level_max_db = -1; % Expected max. level
|
||||
capture_level_min_db = -30; % Expacted min. level
|
||||
|
||||
%% Get device identifier to use
|
||||
if nargin < 1
|
||||
id = 'unknown';
|
||||
end
|
||||
|
||||
if strcmp(id, 'selftest')
|
||||
selftest = 1;
|
||||
% Just some simulated speaker response to use as self test case
|
||||
stb = [ 0.341762453, -0.915611126, 0.482465118, ...
|
||||
1.017612317, -1.722527013, 0.711608745, ...
|
||||
0.630608859, -0.813609935, 0.267690582, ];
|
||||
sta = [ 1.000000000, -3.931695128, 6.630812276, ...
|
||||
-6.339248735, 3.800407709, -1.559376698, ...
|
||||
0.619626250, -0.317702349, 0.097453451, ];
|
||||
else
|
||||
selftest = 0;
|
||||
end
|
||||
measfn = sprintf('mls-%s.wav', id);
|
||||
csvfn = sprintf('mls-%s.txt', id);
|
||||
|
||||
%% Paths
|
||||
addpath('../../test/test_utils');
|
||||
|
||||
%% MLS
|
||||
n_mls = round(fs*t_mls_s);
|
||||
mls = 10^(a_mls_db/20) * (2 * mlsp12(1, n_mls) - 1);
|
||||
mlsfn = 'mls-ref.wav';
|
||||
audiowrite(mlsfn, mls, fs);
|
||||
|
||||
%% Chip markers and parameters for find sync
|
||||
[x1, m1] = sync_chirp(fs, 'up');
|
||||
[x2, m2] = sync_chirp(fs, 'down');
|
||||
fnd.fs = fs; % Sample rate
|
||||
fnd.sm = 3; % Max seek from start
|
||||
fnd.em = 3; % Max seek from end
|
||||
fnd.idle_t = 2; % max idle in start or end
|
||||
fnd.mark_t = m1.t; % Marker length
|
||||
fnd.nf = 1; % One signal (amplitude)
|
||||
fnd.na = 1; % One signal (frequency)
|
||||
fnd.tl = t_mls_s; % Length of signal
|
||||
fnd.mt = 0.1; % Threshold length to issue error
|
||||
fnd.is = 0; % Ignore from start
|
||||
fnd.ie = 0; % Ignore from end
|
||||
|
||||
%% Merge markers and MLS
|
||||
z = zeros(n_mls + m1.n + m2.n, 1);
|
||||
i1 = m1.n + 1;
|
||||
i2 = m1.n + n_mls;
|
||||
z(1:i1 - 1) = x1;
|
||||
z(i1:i2) = mls;
|
||||
z(i2 + 1:end) = x2;
|
||||
|
||||
%% Get config
|
||||
rec_cfg = meas_remote_rec_config(fs, fmt);
|
||||
play_cfg = meas_remote_play_config;
|
||||
|
||||
%% Capture MLS from all playback channel at time
|
||||
mixfn = 'mlsmix.wav';
|
||||
recfn = 'recch.wav';
|
||||
y = [];
|
||||
for i=1:play_cfg.nch
|
||||
fprintf('\n');
|
||||
fprintf('Measure playback channel %d\n', i);
|
||||
if selftest
|
||||
tz =zeros(2*fs+length(z), 1); % Pad 2s
|
||||
tz(fs:fs+length(z)-1) = z;
|
||||
t = filter(stb, sta, tz); % Filter with test response
|
||||
r = t(:) * ones(1, rec_cfg.nch); % Copy to all channels
|
||||
else
|
||||
x = zeros(length(z), play_cfg.nch);
|
||||
x(:,i) = z;
|
||||
mixdfn = sprintf('%s/%s', dir, mixfn);
|
||||
audiowrite(mixdfn, x, fs, 'BitsPerSample', bits);
|
||||
copy_playback(mixdfn, play_cfg);
|
||||
tcap = floor(3 + t_mls_s);
|
||||
remote_capture(recfn, rec_cfg, tcap);
|
||||
pause(1);
|
||||
remote_play(mixfn, play_cfg);
|
||||
pause(3);
|
||||
r = get_recording(recfn, rec_cfg);
|
||||
end
|
||||
[d, nt] = find_test_signal(r(:,1), fnd);
|
||||
for j = 1:rec_cfg.nch
|
||||
y(:, rec_cfg.nch*(i-1) + j) = r(d:d + nt -1, j);
|
||||
end
|
||||
m = 20*log10(max(abs(r)));
|
||||
fprintf('Peak levels for capture channels (dB):');
|
||||
for j = 1:rec_cfg.nch
|
||||
fprintf(' %5.1f', m(j));
|
||||
end
|
||||
fprintf('\n');
|
||||
if max(m) > capture_level_max_db
|
||||
fprintf('Warning: The recording level is too loud.\n');
|
||||
end
|
||||
if min(m) < capture_level_min_db
|
||||
fprintf('Warning: The recording level is too silent.\n');
|
||||
end
|
||||
|
||||
end
|
||||
audiowrite(measfn, y, fs, 'BitsPerSample', bits);
|
||||
fprintf('\n');
|
||||
fprintf('Done.\n');
|
||||
|
||||
[f, m_db, b] = mls_calc_resp(csvfn, mlsfn, measfn, t_tot, np, f_lo, f_hi);
|
||||
|
||||
figure
|
||||
idx = find(f>1e3, 1, 'first') - 1;
|
||||
m_db_align = m_db - m_db(idx);
|
||||
semilogx(f, m_db_align);
|
||||
ax=axis(); axis([f_lo f_hi ax(3:4)]);
|
||||
grid on;
|
||||
xlabel('Frequency (Hz)');
|
||||
ylabel('Magnitude (dB)');
|
||||
|
||||
if selftest
|
||||
title('Measured vs. reference response');
|
||||
h = freqz(stb, sta, f, fs);
|
||||
ref_db = 20*log10(abs(h));
|
||||
ref_db_align = ref_db - ref_db(idx);
|
||||
hold on;
|
||||
plot(f, ref_db_align, 'r--');
|
||||
hold off;
|
||||
idx = find(f < f_hi);
|
||||
idx = find(f(idx) > f_lo);
|
||||
m_lin = 10.^(m_db_align(idx)/20);
|
||||
ref_lin = 10.^(ref_db_align(idx)/20);
|
||||
err_lin = m_lin - ref_lin;
|
||||
e_rms = sqrt(mean(err_lin.^2));
|
||||
e_db = 20*log10(e_rms);
|
||||
figure;
|
||||
semilogx(f(idx), 20*log10(abs(err_lin)))
|
||||
grid on;
|
||||
xlabel('Frequency (Hz)');
|
||||
ylabel('Magnitude (dB)');
|
||||
title('Observed Error in self test');
|
||||
if e_db < -30
|
||||
fprintf('Passed self test. ');
|
||||
else
|
||||
fprintf('Failed self test. ');
|
||||
end
|
||||
fprintf('Response RMS error is %4.1f dB.\n', e_db);
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
function copy_playback(fn, cfg)
|
||||
if cfg.ssh
|
||||
cmd = sprintf('scp %s %s:%s/', fn, cfg.user, cfg.dir);
|
||||
fprintf('Remote copy: %s\n', cmd);
|
||||
system(cmd);
|
||||
else
|
||||
%cmd = sprintf('cp %s %s/', fn, cfg.dir);
|
||||
%fprintf('Local copy: %s\n', cmd);
|
||||
end
|
||||
end
|
||||
|
||||
function y = get_recording(fn, cfg)
|
||||
if cfg.ssh
|
||||
cmd = sprintf('scp %s:%s/%s %s', cfg.user, cfg.dir, fn, fn);
|
||||
fprintf('Remote copy: %s\n', cmd);
|
||||
else
|
||||
cmd = sprintf('cp %s/%s %s', cfg.dir, fn, fn);
|
||||
fprintf('Local copy: %s\n', cmd);
|
||||
end
|
||||
system(cmd);
|
||||
y = audioread(fn);
|
||||
delete(fn);
|
||||
end
|
||||
|
||||
function remote_play(fn, cfg)
|
||||
if cfg.ssh
|
||||
cmd = sprintf('ssh %s aplay -D%s %s/%s', cfg.user, cfg.dev, cfg.dir, fn);
|
||||
fprintf('Remote play: %s\n', cmd);
|
||||
else
|
||||
cmd = sprintf('aplay -D%s %s/%s', cfg.dev, cfg.dir, fn);
|
||||
fprintf('Local play: %s\n', cmd);
|
||||
end
|
||||
system(cmd);
|
||||
end
|
||||
|
||||
function remote_capture(fn, cfg, t)
|
||||
if cfg.ssh
|
||||
cmd = sprintf('ssh %s arecord -q -D%s %s -d %d %s/%s &', ...
|
||||
cfg.user, cfg.dev, cfg.fmt, t, cfg.dir, fn);
|
||||
fprintf('Remote capture: %s\n', cmd);
|
||||
else
|
||||
cmd = sprintf('arecord -q -D%s %s -d %d %s/%s &', ...
|
||||
cfg.dev, cfg.fmt, t, cfg.dir, fn);
|
||||
fprintf('Local capture: %s\n', cmd);
|
||||
end
|
||||
system(cmd);
|
||||
end
|
||||
|
||||
function play = meas_remote_play_config()
|
||||
source mls_play_config.txt;
|
||||
fprintf('\nThe setttings for remote playback are\n');
|
||||
fprintf('Use ssh : %d\n', play.ssh);
|
||||
fprintf('User : %s\n', play.user);
|
||||
fprintf('Directory : %s\n', play.dir);
|
||||
fprintf('Device : %s\n', play.dev);
|
||||
fprintf('Channels : %d\n', play.nch);
|
||||
end
|
||||
|
||||
function rec = meas_remote_rec_config(fs, fmt)
|
||||
source mls_rec_config.txt;
|
||||
rec.fmt = sprintf('-t wav -c %d -f %s -r %d', ...
|
||||
rec.nch, fmt, fs);
|
||||
|
||||
fprintf('\nThe setttings for remote capture are\n');
|
||||
fprintf('Use ssh : %d\n', rec.ssh);
|
||||
fprintf('User : %s\n', rec.user);
|
||||
fprintf('Directory : %s\n', rec.dir);
|
||||
fprintf('Device : %s\n', rec.dev);
|
||||
fprintf('format : %s\n', rec.fmt);
|
||||
fprintf('Channels : %d\n', rec.nch);
|
||||
end
|
||||
|
||||
function [x, seed] = mlsp12(seed, n)
|
||||
|
||||
% Based on book Numerical Recipes in C, chapter 7.4 Generation of
|
||||
% random bits method II p. 298-299 example and 12 bit primitive
|
||||
% polynomial (12, 6, 4, 1, 0)
|
||||
|
||||
x = zeros(1,n);
|
||||
ib1 = 2^(1-1);
|
||||
ib4 = 2^(4-1);
|
||||
ib6 = 2^(6-1);
|
||||
ib12 = 2^(12-1);
|
||||
mask = ib1 + ib4 + ib6 + ib12;
|
||||
|
||||
for i = 1:n
|
||||
if bitand(seed, ib12)
|
||||
seed = bitor(bitxor(seed, mask) * 2 , ib1);
|
||||
x(i) = 1;
|
||||
else
|
||||
seed = seed * 2;
|
||||
x(i) = 0;
|
||||
end
|
||||
end
|
||||
end
|
||||
|
|
@ -0,0 +1,5 @@
|
|||
play.ssh = 0; % Set to use remote ssh commands
|
||||
play.user = ''; % Set user@domain for ssh
|
||||
play.dir = '/tmp'; % directory for temporary files
|
||||
play.dev = 'hw:0,0'; % Audio device for playback
|
||||
play.nch = 2; % Number of playback channels to test
|
||||
|
|
@ -0,0 +1,5 @@
|
|||
rec.ssh = 0; % Set to 1 for remote capture
|
||||
rec.user = ''; % Set to user@domain for ssh
|
||||
rec.dir = '/tmp'; % Directory for temporary files
|
||||
rec.dev = 'hw:0,0'; % Audio capture device
|
||||
rec.nch = 2; % Number audio capture channels to use
|
||||
Loading…
Reference in New Issue