diff --git a/CHANGELOG.md b/CHANGELOG.md
index daf1f7996311f9b079b44424ffbe7c2ce521fbca..88832390ef4e8cb3fd7223f4a812d5ea8e46f47d 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,3 +1,21 @@
+### v0.2.3 ###
+**Virtual Sensing**\
+\+ Added function to estimate the observation filters in the time-domain.\
+\+ Added function to perform estimation with observation filters in the time-domain.
+
+**Signal Processing - Generic**\
+\* Fixed bugs, improved calculations and removed output arguments for the windowed sinc FIR fractional delay filters
+\+ Added four window functions in the calculation of windowed sinc FIR fractional delay filters.\
+
+**Sound Fields**\
+\+ Added function to calculate signals generated by ideal point sources in the time-domain.
+
+
+### v0.2.2 ###
+**Virtual Sensing**\
+\* Update the observation filter and estimation with observation filter functions with better input argument checking and changed the name of variables to match those in the literature.
+
+
### v0.2.1 ###
**Virtual Sensing**\
\* Fix a bug where noise was added to the power spectral density matrix for the optimal observation filter calculations in the noiseless case.
diff --git a/README.md b/README.md
index fd23c6f121163294a72a490cbedd02ed92788c59..ed4f39a3f1f2a6ce28c3d921acb7c445e7ced3a3 100644
--- a/README.md
+++ b/README.md
@@ -3,6 +3,7 @@ This repository holds the open-access codebase used within the [IN-NOVA](https:/
The codebase under each section may or may not include code implementations in different languages but no plan exists to port each implementation in all available programming languages. More information about the codebase, the available implementations and topics/fields can be found in the [Documentation](https://gitlab.com/in-nova/in-nova/-/tree/main/Documentation/latex?ref_type=heads) of the project and its [Wiki](https://gitlab.com/in-nova/in-nova/-/wikis/home).
+
## Current structure
The current structure of the codebase is shown below. More information can be found in the [Wiki](https://gitlab.com/in-nova/in-nova/-/wikis/home) with separate pages for each category, and implementation.
@@ -17,6 +18,8 @@ The current structure of the codebase is shown below. More information can be fo
- Frequency domain
- [Virtual Sensing](https://gitlab.com/in-nova/in-nova/-/tree/main/Virtual%20Sensing)
- Remote Microphone Technique
+ - Frequency domain
+ - Time domain
- [Optimisation](https://gitlab.com/in-nova/in-nova/-/tree/main/Optimisation/Memetic%20Algorithm)
- Memetic Algorithm
- [Signal Processing](https://gitlab.com/in-nova/in-nova/-/tree/main/Signal%20Processing)
@@ -24,17 +27,28 @@ The current structure of the codebase is shown below. More information can be fo
- Generic
- [Sound fields](https://gitlab.com/in-nova/in-nova/-/tree/main/Sound%20Fields)
- Generation
+ - Frequency domain
+ - Time domain
+ - Wave domain
+ - Analysis
+ - Wave domain Discrete Fourier Transforms
+
## Dependencies
Effort has been put towards minimising dependencies external to this repository. Even "inter-repository" dependencies (ex. a function under some field uses a general utility function) will be clearly stated in the relevant code documentation.
+
## 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 this stage, the repository is *under construction* and the main focus lies on uploading the available code, updating documentation and structuring the GitLab project to contain all information related to the codebase.
+Although the repository is not considered to be *under construction* anymore, there are significant gaps in its structure. At the moment, the main focus is placed on updating code with new features and bug fixes and documenting ideas and useful features for the future (mainly as issues).
+
+As a parallel slow-moving project, porting the codebase to [Julia](https://julialang.org/) has also started, aiming to provide full functionality in the future. Although Julia and MATLAB were designed to tackle similar (if not the same) problems, they are quite different. The facilities and intricacies of Julia dictate the need to dig deeper into the language in order to provide an efficient, maintainable and intuitive codebase that will parallel the current MATLAB implementation in functionality. More time must be spent in the design phase for the Julia implementations, which will, most probably, keep the progress to a slow pace. Although it is a long shot, the target is for the Julia codebase to reach the MATLAB implementation's state and then move forward in tandem.
+
## Contributing
If you would like to contribute to the project you could do that in various ways, as listed below:
@@ -46,6 +60,7 @@ If you would like to contribute to the project you could do that in various ways
- 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.
@@ -70,9 +85,8 @@ Part of the `ptsOnSphere()` MATLAB/Octave function is an adaptation of a functio
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/in-nova/-/blob/main/LICENSE) file of the project.
-
## Versioning ##
-The project uses [semantic versioning](https://semver.org/) and the current version is **v0.2.2**.
+The project uses [semantic versioning](https://semver.org/) and the current version is **v0.2.3**.
#### **Important**
diff --git a/Signal Processing/Generic/MATLAB/Functions/winSincFracDel.m b/Signal Processing/Generic/MATLAB/Functions/winSincFracDel.m
index c2e5f3b29052731b46157ca2bb9a57b3d78f074c..3fcfcd620a036056f0748f89cec189170cffd8a4 100644
--- a/Signal Processing/Generic/MATLAB/Functions/winSincFracDel.m
+++ b/Signal Processing/Generic/MATLAB/Functions/winSincFracDel.m
@@ -3,7 +3,7 @@
% Author: Achilles Kappis
% e-mail: axilleaz@protonmail.com
%
-% Date: 17/04/2024 (DD/MM/YYYY)
+% Date: 13/08/2024 (DD/MM/YYYY)
%
% Copyright: MIT
% --------------------------------------------------
@@ -12,8 +12,8 @@
% --------------------------------------------------
% Input
%
-% delay [numeric]: The delay in samples. This must be a non-negative real
-% scalar.
+% del [numeric]: The delay in samples. This must be a non-negative real
+% scalar.
%
% len [numeric] (Optional): This is the length of the filter. This must be
% a positive real integer. If this values is less
@@ -21,18 +21,16 @@
% equals the "delay" value rounded up.
% [Default: ceil(delay)].
%
-% win [numeric] (Optional): A window to apply to the FIR filter. This can
-% be a vector in which case the same window will
-% be applied to the filter for all signals, or it
-% can be a matrix in which case each column will
-% contain the window corresponding to the filter
-% for each signal. The number of rows must be
-% equal to "len" and the number of columns to the
-% number of columns of "sig". If left empty, no
-% window is applied. It has to be noted that the
-% provided window has to be shifted by "delay"
-% samples.
-% [Default: []].
+% winFun [string/char, numeric] (Optional): Windowing functions to be
+% applied to the truncated sinc
+% function. The available options
+% are "Rectangular", "Hann" and
+% "Hamming" and are not
+% case-sensitive. Additionally, a
+% non-negative scalar can be
+% given to be used as the
+% parameter for a Kaiser window.
+% [Default: "Rectangular"].
%
% sig [numeric] (Optional): The signal(s) to be delayed in the time-domain.
% This can be a matrix with each column being a
@@ -52,74 +50,75 @@
% --------------------------------------------------
% Output
%
-% delFilt [numeric]: The fractional delay windowed-sinc FIR filter(s). If
-% multiple windows are provided then a matrix is
-% returned with each column corresponding to the filter
-% for each signal. Otherwise, this argument is a vector
-% with the FIR filter.
+% sincFilt [numeric]: The fractional delay windowed-sinc FIR filter(s).This
+% argument is a vector with the FIR filter. This is the
+% causal version of the filter. If the non-causal
+% version is required, shift to the left by "causDel".
%
-% intDel [numeric]: The causal latency of the filter. This value
-% corresponds to half the length of the filter.
+% causDel [numeric]: The delay required to make the filter causal. This
+% value corresponds to half the length of the filter.
%
-% causalFilt [numeric]: The filter after left shifting "intDel" samples to
-% the left.
-%
-% delSig [numeric]: The delayed signal(s). The length of each signal is
-% equal to that of the original. If the delay is
-% longer than the original signal vector, a zeroed vector
-% is returned (with possible artefacts from the filtering
-% process). Only returned if a signal is provided,
-% otherwise an empty array is returned.
-%
-% causalSig [numeric]: The delayed signal(s) shifted to the left by the
-% causal delay of the filter.
+% dSig [numeric]: The delayed signal(s). The length of each signal is equal
+% to that of the original. If the delay is longer than the
+% original signal vector, a zeroed vector is returned (with
+% possible artefacts from the filtering process). Only
+% returned if a signal is provided, otherwise an empty
+% array is returned.
%
% --------------------------------------------------
% Notes
%
% - The function uses time-domain convolution, which can be quite slow for
-% long signals. If the "delSig" is not asked for this calculation is
+% long signals. If the "dSig" is not required this calculation is
% omitted so, when the signals are not needed (or want to perform
% frequency-domain convolution manually), use only one output argument.
+%
% --------------------------------------------------
-function [delFilt, intDel, causalFilt, delSig, causalSig] = winSincFracDel(delay, len, win, sig, sigLen)
+function [sincFilt, causDel, dSig] = winSincFracDel(del, len, winFun, sig, sigLen)
% ====================================================
% Check for number of arguments
% ====================================================
narginchk(1, 5);
- nargoutchk(0, 5);
+ nargoutchk(0, 3);
% ====================================================
% Validate input arguments
% ====================================================
% Validate mandatory arguments
- validateattributes(delay, "numeric", {'scalar', 'real', 'nonnegative', 'nonnan', 'nonempty', 'finite'}, mfilename, "Delay in samples", 1);
+ validateattributes(del, "numeric", {'scalar', 'real', 'nonnegative', 'nonnan', 'nonempty', 'finite'}, mfilename, "Delay in samples", 1);
% Validate optional arguments
if nargin > 1 && ~isempty(len)
validateattributes(len, "numeric", {'scalar', 'integer', 'real', 'finite', 'nonnan', 'nonnegative'}, mfilename, "Filter length", 2);
- if len < delay
+ if len < del
warning("winSincFracDel(): The length of the filter cannot be less than the delay, setting equal to 'delay',");
- len = ceil(delay);
+ len = ceil(del);
end
else
- len = ceil(delay);
+ len = ceil(del);
end
- if nargin > 2 && ~isempty(win)
- validateattributes(win, "numeric", {'nonnan', '2d', 'real', 'finite', 'nrows', len}, mfilename, "Window(s) to be applied", 3);
+ if nargin > 2 && ~isempty(winFun)
+ validateattributes(winFun, {'char', 'string', 'numeric'}, {'nonempty'}, mfilename, "Window to be applied to the truncated sinc function (or the alpha parameter for a Kaiser window)", 3);
+
+ if ischar(winFun) || isStringScalar(winFun)
+ validatestring(winFun, ["Rectangular", "Hann", "Hamming"], mfilename, "Window to be applied to the truncated sinc function", 3);
+ winFun = string(winFun);
+ elseif isscalar(winFun)
+ validateattributes(winFun, {'numeric'}, {'nonempty', 'nonnan', 'finite', 'nonnegative', 'real', 'scalar'}, mfilename, "Alpha parameter for a Kaiser window to be applied to the truncated sinc function", 3);
+ end
else
- win = ones(len, 1);
+ winFun = "rectangular";
end
if nargin > 3 && ~isempty(sig)
validateattributes(sig, "numeric", {'2d', 'nonempty', 'real'}, mfilename, "Signal(s) to be delayed", 4);
- if size(win, 2) == 1
- win = repmat(win, 1, size(sig, 2));
+ if size(winFun, 2) == 1
+ winFun = repmat(winFun, 1, size(sig, 2));
else
- if size(win, 2) ~= size(sig, 2)
+ if size(winFun, 2) ~= size(sig, 2)
error("Number of windows must match the number of signals provided");
end
end
@@ -134,47 +133,54 @@ function [delFilt, intDel, causalFilt, delSig, causalSig] = winSincFracDel(delay
sigLen = "Full";
end
+
% ====================================================
- % Generate FIR filter
+ % Calculate some parameters
% ====================================================
- delFilt = -len + 1:len; delFilt = delFilt.'; % Filter sample indices
-
- intDel = -floor(delFilt(1)/2); % The causal delay of the filter
- delFilt = sinc(delFilt - delay); % Calculate the filter
+ fracDel = rem(del, 1); % Fractional part of the delay
- % Keep only the central part of the filter to match length asked for
- delFilt = delFilt(intDel + 1:intDel + len);
- % Replicate filter and apply window function
- delFilt = repmat(delFilt, 1, size(sig, 2)) .* win;
+ % ====================================================
+ % Generate FIR filter for fractional delay
+ % ====================================================
+ idx = floor(-len/2 + 1):floor(len/2); % Filter sample indices
+ causDel = -idx(1); % The causal delay of the filter
+ sincFilt = sinc(idx - fracDel); % Calculate the filter
+
% ====================================================
- % Filter the signals
+ % Apply windowing
% ====================================================
- % Return the shifted filter
- if nargout > 2
- causalFilt = delFilt(intDel:end);
- end
-
- if nargout > 3
- for idx = size(sig, 2):-1:1
- delSig(:, idx) = conv(sig, delFilt, "full");
+ % Generate window
+ if isstring(winFun)
+ switch lower(winFun)
+ case "rectangular"
+ win = ones(1, length(idx));
+ case "hann"
+ win = sin(pi * (idx - fracDel + causDel + 1)/len).^2;
+ case "hamming"
+ alpha = 25/46;
+ win = alpha - (1 - alpha) * cos(2 * pi * (idx - fracDel + causDel + 1)/len);
end
+ else
+ % Kaiser window
+ win = besseli(0, pi * winFun * sqrt(1 - (2 * (idx - fracDel)/len).^2))/besseli(0, pi * winFun);
end
+
+ % Apply window to the filter
+ sincFilt = sincFilt .* win;
- % Shift signal to the left by the amount of the "causal" delay of the filter
- if nargout > 4
- causalSig = delSig(intDel:end, :);
- end
+ % Add the integral delay
+ sincFilt = [zeros(1, fix(del)), sincFilt(1:end - fix(del))];
- % Return the length of the signal asked for
- if strcmpi(sigLen, "Same")
- if nargout > 3
- delSig = delSig(1:size(sig, 1), :);
- end
- if nargout > 4
- causalSig = causalSig(1:size(sig, 1), :);
+
+ % ====================================================
+ % Filter the signals
+ % ====================================================
+ if nargout > 2
+ for idx = size(sig, 2):-1:1
+ dSig(:, idx) = conv(sig, sincFilt.', lower(sigLen));
end
end
end
\ No newline at end of file
diff --git a/Sound Fields/MATLAB/Functions/ptSrcFieldTD.m b/Sound Fields/MATLAB/Functions/ptSrcFieldTD.m
new file mode 100644
index 0000000000000000000000000000000000000000..01fa2c62872193b9362ebd5ccb9fd37ef4134279
--- /dev/null
+++ b/Sound Fields/MATLAB/Functions/ptSrcFieldTD.m
@@ -0,0 +1,187 @@
+%% Signal generated by point sources in the time-domain
+% --------------------------------------------------
+% Author: Achilles Kappis
+% e-mail: axilleaz@protonmail.com
+%
+% Date: 14/09/2024 (DD/MM/YYYY)
+%
+% Copyright: MIT
+% --------------------------------------------------
+% Functionality: Calculate the signals generated by ideal point sources in
+% the time-domain.
+% --------------------------------------------------
+% Input
+%
+% sPos [numeric]: The Cartesian coordinates of the source positions. This
+% must be an Nx3 matrix where N is the number of sources.
+%
+% rPos [numeric]: The Cartesian coordinates of the receiver positions. This
+% must be an Mx3 matrix where M is the number of receivers.
+%
+% fs [numeric]: The sampling frequency. This is required in order to
+% calculate the time-of-flight in samples and must be a
+% positive real scalar.
+%
+% Q [numeric] (Optional): The source signals. This can be either an IxNxJ
+% array, where I is the length of the source
+% signals in samples and J is the number of
+% trials/sound field realisations, a vector of
+% positive real values holding the source strengths
+% (this is just an amplitude scaling factor), or a
+% real positive integer denoting the common
+% strength of all the source signals. [Default: 1].
+%
+% sigLen [numeric] (Optional): The length of the source signals. This must
+% be real positive integer. The generated
+% signals are zero-mean, (approximately)
+% uniformly distributed and normalised to
+% unity. If "Q" is an array, "sigLen" is used
+% to set the length of the signals. If
+% "sigLen" is smaller than I, the signals are
+% truncated and if it is larger the signals
+% are padded with zeros. Leave empty to have
+% the signals unchanged.
+% [Default: 128 or if Q is matrix size(Q, 1)].
+%
+% nTrials [numeric] (Optional): This is the number of trials/sound field
+% realisations that will be generated. If "Q"
+% is provided, this argument is ignored. It
+% must be a positive real scalar.
+% [Default: 1]
+%
+% c [numeric] (Optional): The speed of sound. [Default: 343].
+%
+% --------------------------------------------------
+% Output
+%
+% rSig [numeric]: The signals at the receiver positions for all
+% trials/sournd field realisations. This is an IxMxJ array.
+%
+% rSigMean [numeric]: The signals at the receiver positions averaged over
+% all trials/sound field realisations. This is an IxM
+% matrix.
+%
+% rSigMtx [numeric]: The signals at each receiver position due to each
+% source for each trial/sound field realisation. This is
+% an IxMxNxJ array.
+%
+% rSigMtxMean [numeric]: The signals at each receiver position due to each
+% source averaged over the trials/sound field
+% realisations. This is an IxMxN array.
+%
+% Q [numeric]: The source signals. If Q is provided, the same variable is
+% returned here. If Q is generated internally, the signals are
+% returned.
+%
+% --------------------------------------------------
+% Notes
+%
+% Dependencies: - twoPtDist(): To calculate the distances between sources
+% and receivers.
+%
+% - It is the responsibility of the user to pick (or provide) adequately
+% long source signals
+%
+% --------------------------------------------------
+function [rSig, rSigMean, rSigMtx, rSigMtxMean, Q] = ptSrcFieldTD(sPos, rPos, fs, Q, sigLen, nTrials, c)
+ % ====================================================
+ % Check for number of arguments
+ % ====================================================
+ narginchk(3, 7);
+ nargoutchk(0, 5);
+
+ % ====================================================
+ % Validate input arguments
+ % ====================================================
+ % Validate mandatory arguments
+ validateattributes(sPos, "numeric", {'2d', 'real', 'nonnan', 'nonempty', 'finite', 'ncols', 3}, mfilename, "Cartesian coordinates of the source positions", 1);
+ validateattributes(rPos, "numeric", {'2d', 'real', 'nonnan', 'nonempty', 'finite', 'ncols', 3}, mfilename, "Cartesian coordinates of the receiver positions", 2);
+ validateattributes(fs, "numeric", {'scalar', 'real', 'nonnan', 'nonempty', 'finite', 'positive'}, mfilename, "The sampling frequency", 3);
+
+ % Validate optional arguments
+ if nargin > 3 && ~isempty(Q)
+ if isscalar(Q)
+ validateattributes(Q, "numeric", {'scalar', 'real', 'nonnan', 'nonempty', 'finite', 'positive', 'integer'}, mfilename, "Length of source signals in samples", 4);
+ elseif isvector(Q)
+ validateattributes(Q, "numeric", {'vector', 'real', 'nonnan', 'nonempty', 'finite', 'numel', size(sPos, 1)}, mfilename, "Source signals", 4);
+ else
+ validateattributes(Q, "numeric", {'3d', 'real', 'nonnan', 'nonempty', 'finite', 'ncols', size(sPos, 1)}, mfilename, "Source signals", 4);
+ end
+ else
+ Q = 1;
+ end
+
+ if nargin > 4 && ~isempty(sigLen)
+ validateattributes(sigLen, "numeric", {'scalar', 'nonempty', 'nonnan', 'finite', 'real', 'positive', 'integer'}, mfilename, "The length of the source signals", 5)
+ else
+ if ~isvector(Q)
+ sigLen = size(Q, 1);
+ else
+ sigLen = 128;
+ end
+ end
+
+
+ if ~ismatrix(Q)
+ nTrials = size(Q, 3);
+ elseif nargin > 5 && ~isempty(nTrials)
+ validateattributes(nTrials, "numeric", {'scalar', 'positive', 'nonnan', 'nonempty', 'finite', 'integer'}, mfilename, "Number of trials/sound field realisations", 6);
+ else
+ nTrials = 1;
+ end
+
+ if nargin > 6 && ~isempty(c)
+ validateattributes(c, "numeric", {'scalar', 'nonempty', 'nonnan', 'finite', 'real', 'positive'}, mfilename, "The speed of sound", 7);
+ else
+ c = 343;
+ end
+
+
+ % ====================================================
+ % Pre-process data
+ % ====================================================
+ % Generate source signals
+ if isvector(Q)
+ tmp = rand(sigLen, size(sPos, 1), nTrials);
+ tmp = tmp - mean(tmp);
+ tmp = tmp./max(abs(tmp));
+ Q = Q(:).' .* tmp;
+ elseif sigLen ~= size(Q, 1)
+ if sigLen > size(Q, 1)
+ Q = cat(1, Q, zeros(sigLen - size(Q, 1), size(Q, 2), size(Q, 3)));
+ else
+ Q = Q(1:sigLen, :, :);
+ end
+ end
+
+ % Calculate source-receiver distances
+ dist = twoPtDist(sPos, rPos); % Distances
+ del = dist/c; % Delays in seconds
+
+ % ====================================================
+ % Calculate signals
+ % ====================================================
+ % Go through the trials/sound field realisations
+ for jIdx = nTrials:-1:1
+ % Go through the sources (calculate for all receivers in one go)
+ for sIdx = size(dist, 1):-1:1
+ rSigMtx(:, :, sIdx, jIdx) = (1./dist(sIdx, :)) .* delayseq(Q(:, sIdx, jIdx), del(sIdx, :), fs);
+ end
+ end
+
+ % Sum source signals at each receiver position
+ rSig = squeeze(sum(rSigMtx, 3));
+
+ % ====================================================
+ % Calculate output arguments
+ % ====================================================
+ % Mean receiver signal over all trials/sound field realisations
+ if nargout > 1
+ rSigMean = squeeze(mean(rSig, 3));
+ end
+
+ % Mean receiver signal due to each source over all trials/sound field realisations
+ if nargout > 3
+ rSigMtxMean = squeeze(mean(rSigMtx, 4));
+ end
+end
\ No newline at end of file
diff --git a/Sound Fields/README.md b/Sound Fields/README.md
index 42d3fc0aa29707cf646ceeebe8b2a7207ce163d8..b843dbfea533f78d8725960d91ebfc359adb35c6 100644
--- a/Sound Fields/README.md
+++ b/Sound Fields/README.md
@@ -9,6 +9,7 @@ The current structure of the section is shown below.
- Generation
- Point source
- Frequency domain
+ - Time domain
- Wave domain (spherical harmonics)
- Plane wave
- Frequency domain
@@ -34,5 +35,6 @@ The current structure of the section is shown below.
Some implementations use dependencies from [Utilities/Generic](https://gitlab.com/in-nova/in-nova/-/tree/main/Utilities/Generic?ref_type=heads). The necessary files, as well as the relevant functions are listed below:
- MATLAB
- - `ptSrcField()` requires `twoPtDist()` from [Utilities/Generic](https://gitlab.com/in-nova/in-nova/-/tree/main/Utilities/Generic?ref_type=heads)
- - `planeWaveSH()` requires `sphHarm()` and `idsft()` from [Sound Fields](https://gitlab.com/in-nova/in-nova/-/tree/main/Sound%20Fields) (this folder)
\ No newline at end of file
+ - `ptSrcField()` requires `twoPtDist()` from [Utilities/Generic](https://gitlab.com/in-nova/in-nova/-/tree/main/Utilities/Generic?ref_type=heads).
+ - `planeWaveSH()` requires `sphHarm()` and `idsft()` from [Sound Fields](https://gitlab.com/in-nova/in-nova/-/tree/main/Sound%20Fields) (this folder).
+ - `ptSrcFieldTD()` requires `twoPtDist()` from from [Utilities/Generic](https://gitlab.com/in-nova/in-nova/-/tree/main/Utilities/Generic?ref_type=heads).
diff --git a/Virtual Sensing/README.md b/Virtual Sensing/README.md
index 973f143b5cfa557145038cd44a3427bbce70b40e..fc21835acb4ddb4ce8ac5739e5a18fd2dae14c93 100644
--- a/Virtual Sensing/README.md
+++ b/Virtual Sensing/README.md
@@ -24,6 +24,7 @@ There are no functional dependencies, however, the implementations are based on
- `virtMicGeo()`
- Sound Fields
- `ptSrcField()`
+ - `ptSrcFieldTD()`
- `planeWave()`
## Examples
diff --git a/Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/obsFiltEstTD.m b/Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/obsFiltEstTD.m
new file mode 100644
index 0000000000000000000000000000000000000000..d8827b8322fc1ae625d62cf36626692ae28e9b19
--- /dev/null
+++ b/Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/obsFiltEstTD.m
@@ -0,0 +1,131 @@
+%% Perform estimation with observation filters in the Time-Domain
+% --------------------------------------------------
+% Author: Achilles Kappis
+% e-mail: axilleaz@protonmail.com
+%
+% Date: 14/09/2024 (DD/MM/YYYY)
+%
+% Copyright: MIT
+% --------------------------------------------------
+% Functionality: Perform estimation with observation filters in the
+% Time-Domain.
+% --------------------------------------------------
+% Input
+%
+% 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.
+%
+% O [numeric]: The observation filters. This must be an NxIxMxJ array,
+% where N is the number of virtual microphones and I the
+% length of the observation filters (in samples). Here J can
+% either be equal to the third dimension of "m", or 1 in
+% which case a single set of observation filters will be
+% applied to all trials/realisations of the sound field.
+%
+%
+% e [numeric] (Optional): The measurement(s) at the virtual microphone(s)
+% position(s). This must be a KxNxJ array. If this
+% argument is not provided, the output argument
+% "err" is an empty array.
+%
+% --------------------------------------------------
+% Output
+%
+% estPerMic [numeric]: This array holds the monitoring microphone signals
+% filtered by the corresponding observation filters.
+% It is a KxNxMxJ array.
+%
+% est [numeric]: The estimated signals at the virtual microphone positions.
+% This is an KxNxJ array and corresponds to the sum of
+% "estPerMic" output argument over its third dimension.
+%
+% err [numeric]: The error signals at the virtual microphone positions.
+% This corresponds to the difference between the true, "e"
+% and the estimated, "est", signals. If "e" is not provided
+% this argument is an empty array.
+%
+% estMean [numeric]: These are the estimated virtual microphone signals
+% averaged over the trials/sound field realisations.
+%
+% errMean [numeric]: These are the error signals averaged over the
+% trials/sound field realisations. If "e" is not
+% provided this argument is an empty array.
+%
+% --------------------------------------------------
+% Notes
+%
+% --------------------------------------------------
+function [estPerMic, est, err, estMean, errMean] = obsFiltEstTD(m, O, e)
+ % ====================================================
+ % Check for number of arguments
+ % ====================================================
+ narginchk(2, 3);
+ nargoutchk(0, 5);
+
+ % ====================================================
+ % Validate input arguments
+ % ====================================================
+ % Validate mandatory arguments
+ validateattributes(m, "numeric", {'3d', 'real', 'nonnan', 'nonempty', 'finite'}, mfilename, "Monitoring microphone measurements", 1);
+ validateattributes(O, "numeric", {'real', 'nonnan', 'nonempty', 'finite'}, mfilename, "Observation filters", 2);
+
+ % Check for observation filter size
+ if ndims(O) > 4
+ error("Observation filter array cannot have more than four dimensions. See documentation (at the top of the function file) for more information.");
+ elseif size(O, 4) ~= size(m, 3) && size(O, 4) ~= 1
+ error("The number of observation filter sets must be either equal to the trials/sound field realisations, or 1.")
+ end
+
+ % Validate optional arguments
+ if nargin > 2 && ~isempty(e)
+ validateattributes(e, "numeric", {'3d', 'real', 'nonnan', 'nonempty', 'finite', 'size', [size(m, 1), NaN, size(m, 3)]}, mfilename, "Virtual microphone measurements", 3);
+ else
+ err = [];
+ end
+
+
+ % ====================================================
+ % Pre-process the data
+ % ====================================================
+ % Replicate observation filters if only one set is provided
+ if size(O, 4) ~= size(m, 3)
+ O = repmat(O, 1, 1, 1, size(m, 3));
+ end
+
+ % ====================================================
+ % Filter signals to perform estimation
+ % ====================================================
+ % Go through the trials/sound field realisations
+ for jIdx = size(m, 3):-1:1
+ % Go through the virtual microphone positions
+ for eIdx = size(O, 1):-1:1
+ % Go through the monitoring microphone positions
+ for mIdx = size(O, 3):-1:1
+ estPerMic(:, eIdx, mIdx, jIdx) = filter(O(eIdx, :, mIdx, jIdx), 1, m(:, mIdx, jIdx));
+ end
+ end
+ end
+
+ % Sum the estimates of each monitoring microphone to get the estimated virtual microphone signals
+ if nargout > 1
+ est = squeeze(sum(estPerMic, 3));
+ end
+
+ % Calculate the error signals
+ if nargout > 2 && ~exist("err", "var")
+ err = e - est;
+ end
+
+ % Average of estimated signals
+ if nargout > 3
+ estMean = mean(est, 3);
+ end
+
+ % Mean error signal
+ if nargout > 4
+ errMean = mean(err, 3);
+ end
+end
\ No newline at end of file
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..49b3b4c59d4b4bb76428fcd40facba10fdd1c029
--- /dev/null
+++ b/Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/obsFiltTD.m
@@ -0,0 +1,283 @@
+%% Calculate optimal observation filters in the time domain
+% --------------------------------------------------
+% Author: Achilles Kappis
+% e-mail: axilleaz@protonmail.com
+%
+% Date: 14/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
+% can be either a scalar, in which case the
+% same delay is applied to all virtual
+% microphone signals, or a vector with each
+% element holding the delay for each virtual
+% microphone signal respectively. The delays
+% must be real non-negative values less than
+% or equal to I, the length of the virtual
+% microphone signals. If "fs" is provided,
+% these values correspond to seconds otherwise
+% they are in samples and must be integer.
+% [Default: 0].
+%
+% fs [numeric] (Optional): The sampling frequency. If this value is given,
+% the delay is treated as being in seconds,
+% otherwise it is in samples. This argument must
+% be a positive real integer.
+%
+% --------------------------------------------------
+% Output
+%
+% O [numeric]: The observation filters. This is an NxfirLenxMxJ 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 calculated with the correlation
+% matrices RmeMtx and RmmMtx averaged over the
+% trials/sound field realisations (these are the
+% "RmeMtxMean" and "RmmMtxMean" output arguments). This is
+% an NxfirLenxM 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, calculated with the correlation matrices RmeMtx
+% and RmmMtx averaged over all trials/sound field
+% realisations (the "RmeMtxMean" and "RmmMtxMean" output
+% arguments). The dimensions of this matrix are
+% NxTotFiltLen.
+%
+% 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, fs)
+ % ====================================================
+ % Check for number of arguments
+ % ====================================================
+ narginchk(2, 6);
+ 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', 'real'}, mfilename, "Regularisation factor", 3);
+ else
+ beta = 0;
+ end
+
+ if nargin > 3 && ~isempty(filtLen)
+ validateattributes(filtLen, "numeric", {'scalar', 'positive', 'nonnan', 'real', '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", {'vector', 'nonnegative', 'nonnan', 'finite', 'real', 'nonempty', '<=', size(e, 1)}, mfilename, "Delay to be added to virtual microphone signals to implement the delayed Remote Microphone Technique", 5);
+
+ if isscalar(delay)
+ delay = delay * ones(size(e, 2), 1);
+ elseif length(delay) ~= size(e, 2)
+ error("The 'delay' argument must be either a scalar, or a vector with number of elements equal to the number of virtual microphone signals.");
+ end
+ else
+ delay = zeros(size(e, 2), 1);
+ end
+
+ if nargin > 5 && ~isempty(fs)
+ validateattributes(fs, "numeric", {'scalar', 'integer', 'positive', 'nonnan', 'nonempty', 'finite', 'real'}, mfilename, "The sampling frequency", 6);
+ else
+ fs = [];
+ end
+
+
+ % ====================================================
+ % Pre-process data
+ % ====================================================
+ % Delay the virtual microphone signals
+ if isempty(fs)
+ for jIdx = size(e, 3):-1:1
+ e(:, :, jIdx) = delayseq(e(:, :, jIdx), delay);
+ end
+ else
+ for jIdx = size(e, 3):-1:1
+ e(:, :, jIdx) = delayseq(e(:, :, jIdx), delay, fs);
+ end
+ 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(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
+
+ % Observation filter calculated with the mean Rme and Rmm over the trials/sound field realisation
+ if nargout > 7
+ % Average the RmeMtx and RmmMtx matrices
+ RmeMtxMean = mean(RmeMtx, 3);
+ RmmMtxMean = mean(RmmMtx, 3);
+
+ % Calculate the observation filter
+ Oopt = RmeMtxMean.'/(RmmMtxMean + beta * eye(size(RmmMtxMean)));
+
+ % Reshape
+ Omean = permute(reshape(Oopt.', filtLen, size(m, 2), size(e, 2)), [3, 1, 2]);
+ 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
+end
\ No newline at end of file