diff --git a/Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/obsFiltTD.m b/Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/obsFiltTD.m
new file mode 100644
index 0000000000000000000000000000000000000000..09a53a42a368f7482eadcaa05b4ac3afa0bad836
--- /dev/null
+++ b/Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/obsFiltTD.m
@@ -0,0 +1,247 @@
+%% Calculate optimal observation filters in the time domain
+% --------------------------------------------------
+% Author: Achilles Kappis
+% e-mail: axilleaz@protonmail.com
+%
+% Date: 12/09/2024 (DD/MM/YYYY)
+%
+% Copyright: MIT
+% --------------------------------------------------
+% Functionality: Calculate the optimal, in the least squares sense,
+% observation filters for the Remote Microphone Technique in
+% the time domain.
+% --------------------------------------------------
+% Input
+%
+% e [numeric]: The measurement(s) at the virtual microphone position(s).
+% This must be an IxNxJ matrix with I the length of the
+% measurements in samples, N the number of virtual microphones
+% and J the number of trials/sound field realisations.
+%
+% m [numeric]: The measurement(s) at the monitoring microphone position(s).
+% This must be an KxMxJ matrix with K representing the length
+% of the measurements in samples, M the number of monitoring
+% microphones and J the number of trials/sound field
+% realisations. The number of trials must be the same as in
+% the argument "e".
+%
+% beta [numeric] (Optional): The regularisation factor. This must be a real
+% non-negative scalar. [Default: 0].
+%
+% filtLen [numeric] (Optional): The length of the observation filters in
+% samples. This must be a positive scalar, at
+% less than or equal to K, the number of
+% samples in the monitoring microphone
+% signals. [Default: size(m, 2)].
+%
+% delay [numeric] (Optional): The delay to be added to the virtual
+% microphone signals, effectively implementing
+% the delayed Remote Microphone Technique. This
+% must be a real non-negative scalar less than
+% or equal to I, the length of the virtual
+% microphone signals. [Default: 0].
+%
+% --------------------------------------------------
+% Output
+%
+% O [numeric]: The observation filters. This is an firLenxNxMxJ array.
+%
+% Rme [numeric]: The cross-correlation matrices between the monitoring and
+% virtual microphone signals. This is an filtLenxMxNxJ
+% array.
+%
+% Rmm [numeric]: The cross-correlation matrices of the monitoring
+% microphone signals. This is an filtLenxfiltLenxMxMxJ
+% array.
+%
+% Ovec [numeric]: This is a matrix with the observation filters
+% stacked/vectorised, so that they can be applied to a
+% stacked/vectorised monitoring microphone measurements
+% vector. The dimensions of this matrix are NxTotFiltLenxJ,
+% where TotFiltLen is equal to M times filtLen. This way,
+% stacking the monitoring microphone measured signals one
+% can perform Ovec(:, :, jIdx) * dm (dm is the vectorised
+% monitoring microphone signals vector) to calculate the
+% estimated measurements at the N virtual microphone
+% positions.
+%
+% RmeMtx [numeric]: This is the matrix with the cross-correlation vectors
+% between the monitoring and virtual microphones
+% stacked/vectorised. The dimensions of the matrix are
+% NxTotFiltLenxJ.
+%
+% RmmMtx [numeric]: This is the matrix with the cross-correlation matrices
+% of the monitoring microphone signals stacked together.
+% The dimensions are TotFiltLenxTotFiltLenxJ.
+%
+% mMtx [numeric]: The monitoring microphone signals vectorised so that they
+% can be used directly with either "Ovec" or "Oopt"
+% arguments to perform the filtering/estimation.
+%
+% Omean [numeric]: The observation filters averaged over the trials/sound
+% field realisations. This is an NxMxfirLen array.
+%
+% Rmemean [numeric]: The cross-correlation matrices between the monitoring
+% and virtual microphone signals averaged over the
+% trials/sound field realisations. This is an
+% NxMxfirLen array.
+%
+% RmmMean [numeric]: The cross-correlation matrices of the monitoring
+% microphone signals averaged over the trials/sound
+% field realisations. This is an MxMxfirLenxfirLen
+% array.
+%
+% Oopt [numeric]: This is a matrix with the observation filters
+% stacked/vectorised, so that they can be applied to a
+% stacked/vectorised monitoring microphone measurements
+% vector, averaged over the trials/sound field
+% realisations. The dimensions of this matrix are
+% NxTotFiltLen. These filters are denoted "optimal" since
+% they are averaged over many realisations (if J > 1).
+%
+% RmeMtxMean [numeric]: This is the matrix with the cross-correlation
+% vectors between the monitoring and virtual
+% microphones stacked/vectorised, averaged over the
+% trials/sound field realisations. The dimensions of
+% the matrix are NxTotFiltLen.
+%
+% RmmMtxMean [numeric]: This is the matrix with the cross-correlation
+% matrices of the monitoring microphone signals
+% stacked together and averaged over the trials/sound
+% field realisations. The dimensions are
+% TotFiltLenxTotFiltLen.
+%
+% --------------------------------------------------
+% Notes
+%
+% - The calculations are based on the paper: "Modeling local active sound
+% control with remote sensors in spatially random pressure fields" by
+% Stephen J. Elliott and Jordan Cheer.
+%
+% --------------------------------------------------
+function [O, Rme, Rmm, Ovec, RmeMtx, RmmMtx, mMtx, Omean, RmeMean, RmmMean, Oopt, RmeMtxMean, RmmMtxMean] = obsFiltTD(e, m, beta, filtLen, delay)
+ % ====================================================
+ % Check for number of arguments
+ % ====================================================
+ narginchk(2, 5);
+ nargoutchk(0, 13);
+
+ % ====================================================
+ % Validate input arguments
+ % ====================================================
+ % Validate mandatory arguments
+ validateattributes(e, "numeric", {'3d', 'real', 'nonnan', 'nonempty', 'finite'}, mfilename, "Virtual microphone measurements", 1);
+ validateattributes(m, "numeric", {'3d', 'real', 'nonnan', 'nonempty', 'finite', 'size', [NaN, NaN, size(e, 3)]}, mfilename, "Monitoring microphone measurements", 2);
+
+ % Validate optional arguments
+ if nargin > 2 && ~isempty(beta)
+ validateattributes(beta, "numeric", {'scalar', 'nonnegative', 'nonnan', 'nonempty', 'finite'}, mfilename, "Regularisation factor", 3);
+ else
+ beta = 0;
+ end
+
+ if nargin > 3 && ~isempty(filtLen)
+ validateattributes(filtLen, "numeric", {'scalar', 'positive', 'nonnan', 'nonempty', 'finite', '<=', size(m, 1)}, mfilename, "Length of observation filters", 4);
+ else
+ filtLen = size(m, 2);
+ end
+
+ if nargin > 4 && ~isempty(delay)
+ validateattributes(delay, "numeric", {'scalar', 'nonnegative', 'nonnan', 'finite', 'nonempty', '<=', size(e, 1)}, mfilename, "Delay to be added to virtual microphone signals to implement the delayed Remote Microphone Technique", 5);
+ else
+ delay = 0;
+ end
+
+
+ % ====================================================
+ % Calculate cross- and auto-correlation matrices
+ % ====================================================
+ % Go through the trials/sound field realisations
+ for jIdx = size(m, 3):-1:1
+ % Go through the monitoring microphones
+ 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);
+
+ Rme(:, mIdx, eIdx, jIdx) = corr(delay + lIdx:-1:lIdx-filtLen + 1);
+ 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);
+
+ Rmm(:, :, mIdx, mmIdx, jIdx) = toeplitz(corr(lIdx:-1:lIdx - filtLen + 1));
+ else
+ [corr, lags] = xcorr(m(:, mIdx, jIdx), m(:, mmIdx, jIdx));
+ lIdx = find(lags == 0);
+
+ % Cross-correlation matrices
+ for iIdx = filtLen-1:-1:0
+ Rmm(:, iIdx + 1, mIdx, mmIdx, jIdx) = corr(iIdx + (lIdx:-1:lIdx - filtLen + 1));
+ end
+ Rmm(:, :, mmIdx, mIdx, jIdx) = squeeze(Rmm(:, :, mIdx, mmIdx, jIdx)).';
+ end
+ end
+ end
+ end
+
+ % ====================================================
+ % 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));
+
+ % ====================================================
+ % Calculate observation filters
+ % ====================================================
+ for jIdx = size(RmmMtx, 3):-1:1
+ 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]);
+
+ % ====================================================
+ % Provide additional output arguments
+ % ====================================================
+ % The monitoring microphone signals vectorised per trial/sound field realisation
+ if nargout > 6
+ mMtx = reshape(m, prod(size(m, [1, 2])), size(m, 3));
+ end
+
+ % Mean observation over trials/sound field realisations
+ if nargout > 7
+ Omean = mean(O, 4);
+ end
+
+ % Mean cross-correlations between monitoring and virtual microphones over trials/sound field realisations
+ if nargout > 8
+ RmeMean = mean(Rme, 4);
+ end
+
+ % Mean cross-correlations of monitoring microphones over trials/sound field realisations
+ if nargout > 9
+ RmmMean = mean(Rmm, 5);
+ end
+
+ % Average of optimal observation filters over trials/sound field realisations
+ if nargout > 10
+ Oopt = mean(Ovec, 3);
+ end
+
+ % Mean cross-correlation matrix between monitoring and virtual microphones over trials/sound field realisations
+ if nargout > 11
+ RmeMtxMean = mean(RmeMtx, 3);
+ end
+
+ % Mean cross-correlation matrix of monitoring microphones over trials/sound field realisations
+ if nargout > 12
+ RmmMtxMean = mean(RmmMtx, 3);
+ end
+end
\ No newline at end of file