diff --git a/Active Noise Control/Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/mulCoh.m b/Active Noise Control/Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/mulCoh.m
new file mode 100644
index 0000000000000000000000000000000000000000..d9c482d617521a7392f793da87735827296014cd
--- /dev/null
+++ b/Active Noise Control/Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/mulCoh.m
@@ -0,0 +1,81 @@
+%% Calculate the "multiple coherence"
+% --------------------------------------------------
+% Author: Achilles Kappis
+% e-mail: axilleaz@protonmail.com
+%
+% Date: 20/03/2024 (DD/MM/YYYY)
+%
+% Copyright: MIT
+% --------------------------------------------------
+% Functionality: Calculate "multiple coherence" between positions.
+% --------------------------------------------------
+% Input
+%
+% Px [numeric]: The input matrix.
+%
+% Py [numeric]: The output matrix. The number of rows must match the number
+% of columns of "Px".
+%
+% Sqq [numeric] (Optional): The source power spectral density matrix. This
+% must be a square matrix with each dimension
+% matching the number of columns of "Px".
+%
+% xRegFac [numeric] (Optional): A regularisation factor to be used in the
+% inversion of the cross spectral density
+% matrix of "Px". [Default: 0].
+%
+% --------------------------------------------------
+% Output
+%
+% coh [numeric]: The multiple coherence.
+%
+% --------------------------------------------------
+% Notes
+%
+% - The implementation is based on the paper: "A study of coherence between
+% virtual signal and physical signals in remote acoustic sensing" by
+% P. Zhang, S. Wang, H. Duan, J. Tao, H. Zou and X. Qiu.
+% --------------------------------------------------
+function coh = mulCoh(Px, Py, Sqq, xRegFac)
+ % ====================================================
+ % Check for number of arguments
+ % ====================================================
+ narginchk(2, 4);
+ nargoutchk(0, 2);
+
+ % ====================================================
+ % Validate input arguments
+ % ====================================================
+ % Validate mandatory arguments
+ validateattributes(Px, "numeric", {'nonempty', '2d'}, mfilename, "Input matrix", 1);
+ validateattributes(Py, "numeric", {'nonempty', '2d', 'ncols', size(Px, 2)}, mfilename, "Output matrix", 2);
+
+ % Validate optional arguments
+ if nargin > 2 && ~isempty(Sqq)
+ validateattributes(Sqq, "numeric", {'2d'}, mfilename, "Source strength(s)", 3);
+
+ if isscalar(Sqq) || isvector(Sqq)
+ Sqq = Sqq * eye(size(Px, 2));
+ else
+ validateattributes(Sqq, "numeric", {'2d', 'square', 'ncols', size(Px, 2)}, mfilename, "Source cross spectral density matrix", 3);
+ end
+ else
+ Sqq = eye(size(Px, 2));
+ end
+
+ if nargin > 3 && ~isempty(xRegFac)
+ validateattributes(xRegFac, "numeric", {'scalar', 'real', 'finite', 'nonnegative', 'nonnan'}, mfilename, "'Px' cross spectral density matrix regularisation value", 4);
+
+
+ % ====================================================
+ % Validate input arguments
+ % ====================================================
+ % Calculate power spectral density matrices
+ Sxx = Px * Sqq * Px';
+ Syy = Py * Sqq * Py';
+ Sxy = Px * Sqq * Py';
+
+ % Calculate multiple coherence
+ coh = diag(Sxy'/(Sxx + xRegFac * eye(size(Sxx, 1))) * Sxy)./diag(Syy);
+ coh = real(coh); % Make sure to get rid of residual imaginary parts
+end
\ No newline at end of file
diff --git a/Active Noise Control/Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/obsFilt.m b/Active Noise Control/Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/obsFilt.m
new file mode 100644
index 0000000000000000000000000000000000000000..913849a96b693e34fae409aceb28edd59ac55ae2
--- /dev/null
+++ b/Active Noise Control/Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/obsFilt.m
@@ -0,0 +1,163 @@
+%% Calculate optimal observation filters
+% --------------------------------------------------
+% Author: Achilles Kappis
+% e-mail: axilleaz@protonmail.com
+%
+% Date: 13/02/2024 (DD/MM/YYYY)
+%
+% Copyright: MIT
+% --------------------------------------------------
+% Functionality: Calculate the optimal, in the least squares sense,
+% observation filters for the Remote Microphone Technique.
+% --------------------------------------------------
+% Input
+%
+% Pe [numeric]: The transfer function from sources to the virtual
+% microphones. The dimensions of the matrix must be NxM,
+% where N is the number of virtual microphones and M the
+% number of sources.
+%
+% Pm [numeric]: The transfer function from sources to the monitoring
+% microphones. The dimensions of the matrix must be KxM,
+% where K is the number of monitoring microphones and M is
+% the number of sources.
+%
+% srcCsd [numeric] (Optional): The source cross spectral density matrix.
+% This must be a square (MxM) symmetric matrix
+% with the cross power spectral density of the
+% sources. [Default: eye(M)]
+%
+% regFacs [numeric] (Optional): The regularisation factors used for the
+% inversion of the cross-spectra of the
+% monitoring microphones. Can be a vector
+% with number of elements equal to K (number
+% of monitoring microphones) or a scalar
+% which will result in the same
+% regularisation factor for all microphones.
+% [Default: 0]
+%
+% --------------------------------------------------
+% Output
+%
+% oOpt [numeric]: The optimal observation filters for each virtual
+% microphone position. The dimensions of the matrix are
+% NxM, where N is the number of virtual positions and M the
+% number of monitoring microphones.
+%
+% Sme[numeric]: The cross-spectra between monitoring and virtual
+% microphones. The matrix has dimensions NxK where N is the
+% number of virtual microphones and M the number of
+% monitoring microphones.
+%
+% Smm[numeric]: The cross-spectra between monitoring microphones. The
+% matrix has dimensions KxK, where K is the number of
+% monitoring microphones.
+%
+% est [numeric]: The estimated complex pressure(s) at the virtual
+% microphone position(s).
+%
+% err [numeric]: The error of the estimated complex pressure(s).
+%
+% errSqr [numeric]: The sum of the squared errors of the estimated
+% complex pressure(s).
+%
+% normErrSqr [numeric]: The sum of the squared errors of the estimated
+% complex pressure(s) normalised to the squared
+% pressure(s) of the virtual microphones.
+%
+% See [numeric]: The normalisation factor (cross spectral density of the
+% error microphones) for each virtual microphone position.
+% "errSqr"./"See" gives "normErrSqr".
+%
+% condNum [numeric]: The condition number of the regularised monitoring
+% microphone cross spectrum matrix.
+%
+% --------------------------------------------------
+% Notes
+%
+% --------------------------------------------------
+function [oOpt, Sme, Smm, est, err, errSqr, normErrSqr, See, condNum] = obsFilt(Pe, Pm, srcCsd, regFacs)
+ % ====================================================
+ % Check for number of arguments
+ % ====================================================
+ narginchk(2, 4);
+ nargoutchk(0, 9);
+
+ % ====================================================
+ % Validate input arguments
+ % ====================================================
+ % Validate mandatory arguments
+ validateattributes(Pe, "numeric", {'nonnan', 'nonempty', 'finite'}, mfilename, "Virtual microphone pressure", 1);
+ validateattributes(Pm, "numeric", {'nonnan', 'nonempty', 'finite'}, mfilename, "Monitoring microphone pressure", 2);
+
+ % Validate optional arguments
+ if nargin > 2 && ~isempty(srcCsd)
+ validateattributes(srcCsd, "numeric", {'2d', 'nonnegative', 'nonnan', 'finite'}, mfilename, "Source cross spectral density matrix", 3)
+ else
+ srcCsd = eye(size(Pe, 2));
+ end
+
+ if nargin > 3 && ~isempty(regFacs)
+ validateattributes(regFacs, "numeric", {'vector', 'nonnegative', 'nonnan', 'nonempty', 'finite'}, mfilename, "Regularisation factors", 4);
+
+ % Make sure regFacs has correct dimensions
+ if ~isscalar(regFacs) && numel(regFacs) ~= size(Pm, 1)
+ error("Regularisation factors must be either a scalar or its length must match the number of monitoring microphones.");
+ end
+ else
+ regFacs = 0;
+ end
+
+
+ % ====================================================
+ % Calculate optimal filters
+ % ====================================================
+ % Calculate needed quantities
+ tmpVal = srcCsd * Pm';
+
+ Sme = Pe * tmpVal; % Virtual-Monitor mics cross-spectra
+ Smm = Pm * tmpVal; % Monitor-Monitor mics cross-spectra
+
+ if isscalar(regFacs)
+ regMat = diag(ones(size(Smm, 1), 1)) .* regFacs; % Regularisation matrix
+ else
+ regMat = diag(regFacs);
+ end
+
+ % Calcualte observation filters
+ invQty = Smm + regMat;
+ oOpt = Sme/invQty;
+
+ % Calculate pressure estimates
+ if nargout > 3
+ est = oOpt * Pm;
+ else
+ return;
+ end
+
+ % Calculate errors
+ if nargout > 4
+ err = Pe - est; % Estimation error
+ else
+ return;
+ end
+
+ if nargout > 5
+ errSqr = diag(err * srcCsd * err'); % Sum of squared estimation errors
+ else
+ return;
+ end
+
+ % Normalised squared errors
+ if nargout > 6
+ See = diag(Pe * srcCsd * Pe'); % Error microphone cross spectral density
+ normErrSqr = errSqr./See; % Normalised error
+ else
+ return
+ end
+
+ % Condition number of auto-correlation (power spectrum) matrix
+ if nargout > 8
+ condNum = cond(invQty);
+ end
+end
\ No newline at end of file
diff --git a/Active Noise Control/Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/rcvGeo.m b/Active Noise Control/Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/rcvGeo.m
new file mode 100644
index 0000000000000000000000000000000000000000..62874ce6cbf5b235da0627ef2a41e7ff67497489
--- /dev/null
+++ b/Active Noise Control/Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/rcvGeo.m
@@ -0,0 +1,458 @@
+%% 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 * rotMat3d(0, 0, -az(measPosIdx) + pi);
+ tetPos(measPosIdx, :, :) = omniPos(measPosIdx, :) + tempPos * rotMat3d(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/Active Noise Control/Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/srcGeo.m b/Active Noise Control/Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/srcGeo.m
new file mode 100644
index 0000000000000000000000000000000000000000..6677141848827623bba8f69af5a58bb954c1654c
--- /dev/null
+++ b/Active Noise Control/Virtual Sensing/Remote Microphone Technique/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 = rotMat3d(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/Active Noise Control/Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/tonalAnc.m b/Active Noise Control/Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/tonalAnc.m
new file mode 100644
index 0000000000000000000000000000000000000000..af7b80b23d2e3b6ebb20d2017b68ed2e74f728e8
--- /dev/null
+++ b/Active Noise Control/Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/tonalAnc.m
@@ -0,0 +1,328 @@
+%% Tonal control in the frequency domain
+% --------------------------------------------------
+% Author: Achilles Kappis
+% e-mail: axilleaz@protonmail.com
+%
+% Date: 22/02/2024 (DD/MM/YYYY)
+%
+% Copyright: MIT
+% --------------------------------------------------
+% Functionality: Perform tonal active noise control in the frequency
+% domain.
+% --------------------------------------------------
+% Input
+%
+% Pe [numeric]: The NxM matrix holding the tranfer functions from the
+% primary sources to the error/virtual microphone positions.
+% N is the number of error microphones and M the number of
+% primary sources.
+%
+% Pm [numeric]: The NxM matrix holding the transfer functions from the
+% primary sources to the monitoring microphones. N is the
+% number of monitoring microphones and M the number of
+% primary sources, which must match the number of columns
+% of Pe.
+%
+% Px [numeric]: The NxM matrix holding the transfer functions from the
+% primary sources to the reference microphone positions. N is
+% the number of reference microphones and M the number of
+% primary sources, which must match the number of columns of
+% Pe.
+%
+% Ge [numeric]: The NxM matrix holding the transfer functions from the
+% secondary sources to the error/virtual microphone
+% positions. N is the number of error/virtual microphones
+% which must match the number of rows of "Pe" and M the
+% number of secondary sources.
+%
+% Gm [numeric]: The NxM matrix holding the transfer functions from the
+% secondary sources to the monitoring microphone positions.
+% N is the number of monitoring microphones which must match
+% the number of rows of "Pm" and N the number of secondary
+% sources, which must match the number of columns of "Ge".
+%
+% GeEst [numeric] (Optional): The NxM matrix holding the estimated transfer
+% functions from the secondary sources to the
+% error/virtual microphone positions. N is the
+% number of error/virtual microphones which
+% must match the number of rows in Pm and M the
+% number of secondary sources, which must match
+% the number of columns of Ge. [Default: Ge].
+%
+% GmEst [numeric] (Optional): The NxM matrix holding the estimated transfer
+% functions from the secondary sources to the
+% monitoring microphone positions. N is the
+% number of monitoring microphones which must
+% match the number of rows in Pm and M the
+% number of secondary sources, which must match
+% the number of columns in Ge. [Default: Gm].
+%
+% v [numeric] (Optional): The primary source strengths or primary source
+% cross power spectral density matrix. If the
+% argument is a vector, its number of elements
+% must match the number of columns in Pe and
+% corresponds to the source strengths. If it is a
+% matrix it must be square with each dimension
+% matching the number of columns of Pe and will be
+% treated as the cross spectral density matrix of
+% the primary sources. [Default: eye(size(Pe, 2))].
+%
+% O [numeric] (Optional): The observation filter in case the remote
+% microphone technique is to be used. This argument
+% must be a matrix of dimensions NxM, with N being
+% the number of error/virtual microphones matching
+% the number of rows in Pe and M the number of
+% monitoring microphones matching the number of
+% rows of Pm. [Default: []].
+%
+% gRegFac [numeric] (Optional): A scalar regularisation factor used for
+% the inversion of the plant responses
+% matrix (G' * G). [Default: 0].
+%
+% xRegFac [numeric] (Optional): A scalar regularisation factor used for
+% the inversion of the reference signal
+% cross spectral density matrix (Sxx).
+% [Default: 0].
+%
+% filterType [char/string] (Optional): Whether the optimal filter, in the
+% least-squares sense, will be used or
+% the FxLMS implementation. Possible
+% values are "Optimal"/"opt" and
+% "FxLMS". The values are NOT case
+% sensitive. [Default: "Optimal"].
+%
+% Peval [numeric] (Optional): The transfer function matrix between the
+% primary sources and the positions of sound
+% field evaluation. The dimensions must be NxM
+% where N is the number of evaluation positions
+% and M the number of primary sources, which
+% must match the columns of "Pe". If left empty
+% or not provided, the values in "Pe" will be
+% used. [Default: []].
+%
+% Geval [numeric] (Optional): The transfer function matrix between the
+% secondary sources and the positions of sound
+% field evaluation (plant responses). The
+% dimensions must be NxM where N is the number
+% of evaluation positions, which must be equal
+% to the positions provided in "Peval", and M
+% is the number of secondary sources which must
+% be equal to the number of columns of "Ge". If
+% left empty or not provided "Ge" will be used.
+% [Default: []].
+%
+% --------------------------------------------------
+% Output
+%
+% W [numeric]: The NxM filter matrix. N is the number of secondary
+% sources matching the number of columns in Ge and M the
+% number of reference microphones, matching the number of
+% rows of Px.
+%
+% micErr [numeric]: This is a vector with number of elements equal to the
+% number of error/virtual microphones (equal to the
+% number of rows of Pe), holding the control error at
+% each error/virtual microphone.
+%
+% normSqMicErr [numeric]: The normalised squared errors at the
+% error/virtual microphones. This vector holds the
+% squares of "micErr" normalised to the power
+% spectral density of the true sound field at the
+% position of each error/virtual microphone.
+%
+% evalErr [numeric]: The vector holding the error at the evaluation
+% positions. This is empty if "Peval" is not provided.
+%
+% normSqEvalErr [numeric]: The vector holding the normalised squared error
+% at the evaluation positions. This is empty if
+% "Peval" is not provided.
+%
+% gCondNum [numeric]: The inversion condition numbers of the plant
+% response matrices to be inverted
+% [G' * G + regFac * eye()].
+%
+% --------------------------------------------------
+% Notes
+%
+% - If O is not provided or given as an empty array, the remote
+% microphone technique will not be used. In this case Pe, Ge and GeEst
+% are ignored in the control phase and only monitoring microphones are
+% assumed. Thus, only Pm, Gm and GmEst are used for control. However,
+% Pe and Ge are still used for the estimation at the error/virtual
+% microphone positions. This allows the evaluation of control without
+% virtual sensing at remote locations.
+%
+%
+% - The implementation is based on:
+% 1) Robust performance of virtual sensing methods for active noise
+% control by J. Zhang, S. J. Elliott and J. Cheer.
+% 2) Modeling local active sound control with remote sensors in
+% spatially random pressure fields by S.J. Elliott and J. Cheer.
+%
+% --------------------------------------------------
+function [W, micErr, normSqMicErr, evalErr, normSqEvalErr, gCondNum] = tonalAnc(Pe, Pm, Px, Ge, Gm, GeEst, GmEst, v, O, gRegFac, xRegFac, filterType, Peval, Geval)
+ % ====================================================
+ % Check for number of arguments
+ % ====================================================
+ narginchk(5, 14);
+ nargoutchk(0, 6);
+
+ % ====================================================
+ % Validate input arguments
+ % ====================================================
+ % Validate mandatory arguments
+ validateattributes(Pe, "numeric", {'nonempty', 'finite', 'nonnan'}, mfilename, "Tranfer function matrix from primary sources to error/virtual the microphones.", 1);
+ validateattributes(Pm, "numeric", {'nonempty', 'finite', 'nonnan', 'ncols', size(Pe, 2)}, mfilename, "Tranfer function matrix from primary sources to the monitoring microphones.", 2);
+ validateattributes(Px, "numeric", {'nonempty', 'finite', 'nonnan', 'ncols', size(Pe, 2)}, mfilename, "Tranfer function matrix from primary sources to the reference microphones.", 3);
+ validateattributes(Ge, "numeric", {'nonempty', 'finite', 'nonnan', 'nrows', size(Pe, 1)}, mfilename, "Plant response matrix from the secondary sources to the error/virtual microphones", 4);
+ validateattributes(Gm, "numeric", {'nonempty', 'finite', 'nonnan', 'nrows', size(Pm, 1)}, mfilename, "Plant response matrix from the secondary sources to the monitoring microphones", 5);
+
+ % Validate optional arguments
+ if nargin > 5 && ~isempty(GeEst)
+ validateattributes(GeEst, "numeric", {'finite', 'nonnan', 'nrows', size(Pe, 1), 'ncols', size(Ge, 2)}, mfilename, "Estimated plant response matrix from the secondary sources to the error/virtual microphones", 6);
+ else
+ GeEst = Ge;
+ end
+
+ if nargin > 6 && ~isempty(GmEst)
+ validateattributes(GmEst, "numeric", {'finite', 'nonnan', 'nrows', size(Pm, 1), 'ncols', size(Ge, 2)}, mfilename, "Estimated plant response matrix from the secondary sources to the monitoring microphones", 7);
+ else
+ GmEst = Gm;
+ end
+
+ if nargin > 7 && ~isempty(v)
+ if isvector(v)
+ validateattributes(v, "numeric", {'nonnan', 'finite', 'real', 'numel', size(Pe, 2)}, mfilename, "Primary source strengths.", 8);
+ Svv = diag(v.^2);
+ else
+ validateattributes(v, "numeric", {'2d', 'square', 'finite', 'nonnan', 'real', 'ncols', size(Pe, 2)}, mfilename, "Primary source cross spectral density matrix.", 8);
+ Svv = v;
+ end
+ else
+ Svv = eye(size(Pe, 2));
+ end
+
+ if nargin > 8 && ~isempty(O)
+ validateattributes(O, "numeric", {'2d', 'finite', 'nonnan', 'nrows', size(Pe, 1), 'ncols', size(Pm, 1)}, mfilename, "Observation filter.", 9);
+ else
+ O = [];
+ end
+
+ if nargin > 9 && ~isempty(gRegFac)
+ validateattributes(gRegFac, "numeric", {'scalar', 'nonnegative', 'finite', 'nonnan', 'real'}, mfilename, "Plant response regularisation factor.", 10);
+ else
+ gRegFac = 0;
+ end
+
+ if nargin > 10 && ~isempty(xRegFac)
+ validateattributes(gRegFac, "numeric", {'scalar', 'nonnegative', 'finite', 'nonnan', 'real'}, mfilename, "Reference signal cross spectral density matrix regularisation factor.", 11);
+ else
+ xRegFac = 0;
+ end
+
+ if nargin > 11 && ~isempty(filterType)
+ validateattributes(filterType, {'char', 'string'}, {'scalartext'}, mfilename, "Filter type", 12);
+ validatestring(filterType, ["Optimal", "Opt", "FxLMS"], mfilename, "Filter type", 12);
+ else
+ filterType = "optimal";
+ end
+
+ if nargin > 12 && ~isempty(Peval)
+ validateattributes(Peval, {'numeric'}, {'2d', 'ncols', size(Pe, 2)}, mfilename, "Primary path transfer functions to evaluation positions", 13);
+ else
+ Peval = [];
+ end
+
+ if nargin > 13 && ~isempty(Geval)
+ if isempty(Peval)
+ error("'Peval' cannot be empty if 'Geval' is provided.");
+ end
+
+ validateattributes(Geval, {'numeric'}, {'2d', 'ncols', size(Ge, 2)}, mfilename, "Secondary path transfer functions to evaluation positions (plant responses)", 14);
+ else
+ Geval = [];
+ end
+
+
+ % ====================================================
+ % Calculated needed quantities
+ % ====================================================
+ % Calculate cross spectral densities
+ Sxm = Pm * Svv * Px';
+ Sxx = Px * Svv * Px';
+
+ % Calculate plant response matrices
+ if isempty(O)
+ G = GmEst;
+ else
+ G = GeEst + O * (Gm - GmEst);
+ end
+
+
+
+ % ====================================================
+ % Calculate filters
+ % ====================================================
+ if isempty(O)
+ % No remote microphone technique
+ invG = (G' * G);
+ if gRegFac ~= 0
+ invG = invG + gRegFac * eye(size(invG));
+ end
+
+ % Optimal solution and FxLMS are the same
+ W = -invG\G';
+ W = W * Sxm/(Sxx + xRegFac * eye(size(Sxx)));
+ else
+ % With remote microphone technique
+ if strcmpi(filterType, "FxLMS")
+ invG = GeEst' * G;
+ if gRegFac ~= 0
+ invG = invG + gRegFac * eye(size(invG));
+ end
+
+ W = -invG\GeEst' * O * Sxm/(Sxx + xRegFac * eye(size(Sxx)));
+ else
+ invG = G' * G;
+ if gRegFac ~= 0
+ invG = invG + gRegFac * eye(size(invG));
+ end
+
+ W = -invG\G' * O * Sxm/(Sxx + xRegFac * eye(size(Sxx)));
+ end
+ end
+
+
+
+ % ====================================================
+ % Calculate errors
+ % ====================================================
+ % Control error
+ if nargout > 1
+ micErr = Pe + Ge * W * Px;
+ end
+
+ % Normalised squared control error
+ if nargout > 2
+ normSqMicErr = diag(micErr * Svv * micErr')./diag(Pe * Svv * Pe');
+ end
+
+ % Error at evaluation positions
+ if nargout > 3 && ~isempty(Peval)
+ evalErr = Peval + Geval * W * Px;
+ else
+ evalErr = [];
+ end
+
+ % Normalised squared error at evaluation positions
+ if nargout > 4 && ~isempty(evalErr)
+ normSqEvalErr = diag(evalErr * Svv * evalErr')./diag(Peval * Svv * Peval');
+ else
+ normSqEvalErr = [];
+ end
+
+ % Condition number of inverse of (G' * G + regFac * eye())
+ if nargout > 5
+ gCondNum = cond(invG);
+ end
+end
\ No newline at end of file
diff --git a/Active Noise Control/Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/virtMicGeo.m b/Active Noise Control/Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/virtMicGeo.m
new file mode 100644
index 0000000000000000000000000000000000000000..b5086ce7841e3568a7b733e8d52362ff91c1d857
--- /dev/null
+++ b/Active Noise Control/Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/virtMicGeo.m
@@ -0,0 +1,225 @@
+%% Calculate virtual microphone positions for specific geometries
+% --------------------------------------------------
+% Author: Achilles Kappis
+% e-mail: axilleaz@protonmail.com
+%
+% Date: 25/03/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)
+ 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(strmpci(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
diff --git a/Active Noise Control/Virtual Sensing/Remote Microphone Technique/MATLAB/README.md b/Active Noise Control/Virtual Sensing/Remote Microphone Technique/MATLAB/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..9b7db07ab447fbb12e40bc6666ad91d101dc1cac
--- /dev/null
+++ b/Active Noise Control/Virtual Sensing/Remote Microphone Technique/MATLAB/README.md
@@ -0,0 +1,89 @@
+# Remote Microphone Technique
+This project holds the codebase used in the sound field estimation simulations with the *Remote Microphone Technique* (RMT) within the *DC12 - Directional microphone arrays for remote microphone virtual sensing* of the [IN-NOVA](https://in-nova-horizon.eu/) project.
+
+Currently, only MATLAB® scripts and functions are included. Updates are done on an as-needed basis, and at the moment, there is no plan for in porting to other programming languages or environments.
+
+The project uses [semantic versioning](https://semver.org/) and the current version is **v0.1.0**.
+
+
+## Brief introduction to the Remote Microphone Technique
+The *Remote Microphone Technique* (RMT) was introduced by Roure and Albarrazin \[[1](#references)\] to control a sound field at locations remote from sensors. The sound field is estimated at a remote position based on the measurements at the physical sensors through an observation filter, constituting the transfer function from the physical sensors to the remote location termed *virtual microphone*. An optimal, in the least-squares sense, observation filter is estimated in a preliminary *identification phase* and consequently used during the control phase.
+
+More information on the formulation used in the development of this project can be found in the [*Documentation*]() file.
+
+
+## Dependencies
+The project uses external dependencies for the implementation of the functions. The dependencies can be found in the [MATLAB® Utilities project](https://gitlab.com/in-nova/utilities/matlab). The necessary files are listed below:
+
+- srcOnSphere.m
+- srcRcvDist.m
+- rotMat3d.m
+
+## Examples
+A condensed example of a very simple estimation case is shown below and more examples can be found in the [*Documentation*]() folder, along with a complete description of all the functions.
+
+```matlab
+% Place primary sources, physical and virtual microphones
+mPos = rcvGeo("UCA", 4, 0.5); % Omni mic positions
+vPos = srcGeo("Diffuse", 0, 3, 10); % Primary source positions
+ePos = virtMicGeo("Grid", 2, 0, 0, 0, 11^2, "xy"); % Virtual/error mic positions on a grid
+
+% Calculate transfer functions
+Pm = ptSrcField(vPos, mPos, 7.5e2); % Primary-to-monitoring transfer function
+Pe = ptSrcField(vPos, ePos, 7.5e2); % Primary-to-virtual transfer function
+
+% Perform estimation and acquire the normalised square estimation error
+[~, ~, ~, ~, ~, ~, normSqrErr] = obsFilt(Pe, Pm, [], 1e2); % Use regularisation factor 100
+```
+
+## Support
+If you find any bugs, other issues and/or inconsistencies you can [e-mail](mailto::axilleaz@protonmail.com) Achilles Kappis at *<axilleaz@protonmail.com>* or *<A.Kappis@soton.ac.uk>*. More contacts will be added if/when new contributors join the project.
+
+## Roadmap
+The repository will be updated on an as-needed basis.
+
+At the moment, it is planned to include modal sound fields in the available primary source fields.
+
+It is expected that in the near future, active noise control will also be implemeneted. Some of the functions related to both virtual sensing and control (such as the source and receiver position functions) will be moved to a project higher in the hierarchy or in the [MATLAB® Utilities project](https://gitlab.com/in-nova/utilities/matlab).
+
+## Contributing
+If you would like to contribute to the project you could do that in various ways, as listed below:
+
+- Create issues for bugs, improvements and/or features.
+- Fix a bug or develop an improvement and/or feature yourself. Clone the repository, create a new branch when done create a merge request.
+- Report typos and/or inconsistencies in the documentation.
+- Contribute to the documentation and/or the Wiki of the project.
+- Create a CI configuration for automated testing, or develop unit tests.
+- Port the codebase to other programming languages and/or environments.
+
+#### Remarks
+
+Since there is no CI set at the moment to perform automated testing, you are highly encouraged to test your contributions as much as possible and even contribute some unit tests if this aligns with your contributed content.
+
+To try and maintain a coherent style in the codebase you are encouraged to use the template available in the [MATLAB® Utilities project](https://gitlab.com/in-nova/utilities/matlab), although this is not necessary.
+Irrespective of style, it is highly advisable to name your variables in a meaningful manner and keep them as short as possible while still communicating their purpose. Comments are highly recommended in your code, especially where they can facilitate a better understanding of the code and the intent of the developer. Over-commented code is better than under-commented code.
+
+
+## Authors and acknowledgment
+
+#### Authors
+
+The main contributor is Achilles Kappis.
+
+
+#### Acknowledgments
+
+None.
+
+
+## License
+The project is licensed under the ***MIT License***, which is a rather unrestrictive license, effectively passing the code to the public scope. For more information see the [MIT License](https://mit-license.org/), or the [LICENSE](https://gitlab.com/in-nova/virtual-sensing/remote-microphone-technique/-/blob/main/LICENSE) file of the project.
+
+#### **Important**
+
+If you are unable, or do not want to contribute under this license please do **not** contribute to this project at all. All content within the project must be under the same license.
+
+
+## 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).
\ No newline at end of file