1function [f, m_db] = mls_freq_resp(id) 2%% Measure frequency response with MLS test signal 3% 4% [f, m] = mls_freq_resp(id) 5% 6% Input parameters 7% id - A string identifier for test case. An id 'selftest' is for special 8% usage. It calculates response of filtered MLS signal and computes 9% the measurement vs. known. Deviation is reported as error. It can 10% be useful if the internal MLS measurement parameters are adjusted. 11% 12% Output parameters 13% f - Frequency vector in Hz 14% m - Measured magnitude responses in dB 15% 16% Configuration (edit these): 17% mls_play_config.txt 18% mls_rec_config.txt 19% 20% The script will return also a text CSV format file with name mls-<id>.txt. 21% 22 23% SPDX-License-Identifier: BSD-3-Clause 24% 25% Copyright (c) 2018-2020, Intel Corporation. All rights reserved. 26% 27% Author: Seppo Ingalsuo <seppo.ingalsuo@linux.intel.com> 28 29%% Settings 30np = 1024; % Number of frequency points to use 31f_lo = 100; % Lower frequency limit for analysis 32f_hi = 20e3; % Upper frequency limit for analysis 33t_tot = 10e-3; % MLS analysis window length in s 34t_mls_s = 1.0; % MLS test signal length in s 35a_mls_db = -10; % MLS test signal amplitude in dB 36fs = 48e3; % Sample rate in Hz 37bits = 16; % Audio format to use (bits) 38fmt = 'S16_LE'; % Audio format to use (ALSA) 39dir = '/tmp'; % Directory for temporary files 40capture_level_max_db = -1; % Expected max. level 41capture_level_min_db = -30; % Expacted min. level 42 43%% Get device identifier to use 44if nargin < 1 45 id = 'unknown'; 46end 47 48if strcmp(id, 'selftest') 49 selftest = 1; 50 % Just some simulated speaker response to use as self test case 51 stb = [ 0.341762453, -0.915611126, 0.482465118, ... 52 1.017612317, -1.722527013, 0.711608745, ... 53 0.630608859, -0.813609935, 0.267690582, ]; 54 sta = [ 1.000000000, -3.931695128, 6.630812276, ... 55 -6.339248735, 3.800407709, -1.559376698, ... 56 0.619626250, -0.317702349, 0.097453451, ]; 57else 58 selftest = 0; 59end 60measfn = sprintf('mls-%s.wav', id); 61csvfn = sprintf('mls-%s.txt', id); 62 63%% Paths 64addpath('../../test/audio/test_utils'); 65 66%% MLS 67n_mls = round(fs*t_mls_s); 68mls = 10^(a_mls_db/20) * (2 * mlsp12(1, n_mls) - 1); 69mlsfn = 'mls-ref.wav'; 70audiowrite(mlsfn, mls, fs); 71 72%% Chip markers and parameters for find sync 73[x1, m1] = sync_chirp(fs, 'up'); 74[x2, m2] = sync_chirp(fs, 'down'); 75fnd.fs = fs; % Sample rate 76fnd.sm = 5; % Max seek from start 77fnd.em = 3; % Max seek from end 78fnd.idle_t = 2; % max idle in start or end 79fnd.mark_t = m1.t; % Marker length 80fnd.nf = 1; % One signal (amplitude) 81fnd.na = 1; % One signal (frequency) 82fnd.tl = t_mls_s; % Length of signal 83fnd.mt = 0.1; % Threshold length to issue error 84fnd.is = 0; % Ignore from start 85fnd.ie = 0; % Ignore from end 86 87%% Merge markers and MLS 88z = zeros(n_mls + m1.n + m2.n, 1); 89i1 = m1.n + 1; 90i2 = m1.n + n_mls; 91z(1:i1 - 1) = x1; 92z(i1:i2) = mls; 93z(i2 + 1:end) = x2; 94 95%% Get config 96rec_cfg = meas_remote_rec_config(fs, fmt); 97play_cfg = meas_remote_play_config; 98 99%% Capture MLS from all playback channel at time 100mixfn = 'mlsmix.wav'; 101recfn = 'recch.wav'; 102y = []; 103for i=1:play_cfg.nch 104 fprintf('\n'); 105 fprintf('Measure playback channel %d\n', i); 106 if selftest 107 tz =zeros(2*fs+length(z), 1); % Pad 2s 108 tz(fs:fs+length(z)-1) = z; 109 t = filter(stb, sta, tz); % Filter with test response 110 r = t(:) * ones(1, rec_cfg.nch); % Copy to all channels 111 else 112 x = zeros(length(z), play_cfg.nch); 113 x(:,i) = z; 114 mixdfn = sprintf('%s/%s', dir, mixfn); 115 audiowrite(mixdfn, x, fs, 'BitsPerSample', bits); 116 copy_playback(mixdfn, play_cfg); 117 tcap = floor(6 + t_mls_s); % Capture for MLS +6s 118 remote_capture(recfn, rec_cfg, tcap); 119 pause(1); 120 remote_play(mixfn, play_cfg); 121 pause(3); 122 r = get_recording(recfn, rec_cfg); 123 end 124 [d, nt] = find_test_signal(r(:,1), fnd); 125 if isempty(d) 126 figure; 127 sr = size(r); 128 ts = (0:sr(1)-1)/fs; 129 plot(ts, r(:,1)); 130 grid on; 131 xlabel('Time (s)'); 132 ylabel('Sample value'); 133 title('Captured audio test waveform'); 134 fprintf(1, 'Error: check the plot for skew in capture/playback.\n'); 135 f = []; 136 m_db = []; 137 return 138 end 139 for j = 1:rec_cfg.nch 140 y(:, rec_cfg.nch*(i-1) + j) = r(d:d + nt -1, j); 141 end 142 m = 20*log10(max(abs(r))); 143 fprintf('Peak levels for capture channels (dB):'); 144 for j = 1:rec_cfg.nch 145 fprintf(' %5.1f', m(j)); 146 end 147 fprintf('\n'); 148 if max(m) > capture_level_max_db 149 fprintf('Warning: The recording level is too loud.\n'); 150 end 151 if min(m) < capture_level_min_db 152 fprintf('Warning: The recording level is too silent.\n'); 153 end 154 155end 156audiowrite(measfn, y, fs, 'BitsPerSample', bits); 157fprintf('\n'); 158fprintf('Done.\n'); 159 160[f, m_db, b] = mls_calc_resp(csvfn, mlsfn, measfn, t_tot, np, f_lo, f_hi); 161 162[f, m_db] = apply_mic_calibration(f, m_db, rec_cfg); 163 164figure 165idx = find(f>1e3, 1, 'first') - 1; 166m_db_align = m_db - m_db(idx); 167semilogx(f, m_db_align); 168ax=axis(); axis([f_lo f_hi ax(3:4)]); 169grid on; 170xlabel('Frequency (Hz)'); 171ylabel('Magnitude (dB)'); 172 173if selftest 174 title('Measured vs. reference response'); 175 h = freqz(stb, sta, f, fs); 176 ref_db = 20*log10(abs(h)); 177 ref_db_align = ref_db - ref_db(idx); 178 hold on; 179 plot(f, ref_db_align, 'r--'); 180 hold off; 181 idx = find(f < f_hi); 182 idx = find(f(idx) > f_lo); 183 m_lin = 10.^(m_db_align(idx)/20); 184 ref_lin = 10.^(ref_db_align(idx)/20); 185 err_lin = m_lin - ref_lin; 186 e_rms = sqrt(mean(err_lin.^2)); 187 e_db = 20*log10(e_rms); 188 figure; 189 semilogx(f(idx), 20*log10(abs(err_lin))) 190 grid on; 191 xlabel('Frequency (Hz)'); 192 ylabel('Magnitude (dB)'); 193 title('Observed Error in self test'); 194 if e_db < -30 195 fprintf('Passed self test. '); 196 else 197 fprintf('Failed self test. '); 198 end 199 fprintf('Response RMS error is %4.1f dB.\n', e_db); 200end 201 202end 203 204function copy_playback(fn, cfg) 205 if cfg.ssh 206 cmd = sprintf('scp %s %s:%s/', fn, cfg.user, cfg.dir); 207 fprintf('Remote copy: %s\n', cmd); 208 system(cmd); 209 else 210 %cmd = sprintf('cp %s %s/', fn, cfg.dir); 211 %fprintf('Local copy: %s\n', cmd); 212 end 213end 214 215function y = get_recording(fn, cfg) 216 if cfg.ssh 217 cmd = sprintf('scp %s:%s/%s %s', cfg.user, cfg.dir, fn, fn); 218 fprintf('Remote copy: %s\n', cmd); 219 else 220 cmd = sprintf('cp %s/%s %s', cfg.dir, fn, fn); 221 fprintf('Local copy: %s\n', cmd); 222 end 223 system(cmd); 224 y = audioread(fn); 225 delete(fn); 226end 227 228function remote_play(fn, cfg) 229 if cfg.ssh 230 cmd = sprintf('ssh %s aplay -D%s %s/%s', cfg.user, cfg.dev, cfg.dir, fn); 231 fprintf('Remote play: %s\n', cmd); 232 else 233 cmd = sprintf('aplay -D%s %s/%s', cfg.dev, cfg.dir, fn); 234 fprintf('Local play: %s\n', cmd); 235 end 236 system(cmd); 237end 238 239function remote_capture(fn, cfg, t) 240 if cfg.ssh 241 cmd = sprintf('ssh %s arecord -q -D%s %s -d %d %s/%s &', ... 242 cfg.user, cfg.dev, cfg.fmt, t, cfg.dir, fn); 243 fprintf('Remote capture: %s\n', cmd); 244 else 245 cmd = sprintf('arecord -q -D%s %s -d %d %s/%s &', ... 246 cfg.dev, cfg.fmt, t, cfg.dir, fn); 247 fprintf('Local capture: %s\n', cmd); 248 end 249 system(cmd); 250end 251 252function play = meas_remote_play_config() 253 source mls_play_config.txt; 254 fprintf('\nThe setttings for remote playback are\n'); 255 fprintf('Use ssh : %d\n', play.ssh); 256 fprintf('User : %s\n', play.user); 257 fprintf('Directory : %s\n', play.dir); 258 fprintf('Device : %s\n', play.dev); 259 fprintf('Channels : %d\n', play.nch); 260end 261 262function rec = meas_remote_rec_config(fs, fmt) 263 source mls_rec_config.txt; 264 rec.fmt = sprintf('-t wav -c %d -f %s -r %d', ... 265 rec.nch, fmt, fs); 266 267 fprintf('\nThe setttings for remote capture are\n'); 268 fprintf('Use ssh : %d\n', rec.ssh); 269 fprintf('User : %s\n', rec.user); 270 fprintf('Directory : %s\n', rec.dir); 271 fprintf('Device : %s\n', rec.dev); 272 fprintf('format : %s\n', rec.fmt); 273 fprintf('Channels : %d\n', rec.nch); 274 fprintf('Calibration Data : %s\n', rec.cal); 275 276 if ~isempty(rec.cal) 277 if exist(rec.cal, 'file') 278 [rec.cf, rec.cm, rec.sens, rec.cs] = get_calibration(rec.cal); 279 else 280 error('The calibration file does not exist'); 281 end 282 else 283 rec.cf = []; 284 rec.cm = []; 285 rec.sens = []; 286 rec.cs = ''; 287 end 288end 289 290% The syntax of ASCII text calibration data is such that first line is 291% text string sometimes within "" or text line about Sensitivity. Such 292% lines are read and stored to description text. The sensitivty is extracted 293% and printed if success but currently it is not utilized by the code. 294% 295% The next lines are <frequency Hz> <magnitude decibels> for 296% measurement data for microphone. If there are more than two numbers 297% per line the other than two first columns are ignored. This code 298% applies inverse of calibration data for measured response calibration. 299 300function [f, m, sens, desc] = get_calibration(fn) 301 fh = fopen(fn, 'r'); 302 if fh < 0 303 error('Cannot open calibration data file'); 304 end 305 n = 1; 306 f = []; 307 m = []; 308 sens =[]; 309 desc = ''; 310 str = fgets(fh); 311 idx = findstr(str, '"'); 312 while length(idx) > 0 313 line = str(idx(1)+1:idx(2)-1); 314 desc = sprintf('%s%s ', desc, line); 315 str = fgets(fh); 316 idx = findstr(str, '"'); 317 end 318 if length(strfind(str, 'Sens')) 319 desc = str; 320 str = fgets(fh); 321 end 322 while str ~= -1 323 d=sscanf(str,'%f'); 324 f(n) = d(1); 325 m(n) = d(2); 326 n = n + 1; 327 str = fgets(fh); 328 end 329 330 % Strip possible linefeed from description end 331 if double(desc(end)) == 10 332 desc = desc(1:end-1); 333 end 334 fprintf('Calibration Info : %s\n', desc); 335 i1 = strfind(desc, 'Sens Factor ='); 336 i2 = strfind(desc, 'dB'); 337 if length(i1) == 1 && length(i2) == 1 338 sens = sscanf(desc(i1+13:i2-1), '%f'); 339 fprintf('Calibration Sens : %6.2f dB\n', sens); 340 end 341 342 fprintf('Calibration range: %.2f .. %.2f Hz\n', min(f), max(f)); 343 fprintf('Calibration range: %.2f .. %.2f dB\n', min(m), max(m)); 344 345 figure; 346 semilogx(f, m); 347 grid on; 348 title(desc); 349 xlabel('Frequency (Hz)'); 350 ylabel('Magnitude (dB)'); 351end 352 353function [x, seed] = mlsp12(seed, n) 354 355% Based on book Numerical Recipes in C, chapter 7.4 Generation of 356% random bits method II p. 298-299 example and 12 bit primitive 357% polynomial (12, 6, 4, 1, 0) 358 359 x = zeros(1,n); 360 ib1 = 2^(1-1); 361 ib4 = 2^(4-1); 362 ib6 = 2^(6-1); 363 ib12 = 2^(12-1); 364 mask = ib1 + ib4 + ib6 + ib12; 365 366 for i = 1:n 367 if bitand(seed, ib12) 368 seed = bitor(bitxor(seed, mask) * 2 , ib1); 369 x(i) = 1; 370 else 371 seed = seed * 2; 372 x(i) = 0; 373 end 374 end 375end 376 377%% Calibration apply function 378% Resample microphone calibration data into used grid and 379% then subtract calibration response from measured 380% response. 381 382function [cal_f, cal_m_db] = apply_mic_calibration(f, m_db, rec) 383 384 if length(rec.cm) > 0 385 if ~isvector(rec.cm) 386 error('Calibration can be for one channel only'); 387 end 388 mic_m_db = interp1(rec.cf, rec.cm, f, 'linear'); 389 nans = isnan(mic_m_db); 390 idx = find(nans == 0); 391 cal_f = f(idx); 392 cal_db = mic_m_db(idx); 393 s = size(m_db); 394 for i = 1:s(2) 395 cal_m_db(:,i) = m_db(idx,i) - cal_db(:); 396 end 397 else 398 cal_m_db = m_db; 399 cal_f = f; 400 end 401 402end 403
