Skip to content
Snippets Groups Projects
Commit 818bfacf authored by Achilles Kappis's avatar Achilles Kappis
Browse files

Add function to calculate the time domain equivalent of the array manifold

parent b3ccf5ce
No related branches found
No related tags found
1 merge request!10Update 3D rotations function interface
%% Time-domain array manifold (steering vector)
% --------------------------------------------------
% Author: Achilles Kappis
% e-mail: axilleaz@protonmail.com
%
% Date: 21/10/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.
%
% nRows [numeric]: The number of rows in the final "array manifold" matrix
% for each direction (see notes).
%
% 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 "nRows", K
% is equal to nRows + 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.
%
% 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 "nRows" argument represents the length of the DMA filter if this
% function is used in the context of time-domain DMA filter design.
%
% - 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, amMtx, maxDel] = arrManTD(mPos, dirs, nRows, 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(nRows, "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
% =============================================
w = warning("off", "MATLAB:toeplitz:DiagonalConflict");
for angIdx = size(h, 3):-1:1
for mIdx = size(h, 2):-1:1
am(:, :, mIdx, angIdx) = toeplitz([h(:, mIdx, angIdx); zeros(nRows - 1, 1)], zeros(nRows, 1)).';
end
end
warning(w.state, w.identifier);
% =============================================
% Return additional output arguments
% =============================================
if nargout > 1
amMtx = reshape(permute(am, [1, 3, 2, 4]), [], size(am, 2), size(am, 4));
end
if nargout > 2
maxDel = ceil(maxDel);
end
end
\ No newline at end of file
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment