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

Bring back arrManTd.m after renaming

parent ce4f60a7
No related branches found
No related tags found
1 merge request!12Implementation of first order differential arrays in the time domain
%% 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 "On the Design of Time-Domain
% Differential Microphone Arrays" 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
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment