diff --git a/Utilities/Geometries/MATLAB/Functions/rcvGeo.m b/Utilities/Geometries/MATLAB/Functions/rcvGeo.m new file mode 100644 index 0000000000000000000000000000000000000000..e8fbc59fd159f48b0b5b076b6119014bcf5af02e --- /dev/null +++ b/Utilities/Geometries/MATLAB/Functions/rcvGeo.m @@ -0,0 +1,457 @@ +%% Calculate positions for specific receiver geometries +% -------------------------------------------------- +% Author: Achilles Kappis +% e-mail: axilleaz@protonmail.com +% +% Date: 10/03/2024 (DD/MM/YYYY) +% +% Copyright: MIT +% -------------------------------------------------- +% Functionality: Calculate receiver positions for specific receiver +% geometries. +% -------------------------------------------------- +% Input +% +% gType [char/string]: The type of the receiver geometry. At the moment the +% available options are "Single", "ULA", "UCA" and +% "Log". +% +% nSens [numeric] (Optional): The number of sensors. This value is +% ignored for the "Single" geometry. +% [Default: 4]. +% +% elemDist [numeric] (Optional): The distance between the elements of +% the array configuration. For "Single" +% it is ignored, for "ULA" corresponds +% to the distance between successive +% elements, for "UCA" it corresponds to +% the radius of the circle on which the +% sensors are placed. [Default: 0.1]. +% +% dmaElemDist [numeric] (Optional): The distance between the elements +% of the directional configurations. +% The available configurations are +% "Figure-of-Eight", "Triangle" and +% "Box", where the first comprises of +% two elements with their axis +% parallel to the y-axis, the +% "Triangle" has an addition element +% with positive z-coordinate forming +% an equilateral triangle and the +% "Box" configuration consists of +% three "Figure-of-Eights", one on +% each axis. If the given value is +% not positive or left empty the +% respective output variables are +% empty. [Default: 0.05]. +% +% xOff [numeric] (Optional): An offset along the x-axis of the microphone +% configuration. [Default: 0]. +% +% yOff [numeric] (Optional): An offset along the y-axis of the microphone +% configuration. [Default: 0]. +% +% zOff [numeric] (Optional): An offset along the z-axis of the microphone +% configuration. [Default: 0]. +% +% -------------------------------------------------- +% Output +% +% omniPos [numeric]: This matrix contains the coordinates of the +% omnidirectional receivers and has dimensions Nx3, +% where N is the number of receivers and the other +% dimension corresponds to the x, y and z Cartesian +% coordinates. +% +% fig8Pos [numeric] (Optional): This matrix contains the coordinates of the +% Figure-of-Eight elements (two sensors that +% can be used to create a figure-of-eight +% pattern when combined as a first order +% DMA). The dimensions of the matrix are +% Nx3x2 where N is the number of sensors, the +% second dimension corresponds to the x, y +% and z Cartesian coordinates and the last +% dimension corresponds to the +% "Figure-of-Eight" elements. For the +% "Single" and "ULA" geometries, the elements +% are positioned so that their axis is +% parallel to the y-axis and for the "UCA" +% geometry, the elements axis passes through +% the origin of the geometry. +% +% triPos [numeric] (Optional): This matrix contains the coordinates of the +% "Triangle" geometry where in addition to the +% "figure-of-eight" elements, one more is +% placed in the middle of ther axis and higher +% in order to create an equilateral triangle. +% The matrix has dimensions Nx3x3 where N is +% the number of microphone positions, the +% second dimension corresponds to the x, y and +% z Cartesian coordinates and the third +% dimension to individual elements. For the +% "Single" and "ULA geometries, the axis of +% the "Figure-of-Eight" elements on the x-y +% plane is parallel to the y-axis. For the +% "UCA" geometry the axis passes through the +% origin of the geometry. For "Single" and +% "ULA", the first element is the one with +% positive y-coordinate and for "UCA" the +% furthest from the origin. For all +% configurations the element with postitive +% z-coordinate is the last. +% +% boxPos [numeric] (Optional): This matrix contains the coordinates of +% three "Figure-of-Eight" sensors (two sensors +% which can be combined as a first order DMA +% to create a figure-of-eight sensor). The +% dimensions of the matrix are Nx3x2x3 where N +% is the number of sensor positions, the +% second dimension corresponds to the x, y and +% z Cartesian coordinates. The third dimension +% corresponds to the elements of each +% figure-of-eight sensor with the first +% element being the one having a positive +% coordinate value. The last dimension +% corresponds to the figure-of-eight sensors +% with the first being parallel to the x-axis, +% the second to the y-axis and the third to +% the z-axis. +% +% box2DPos [numeric] (Optional): This matrix contains the coordinates of +% the elements of "boxPos" lying on the x-y +% plane (the z-axis elements are discarded). +% The dimensions of the matrix are Nx3x2x2. +% +% tetPos [numeric]: This matrix contains the positions of a normal +% tetrahedron, with four vertices (points) non-parallel +% to the Cartesian axes. The dimensions of the matrix are +% Nx3x4 where N is the number of measurement positions, +% the second dimension is the x, y, and z Cartesian +% coordinates and the last dimension corresponds to the +% elements of the configuration. The element indices are: +% 1) Positive x, y and z, 2) negative x, positive y, +% zero z, 3) negative x and y, zero z and 4) negative x, +% y and z. +% +% fig8Vec [numeric]: This is an Nx3 matrix containing the positions of +% the figure-of-eight elements like: +% Pos_1_Elem_1_x, Pos_1_Elem_2_y, Pos_1_Elem_1_z +% Pos_1_Elem_2_x, Pos_1_Elem_2_y, Pos_1_Elem_2_z +% Pos_2_Elem_1_x, Pos_2_Elem_2_y, Pos_2_Elem_1_z +% Pos_2_Elem_2_x, Pos_2_Elem_2_y, Pos_2_Elem_2_z +% . . . +% . . . +% . . . +% Pos_M_Elem_1_x, Pos_M_Elem_2_y, Pos_M_Elem_1_z +% Pos_M_Elem_2_x, Pos_M_Elem_2_y, Pos_M_Elem_2_z +% +% triVec [numeric]: This is an Nx3 matrix containing the positions of +% the triangle elements like: +% Pos_1_Elem_1_x, Pos_1_Elem_2_y, Pos_1_Elem_1_z +% Pos_1_Elem_2_x, Pos_1_Elem_2_y, Pos_1_Elem_2_z +% Pos_1_Elem_3_x, Pos_1_Elem_3_y, Pos_1_Elem_3_z +% Pos_1_Elem_1_x, Pos_1_Elem_2_y, Pos_1_Elem_1_z +% Pos_1_Elem_2_x, Pos_1_Elem_2_y, Pos_1_Elem_2_z +% Pos_1_Elem_3_x, Pos_1_Elem_3_y, Pos_1_Elem_3_z +% . . . +% . . . +% . . . +% Pos_M_Elem_1_x, Pos_M_Elem_2_y, Pos_M_Elem_1_z +% Pos_M_Elem_2_x, Pos_M_Elem_2_y, Pos_M_Elem_2_z +% Pos_M_Elem_3_x, Pos_M_Elem_3_y, Pos_M_Elem_3_z +% +% boxVec [numeric]: This is an Nx3 matrix containing the positions of +% the box configuration elements like: +% Pos_1_X_Elem_1_x, Pos_1_X_Elem_2_y, Pos_1_X_Elem_1_z +% Pos_1_X_Elem_2_x, Pos_1_X_Elem_2_y, Pos_1_X_Elem_2_z +% Pos_1_Y_Elem_1_x, Pos_1_Y_Elem_2_y, Pos_1_Y_Elem_1_z +% Pos_1_Y_Elem_2_x, Pos_1_Y_Elem_2_y, Pos_1_Y_Elem_2_z +% Pos_1_Z_Elem_1_x, Pos_1_Z_Elem_2_y, Pos_1_Z_Elem_1_z +% Pos_1_Z_Elem_2_x, Pos_1_Z_Elem_2_y, Pos_1_Z_Elem_2_z +% Pos_2_X_Elem_1_x, Pos_2_X_Elem_2_y, Pos_2_X_Elem_1_z +% Pos_2_X_Elem_2_x, Pos_2_X_Elem_2_y, Pos_2_X_Elem_2_z +% Pos_2_Y_Elem_1_x, Pos_2_Y_Elem_2_y, Pos_2_Y_Elem_1_z +% Pos_2_Y_Elem_2_x, Pos_2_Y_Elem_2_y, Pos_2_Y_Elem_2_z +% Pos_2_Z_Elem_1_x, Pos_2_Z_Elem_2_y, Pos_2_Z_Elem_1_z +% Pos_2_Z_Elem_2_x, Pos_2_Z_Elem_2_y, Pos_2_Z_Elem_2_z +% . . . +% . . . +% . . . +% Pos_N_X_Elem_1_x, Pos_N_X_Elem_2_y, Pos_N_X_Elem_1_z +% Pos_N_X_Elem_2_x, Pos_N_X_Elem_2_y, Pos_N_X_Elem_2_z +% Pos_N_Y_Elem_1_x, Pos_N_Y_Elem_2_y, Pos_N_Y_Elem_1_z +% Pos_N_Y_Elem_2_x, Pos_N_Y_Elem_2_y, Pos_N_Y_Elem_2_z +% Pos_N_Z_Elem_1_x, Pos_N_Z_Elem_2_y, Pos_N_Z_Elem_1_z +% Pos_N_Z_Elem_2_x, Pos_N_Z_Elem_2_y, Pos_N_Z_Elem_2_z +% +% box2DVec [numeric]: This is an Nx3 matrix containing the positions of +% the 2D box configuration elements like: +% Pos_1_X_Elem_1_x, Pos_1_X_Elem_2_y, Pos_1_X_Elem_1_z +% Pos_1_X_Elem_2_x, Pos_1_X_Elem_2_y, Pos_1_X_Elem_2_z +% Pos_1_Y_Elem_1_x, Pos_1_Y_Elem_2_y, Pos_1_Y_Elem_1_z +% Pos_1_Y_Elem_2_x, Pos_1_Y_Elem_2_y, Pos_1_Y_Elem_2_z +% Pos_2_X_Elem_1_x, Pos_2_X_Elem_2_y, Pos_2_X_Elem_1_z +% Pos_2_X_Elem_2_x, Pos_2_X_Elem_2_y, Pos_2_X_Elem_2_z +% Pos_2_Y_Elem_1_x, Pos_2_Y_Elem_2_y, Pos_2_Y_Elem_1_z +% Pos_2_Y_Elem_2_x, Pos_2_Y_Elem_2_y, Pos_2_Y_Elem_2_z +% . . . +% . . . +% . . . +% Pos_N_X_Elem_1_x, Pos_N_X_Elem_2_y, Pos_N_X_Elem_1_z +% Pos_N_X_Elem_2_x, Pos_N_X_Elem_2_y, Pos_N_X_Elem_2_z +% Pos_N_Y_Elem_1_x, Pos_N_Y_Elem_2_y, Pos_N_Y_Elem_1_z +% Pos_N_Y_Elem_2_x, Pos_N_Y_Elem_2_y, Pos_N_Y_Elem_2_z +% +% tetVec [numeric]: This is an Nx3 matrix containg the positions of the +% normal tetrahedron configuration like: +% Pos_1_Elem_1_x, Pos_1_Elem_1_y, Pos_1_Elem_1_z +% Pos_1_Elem_2_x, Pos_1_Elem_2_y, Pos_1_Elem_2_z +% Pos_1_Elem_3_x, Pos_1_Elem_3_y, Pos_1_Elem_3_z +% Pos_1_Elem_4_x, Pos_1_Elem_4_y, Pos_1_Elem_4_z +% . . . +% . . . +% . . . +% Pos_N_Elem_1_x, Pos_1_Elem_1_y, Pos_1_Elem_1_z +% Pos_N_Elem_2_x, Pos_1_Elem_2_y, Pos_1_Elem_2_z +% Pos_N_Elem_3_x, Pos_1_Elem_3_y, Pos_1_Elem_3_z +% Pos_N_Elem_4_x, Pos_1_Elem_4_y, Pos_1_Elem_4_z +% +% -------------------------------------------------- +% Notes +% +% -------------------------------------------------- +function [omniPos, fig8Pos, triPos, boxPos, box2DPos, tetPos, fig8Vec, triangleVec, boxVec, box2DVec, tetVec] = rcvGeo(gType, nSens, elemDist, dmaElemDist, xOff, yOff, zOff) + % ==================================================== + % Check for number of arguments + % ==================================================== + narginchk(1, 7); + nargoutchk(0, 11); + + % ==================================================== + % Validate input arguments + % ==================================================== + % Validate mandatory arguments + validateattributes(gType, {'char', 'string'}, {'scalartext', 'nonempty'}, mfilename, "Geometry type", 1); + validatestring(gType, ["Single", "ULA", "UCA", "Log"], mfilename, "Geometry type", 1); + + % Validate optional arguments + if nargin > 6 && ~isempty(zOff) + validateattributes(zOff, "numeric", {'scalar', 'real', 'nonnan', 'finite'}, mfilename, "Z offset of the microphone setup", 7); + else + zOff = 0; + end + + if nargin > 5 && ~isempty(yOff) + validateattributes(yOff, "numeric", {'scalar', 'real', 'nonnan', 'finite'}, mfilename, "Y offset of the microphone setup", 6); + else + yOff = 0; + end + + if nargin > 4 && ~isempty(xOff) + validateattributes(xOff, "numeric", {'scalar', 'real', 'nonnan', 'finite'}, mfilename, "X offset of the microphone setup", 5); + else + xOff = 0; + end + + if nargin > 3 && ~isempty(dmaElemDist) + validateattributes(dmaElemDist, "numeric", {'scalar', 'real', 'nonnan', 'finite'}, mfilename, "Directional configurations inter-element distance", 4); + else + dmaElemDist = 0.05; + end + + if nargin > 2 && ~isempty(elemDist) + validateattributes(elemDist, "numeric", {'scalar', 'real', 'positive', 'nonnan', 'finite'}, mfilename, "Array inter-element distance", 3); + else + elemDist = 0.1; + end + + if nargin > 1 && ~isempty(nSens) + validateattributes(nSens, "numeric", {'scalar', 'positive', 'real', 'nonnan', 'finite'}, mfilename, "Number of sensors", 2); + else + nSens = 4; + end + + + % ==================================================== + % Calculate omnidirectional sensor positions + % ==================================================== + % Omni sensor coordinates + switch lower(gType) + case "single" + omniPos = [0, 0, 0]; + case "ula" + omniPos = elemDist * ((nSens - 1)/2) * linspace(-1, 1, nSens); + omniPos = [omniPos; zeros(2, nSens)]'; + case "uca" + angsOfRot = (0:1/nSens:(1 - 1/nSens)) * 2 * pi; + [x, y] = pol2cart(angsOfRot, elemDist); + omniPos = [x; y; zeros(size(x))].'; + end + + % ==================================================== + % Calculate directional configurations sensor positions + % ==================================================== + % Check whether we have to calculate the directional configuration sensor position + if dmaElemDist <= 0 || isempty(dmaElemDist) + fig8Pos = []; fig8Vec = []; + triPos = []; triangleVec = []; + boxPos = []; boxVec = []; + box2DPos = []; box2DVec = []; + tetPos = []; tetVec = []; + + % Translate + if nargout > 0 + omniPos = omniPos + [xOff, yOff, zOff]; + end + return; + end + + % Figure-of-Eight positions + if nargout > 1 + switch lower(gType) + case {"single", "ula"} + fig8Pos = omniPos + [0, dmaElemDist/2, 0]; + fig8Pos = cat(3, fig8Pos, omniPos - [0, dmaElemDist/2, 0]); + case "uca" + [x, y] = pol2cart(angsOfRot, elemDist + dmaElemDist/2); + fig8Pos = [x; y; zeros(size(x))].'; + [x, y] = pol2cart(angsOfRot, elemDist - dmaElemDist/2); + fig8Pos = cat(3, fig8Pos, [x; y; zeros(size(x))].'); + end + end + + % Triangle positions + if nargout > 2 + triPos = cat(3, fig8Pos, omniPos + [0, 0, dmaElemDist * sqrt(3)/2]); % Add an element to create an equilateral triangle + triPos = triPos - [0, 0, dmaElemDist * tand(30)/2]; % Move it down so that the array centre will be at the coordinates of the omni microphone + end + + % Box positions + if nargout > 3 + switch lower(gType) + case {"single", "ula"} + % X-axis figure-of-eight elements + tempBox = omniPos + [dmaElemDist/2, 0, 0]; + tempBox = cat(3, tempBox, omniPos - [dmaElemDist/2, 0, 0]); + boxPos = cat(4, tempBox, fig8Pos); + case "uca" + % X-Y plane perpendicular to the radial direction + % figure-of-eight elements (the "X-axis figure-of-eight") + [x, y] = pol2cart(angsOfRot + deg2rad(90), dmaElemDist/2); + tempBox = [x; y; zeros(size(x))].'; + [x, y] = pol2cart(angsOfRot - deg2rad(90), dmaElemDist/2); + tempBox = cat(3, tempBox, [x; y; zeros(size(x))].'); + tempBox = omniPos + tempBox; + boxPos = cat(4, tempBox, fig8Pos); + end + + % Z-axis figure-of-eight elements (common to all geometries) + tempBox = omniPos + [0, 0, dmaElemDist/2]; + tempBox = cat(3, tempBox, omniPos - [0, 0, dmaElemDist/2]); + boxPos = cat(4, boxPos, tempBox); + end + + % 2D box positions + if nargout > 4 + box2DPos = boxPos; % Get all elements from box + box2DPos(:, :, :, 3:3:end) = []; % Get rid of the z-axis elements + end + + % Normal tetrahedron positions + if nargout > 5 + % Generate normal tetrahedron with lower face parallel to x-y plane + % and edge length sqrt(8/3) [from: https://en.wikipedia.org/wiki/Tetrahedron] + tempPos = [sqrt(8/9), 0, -1/3;... + -sqrt(2/9), sqrt(2/3), -1/3;... + -sqrt(2/9), -sqrt(2/3), -1/3;... + 0, 0, 1]; + + % Set edge length equal to the given DMA inter-element distance + tempPos = dmaElemDist * tempPos/(vecnorm(tempPos(1, :) - tempPos(2, :))); + + switch lower(gType) + case "ula" + % Rotate the tetrahedrals on the positive side of they-axis + rotIdx = sum(omniPos(:, 1) > 0); % Calculate how many tetrahedrals we have to rotate + + % Go through the tetrahedrals + for measPosIdx = size(omniPos, 1):-1:1 + % Calculate rotation angle + if measPosIdx > rotIdx + rot = 60; + else + rot = 0; + end + + % Rotate and position + tetPos(measPosIdx, :, :) = omniPos(measPosIdx, :) + tempPos * rotMat3dDeg(0, 0, rot); + end + case "uca" + % Get the angle of the positions on the circle + az = cart2sph(omniPos(:, 1), omniPos(:, 2), omniPos(:, 3)); + + % Rotate and position tetrahedrals + for measPosIdx = numel(az):-1:1 + tetPos(measPosIdx, :, :) = omniPos(measPosIdx, :) + tempPos * rotMat3dDeg(0, 0, -az(measPosIdx)); + end + case "single" + tetPos = omniPos + tempPos * rotMat3dDeg(30, 30, 0); + end + + if ~strcmpi(gType, "Single") + tetPos = permute(tetPos, [1, 3, 2]); + end + end + + % ==================================================== + % Translate geometries + % ==================================================== + if nargout > 0 + omniPos = omniPos + [xOff, yOff, zOff]; + end + + if nargout > 1 + fig8Pos = fig8Pos + [xOff, yOff, zOff]; + end + + if nargout > 2 + triPos = triPos + [xOff, yOff, zOff]; + end + + if nargout > 3 + boxPos = boxPos + [xOff, yOff, zOff]; + end + + if nargout > 4 + box2DPos = box2DPos + [xOff, yOff, zOff]; + end + + if nargout > 5 + tetPos = tetPos + [xOff, yOff, zOff]; + end + + % Provide a Nx3 "Figure-of-Eight" matrix + if nargout > 6 + fig8Vec = reshape(permute(fig8Pos, [3, 1, 2]), [], 3); + end + + % Provide a Nx3 "Triangle" matrix + if nargout > 7 + triangleVec = reshape(permute(triPos, [3, 1, 2]), [], 3); + end + + % Provide a Nx3 "Box" matrix + if nargout > 8 + boxVec = reshape(permute(boxPos, [3, 4, 1, 2]), [], 3); + end + + % Provide a Nx3 "2DBox" matrix + if nargout > 9 + box2DVec = reshape(permute(box2DPos, [3, 4, 1, 2]), [], 3); + end + + % Provide a Nx3 Tetrahedron matrix + if nargout > 10 + tetVec = reshape(permute(tetPos, [3, 1, 2]), [], 3); + end +end \ No newline at end of file diff --git a/Utilities/Geometries/MATLAB/Functions/srcGeo.m b/Utilities/Geometries/MATLAB/Functions/srcGeo.m new file mode 100644 index 0000000000000000000000000000000000000000..49dc22a77ad7989714d923edbe105b0a365302f6 --- /dev/null +++ b/Utilities/Geometries/MATLAB/Functions/srcGeo.m @@ -0,0 +1,266 @@ +%% Calculate positions for specific source geometries +% -------------------------------------------------- +% Author: Achilles Kappis +% e-mail: axilleaz@protonmail.com +% +% Date: 24/02/2024 (DD/MM/YYYY) +% +% Copyright: MIT +% -------------------------------------------------- +% Functionality: Calculate source positions for specific source geometries. +% -------------------------------------------------- +% Input +% +% gType [char/string]: The type of the source geometry. At the moment the +% available options are "Single", "Causal", +% "Anti-Causal", "Side", "Diagonal", "Diffuse" and +% "Directional", "Circle". +% +% srcLen [numeric]: It is used for all but the "Diffuse" and "Circle" +% geometries. For the "Single" geometry it specifies +% the x coordinate of the source and for all other +% geometries the length of the source array. +% +% originDist [numeric]: For the "Diffuse" and "Circle" geometries, this is +% the radius of the sphere on whose surface the +% sources are placed. For the rest of the geometries, +% it specifies the distance between the origin and +% the centre of the source array. +% +% ang [numeric] (Optional): It is used for the creation of "Diagonal" and +% "Diffuse" geometries (for the rest is ignored +% and can even be set to []). For "Diagonal" is +% the angle for which the source array is rotated +% (clockwise on the x-y plane) and for the +% "Diffuse" field it is the angle between sources +% on the sphere. Currently, this is the ONLY way +% to create a "Diffuse" field and so it is +% necessary for this geometry. The value is in +% degrees. +% [Diagonal default: 45, Diffuse default: 10] +% +% nSrc [numeric] (Optional): Used for all but the "Diffuse" and +% "Single" geometries and dictates the +% number of sources that will comprise the +% source array. It is ignored if "Diffuse" +% or "Single" is chosen as the gType. For +% the "directional" configuration, it +% defines the number of neighbouring sources +% picked for the dominant sources. +% [Default: 5] +% +% azim [numeric] (Optional): Used for the directional source configuration. +% It is a vector with the azimuthal angles, in +% degrees, of the "central" dominant sources. +% [Default: 90] +% +% elev [numeric] (Optional): Used for the directional source configuration. +% It is a vector with the elevation angles of +% the "central" dominant sources. It must have +% the same number of elements as the "azim" +% parameter. [Default: 0] +% +% SNR [numeric] (Optional): Used for the directional source configuration. +% It scales the dominant and background source +% strengths so that the resulting SNR at the +% origin is equal to the given value. +% [Default: Inf] +% +% -------------------------------------------------- +% Output +% +% sPos [numeric]: The matrix with the source coordinates in the Cartesian +% system. The matrix has dimensions Nx3, where N is the +% number of sources and the other dimension correspond to +% the x, y and z coodrinates. +% +% Q [numeric] (Optional): The source strengths for the directional source +% configuration. For all other configurations this +% parameter is empty. +% +% domIdx [numeric] (Optional): The dominant source indices in the +% directional source configuration. This may +% be used to extract the Cartesian coordinates +% of the dominant sources as well as their +% strenghts. +% +% +% -------------------------------------------------- +% Notes +% +% - The "Directional" source configuration creates a diffuse field +% configuration and from those sources, the ones that fall closer to the +% given azimuthal and elevation angles are picked to be the "dominant" +% ones. +% The "nSrc" parameter defines how many neighbouring sources will be +% used for the generation of the "directional" field. For example if one +% azimuthal and elevation angle is chosen and sourceLen is equal to 3, +% the source with azimuthal and elevation angles closer to the given ones +% will be chosen and two neighbouring (on the horizontal plane) will be +% chosen to act as "dominant" sources too. This allows for "distributed" +% sources to be used and decrease the coherency of a "directional" +% source. +% +% Dependencies: - ptsOnSphere(): For the calculation of the "Diffuse" and +% "Directional" geometry source positions. +% - rotMat3dDeg(): To rotate the source configurations. +% +% -------------------------------------------------- +function [sPos, Q, domIdx] = srcGeo(gType, srcLen, originDist, ang, nSrc, azim, elev, SNR) + % ==================================================== + % Check for number of arguments + % ==================================================== + narginchk(3, 8); + nargoutchk(0, 3); + + + % ==================================================== + % Validate input arguments + % ==================================================== + % Validate mandatory arguments + validateattributes(gType, {'char', 'string'}, {'scalartext', 'nonempty'}, mfilename, "Geometry type", 1); + validatestring(gType, ["Single", "Causal", "Anti-Causal", "AntiCausal", "Side", "Diagonal", "Diffuse", "Directional", "Circle"], mfilename, "Geometry type", 1); + + validateattributes(srcLen, "numeric", {'scalar', 'real', 'nonnan', 'nonnegative', 'finite', 'nonempty'}, mfilename, "Source length", 2); + validateattributes(originDist, "numeric", {'scalar', 'nonnegative', 'real', 'nonnan', 'finite', 'nonempty'}, mfilename, "Distance from the origin", 3); + + % Validate optional arguments + if nargin > 3 && ~isempty(ang) + validateattributes(ang, "numeric", {'scalar', 'real', 'nonnan', 'finite', 'nonempty'}, mfilename, "Angle value", 4); + + if strcmpi(gType, "Diffuse") + if ang == 0 + error("For the Diffuse geometry a non-zero value must be given for angle."); + elseif ang < 0 + warning("For the Diffuse geometry the angle must be positive. A negative value is given so its absolute value will be used."); + end + end + else + if strcmpi(gType, "Diffuse") || strcmpi(gType, "Directional") + ang = 10; + else + ang = 45; + end + end + + if nargin > 4 && ~isempty(nSrc) + validateattributes(nSrc, "numeric", {'scalar', 'real', 'nonnegative', 'nonnan', 'finite'}, mfilename, "Number of sources", 5); + else + nSrc = 5; + end + + if nargin > 5 && ~isempty(azim) + validateattributes(azim, "numeric", {'vector', 'real', 'nonnan', 'finite'}, mfilename, "Directional field azimuthal angles", 6); + + azim(azim > 180) = 360 - azim(azim > 180); + azim(azim < -180) = 360 + azim(azim < -180); + azim = deg2rad(azim); + else + azim = 0; + end + + if nargin > 6 && ~isempty(elev) + validateattributes(elev, "numeric", {'vector', 'real', 'nonnan', 'finite'}, mfilename, "Directional field elevation angles", 7); + + if numel(elev) ~= numel(azim) + error("Elevation angles vector must have the same number of elements as the azimuth angles vector parameter."); + end + + elev(elev > 90) = 90 - elev(elev > 90); + elev(elev < -90) = 90 + elev(elev < -90); + elev = deg2rad(elev); + else + elev = zeros(length(azim), 1); + end + + if nargin > 7 && ~isempty(SNR) + validateattributes(SNR, "numeric", {'scalar', 'real', 'nonnan'}, mfilename, "Directional field SNR", 8); + else + SNR = Inf; + end + + + % ==================================================== + % Calculate source positions + % ==================================================== + if strcmpi(gType, "Single") + sPos = [srcLen, originDist, 0]; + domIdx = []; + elseif sum(strcmpi(gType, ["Causal", "Anti-Causal", "AntiCausal", "Side", "Diagonal"])) > 0 + if srcLen ~= 0 + sPos = [linspace(-srcLen/2, srcLen/2, nSrc); originDist * ones(1, nSrc); zeros(1, nSrc)].'; + else + sPos = [srcLen, originDist, 0]; + end + + if nSrc ~= 0 + switch lower(gType) + case {"anti-causal", "anticausal"} + sPos = sPos * rotMat3dDeg(0, 0, 180); % Rotate 180 degrees (clockwise) + case "side" + sPos = sPos * rotMat3dDeg(0, 0, 90); % Rotate 90 degrees (clockwise) + case "diagonal" + sPos = rotMat3dDeg(0, 0, -90 + ang) * sPos.'; % Rotate specified degrees (clockwise) + sPos = sPos.'; + end + end + domIdx = []; + elseif sum(strcmpi(gType, ["Diffuse", "Directional"]) > 0) + sPos = ptsOnSphere(originDist, ang); % Calculate the source positions on the surface of the sphere + + if strcmpi(gType, "Directional") && nSrc ~= 0 + % Get angles of sources + [az, el] = cart2sph(sPos(:, 1), sPos(:, 2), sPos(:, 3)); + + domIdx = []; + % Go through the angles + for angIdx = length(azim):-1:1 + closestVal = abs(elev(angIdx) - el); + elIdx = find((closestVal - 1e-5) <= min(closestVal)); + + closestVal = abs(azim(angIdx) - az(elIdx)); + closestIdx = elIdx((closestVal - 1e-5) <= min(closestVal)); + closestIdx = closestIdx(1); % Make sure to get only one index if there are two closest value candidates + + % Take care of "distributed" sources + if nSrc == 1 + domIdx = cat(1, domIdx, closestIdx); + else + for distSrcIdx = -ceil(nSrc/2)+1:floor(nSrc/2) + tempIdx = closestIdx + distSrcIdx; + + if tempIdx > max(elIdx) + tempIdx = min(elIdx) + mod(tempIdx, max(elIdx)) - 1; + elseif tempIdx < min(elIdx) + tempIdx = max(elIdx) + mod(tempIdx, -min(elIdx)) + 1; + end + domIdx = cat(1, domIdx, tempIdx); + end + end + end + domIdx = unique(domIdx); % Make sure there's no duplicates + else + domIdx = []; + end + elseif strcmpi(gType, "circle") + sPos = (0:2 * pi/nSrc:2 * pi - 1/nSrc).'; + [sPosX, sPosY] = sph2cart(sPos, 0, originDist); + sPos = [sPosX, sPosY, zeros(size(sPosX))]; + domIdx = []; + end + + % Calculate source strengths + if nargout > 1 + % Set the amplitude of the "weak" sources + Q = (10^(-SNR/10))/sum(~ismember(1:size(sPos, 1), domIdx)) * ones(size(sPos, 1), 1); + + % Set the amplitude of the dominant sources + if ~isinf(SNR) + Q(domIdx) = 1./length(domIdx); + else + Q(domIdx) = 1; + end + else + Q = []; + end +end \ No newline at end of file diff --git a/Utilities/Geometries/MATLAB/Functions/virtMicGeo.m b/Utilities/Geometries/MATLAB/Functions/virtMicGeo.m new file mode 100644 index 0000000000000000000000000000000000000000..204fd87e4442043d07a1a5e2ebe1ff0930e47770 --- /dev/null +++ b/Utilities/Geometries/MATLAB/Functions/virtMicGeo.m @@ -0,0 +1,225 @@ +%% Calculate virtual microphone positions for specific geometries +% -------------------------------------------------- +% Author: Achilles Kappis +% e-mail: axilleaz@protonmail.com +% +% Date: 12/07/2024 (DD/MM/YYYY) +% +% Copyright: MIT +% -------------------------------------------------- +% Functionality: Calculate virtual microphone positions for specific +% geometries. +% -------------------------------------------------- +% Input +% +% gType [char/string]: The type of the source geometry. At the moment the +% available options are "Single", "Array", "Grid", +% "Cube" and "Dual". The last is an arrangement with 2 +% sensors "xLen" apart that can be translated and +% rotated. +% +% +% xLen [numeric] (Optional): The width of the "Array" geometry, the length +% of the side of the "Grid" and "Cube geometries +% and the distance between positions in "Dual" +% geometry. +% [Default: Array/Grid/Cube - 1, Dual - 0.2]. +% +% xOff [numeric] (Optional): X-axis offset. [Default: 0]. +% +% yOff [numeric] (Optional): Y-axis offset. [Default: 0]. +% +% zOff [numeric] (Optional): Z-axis offset. [Default: 0]. +% +% nSens [numeric] (Optional): The number of sensors in the array. +% For the "Grid" geometry the square root of +% "nSens" must be an integer value. For the +% "Cube" geometry the cubic root of "nSens" +% must be an integer value. This parameter is +% ignored for the "Single" geometry. +% [Default: Array/Grid - 100, Cube - 1000]. +% +% orient [char/string/numeric] (Optional): The orientation of the +% geometry. For the "Single" and +% "Cube" geometries this value is +% ignored. For the "Dual" +% configuration only numerical +% values are used to declare the +% rotation of the configuration +% with respect to the local z-axis +% in degrees. For the rest of the +% configurations, the string +% values allowed are "XY", "XZ" +% and "YZ" and declare the plane +% on which the configuration will +% be placed. [Default: "XY" | 0]. +% +% -------------------------------------------------- +% Output +% +% vPos [numeric]: The matrix with the source coordinates in the Cartesian +% system. The matrix has dimensions Nx3, where N is the +% number of sources and the other dimension correspond to +% the x, y and z coodrinates. +% +% vPosMesh [numeric]: The matrix has dimensions sqrt(N)xsqrt(N)x3 where N +% is the number of sources. It contains the x, y and z +% Cartesian coordinates for each position in a +% rectangular grid with dimensions sqrt(N)xsqrt(N). +% This is used only with the "Grid" geometry, otherwise +% it is empty. +% +% -------------------------------------------------- +% Notes +% +% - Dependencies: rotMat3dDeg(): To rotate sources. +% +% - The position of the geometry is calculated like this: First the +% geometry is placed on the x-y plane, then rotated to match the +% orientation defined by the "orient" argument and finally, the offsets +% on each axis are applied. +% -------------------------------------------------- +function [vPos, vPosMesh] = virtMicGeo(gType, xLen, xOff, yOff, zOff, nSens, orient) + % ==================================================== + % Check for number of arguments + % ==================================================== + narginchk(1, 7); + nargoutchk(0, 2); + + % ==================================================== + % Validate input arguments + % ==================================================== + % Validate mandatory arguments + validateattributes(gType, {'char', 'string'}, {'scalartext', 'nonempty'}, mfilename, "Geometry type", 1); + validatestring(gType, ["Single", "Dual", "Array", "Grid", "Cube"], mfilename, "Geometry type", 1); + + % Validate optional arguments + if nargin > 1 && ~isempty(xLen) && ~strcmpi(gType, "Single") + validateattributes(xLen, "numeric", {'scalar', 'real', 'nonnan', 'finite', 'positive'}, mfilename, "Width of the geometry, or distance between positions in Dual geometry", 2); + else + if strcmpi(gType, "Dual") + xLen = 0.2; + else + xLen = 1; + end + end + + if nargin > 2 && ~isempty(xOff) + validateattributes(xOff, "numeric", {'scalar', 'real', 'nonnan', 'finite'}, mfilename, "X-offset", 3); + else + xOff = 0; + end + + if nargin > 3 && ~isempty(yOff) + validateattributes(yOff, "numeric", {'scalar', 'real', 'nonnan', 'finite'}, mfilename, "Y-offset", 4); + else + yOff = 0; + end + + if nargin > 4 && ~isempty(zOff) + validateattributes(zOff, "numeric", {'scalar', 'real', 'nonnan', 'finite'}, mfilename, "Z-offset", 5); + else + zOff = 0; + end + + if nargin > 5 && ~isempty(nSens) + if ~strcmpi(gType, "Single") + validateattributes(nSens, "numeric", {'scalar', 'real', 'nonnegative', 'finite', 'nonnan'}, mfilename, "Number of sensors", 6); + end + else + if sum(strcmpi(gType, ["Array", "Grid"])) > 0 + nSens = 1e2; + elseif strcmpi(gType, "Cube") + nSens = 1e3; + else + nSens = 0; + end + end + + if nargin > 6 && ~isempty(orient) + if ischar(orient) || isstring(orient) + validateattributes(orient, {'char', 'string'}, {'scalartext', 'nonempty'}, mfilename, "Geometry orientation", 7); + validatestring(orient, ["XY", "XZ", "YZ"], mfilename, "Geometry orientation", 7); + else + if ~strcmpi(gType, "Dual") + error("Planes can be provided for the orientation of only the Array and Grid geometries."); + end + + validateattributes(orient, {'numeric'}, {'scalar', 'nonempty', 'nonnan', 'finite', 'real'}, mfilename, "Geometry orientation", 7); + end + else + if strcmpi(gType, "Dual") + orient = 0; + else + orient = "XY"; + end + end + + % ==================================================== + % Check we have all needed parameters and they have + % acceptable values + % ==================================================== + + % For "Grid" geometry the square root of the number of sensors must be an integral value + if strcmpi(gType, "Grid") && mod(sqrt(nSens), 1) ~= 0 + error("For the Grid geometry, the square root of the number of sensors must be an integer."); + elseif strcmpi(gType, "Cube") && mod(nthroot(nSens, 3), 1) ~= 0 + error("For the Cube geometry, the cubic root of the number of sensors must be an integer (a tolerance of 1e-6 has been used for the calculation)."); + end + + % ==================================================== + % Calculate virtual microphone positions + % ==================================================== + switch lower(gType) + case "single" + % Create a single point at offset coordinates and return + vPos = [xOff, yOff, zOff]; + vPosMesh = []; + return; + case "dual" + vPos = [-xLen/2, 0, 0; ... + xLen/2, 0, 0]; + case "array" + vPos = linspace(-xLen/2, xLen/2, nSens); + vPos = [vPos(:), zeros(numel(vPos), 2)]; + case "grid" + vPos = linspace(-xLen/2, xLen/2, sqrt(nSens)); + [x, y] = meshgrid(vPos, vPos); + vPos = [x(:), y(:), zeros(numel(x), 1)]; + case "cube" + vPos = linspace(-xLen/2, xLen/2, round(nSens^(1/3))); + [x, y, z] = meshgrid(vPos, vPos, vPos); + vPos = [x(:), y(:), z(:)]; + otherwise + error("Oops... something went wrong here... not known geometry...!!!"); + end + + % Rotate + if ~strcmpi(gType, "Cube") + if isnumeric(orient) + vPos = vPos * rotMat3dDeg(0, 0, -orient); + elseif strcmpi(orient, "XZ") + if strcmpi(gType, "array") + vPos = vPos * rotMat3dDeg(0, 0, 90); + else + vPos = vPos * rotMat3dDeg(90, 0, 0); + end + elseif strcmpi(orient, "YZ") + vPos = vPos * rotMat3dDeg(0, 90, 0); + end + end + + % Translate + if ~strcmpi(gType, "Cube") + vPos = vPos + [xOff, yOff, zOff]; + end + + % For "Grid" and "Cube" geometry provide the coordinates in a "mesh format" + if strcmpi(gType, "Grid") + vPosMesh = reshape(vPos, sqrt(nSens), [], 3); + elseif strcmpi(gType, "Cube") + vPosMesh = reshape(vPos, round(nSens^(1/3)), round(nSens^(1/3)), [], 3); + else + vPosMesh = []; + end +end \ No newline at end of file