diff --git a/CHANGELOG.md b/CHANGELOG.md index 88832390ef4e8cb3fd7223f4a812d5ea8e46f47d..8c36569c5adf3bd2c8601efb7d2b2c1cca1fee78 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,3 +1,16 @@ +### v0.2.4 ### +**Virtual Sensing**\ +\* Corrected time-domain observation filter calculations.\ +\* Removed scaling from the cross-correlation calculations in the time-domain optimal observation filters.\ +\* Fixed bug in time-domain observation filter that screwed up dimensions of the filters when only one source was present. + +**Sound Fields**\ +\* Fixed bug in time-domain point source signal calculations that screwed up dimensions of the resulting arrays when only one source was present. + +**Utilities - Geometries**\ +\* The virtual microphone geometry generation function has been heavily modified and its interface has changed in a **backward incompatible** way. It is now more generic and allows for "arbitrary" size of each Cartesian dimension of the geometries as well as different number of sensors along each dimension. + + ### v0.2.3 ### **Virtual Sensing**\ \+ Added function to estimate the observation filters in the time-domain.\ diff --git a/README.md b/README.md index ed4f39a3f1f2a6ce28c3d921acb7c445e7ced3a3..8a3cb464712ce6e3f6322573fa503d59f17db195 100644 --- a/README.md +++ b/README.md @@ -86,7 +86,7 @@ The project is licensed under the ***MIT License***, which is a rather unrestric ## Versioning ## -The project uses [semantic versioning](https://semver.org/) and the current version is **v0.2.3**. +The project uses [semantic versioning](https://semver.org/) and the current version is **v0.2.4**. #### **Important** diff --git a/Sound Fields/MATLAB/Functions/ptSrcFieldTD.m b/Sound Fields/MATLAB/Functions/ptSrcFieldTD.m index 422cb6ead77fc2e3aae39345dc04881c21a7b164..ddfa6a6e4adebcdd80f61a9f53514813ef2d875a 100644 --- a/Sound Fields/MATLAB/Functions/ptSrcFieldTD.m +++ b/Sound Fields/MATLAB/Functions/ptSrcFieldTD.m @@ -3,7 +3,7 @@ % Author: Achilles Kappis % e-mail: axilleaz@protonmail.com % -% Date: 15/09/2024 (DD/MM/YYYY) +% Date: 17/09/2024 (DD/MM/YYYY) % % Copyright: MIT % -------------------------------------------------- @@ -170,7 +170,7 @@ function [rSig, rSigMean, rSigMtx, rSigMtxMean, Q] = ptSrcFieldTD(sPos, rPos, fs end % Sum source signals at each receiver position - rSig = squeeze(sum(rSigMtx, 3)); + rSig = reshape(sum(rSigMtx, 3), size(rSigMtx, 1), size(rSigMtx, 2), size(rSigMtx, 4)); % ==================================================== % Calculate output arguments diff --git a/Utilities/Geometries/MATLAB/Functions/virtMicGeo.m b/Utilities/Geometries/MATLAB/Functions/virtMicGeo.m index 564dfa39f7ed587750f32163456f71430b73b86b..b6947bea4181865388d190d9d660bc204e60f0d9 100644 --- a/Utilities/Geometries/MATLAB/Functions/virtMicGeo.m +++ b/Utilities/Geometries/MATLAB/Functions/virtMicGeo.m @@ -3,7 +3,7 @@ % Author: Achilles Kappis % e-mail: axilleaz@protonmail.com % -% Date: 16/07/2024 (DD/MM/YYYY) +% Date: 18/09/2024 (DD/MM/YYYY) % % Copyright: MIT % -------------------------------------------------- @@ -13,17 +13,35 @@ % 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. +% available options are "Single", "Array", "Cube" and +% "Dual". The last is an arrangement a distance apart +% that can be translated and rotated. % +% geoDim [numeric] (Optional): The dimensions of the geometry. This is a +% vector with three elements holding the +% dimensions of the setup in each Cartesian +% direction. Based on the "gType" some or all +% of them are used. For "gType" set either as +% "Array" or "Dual", the first value is used +% as the length of the "Array" or distance +% between the two position in "Dual". When +% "gType" is "Single" this argument is +% ignored. [Default: ones(3, 1)]. % -% 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]. +% nSens [numeric] (Optional): This can be either a real integer value +% denoting the number of sensors in the array +% or a vector with three elements denoting +% the number of sensor along each Cartesian +% dimension. If a dimension is equal to zero +% the corresponding value is ignored. If this +% is a scalar value the used must ensure that +% it is correctly divisible over the non-zero +% dimensions of the geometry uniformly (see +% Notes for info). For the "Array" geometry, +% only the first value (if "nSens" is a +% vector) is used and for the "Single" +% geometry the argument is ignored. +% [Default: 10 * ones(3, 1)]. % % xOff [numeric] (Optional): X-axis offset. [Default: 0]. % @@ -31,28 +49,10 @@ % % 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]. +% orient [numeric] (Optional): This is a real vector holding the rotation +% for the geometry in degrees, along each +% Cartesian axis. The rotations are performed +% clockwise. [Default: zeros(3, 1)]. % % -------------------------------------------------- % Output @@ -62,24 +62,26 @@ % 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. +% vPosMesh [numeric]: The matrix has dimensions nSens(1)xnSens(2)x3. It +% holds the x, y and z Cartesian coordinates for each +% position of the "Grid" geometry. If the geometry is +% otherwise not "Grid", this is an empty array. % % -------------------------------------------------- % Notes % -% - Dependencies: rotMat3d(): To rotate sources. +% - Dependencies: rotMat3d(): To rotate the geometry. +% +% - If the number of sources are provided as a single integer, the user +% must make sure that they can be uniformly distributed in the non-zero +% dimensions of the geometry. For example, for the "Cube" geometry with +% two non-zero dimensions, the square root of the number of sensors must +% be an integer and for three non-zero dimensions its third root must be +% an integer. For the "Array" configuration, there is not an issue as the +% setup is one dimensional. % -% - 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) +function [vPos, vPosMesh] = virtMicGeo(gType, geoDim, nSens, xOff, yOff, zOff, orient) % ==================================================== % Check for number of arguments % ==================================================== @@ -91,81 +93,67 @@ function [vPos, vPosMesh] = virtMicGeo(gType, xLen, xOff, yOff, zOff, nSens, ori % ==================================================== % Validate mandatory arguments validateattributes(gType, {'char', 'string'}, {'scalartext', 'nonempty'}, mfilename, "Geometry type", 1); - validatestring(gType, ["Single", "Dual", "Array", "Grid", "Cube"], mfilename, "Geometry type", 1); + validatestring(gType, ["Single", "Dual", "Array", "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); + if nargin > 1 && ~isempty(geoDim) && ~strcmpi(gType, "Single") + validateattributes(geoDim, "numeric", {'vector', 'real', 'nonnan', 'finite', 'nonempty', 'positive'}, mfilename, "Geometry dimensions", 2); + + if numel(geoDim) > 3 || numel(geoDim) == 2 + error("The dimensions argument must be either a scalar or a vector with three elements"); + end + + if isscalar(geoDim) + geoDim = geoDim * ones(3, 1); + end else - if strcmpi(gType, "Dual") - xLen = 0.2; - else - xLen = 1; + geoDim = ones(3, 1); + end + + if nargin > 2 && ~isempty(nSens) && sum(strcmpi(gType, ["Single", "Dual"])) == 0 + validateattributes(nSens, "numeric", {'vector', 'real', 'nonempty', 'nonnan', 'nonnegative', 'integer', 'finite'}, mfilename, "Number of sensors in the geometry", 3); + + if numel(nSens) > 3 || numel(nSens) == 2 + error("The number of sensors must be either a scalar or a vector with three elements"); + end + + if isscalar(nSens) && strcmpi(gType, "Cube") + if mod(nthroot(nSens, sum(geoDim ~= 0)), 1) ~= 0 + error("The number of sources is cannot be divided uniformly over the non-zero dimensions"); + else + nSens = nthroot(nSens, sum(geoDim ~= 0)) * double(geoDim ~= 0); + end end + elseif strcmpi(gType, "Dual") + nSens = [2, 0, 0]; + else + nSens = 10 * ones(3, 1); end - if nargin > 2 && ~isempty(xOff) - validateattributes(xOff, "numeric", {'scalar', 'real', 'nonnan', 'finite'}, mfilename, "X-offset", 3); + if nargin > 3 && ~isempty(xOff) + validateattributes(xOff, "numeric", {'scalar', 'real', 'nonnan', 'finite'}, mfilename, "X-offset", 4); else xOff = 0; end - if nargin > 3 && ~isempty(yOff) - validateattributes(yOff, "numeric", {'scalar', 'real', 'nonnan', 'finite'}, mfilename, "Y-offset", 4); + if nargin > 4 && ~isempty(yOff) + validateattributes(yOff, "numeric", {'scalar', 'real', 'nonnan', 'finite'}, mfilename, "Y-offset", 5); else yOff = 0; end - if nargin > 4 && ~isempty(zOff) - validateattributes(zOff, "numeric", {'scalar', 'real', 'nonnan', 'finite'}, mfilename, "Z-offset", 5); + if nargin > 5 && ~isempty(zOff) + validateattributes(zOff, "numeric", {'scalar', 'real', 'nonnan', 'finite'}, mfilename, "Z-offset", 6); 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 + if nargin > 6 && ~isempty(orient) && ~strcmpi(gType, "Single") + validateattributes(orient, {'numeric'}, {'vector', 'nonempty', 'nonnan', 'finite', 'real', 'numel', 3}, mfilename, "Geometry rotation around the Cartesian axes", 7); else - if strcmpi(gType, "Dual") - orient = 0; - else - orient = "XY"; - end + orient = zeros(3, 1); 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 @@ -173,53 +161,49 @@ function [vPos, vPosMesh] = virtMicGeo(gType, xLen, xOff, yOff, zOff, nSens, ori switch lower(gType) case "single" % Create a single point at offset coordinates and return - vPos = [xOff, yOff, zOff]; - vPosMesh = []; - return; + vPos = zeros(1, 3); case "dual" - vPos = [-xLen/2, 0, 0; ... - xLen/2, 0, 0]; + vPos = [-geoDim(1)/2, 0, 0; ... + geoDim(1)/2, 0, 0]; case "array" - vPos = linspace(-xLen/2, xLen/2, nSens); + vPos = linspace(-geoDim(1)/2, geoDim(1)/2, nSens(1)); 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); + % Make we don't get empty arrays when dimensions are 0 + if nSens(1) == 0 + x = 0; + else + x = linspace(-geoDim(1)/2, geoDim(1)/2, nSens(1)); + end + + if nSens(2) == 0 + y = 0; + else + y = linspace(-geoDim(2)/2, geoDim(2)/2, nSens(2)); + end + + if nSens(3) == 0 + z = 0; + else + z = linspace(-geoDim(3)/2, geoDim(3)/2, nSens(3)); + end + + [x, y, z] = meshgrid(x, y, z); 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 * rotMat3d(0, 0, -orient, "Degs"); - elseif strcmpi(orient, "XZ") - if strcmpi(gType, "array") - vPos = vPos * rotMat3d(0, 0, 90, "Degs"); - else - vPos = vPos * rotMat3d(90, 0, 0, "Degs"); - end - elseif strcmpi(orient, "YZ") - vPos = vPos * rotMat3d(0, 90, 0, "Degs"); - end - end + vPos = vPos * rotMat3d(-orient(1), -orient(2), -orient(3), "Degs"); % Translate - if ~strcmpi(gType, "Cube") - vPos = vPos + [xOff, yOff, zOff]; - end + vPos = vPos + [xOff, yOff, zOff]; - % 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 + % For "Cube" geometry provide the coordinates in a "mesh format" + if nargout > 1 && strcmpi(gType, "Cube") + vPosMesh = cat(3, x, y, z); + elseif nargout > 1 vPosMesh = []; end end \ No newline at end of file diff --git a/Virtual Sensing/README.md b/Virtual Sensing/README.md index fc21835acb4ddb4ce8ac5739e5a18fd2dae14c93..953cc48bea4b63aceb1ba013bef85963907cc09c 100644 --- a/Virtual Sensing/README.md +++ b/Virtual Sensing/README.md @@ -11,6 +11,7 @@ The current structure of the section is shown below - Remote Microphone Technique - Observation filter calculation - Frequency domain (tonal) + - Time domain (broadband and causal) ## Dependencies @@ -51,4 +52,5 @@ None. ## References -[1] A. Roure, A. Albarrazin, *"The remote microphone technique for active noise control"*, Inter-Noise and Noise-Con Congress and Conference Proceedings, Active99, Fort Lauderdale FL, pp. 1233-1244(12). +[1] A. Roure, A. Albarrazin, *"The remote microphone technique for active noise control"*, Inter-Noise and Noise-Con Congress and Conference Proceedings, Active99, Fort Lauderdale FL, pp. 1233-1244(12). +[2] S. J. Elliott, J. Cheer, *"Modeling local active sound control with remote sensors in spatially random pressure fields", Journal of the Acoustical Society of America, 137, pp. 1936-1946, (2015). diff --git a/Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/obsFiltEstTD.m b/Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/obsFiltEstTD.m index d8827b8322fc1ae625d62cf36626692ae28e9b19..4ac8d107c77a9bf587675003844dfca50b7f11fc 100644 --- a/Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/obsFiltEstTD.m +++ b/Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/obsFiltEstTD.m @@ -3,7 +3,7 @@ % Author: Achilles Kappis % e-mail: axilleaz@protonmail.com % -% Date: 14/09/2024 (DD/MM/YYYY) +% Date: 17/09/2024 (DD/MM/YYYY) % % Copyright: MIT % -------------------------------------------------- @@ -111,7 +111,7 @@ function [estPerMic, est, err, estMean, errMean] = obsFiltEstTD(m, O, e) % Sum the estimates of each monitoring microphone to get the estimated virtual microphone signals if nargout > 1 - est = squeeze(sum(estPerMic, 3)); + est = reshape(sum(estPerMic, 3), size(estPerMic, 1), size(estPerMic, 2), size(estPerMic, 4)); end % Calculate the error signals diff --git a/Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/obsFiltTD.m b/Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/obsFiltTD.m index 77ecd8770ab075a325f4e3b7291ba51b5a1a191b..43cf7fb9b69751d65d3dcf0eb7530fd916bc3463 100644 --- a/Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/obsFiltTD.m +++ b/Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/obsFiltTD.m @@ -3,7 +3,7 @@ % Author: Achilles Kappis % e-mail: axilleaz@protonmail.com % -% Date: 15/09/2024 (DD/MM/YYYY) +% Date: 17/09/2024 (DD/MM/YYYY) % % Copyright: MIT % -------------------------------------------------- @@ -211,27 +211,24 @@ function [O, Rme, Rmm, Ovec, RmeMtx, RmmMtx, condNum, mMtx, Omean, RmeMean, RmmM for mIdx = size(m, 2):-1:1 % Calculate the cross-correlations between virtual and monitoring microphones for eIdx = size(e, 2):-1:1 - [corr, lags] = xcorr(m(:, mIdx, jIdx), e(:, eIdx, jIdx)); - lIdx = find(lags == 0); + corr = xcorr(m(:, mIdx, jIdx), e(:, eIdx, jIdx), filtLen); - Rme(:, mIdx, eIdx, jIdx) = corr(lIdx:-1:lIdx-filtLen + 1); + Rme(:, mIdx, eIdx, jIdx) = corr(filtLen + 1:-1:2); end % Go through the monitoring microphones to calculate the monitoring microphone correlation matrices for mmIdx = mIdx:-1:1 % Auto-correlation matrices are Toeplitz symmetric if mIdx == mmIdx - [corr, lags] = xcorr(m(:, mmIdx, jIdx), m(:, mmIdx, jIdx)); - lIdx = find(lags == 0); + corr = xcorr(m(:, mmIdx, jIdx), m(:, mmIdx, jIdx), filtLen); - Rmm(:, :, mIdx, mmIdx, jIdx) = toeplitz(corr(lIdx:-1:lIdx - filtLen + 1)); + Rmm(:, :, mIdx, mmIdx, jIdx) = toeplitz(corr(filtLen + 1:-1:2)); else - [corr, lags] = xcorr(m(:, mIdx, jIdx), m(:, mmIdx, jIdx)); - lIdx = find(lags == 0); + corr = xcorr(m(:, mIdx, jIdx), m(:, mmIdx, jIdx), filtLen); % Cross-correlation matrices for iIdx = filtLen-1:-1:0 - Rmm(:, iIdx + 1, mIdx, mmIdx, jIdx) = corr(iIdx + (lIdx:-1:lIdx - filtLen + 1)); + Rmm(:, iIdx + 1, mIdx, mmIdx, jIdx) = corr(iIdx + (filtLen + 1:-1:2)); end Rmm(:, :, mmIdx, mIdx, jIdx) = squeeze(Rmm(:, :, mIdx, mmIdx, jIdx)).'; end @@ -243,8 +240,8 @@ function [O, Rme, Rmm, Ovec, RmeMtx, RmmMtx, condNum, mMtx, Omean, RmeMean, RmmM % Post-process cross- and auto-correlation matrices % ==================================================== % "Reshape" the data - RmeMtx = reshape(Rme, prod(size(Rme, [1, 2])), size(Rme, 3), size(Rme, 4)); - RmmMtx = reshape(permute(Rmm, [1, 3, 2, 4, 5]), prod(size(Rmm, [1, 3])), prod(size(Rmm, [2, 4])), size(Rmm, 5)); + RmeMtx = reshape(permute(Rme, [2, 1, 3, 4]), prod(size(Rme, [1, 2])), size(Rme, 3), size(Rme, 4)); + RmmMtx = reshape(permute(Rmm, [3, 1, 4, 2, 5]), prod(size(Rmm, [1, 3])), prod(size(Rmm, [2, 4])), size(Rmm, 5)); % ==================================================== % Calculate observation filters @@ -253,8 +250,8 @@ function [O, Rme, Rmm, Ovec, RmeMtx, RmmMtx, condNum, mMtx, Omean, RmeMean, RmmM Ovec(:, :, jIdx) = RmeMtx(:, :, jIdx).'/(RmmMtx(:, :, jIdx) + beta * eye(size(RmmMtx, 1))); end - % "Split" observation filter vector to observation filters per monitoring and virtual microphone - O = permute(reshape(permute(Ovec, [2, 1, 3]), filtLen, size(m, 2), size(e, 2), size(m, 3)), [3, 1, 2, 4]); + % "Split" observation filter vector to observation filters per monitoring and virtual microphone + O = permute(reshape(Ovec, size(Ovec, 1), size(m, 2), filtLen, size(Ovec, 3)), [1, 3, 2, 4]); % ==================================================== % Provide additional output arguments @@ -281,7 +278,7 @@ function [O, Rme, Rmm, Ovec, RmeMtx, RmmMtx, condNum, mMtx, Omean, RmeMean, RmmM Oopt = RmeMtxMean.'/(RmmMtxMean + beta * eye(size(RmmMtxMean))); % Reshape - Omean = permute(reshape(Oopt.', filtLen, size(m, 2), size(e, 2)), [3, 1, 2]); + Omean = permute(reshape(Oopt, size(Oopt, 1), size(m, 2), filtLen), [1, 3, 2]); end % Mean cross-correlations between monitoring and virtual microphones over trials/sound field realisations