1% bf = bf_design(bf) 2% 3% This script calculates beamformer filters with superdirective design 4% criteria. 5 6% SPDX-License-Identifier: BSD-3-Clause 7% 8% Copyright (c) 2020, Intel Corporation. All rights reserved. 9% 10% Author: Seppo Ingalsuo <seppo.ingalsuo@linux.intel.com> 11 12function bf = bf_design(bf) 13 14addpath('../../test/audio/test_utils'); 15addpath('../../test/audio/std_utils'); 16mkdir_check('plots'); 17mkdir_check('data'); 18 19if length(bf.steer_az) ~= length(bf.steer_el) 20 error('The steer_az and steer_el lengths need to be equal.'); 21end 22 23if max(bf.steer_az) > 180 || min(bf.steer_az) < -180 24 error('The steer_az angles need to be -180 to +180 degrees'); 25end 26 27if max(bf.steer_el) > 90 || min(bf.steer_el) < -90 28 error('The steer_el angles need to be -90 to +90 degrees'); 29end 30 31switch lower(bf.array) 32 case 'line' 33 bf = bf_array_line(bf); 34 case 'circular' 35 bf = bf_array_circ(bf); 36 case 'rectangle' 37 bf = bf_array_rect(bf); 38 case 'lshape' 39 bf = bf_array_lshape(bf); 40 case 'xyz' 41 bf = bf_array_xyz(bf); 42 otherwise 43 error('Invalid array type') 44end 45 46if isempty(bf.num_filters) 47 if isempty(bf.input_channel_select) 48 bf.num_filters = bf.mic_n; 49 else 50 bf.num_filters = length(bf.input_channel_select); 51 end 52end 53 54bf = bf_array_rot(bf); 55 56% The design function handles only single (az, el) value, so need to 57% loop every steer angle. 58bf.num_angles = length(bf.steer_az); 59all_az = bf.steer_az; 60all_el = bf.steer_el; 61all_array_id = bf.array_id; 62all_noiserot_fn = bf.noiserot_fn; 63all_sinerot_fn = bf.sinerot_fn; 64all_diffuse_fn = bf.diffuse_fn; 65all_random_fn = bf.random_fn; 66all_mat_fn = bf.mat_fn; 67 68for n = 1:bf.num_angles 69 bf.steer_az = all_az(n); 70 bf.steer_el = all_el(n); 71 bf.array_id = all_array_id{n}; 72 bf.noiserot_fn = all_noiserot_fn{n}; 73 bf.sinerot_fn = all_sinerot_fn{n}; 74 bf.diffuse_fn = all_diffuse_fn{n}; 75 bf.random_fn = all_random_fn{n}; 76 bf.mat_fn = all_mat_fn{n}; 77 bf = bf_one_design(bf); 78 w_all(:,:,n) = bf.w; 79end 80 81bf.steer_az = all_az; 82bf.steer_el = all_el; 83bf.array_id = all_array_id; 84bf.noiserot_fn = all_noiserot_fn; 85bf.sinerot_fn = all_sinerot_fn; 86bf.diffuse_fn = all_diffuse_fn; 87bf.random_fn = all_random_fn; 88bf.mat_fn = all_mat_fn; 89bf.w = w_all; 90 91end 92 93function bf = bf_one_design(bf) 94 95 96%% Defaults 97j = complex(0,-1); 98fs = bf.fs; 99N = 1024; 100N_half = N/2+1; 101f = (0:N/2)*fs/N'; 102phi_rad = (-180:180)*pi/180; 103phi_rad = phi_rad(1:end-1); 104n_phi = length(phi_rad); 105steer_az = bf.steer_az*pi/180; 106steer_el = bf.steer_el*pi/180; 107mu = ones(1,N_half) * 10^(bf.mu_db/20); 108%% Source at distance r 109[src_x, src_y, src_z] = source_xyz(bf.steer_r, steer_az, steer_el); 110 111%% Default frequency domain weights 112W = zeros(N_half, bf.num_filters); 113 114%% Coherence matrix, diffuse field 115% Equation 2.11 116Gamma_vv = zeros(N_half, bf.num_filters, bf.num_filters); 117for n=1:bf.num_filters 118 for m=1:bf.num_filters 119 % Equation 2.17 120 lnm = sqrt( (bf.mic_x(n) - bf.mic_x(m))^2 ... 121 +(bf.mic_y(n) - bf.mic_y(m))^2 ... 122 +(bf.mic_z(n) - bf.mic_z(m))^2); 123 % singular value decomposition 124 [U,S,V] = svd(sinc(2*pi*f*lnm/bf.c), 'econ'); 125 Gamma_vv(:,n,m) = U*S*V'; 126 end 127end 128 129%% Delays from source to each mic 130dt = delay_from_source(bf, src_x, src_y, src_z); 131dt = dt-min(dt); 132 133%% Create array vector 134tau0 = zeros(n_phi, bf.num_filters); 135A = zeros(N_half, n_phi, bf.num_filters); 136d = zeros(N_half, bf.num_filters); 137for n=1:bf.num_filters 138 % Equation 2.27 139 d(:,n) = exp(-j*2*pi*f*dt(n)); % Delays to steer direction 140 for ip = 1:n_phi; 141 phi = phi_rad(ip); 142 x_phi = bf.steer_r*cos(phi)*cos(steer_el); 143 y_phi = bf.steer_r*sin(phi)*cos(steer_el); 144 z_phi = bf.steer_r*sin(steer_el); 145 tau0(ip, n) = sqrt((x_phi-bf.mic_x(n))^2 ... 146 + (y_phi-bf.mic_y(n))^2 ... 147 + (z_phi-bf.mic_z(n))^2)/bf.c; 148 end 149end 150tau = tau0-min(min(tau0)); 151 152for n=1:bf.num_filters 153 for ip = 1:n_phi 154 % N_half x n_phi x Nm 155 A(:, ip, n) = exp(-j*2*pi*f*tau(ip, n)); 156 end 157end 158 159switch lower(bf.type) 160 case 'sdb' 161 % Superdirective 162 for iw = 1:N_half 163 % Equation 2.33 164 I = eye(bf.num_filters, bf.num_filters); 165 d_w = d(iw,:).'; 166 Gamma_vv_w = squeeze(Gamma_vv(iw, :, :)); 167 Gamma_vv_w_diagload = Gamma_vv_w + mu(iw)*I; 168 Gamma_vv_w_inv = inv(Gamma_vv_w_diagload); 169 num = Gamma_vv_w_inv * d_w; 170 denom1 = d_w' * Gamma_vv_w_inv; 171 denom2 = denom1 * d_w; 172 W_w = num / denom2; 173 W(iw, :) = W_w.'; 174 end 175 case 'dsb' 176 % Delay and sum 177 for iw = 1:N_half 178 % Equation 2.31 179 % W = 1/N*d 180 d_w = d(iw, :); 181 W_w = 1/bf.num_filters * d_w; 182 W(iw, :) = W_w.'; 183 end 184 otherwise 185 error('Invalid type, use SDB or DSB'); 186end 187 188%% Convert w to time domain 189W_full = zeros(N, bf.num_filters); 190W_full = W(1:N_half, :); 191for i=N_half+1:N 192 W_full(i,:) = conj(W(N_half-(i-N_half),:)); 193end 194bf.w = zeros(bf.fir_length, bf.num_filters); 195w_tmp = zeros(N, bf.num_filters); 196idx_max = zeros(bf.num_filters, 1); 197for i=1:bf.num_filters 198 h = real(fftshift(ifft(W_full(:,i)))); 199 w_tmp(:,i) = h; 200 idx_max(i) = find(h == max(h)); 201end 202 203center_idx = round(mean(idx_max)); 204start = center_idx - floor(bf.fir_length/2); 205switch lower(bf.type_filt) 206 case 'hann' 207 win0 = hann (bf.fir_length); 208 case 'hamming' 209 win0 = hamming(bf.fir_length); 210 case 'taylorwin' 211 win0 = taylorwin(bf.fir_length, bf.taylorwin_nbar, bf.taylorwin_sidelobe); 212 case 'chebwin' 213 win0 = chebwin(bf.fir_length, bf.chebwin_sidelobe); 214 case 'kaiser' 215 win0 = kaiser(bf.fir_length,bf.kaiser_beta); 216end 217for i=1:bf.num_filters 218 win = zeros(bf.fir_length, 1); 219 win_shift = center_idx - idx_max(i) - 1; 220 if (win_shift >= 0) 221 win(1:end-win_shift) = win0(win_shift+1:end); 222 else 223 win(-win_shift:end) = win0(1:end+win_shift+1); 224 end 225 bf.w(:,i) = w_tmp(start:start + bf.fir_length - 1, i) .* win; 226end 227 228%% Back to frequency domain to check spatial response 229W2_full = zeros(N, bf.num_filters); 230for i=1:bf.num_filters 231 % Zero pad 232 h2 = zeros(1,N); 233 h2(start:start + bf.fir_length - 1) = bf.w(:,i); 234 W2_full(:,i) = fft(h2); 235end 236W2 = W2_full(1:N_half, :); 237B2 = zeros(N_half, n_phi); 238for iw = 1:N_half 239 WT = (W2(iw,:)').'; 240 AS = squeeze(A(iw,:,:)).'; 241 WA = WT * AS; 242 B2(iw,:) = WA; 243end 244bf.resp_fa = B2'; 245bf.resp_angle = phi_rad * 180/pi; 246 247%% Directivity in diffuse field 248% Equation 2.18 249% DI(exp(j Omega) = 10*log10( abs(W^H d)^2 / (W^H Gamma_vv W)) 250bf.f = f; 251bf.di_db = zeros(1, N_half); 252for iw = 1:N_half 253 W_w = W2(iw,:).'; 254 d_w = d(iw,:).'; 255 Gamma_vv_w = squeeze(Gamma_vv(iw,:,:)); 256 W_wh = W_w'; 257 num = abs(W_wh * d_w)^2; 258 denom1 = W_wh * Gamma_vv_w; 259 denom2 = denom1 * W_w; 260 di = num / denom2; 261 bf.di_db(iw) = 10*log10(abs(di)); 262end 263 264 265%% White noise gain 266for iw = 1:N_half 267 % WNG = abs(^w^H d)^2/(w^H w); 268 W_w = W2(iw,:).'; 269 d_w = d(iw,:).'; 270 W_wh = W_w'; 271 num = abs(W_wh * d_w)^2; 272 denom = W_wh * W_w; 273 wng = num / denom2; 274 wng_db(iw) = 10*log10(abs(wng)); 275end 276bf.wng_db = wng_db; 277 278if bf.do_plots 279 %% Array 280 fh(1) = figure(bf.fn); 281 plot3(bf.mic_x(1), bf.mic_y(1), bf.mic_z(1), 'ro'); 282 hold on; 283 plot3(bf.mic_x(2:end), bf.mic_y(2:end), bf.mic_z(2:end), 'bo'); 284 plot3(src_x, src_y, src_z, 'gx'); 285 plot3([0 src_x],[0 src_y],[0 src_z],'c--') 286 for n=1:bf.num_filters 287 text(bf.mic_x(n), bf.mic_y(n), bf.mic_z(n) + 20e-3, ... 288 num2str(n)); 289 end 290 hold off 291 pb2 = bf.plot_box / 2; 292 axis([-pb2 pb2 -pb2 pb2 -pb2 pb2]); 293 axis('square'); 294 grid on; 295 xlabel('x (m)'); ylabel('y (m)'); zlabel('z (m)'); 296 view(-50, 30); 297 title(['Geometry ' bf.array_id], 'Interpreter','none'); 298 299 %% Coef 300 fh(2) = figure(bf.fn + 1); 301 plot(bf.w) 302 grid on; 303 xlabel('FIR coefficient'); ylabel('Tap value'); 304 title(['FIR filters ' bf.array_id], 'Interpreter','none'); 305 ch_legend(bf.num_filters); 306 307 %% Frequency responses 308 fh(3) = figure(bf.fn + 2); 309 f = linspace(0, bf.fs/2, 256); 310 h = zeros(256, bf.num_filters); 311 for i = 1:bf.num_filters 312 h(:,i) = freqz(bf.w(:,i), 1, f, bf.fs); 313 end 314 plot(f, 20*log10(abs(h))); 315 grid on 316 ylabel('Magnitude (dB)'); 317 xlabel('Frequency (Hz)'); 318 title(['FIR magnitude responses ' bf.array_id], 'Interpreter','none'); 319 ch_legend(bf.num_filters); 320 321 %% Group delays 322 fh(4) = figure(bf.fn + 3); 323 gd = zeros(256, bf.num_filters); 324 for i = 1:bf.num_filters 325 gd(:,i) = grpdelay(bf.w(:,i), 1, f, bf.fs); 326 end 327 plot(f, gd/bf.fs*1e6); 328 grid on 329 ylabel('Group delay (us)'); 330 xlabel('Frequency (Hz)'); 331 title(['FIR group delays ' bf.array_id], 'Interpreter','none'); 332 ch_legend(bf.num_filters); 333 334 %% DI 335 bf.fh(5) = figure(bf.fn + 4); 336 semilogx(bf.f(2:end), bf.di_db(2:end)) 337 xlabel('Frequency (Hz)'); ylabel('DI (dB)'); grid on; 338 legend('Suppression of diffuse field noise','Location','SouthEast'); 339 title(['Directivity Index ' bf.array_id], 'Interpreter','none'); 340 341 %% WNG 342 bf.fh(6) = figure(bf.fn + 5); 343 semilogx(bf.f(2:end), bf.wng_db(2:end)) 344 xlabel('Frequency (Hz)'); ylabel('WNG (dB)'); grid on; 345 legend('Attenuation of uncorrelated noise','Location','SouthEast'); 346 title(['White noise gain ' bf.array_id], 'Interpreter','none'); 347 drawnow; 348 349 %% 2D 350 bf.fh(7) = figure(bf.fn + 6); 351 colormap(jet); 352 phi_deg = phi_rad*180/pi; 353 imagesc(bf.f, bf.resp_angle, 20*log10(abs(bf.resp_fa)), [-30 0]); 354 set(gca,'YDir','normal') 355 grid on; 356 colorbar; 357 xlabel('Frequency (Hz)'); ylabel('Angle (deg)'); 358 title(['Spatial response ' bf.array_id], 'Interpreter','none'); 359 360 %% Polar 361 bf.fh(8) = figure(bf.fn + 7); 362 flist = [1000 2000 3000 4000]; 363 idx = []; 364 for i = 1:length(flist) 365 idx(i) = find(f > flist(i), 1, 'first'); 366 end 367 bf.resp_polar = abs(B2(idx,:)); 368 if exist('OCTAVE_VERSION', 'builtin') 369 polar(phi_rad, bf.resp_polar); 370 else 371 polarplot(phi_rad, bf.resp_polar); 372 end 373 legend('1 kHz','2 kHz','3 kHz','4 kHz'); 374 title(['Polar response ' bf.array_id], 'Interpreter','none'); 375end 376 377%% Create data for simulation 1s per angle 378 379if isempty(bf.sinerot_fn) 380 fprintf(1, 'No file for 360 degree sine source rotate specified\n'); 381else 382 rotate_sound_source(bf, bf.sinerot_fn, 'sine'); 383end 384 385if isempty(bf.noiserot_fn) 386 fprintf(1, 'No file for 360 degree random noise source rotate specified\n'); 387else 388 rotate_sound_source(bf, bf.noiserot_fn, 'noise'); 389end 390 391 392if isempty(bf.diffuse_fn) 393 fprintf(1, 'No file for diffuse noise field specified\n'); 394else 395 fprintf(1, 'Creating diffuse noise field...\n'); 396 fsi = 384e3; % Target interpolated rate 397 p = round(fsi / bf.fs); % Interpolation factor 398 fsi = p * bf.fs; % Recalculate high rate 399 ti = 1/fsi; % period at higher rate 400 t_add = 10.0/bf.c; % Additional signal time for max 20m propagation 401 t0 = bf.diffuse_t + t_add; % Total sine length per angle 402 n0 = floor(bf.fs * t0); 403 nt = floor(bf.fs * bf.diffuse_t); % Number samples output per angle 404 nti = p * nt; % Number samples output per angle at high rate 405 el = 0; 406 for az_deg = -160:20:180 % Azimuth plane only noise with sources 407 az = az_deg * pi/180; 408 [nx, ny, nz] = source_xyz(bf.steer_r, az, el); 409 dt = delay_from_source(bf, nx, ny, nz); 410 dn = round(dt / ti); 411 ns = rand(n0, 1) + rand(n0, 1) - 1; 412 nsi = interp(ns, p); 413 nmi = zeros(nti, bf.mic_n); 414 for j = 1:bf.mic_n 415 nmi(:,j) = nmi(:,j) + nsi(end-dn(j)-nti+1:end-dn(j)); 416 end 417 end 418 nm = nmi(1:p:end, :); 419 nlev = level_dbfs(nm(:,1)); 420 nm = nm * 10^((bf.diffuse_lev - nlev)/20); 421 myaudiowrite(bf.diffuse_fn, nm, bf.fs); 422end 423 424if isempty(bf.random_fn) 425 fprintf(1, 'No file for random noise specified\n'); 426else 427 fprintf(1, 'Creating random noise ...\n'); 428 nt = bf.fs * bf.random_t; 429 rn = rand(nt, bf.num_filters) + rand(nt, bf.num_filters) - 1; 430 431 nlev = level_dbfs(rn(:,1)); 432 rn = rn * 10^((bf.random_lev - nlev)/20); 433 myaudiowrite(bf.random_fn, rn, bf.fs); 434end 435 436if isempty(bf.mat_fn) 437 fprintf(1, 'No file for beam pattern simulation data specified.\n'); 438else 439 fprintf(1, 'Saving design to %s\n', bf.mat_fn); 440 bf_copy = bf; 441 bf.fh = []; % Don't save the large figures, this avoids a warning print too 442 mkdir_check(bf.data_path); 443 save(bf.mat_fn, 'bf'); 444 bf = bf_copy; 445end 446 447fprintf(1, 'Done.\n'); 448 449end 450 451%% Helper functions 452 453function rotate_sound_source(bf, fn, type); 454 fprintf(1, 'Creating 360 degree %s source rotate...\n', type); 455 fsi = 384e3; % Target interpolated rate 456 t_add = 10.0/bf.c; % Additional signal time for max 10m propagation 457 t_ramp = 5e-3; % 5 ms ramp to start and end of each angle tone 458 t_idle = 25e-3; % 25 ms idle to start and end of each angle tone 459 p = round(fsi / bf.fs); % Interpolation factor 460 fsi = p * bf.fs; % Recalculate high rate 461 ti = 1/fsi; % perid at higher rate 462 tt0 = bf.sinerot_t + t_add; % Total sine length per angle 463 nt = bf.fs * bf.sinerot_t; % Number samples output per angle 464 nt0 = floor(bf.fs * tt0); % Number samples output per angle 465 nti = p * nt; % Number samples output per angle at high rate 466 n_ramp = round(t_ramp * fsi); 467 n_idle = round(t_idle * fsi); 468 n_win = round(tt0 * fsi); 469 win = ones(n_win, 1); % Make a kind of hat shape window 470 win(n_idle:(n_idle + n_ramp - 1)) = linspace(0, 1, n_ramp); 471 win((end - n_idle - n_ramp + 1):(end - n_idle)) = linspace(1, 0, n_ramp); 472 win(1:n_idle) = zeros(1, n_idle); 473 win(end - n_idle + 1: end) = zeros(1, n_idle); 474 switch lower(type) 475 case 'sine' 476 si = multitone(fsi, bf.sinerot_f, bf.sinerot_a, tt0); 477 si = si .* win; 478 case 'noise' 479 [b, a] = butter(4, 2*[100 7000]/bf.fs); 480 ns0 = bf.sinerot_a * (2 * rand(round(tt0*bf.fs) + 1, 1) - 1); 481 nsf = filter(b, a, ns0); 482 si = interp(nsf, p); 483 si = si(1:n_win) .* win; 484 otherwise 485 error('no test signal type defined'); 486 end 487 test_az = (bf.sinerot_az_start:bf.sinerot_az_step:bf.sinerot_az_stop) * pi/180; 488 test_n = length(test_az); 489 test_el = zeros(1, test_n); 490 [test_x, test_y, test_z] = source_xyz(bf.steer_r, test_az, test_el); 491 td = zeros(test_n * nt, bf.mic_n); 492 for i = 1:length(test_az) 493 dt = delay_from_source(bf, test_x(i), test_y(i), test_z(i)); 494 dn = round(dt / ti); 495 mi = zeros(nti, bf.num_filters); 496 for j = 1:bf.mic_n 497 mi(:,j) = mi(:,j) + si(end-dn(j)-nti+1:end-dn(j)); 498 end 499 i1 = (i - 1) * nt + 1; 500 i2 = i1 + nt -1; 501 for j = 1:bf.mic_n 502 m = mi(1:p:end, j); 503 td(i1:i2, j) = m; 504 end 505 end 506 myaudiowrite(fn, td, bf.fs); 507end 508 509function [x, y, z] = source_xyz(r, az, el) 510 511x = r * cos(az) .* cos(el); 512y = r * sin(az) .* cos(el); 513z = r * sin(el); 514 515end 516 517function dt = delay_from_source(bf, src_x, src_y, src_z) 518 519dm = zeros(1,bf.num_filters); 520for n=1:bf.num_filters 521 dm(n) = sqrt((src_x - bf.mic_x(n))^2 ... 522 + (src_y - bf.mic_y(n))^2 ... 523 + (src_z - bf.mic_z(n))^2); 524end 525dt = dm/bf.c; 526 527end 528 529function ch_legend(n) 530switch n 531 case 2 532 legend('1', '2'); 533 case 3 534 legend('1', '2', '3'); 535 case 4 536 legend('1', '2', '3', '4'); 537 case 5 538 legend('1', '2', '3', '4', '5'); 539 case 6 540 legend('1', '2', '3', '4', '5', '6'); 541 case 7 542 legend('1', '2', '3', '4', '5', '6', '7'); 543 otherwise 544 legend('1', '2', '3', '4', '5', '6', '7', '8'); 545end 546end 547 548function myaudiowrite(fn, x, fs) 549[~, ~, ext] = fileparts(fn); 550if strcmpi(ext, '.raw') 551 s = size(x); 552 xq = zeros(s(1) * s(2), 1, 'int16'); 553 scale = 2^15; 554 xs = round(x * scale); 555 xs(xs > scale - 1) = scale -1; 556 xs(xs < -scale) = -scale; 557 for i = 1:s(2) 558 xq(i:s(2):end) = xs(:,i); 559 end 560 fh = fopen(fn, 'wb'); 561 fwrite(fh, xq, 'int16'); 562 fclose(fh); 563else 564 audiowrite(fn, x, fs); 565end 566end 567
