Skip to content
Snippets Groups Projects

Update 3D rotations function interface

Merged Imported Update 3D rotations function interface
UpdateRotations into main
Merged Imported Achilles Kappis requested to merge UpdateRotations into main
Compare
and
11 files
+ 271
− 75
Compare changes
  • Side-by-side
  • Inline

Files

Signal Processing/Array Processing/MATLAB/Functions/arrManTD.m 0 → 100644
+ 189
− 0
%% Time-domain array manifold (steering vector)
% --------------------------------------------------
% Author: Achilles Kappis
% e-mail: axilleaz@protonmail.com
%
% Date: 02/11/2024 (DD/MM/YYYY)
%
% Copyright: MIT
% --------------------------------------------------
% Functionality: Calculate the time-domain equivalent of the array manifold
% (steering vector) of an arbitrary array.
% --------------------------------------------------
% Input
%
% mPos [numeric]: The positions of the sensors in Cartesian coordinates.
% This must be an 3xM matrix where M denotes the number of
% sensors and for each sensor the coordinates are in
% [x; y; z] format.
%
% dirs [numeric]: These are the directions for which the array manifold
% will be calculated. This must be an 2xN matrix where N is
% the number of directions and for each direction the
% azimuth and elevation/inclination is given. The values
% must be provided in radians. The directions must be
% pointing towards the incoming plane waves and not towards
% the array/origin.
%
% filtLen [numeric]: This argument represents the length of the DMA filter
% if this function is used in the context of time-domain
% DMA filter design. In general, defines the number of
% rows the output array will have per direction.
%
% delLen [numeric] (Optional): The length of the (fractional) delay
% filters used to calculate the "array
% manifold" in samples. If this value is not
% large enough to accommodate the largest
% delay between the point of reference of
% the array (see notes) and the most distant
% sensor, the length of the filters is set
% automatically to that number (internal to
% the "winSincFracDel()" function - see
% dependencies). [Default: 16].
%
% delWin [char/string/numeric] (Optional): The window that will be applied
% to the (fractional) delay
% filter. This can be
% "Rectangular", "Hann", "Hamming"
% or a positive scalar
% corresponding to the Kaiser
% alpha parameter. The string/char
% values are not case-sensitive.
% [Default: 3].
%
% fs [numeric] (Optional): The sampling rate in Hz. This value is used to
% calculate the final length of the filters in the
% returned "array manifold" matrix (the number of
% columns) and "delLen" in case it is not
% provided. For the number of columns in the
% returned "array manifold" see the output
% arguments. [Default: 8e3].
%
% c [numeric] (Optional): The speed of sound in m/s. [Default: 343].
%
% --------------------------------------------------
% Output
%
% am [numeric]: This is an LxKxMxN array, where L is equal to "filtLen", K
% is equal to filtLen + delLen + ceil(2 * max(abs(del)))
% where "del" is the delay between the sensors and the centre
% of mass of the array. M is the number of sensors and N is
% the number of directions.
%
% maxDel [numeric]: The maximum delay from the centre of mass to the
% sensors.
%
% amMtx [numeric]: This is an LMxKxN array, where LM is equal to L * M.
% This corresponds to the "vectorised" version of "am"
% with the "am" of each sensor stacked under each other.
%
% --------------------------------------------------
% Notes
%
% - The point of reference for the array is its centre of mass and it is
% calculated internally with dimensions equal to the mean of each
% dimension over all the sensor positions.
%
% - The implementation is based on "First-Order Differential Microphone
% Arrays from a Time-Domain Broadband Perspective" by Buchris, Cohen and
% Benesty.
%
% - Dependencies: * winSincFracDel() to calculate the fractional delay
% filters.
% --------------------------------------------------
function [am, maxDel, amMtx] = arrManTD(mPos, dirs, filtLen, delLen, delWin, fs, c)
% ====================================================
% Check for number of arguments
% ====================================================
narginchk(3, 7);
nargoutchk(0, 3);
% ====================================================
% Validate input arguments
% ====================================================
% Validate mandatory arguments
validateattributes(mPos, "numeric", {'2d', 'finite', 'nonnan', 'nonempty', 'real', 'nrows', 3}, mfilename, "Position of sensors", 1);
validateattributes(dirs, "numeric", {'2d', 'finite', 'nonnan', 'nonempty', 'real', 'nrows', 2}, mfilename, "Directions of interest", 2);
validateattributes(filtLen, "numeric", {'scalar', 'finite', 'nonnan', 'nonempty', 'positive', 'real'}, mfilename, "Number of rows of the 'array manifold'", 3);
% Optional arguments
if nargin > 3 && ~isempty(delLen)
validateattributes(delLen, "numeric", {'scalar', 'finite', 'nonnan', 'nonempty', 'positive', 'real'}, mfilename, "The length of the delay filters", 4);
else
delLen = 16;
end
if nargin < 5 || (nargin > 4 && isempty(delWin))
delWin = 3;
end
if nargin > 5 && ~isempty(fs)
validateattributes(fs, "numeric", {'scalar', 'nonempty', 'nonnan', 'finite', 'positive', 'real'}, mfilename, "Sampling rate", 6);
else
fs = 8e3;
end
if nargin > 6 && ~isempty(c)
validateattributes(c, "numeric", {'scalar', 'finite', 'nonempty', 'nonnan', 'positive', 'real'}, mfilename, "Speed of sound", 7);
else
c = 343;
end
% =============================================
% Calculate delays and filters
% =============================================
% Set the centre of mass to the origin
mPos = mPos - mean(mPos, 2);
% Get Cartesian coordinates of unit vectors pointing in angle directions
[x, y, z] = sph2cart(dirs(1, :), dirs(2, :), ones(size(dirs(1, :))));
% Calculate each sensor's delay with respect to the centre of mass in samples
del = -[x; y; z].' * fs * mPos/c;
% Calculate the total length of the filters
maxDel = max(abs(del), [], "all");
delLen = ceil(delLen + 2 * maxDel);
% Get the delay filters
for angIdx = size(del, 1):-1:1
for mIdx = size(del, 2):-1:1
h(:, mIdx, angIdx) = winSincFracDel(del(angIdx, mIdx), delLen, delWin);
end
end
% =============================================
% Form Sylvester matrices
% =============================================
% Turn warnings about the generation of Toeplitz matrices off
w = warning("off", "MATLAB:toeplitz:DiagonalConflict");
try
% Generate Sylvester matrices for each mic and angle
for angIdx = size(h, 3):-1:1
for mIdx = size(h, 2):-1:1
am(:, :, mIdx, angIdx) = toeplitz([h(:, mIdx, angIdx); zeros(filtLen - 1, 1)], zeros(filtLen, 1)).';
end
end
catch
% Make sure to turn the warnings on again if an error occurs
warning(w.state, w.identifier);
end
% Reset the warnings state
warning(w.state, w.identifier);
% =============================================
% Return additional output arguments
% =============================================
if nargout > 1
maxDel = ceil(maxDel);
end
if nargout > 2
amMtx = reshape(permute(am, [1, 3, 2, 4]), [], size(am, 2), size(am, 4));
end
end
\ No newline at end of file