From c14444d8bb2bda1240be4fb7be9bc83364678297 Mon Sep 17 00:00:00 2001
From: ZaellixA <axilleaz@protonmail.com>
Date: Thu, 12 Sep 2024 14:42:44 +0100
Subject: [PATCH 01/24] Add function to calculate the optimal observation
filters in the time domain
---
.../MATLAB/Functions/obsFiltTD.m | 247 ++++++++++++++++++
1 file changed, 247 insertions(+)
create mode 100644 Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/obsFiltTD.m
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 0000000..09a53a4
--- /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
--
GitLab
From dd6cd670120df27f0d6a805a4a4cc4861878132f Mon Sep 17 00:00:00 2001
From: ZaellixA <axilleaz@protonmail.com>
Date: Thu, 12 Sep 2024 14:43:19 +0100
Subject: [PATCH 02/24] Add function to perform estimation with the optimal
observation filters in the time domain
---
.../MATLAB/Functions/obsFiltEstTD.m | 131 ++++++++++++++++++
1 file changed, 131 insertions(+)
create mode 100644 Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/obsFiltEstTD.m
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 0000000..dcd0f6b
--- /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: 12/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 ismatrix(O) || ndims(O) > 4
+ error("Observation filter must be either a three- or four-dimensional array. 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, size(O, 1), size(O, 2), size(O, 3), 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
--
GitLab
From d2280a3efa7b66fe692f01d22db5881babdcf51e Mon Sep 17 00:00:00 2001
From: ZaellixA <axilleaz@protonmail.com>
Date: Thu, 12 Sep 2024 19:18:27 +0100
Subject: [PATCH 03/24] Update winSincFracDel.m to fix errors and bugs and
provide four window functions
---
.../Generic/MATLAB/Functions/winSincFracDel.m | 153 +++++++++---------
1 file changed, 78 insertions(+), 75 deletions(-)
diff --git a/Signal Processing/Generic/MATLAB/Functions/winSincFracDel.m b/Signal Processing/Generic/MATLAB/Functions/winSincFracDel.m
index c2e5f3b..3455345 100644
--- a/Signal Processing/Generic/MATLAB/Functions/winSincFracDel.m
+++ b/Signal Processing/Generic/MATLAB/Functions/winSincFracDel.m
@@ -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
@@ -135,46 +134,50 @@ function [delFilt, intDel, causalFilt, delSig, causalSig] = winSincFracDel(delay
end
% ====================================================
- % Generate FIR filter
+ % Generate FIR filter for fractional delay
% ====================================================
- delFilt = -len + 1:len; delFilt = delFilt.'; % Filter sample indices
+ idx = (-len + 1:len).'; % Filter sample indices
- intDel = -floor(delFilt(1)/2); % The causal delay of the filter
- delFilt = sinc(delFilt - delay); % Calculate the filter
+ causDel = -floor(idx(1)/2); % The causal delay of the filter
+ sincFilt = sinc(idx - del); % Calculate the filter
- % 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;
-
% ====================================================
- % Filter the signals
+ % Apply windowing
% ====================================================
- % Return the shifted filter
- if nargout > 2
- causalFilt = delFilt(intDel:end);
- end
+ fracDel = abs(rem(del, 1)); % Total filter delay
- 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 "hann"
+ win = sin(pi * (idx - fracDel)/len).^2;
+ case "hamming"
+ alpha = 25/46;
+ win = alpha - (1 - alpha) * cos(2 * pi * (idx - fracDel)/len);
+ otherwise % Rectangular
+ win = ones(length(idx), 1);
end
+ else
+ % Kaiser window
+ win = besseli(0, pi * winFun * sqrt(1 - (2 * (idx - del)/len).^2))/besseli(0, pi * winFun);
end
-
- % Shift signal to the left by the amount of the "causal" delay of the filter
- if nargout > 4
- causalSig = delSig(intDel:end, :);
+
+ % Shift window function to centre on sinc's peak
+ if ~isscalar(winFun)
+ win = delayseq(win, floor(del) + causDel);
end
- % 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), :);
+ % Apply window function and keep only the causal part up to the length of the filter
+ sincFilt = sincFilt(causDel + 1:causDel + len) .* win(causDel + 1:causDel + len);
+
+
+ % ====================================================
+ % 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
--
GitLab
From 83443ddc9b9287b7ad928d4b4a3b18572d063680 Mon Sep 17 00:00:00 2001
From: ZaellixA <axilleaz@protonmail.com>
Date: Thu, 12 Sep 2024 19:44:26 +0100
Subject: [PATCH 04/24] Update CHANGELOG.md with latest changes and additions
---
CHANGELOG.md | 15 +++++++++++++++
1 file changed, 15 insertions(+)
diff --git a/CHANGELOG.md b/CHANGELOG.md
index daf1f79..6d603ed 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,3 +1,18 @@
+### 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.\
+
+
+### 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.
--
GitLab
From ef369cf50794d62561d59973d18772f846480658 Mon Sep 17 00:00:00 2001
From: ZaellixA <axilleaz@protonmail.com>
Date: Thu, 12 Sep 2024 19:44:53 +0100
Subject: [PATCH 05/24] Increment project version README.md
---
README.md | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/README.md b/README.md
index fd23c6f..3a2a497 100644
--- a/README.md
+++ b/README.md
@@ -72,7 +72,7 @@ The project is licensed under the ***MIT License***, which is a rather unrestric
## Versioning ##
-The project uses [semantic versioning](https://semver.org/) and the current version is **v0.2.2**.
+The project uses [semantic versioning](https://semver.org/) and the current version is **v0.2.3**.
#### **Important**
--
GitLab
From 35787749159e7e5cf672d6a1999ab6fb99643921 Mon Sep 17 00:00:00 2001
From: ZaellixA <axilleaz@protonmail.com>
Date: Thu, 12 Sep 2024 23:41:08 +0100
Subject: [PATCH 06/24] Update project's content with time domain formulation
of virtual sensing and updated roadmap section in README.md
---
README.md | 8 ++++++--
1 file changed, 6 insertions(+), 2 deletions(-)
diff --git a/README.md b/README.md
index 3a2a497..8619738 100644
--- a/README.md
+++ b/README.md
@@ -1,4 +1,4 @@
-# IN-NOVA
+C# IN-NOVA
This repository holds the open-access codebase used within the [IN-NOVA](https://in-nova-horizon.eu/) project. The codebase is divided into various sections, each relevant to a different topic/problem. The sections will be updated and more sections will be added on an as-needed basis, since the purpose of the project is not to create a software library to tackle the scientific problems.
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).
@@ -17,6 +17,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)
@@ -34,7 +36,9 @@ If you find any bugs, other issues and/or inconsistencies you can [e-mail](mailt
## 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:
--
GitLab
From dbe6477388d5e34f7de761287981f36e9edc5a42 Mon Sep 17 00:00:00 2001
From: ZaellixA <axilleaz@protonmail.com>
Date: Fri, 13 Sep 2024 12:37:48 +0100
Subject: [PATCH 07/24] Improved speed and clarity of code for the
winSincFracDel.m and update date
---
.../Generic/MATLAB/Functions/winSincFracDel.m | 36 +++++++++----------
1 file changed, 18 insertions(+), 18 deletions(-)
diff --git a/Signal Processing/Generic/MATLAB/Functions/winSincFracDel.m b/Signal Processing/Generic/MATLAB/Functions/winSincFracDel.m
index 3455345..5f8b3f4 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
% --------------------------------------------------
@@ -133,43 +133,43 @@ function [sincFilt, causDel, dSig] = winSincFracDel(del, len, winFun, sig, sigLe
sigLen = "Full";
end
+
+ % ====================================================
+ % Calculate some parameters
+ % ====================================================
+ fracDel = del - abs(rem(intDel, 1)); % Fractional part of the delay
+
+
% ====================================================
% Generate FIR filter for fractional delay
% ====================================================
- idx = (-len + 1:len).'; % Filter sample indices
+ idx = floor(-len/2 + 1):floor(len/2); % Filter sample indices
- causDel = -floor(idx(1)/2); % The causal delay of the filter
- sincFilt = sinc(idx - del); % Calculate the filter
+ causDel = -idx(1); % The causal delay of the filter
+ sincFilt = sinc(idx - fracDel); % Calculate the filter
% ====================================================
% Apply windowing
% ====================================================
- fracDel = abs(rem(del, 1)); % Total filter delay
-
% Generate window
if isstring(winFun)
switch lower(winFun)
+ case "rectangular"
+ win = ones(length(idx), 1);
case "hann"
- win = sin(pi * (idx - fracDel)/len).^2;
+ win = sin(pi * (idx - fracDel + causDel + 1)/len).^2;
case "hamming"
alpha = 25/46;
- win = alpha - (1 - alpha) * cos(2 * pi * (idx - fracDel)/len);
- otherwise % Rectangular
- win = ones(length(idx), 1);
+ 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 - del)/len).^2))/besseli(0, pi * winFun);
+ win = besseli(0, pi * winFun * sqrt(1 - (2 * (idx - fracDel)/len).^2))/besseli(0, pi * winFun);
end
- % Shift window function to centre on sinc's peak
- if ~isscalar(winFun)
- win = delayseq(win, floor(del) + causDel);
- end
-
- % Apply window function and keep only the causal part up to the length of the filter
- sincFilt = sincFilt(causDel + 1:causDel + len) .* win(causDel + 1:causDel + len);
+ % Apply window to the filter
+ sincFilt = sincFilt .* win;
% ====================================================
--
GitLab
From 24581519e2ff1a96af1516b5a848821123106918 Mon Sep 17 00:00:00 2001
From: ZaellixA <axilleaz@protonmail.com>
Date: Fri, 13 Sep 2024 16:35:17 +0100
Subject: [PATCH 08/24] Fixed bugs in winSincFracDel.m where no integral delay
was applied and corrected calculation of delayed signals
---
.../Generic/MATLAB/Functions/winSincFracDel.m | 13 +++++++++----
1 file changed, 9 insertions(+), 4 deletions(-)
diff --git a/Signal Processing/Generic/MATLAB/Functions/winSincFracDel.m b/Signal Processing/Generic/MATLAB/Functions/winSincFracDel.m
index 5f8b3f4..dde7fd2 100644
--- a/Signal Processing/Generic/MATLAB/Functions/winSincFracDel.m
+++ b/Signal Processing/Generic/MATLAB/Functions/winSincFracDel.m
@@ -137,9 +137,9 @@ function [sincFilt, causDel, dSig] = winSincFracDel(del, len, winFun, sig, sigLe
% ====================================================
% Calculate some parameters
% ====================================================
- fracDel = del - abs(rem(intDel, 1)); % Fractional part of the delay
+ fracDel = rem(del, 1); % Fractional part of the delay
+
-
% ====================================================
% Generate FIR filter for fractional delay
% ====================================================
@@ -156,7 +156,7 @@ function [sincFilt, causDel, dSig] = winSincFracDel(del, len, winFun, sig, sigLe
if isstring(winFun)
switch lower(winFun)
case "rectangular"
- win = ones(length(idx), 1);
+ win = ones(1, length(idx));
case "hann"
win = sin(pi * (idx - fracDel + causDel + 1)/len).^2;
case "hamming"
@@ -171,13 +171,18 @@ function [sincFilt, causDel, dSig] = winSincFracDel(del, len, winFun, sig, sigLe
% Apply window to the filter
sincFilt = sincFilt .* win;
+ % Add the integral delay
+ intDel = zeros(fix(del) + 1, 1); intDel(end) = 1;
+ sincFilt = conv(sincFilt, intDel);
+ sincFilt = sincFilt(1:len);
+
% ====================================================
% Filter the signals
% ====================================================
if nargout > 2
for idx = size(sig, 2):-1:1
- dSig(:, idx) = conv(sig, sincFilt, lower(sigLen));
+ dSig(:, idx) = conv(sig, sincFilt.', lower(sigLen));
end
end
end
\ No newline at end of file
--
GitLab
From 40586023a923683ea225d82f5818a9694596d2a2 Mon Sep 17 00:00:00 2001
From: ZaellixA <axilleaz@protonmail.com>
Date: Fri, 13 Sep 2024 16:48:05 +0100
Subject: [PATCH 09/24] Add function to calculate signals from ideal point
sources in the time domain
---
Sound Fields/MATLAB/Functions/ptSrcFieldTD.m | 148 +++++++++++++++++++
1 file changed, 148 insertions(+)
create mode 100644 Sound Fields/MATLAB/Functions/ptSrcFieldTD.m
diff --git a/Sound Fields/MATLAB/Functions/ptSrcFieldTD.m b/Sound Fields/MATLAB/Functions/ptSrcFieldTD.m
new file mode 100644
index 0000000..909e780
--- /dev/null
+++ b/Sound Fields/MATLAB/Functions/ptSrcFieldTD.m
@@ -0,0 +1,148 @@
+%% Signal generated by point sources in the time-domain
+% --------------------------------------------------
+% Author: Achilles Kappis
+% e-mail: axilleaz@protonmail.com
+%
+% Date: 13/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, or a real
+% positive integer denoting the length of the
+% source signals. In case the length of the source
+% signals is given, the generated signals are
+% zero-mean, (approximately) uniformly distributed
+% and normalised to unity. [Default: 128].
+%
+% 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.
+%
+% --------------------------------------------------
+% Notes
+%
+% - It is the responsibility of the user to pick (or provide) adequately
+% long source signals
+%
+% --------------------------------------------------
+function [rSig, rSigMean, rSigMtx, rSigMtxMean] = ptSrcFieldTD(sPos, rPos, fs, Q, nTrials, c)
+ % ====================================================
+ % Check for number of arguments
+ % ====================================================
+ narginchk(3, 6);
+ nargoutchk(0, 4);
+
+ % ====================================================
+ % 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);
+ else
+ validateattributes(Q, "numeric", {'3d', 'real', 'nonnan', 'nonempty', 'finite', 'ncols', size(sPos, 1)}, mfilename, "Source signals", 4);
+ nTrials = size(Q, 3);
+ end
+ else
+ Q = 128;
+ end
+
+ if nargin > 4 && ~isempty(nTrials) && isscalar(Q)
+ validateattributes(nTrials, "numeric", {'scalar', 'positive', 'nonnan', 'nonempty', 'finite'}, mfilename, "Number of trials/sound field realisations", 5);
+ elseif nargin > 4 && isempty(nTrials)
+ nTrials = 1;
+ end
+
+ if nargin > 5 && ~isempty(c)
+ validateattributes(c, "numeric", {'scalar', 'nonempty', 'nonnan', 'finite', 'real', 'positive'}, mfilename, "The speed of sound", 6);
+ else
+ c = 343;
+ end
+
+
+ % ====================================================
+ % Pre-process data
+ % ====================================================
+ % Generate source signals
+ if isscalar(Q)
+ Q = rand(Q, size(sPos, 1), nTrials);
+ Q = Q - mean(Q);
+ Q = Q./max(abs(Q));
+ 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 = size(Q, 3):-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
--
GitLab
From 0bfc7c393e8dc100b68a017a5f84bfad589e4484 Mon Sep 17 00:00:00 2001
From: ZaellixA <axilleaz@protonmail.com>
Date: Fri, 13 Sep 2024 16:53:56 +0100
Subject: [PATCH 10/24] Update header comments in ptSrcFieldTD.m to state
dependencies
---
Sound Fields/MATLAB/Functions/ptSrcFieldTD.m | 3 +++
1 file changed, 3 insertions(+)
diff --git a/Sound Fields/MATLAB/Functions/ptSrcFieldTD.m b/Sound Fields/MATLAB/Functions/ptSrcFieldTD.m
index 909e780..8301398 100644
--- a/Sound Fields/MATLAB/Functions/ptSrcFieldTD.m
+++ b/Sound Fields/MATLAB/Functions/ptSrcFieldTD.m
@@ -61,6 +61,9 @@
% --------------------------------------------------
% 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
%
--
GitLab
From 8f27bbf4ab4410cb7e063170a619c177585be7bd Mon Sep 17 00:00:00 2001
From: ZaellixA <axilleaz@protonmail.com>
Date: Fri, 13 Sep 2024 16:56:57 +0100
Subject: [PATCH 11/24] Update Sound Fields README.md to include time domain
function information and dependencies
---
Sound Fields/README.md | 6 ++++--
1 file changed, 4 insertions(+), 2 deletions(-)
diff --git a/Sound Fields/README.md b/Sound Fields/README.md
index 42d3fc0..b843dbf 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).
--
GitLab
From f212a90cfbbcc5622f4dd4bd70c25526e160b61a Mon Sep 17 00:00:00 2001
From: ZaellixA <axilleaz@protonmail.com>
Date: Fri, 13 Sep 2024 16:59:09 +0100
Subject: [PATCH 12/24] Update project README.md to included info on latest
Sound Field additions and some of the old functions added in 0.2.0 and
corrected a typo at the top of the file
---
README.md | 7 ++++++-
1 file changed, 6 insertions(+), 1 deletion(-)
diff --git a/README.md b/README.md
index 8619738..d7631f3 100644
--- a/README.md
+++ b/README.md
@@ -1,4 +1,4 @@
-C# IN-NOVA
+# IN-NOVA
This repository holds the open-access codebase used within the [IN-NOVA](https://in-nova-horizon.eu/) project. The codebase is divided into various sections, each relevant to a different topic/problem. The sections will be updated and more sections will be added on an as-needed basis, since the purpose of the project is not to create a software library to tackle the scientific problems.
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).
@@ -26,6 +26,11 @@ 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.
--
GitLab
From 1dd73eec49ae37db8413037598a13e92d70851b9 Mon Sep 17 00:00:00 2001
From: ZaellixA <axilleaz@protonmail.com>
Date: Fri, 13 Sep 2024 16:59:50 +0100
Subject: [PATCH 13/24] Added time domain point source function info in the
CHANGELOG for the 0.2.3 version
---
CHANGELOG.md | 3 +++
1 file changed, 3 insertions(+)
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 6d603ed..8883239 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -7,6 +7,9 @@
\* 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**\
--
GitLab
From c73cc7eb76740224a6e148eb2112f919aab75ea0 Mon Sep 17 00:00:00 2001
From: ZaellixA <axilleaz@protonmail.com>
Date: Fri, 13 Sep 2024 17:18:51 +0100
Subject: [PATCH 14/24] Add the source signals as output argument for
ptSrcFieldTD.m
---
Sound Fields/MATLAB/Functions/ptSrcFieldTD.m | 8 ++++++--
1 file changed, 6 insertions(+), 2 deletions(-)
diff --git a/Sound Fields/MATLAB/Functions/ptSrcFieldTD.m b/Sound Fields/MATLAB/Functions/ptSrcFieldTD.m
index 8301398..f470b8b 100644
--- a/Sound Fields/MATLAB/Functions/ptSrcFieldTD.m
+++ b/Sound Fields/MATLAB/Functions/ptSrcFieldTD.m
@@ -58,6 +58,10 @@
% 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
%
@@ -68,12 +72,12 @@
% long source signals
%
% --------------------------------------------------
-function [rSig, rSigMean, rSigMtx, rSigMtxMean] = ptSrcFieldTD(sPos, rPos, fs, Q, nTrials, c)
+function [rSig, rSigMean, rSigMtx, rSigMtxMean, Q] = ptSrcFieldTD(sPos, rPos, fs, Q, nTrials, c)
% ====================================================
% Check for number of arguments
% ====================================================
narginchk(3, 6);
- nargoutchk(0, 4);
+ nargoutchk(0, 5);
% ====================================================
% Validate input arguments
--
GitLab
From f54889016b0b3f43b86222871215d56128d073c6 Mon Sep 17 00:00:00 2001
From: ZaellixA <axilleaz@protonmail.com>
Date: Fri, 13 Sep 2024 21:28:59 +0100
Subject: [PATCH 15/24] Fixed bug in input argument checks in ptSrcFieldTD.m
---
Sound Fields/MATLAB/Functions/ptSrcFieldTD.m | 11 ++++++-----
1 file changed, 6 insertions(+), 5 deletions(-)
diff --git a/Sound Fields/MATLAB/Functions/ptSrcFieldTD.m b/Sound Fields/MATLAB/Functions/ptSrcFieldTD.m
index f470b8b..3a1ccfc 100644
--- a/Sound Fields/MATLAB/Functions/ptSrcFieldTD.m
+++ b/Sound Fields/MATLAB/Functions/ptSrcFieldTD.m
@@ -93,15 +93,16 @@ function [rSig, rSigMean, rSigMtx, rSigMtxMean, Q] = ptSrcFieldTD(sPos, rPos, fs
validateattributes(Q, "numeric", {'scalar', 'real', 'nonnan', 'nonempty', 'finite', 'positive', 'integer'}, mfilename, "Length of source signals in samples", 4);
else
validateattributes(Q, "numeric", {'3d', 'real', 'nonnan', 'nonempty', 'finite', 'ncols', size(sPos, 1)}, mfilename, "Source signals", 4);
- nTrials = size(Q, 3);
end
else
Q = 128;
end
- if nargin > 4 && ~isempty(nTrials) && isscalar(Q)
- validateattributes(nTrials, "numeric", {'scalar', 'positive', 'nonnan', 'nonempty', 'finite'}, mfilename, "Number of trials/sound field realisations", 5);
- elseif nargin > 4 && isempty(nTrials)
+ if ~isscalar(Q)
+ nTrials = size(Q, 3);
+ elseif nargin > 4 && ~isempty(nTrials) && isscalar(Q)
+ validateattributes(nTrials, "numeric", {'scalar', 'positive', 'nonnan', 'nonempty', 'finite'}, mfilename, "Number of trials/sound field realisations", 5);
+ else
nTrials = 1;
end
@@ -130,7 +131,7 @@ function [rSig, rSigMean, rSigMtx, rSigMtxMean, Q] = ptSrcFieldTD(sPos, rPos, fs
% Calculate signals
% ====================================================
% Go through the trials/sound field realisations
- for jIdx = size(Q, 3):-1:1
+ 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);
--
GitLab
From f495bca546e0724a29fac1e1d4c12a659771689d Mon Sep 17 00:00:00 2001
From: ZaellixA <axilleaz@protonmail.com>
Date: Fri, 13 Sep 2024 21:32:00 +0100
Subject: [PATCH 16/24] Fixed errors in header documentation, fixed
calculations of mean O to be calculated by averaged Rme and Rmm and
introduced separate (fractional) delays for each virtual microphone in either
samples or seconds
---
.../MATLAB/Functions/obsFiltTD.m | 106 ++++++++++++------
1 file changed, 71 insertions(+), 35 deletions(-)
diff --git a/Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/obsFiltTD.m b/Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/obsFiltTD.m
index 09a53a4..117a978 100644
--- a/Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/obsFiltTD.m
+++ b/Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/obsFiltTD.m
@@ -3,7 +3,7 @@
% Author: Achilles Kappis
% e-mail: axilleaz@protonmail.com
%
-% Date: 12/09/2024 (DD/MM/YYYY)
+% Date: 13/09/2024 (DD/MM/YYYY)
%
% Copyright: MIT
% --------------------------------------------------
@@ -37,14 +37,27 @@
% 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
+% 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. [Default: 0].
+% 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 firLenxNxMxJ array.
+% 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
@@ -78,10 +91,13 @@
% 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.
+% 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
+% 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.
@@ -94,10 +110,11 @@
% 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).
+% 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
@@ -119,11 +136,11 @@
% 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)
+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, 5);
+ narginchk(2, 6);
nargoutchk(0, 13);
% ====================================================
@@ -135,24 +152,50 @@ function [O, Rme, Rmm, Ovec, RmeMtx, RmmMtx, mMtx, Omean, RmeMean, RmmMean, Oopt
% Validate optional arguments
if nargin > 2 && ~isempty(beta)
- validateattributes(beta, "numeric", {'scalar', 'nonnegative', 'nonnan', 'nonempty', 'finite'}, mfilename, "Regularisation factor", 3);
+ 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', 'nonempty', 'finite', '<=', size(m, 1)}, mfilename, "Length of observation filters", 4);
+ 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", {'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);
+ 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 = 0;
+ 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
% ====================================================
@@ -165,7 +208,7 @@ function [O, Rme, Rmm, Ovec, RmeMtx, RmmMtx, mMtx, Omean, RmeMean, RmmMean, Oopt
[corr, lags] = xcorr(m(:, mIdx, jIdx), e(:, eIdx, jIdx));
lIdx = find(lags == 0);
- Rme(:, mIdx, eIdx, jIdx) = corr(delay + lIdx:-1:lIdx-filtLen + 1);
+ Rme(:, mIdx, eIdx, jIdx) = corr(lIdx:-1:lIdx-filtLen + 1);
end
% Go through the monitoring microphones to calculate the monitoring microphone correlation matrices
@@ -215,9 +258,17 @@ function [O, Rme, Rmm, Ovec, RmeMtx, RmmMtx, mMtx, Omean, RmeMean, RmmMean, Oopt
mMtx = reshape(m, prod(size(m, [1, 2])), size(m, 3));
end
- % Mean observation over trials/sound field realisations
+ % Observation filter calculated with the mean Rme and Rmm over the trials/sound field realisation
if nargout > 7
- Omean = mean(O, 4);
+ % 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
@@ -229,19 +280,4 @@ function [O, Rme, Rmm, Ovec, RmeMtx, RmmMtx, mMtx, Omean, RmeMean, RmmMean, Oopt
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
--
GitLab
From 3f28073626dcbee70a22473e364afeb0823c2462 Mon Sep 17 00:00:00 2001
From: ZaellixA <axilleaz@protonmail.com>
Date: Fri, 13 Sep 2024 21:33:03 +0100
Subject: [PATCH 17/24] Fixed error in resising the observation filter to apply
in all sound field realisations
---
.../Remote Microphone Technique/MATLAB/Functions/obsFiltEstTD.m | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/obsFiltEstTD.m b/Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/obsFiltEstTD.m
index dcd0f6b..f5c97da 100644
--- a/Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/obsFiltEstTD.m
+++ b/Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/obsFiltEstTD.m
@@ -92,7 +92,7 @@ function [estPerMic, est, err, estMean, errMean] = obsFiltEstTD(m, O, e)
% ====================================================
% Replicate observation filters if only one set is provided
if size(O, 4) ~= size(m, 3)
- O = repmat(O, size(O, 1), size(O, 2), size(O, 3), size(m, 3));
+ O = repmat(O, 1, 1, 1, size(m, 3));
end
% ====================================================
--
GitLab
From e655b9f458ab8ee1853b0b97a3f7d497283c8798 Mon Sep 17 00:00:00 2001
From: ZaellixA <axilleaz@protonmail.com>
Date: Fri, 13 Sep 2024 21:35:00 +0100
Subject: [PATCH 18/24] Update date in obsFiltEstTD.m
---
.../Remote Microphone Technique/MATLAB/Functions/obsFiltEstTD.m | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/obsFiltEstTD.m b/Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/obsFiltEstTD.m
index f5c97da..57543de 100644
--- a/Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/obsFiltEstTD.m
+++ b/Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/obsFiltEstTD.m
@@ -3,7 +3,7 @@
% Author: Achilles Kappis
% e-mail: axilleaz@protonmail.com
%
-% Date: 12/09/2024 (DD/MM/YYYY)
+% Date: 13/09/2024 (DD/MM/YYYY)
%
% Copyright: MIT
% --------------------------------------------------
--
GitLab
From 77422b4892a5438a0f434eeccee7fe11564b5b58 Mon Sep 17 00:00:00 2001
From: ZaellixA <axilleaz@protonmail.com>
Date: Sat, 14 Sep 2024 09:19:29 +0100
Subject: [PATCH 19/24] Simplify the addition of the integral part of the delay
in winSincFracDel.m
---
Signal Processing/Generic/MATLAB/Functions/winSincFracDel.m | 4 +---
1 file changed, 1 insertion(+), 3 deletions(-)
diff --git a/Signal Processing/Generic/MATLAB/Functions/winSincFracDel.m b/Signal Processing/Generic/MATLAB/Functions/winSincFracDel.m
index dde7fd2..3fcfcd6 100644
--- a/Signal Processing/Generic/MATLAB/Functions/winSincFracDel.m
+++ b/Signal Processing/Generic/MATLAB/Functions/winSincFracDel.m
@@ -172,9 +172,7 @@ function [sincFilt, causDel, dSig] = winSincFracDel(del, len, winFun, sig, sigLe
sincFilt = sincFilt .* win;
% Add the integral delay
- intDel = zeros(fix(del) + 1, 1); intDel(end) = 1;
- sincFilt = conv(sincFilt, intDel);
- sincFilt = sincFilt(1:len);
+ sincFilt = [zeros(1, fix(del)), sincFilt(1:end - fix(del))];
% ====================================================
--
GitLab
From daa7c0d1c6a5d21c60ad0c28864cc2d54b5e7480 Mon Sep 17 00:00:00 2001
From: ZaellixA <axilleaz@protonmail.com>
Date: Sat, 14 Sep 2024 09:20:38 +0100
Subject: [PATCH 20/24] Add the ptSrcFieldTD as a 'functional' dependency for
the virtual sensing topic in its README.md
---
Virtual Sensing/README.md | 1 +
1 file changed, 1 insertion(+)
diff --git a/Virtual Sensing/README.md b/Virtual Sensing/README.md
index 973f143..fc21835 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
--
GitLab
From 52f422b6f86d4f90a549bae7e4aa256de4c5f3c1 Mon Sep 17 00:00:00 2001
From: ZaellixA <axilleaz@protonmail.com>
Date: Sat, 14 Sep 2024 09:21:22 +0100
Subject: [PATCH 21/24] Slightly reformatted project README.md
---
README.md | 7 ++++++-
1 file changed, 6 insertions(+), 1 deletion(-)
diff --git a/README.md b/README.md
index d7631f3..ed4f39a 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.
@@ -32,12 +33,15 @@ The current structure of the codebase is shown below. More information can be fo
- 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.
@@ -45,6 +49,7 @@ Although the repository is not considered to be *under construction* anymore, th
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:
@@ -55,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.
@@ -79,7 +85,6 @@ 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.3**.
--
GitLab
From 0d87085ab23fb9dfc73cb5e81d435f5b38aac3f8 Mon Sep 17 00:00:00 2001
From: ZaellixA <axilleaz@protonmail.com>
Date: Sat, 14 Sep 2024 10:36:18 +0100
Subject: [PATCH 22/24] Update date in obsFiltTD.m
---
.../Remote Microphone Technique/MATLAB/Functions/obsFiltTD.m | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/obsFiltTD.m b/Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/obsFiltTD.m
index 117a978..49b3b4c 100644
--- a/Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/obsFiltTD.m
+++ b/Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/obsFiltTD.m
@@ -3,7 +3,7 @@
% Author: Achilles Kappis
% e-mail: axilleaz@protonmail.com
%
-% Date: 13/09/2024 (DD/MM/YYYY)
+% Date: 14/09/2024 (DD/MM/YYYY)
%
% Copyright: MIT
% --------------------------------------------------
--
GitLab
From 50f33059912ac45be540031ecea9906b24cf7241 Mon Sep 17 00:00:00 2001
From: ZaellixA <axilleaz@protonmail.com>
Date: Sat, 14 Sep 2024 10:37:28 +0100
Subject: [PATCH 23/24] Fix input argument checking to allow single monitoring
microphone and source calculations
---
.../MATLAB/Functions/obsFiltEstTD.m | 6 +++---
1 file changed, 3 insertions(+), 3 deletions(-)
diff --git a/Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/obsFiltEstTD.m b/Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/obsFiltEstTD.m
index 57543de..d8827b8 100644
--- a/Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/obsFiltEstTD.m
+++ b/Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/obsFiltEstTD.m
@@ -3,7 +3,7 @@
% Author: Achilles Kappis
% e-mail: axilleaz@protonmail.com
%
-% Date: 13/09/2024 (DD/MM/YYYY)
+% Date: 14/09/2024 (DD/MM/YYYY)
%
% Copyright: MIT
% --------------------------------------------------
@@ -73,8 +73,8 @@ function [estPerMic, est, err, estMean, errMean] = obsFiltEstTD(m, O, e)
validateattributes(O, "numeric", {'real', 'nonnan', 'nonempty', 'finite'}, mfilename, "Observation filters", 2);
% Check for observation filter size
- if ismatrix(O) || ndims(O) > 4
- error("Observation filter must be either a three- or four-dimensional array. See documentation (at the top of the function file) for more information.");
+ 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
--
GitLab
From 89e91a97ae41941ed0268012c4576cdaccdbef2a Mon Sep 17 00:00:00 2001
From: ZaellixA <axilleaz@protonmail.com>
Date: Sat, 14 Sep 2024 22:49:03 +0100
Subject: [PATCH 24/24] Separate the source strengths and signal length input
arguments, updated input argument checks and source signal calculations
---
Sound Fields/MATLAB/Functions/ptSrcFieldTD.m | 69 ++++++++++++++------
1 file changed, 50 insertions(+), 19 deletions(-)
diff --git a/Sound Fields/MATLAB/Functions/ptSrcFieldTD.m b/Sound Fields/MATLAB/Functions/ptSrcFieldTD.m
index 3a1ccfc..01fa2c6 100644
--- a/Sound Fields/MATLAB/Functions/ptSrcFieldTD.m
+++ b/Sound Fields/MATLAB/Functions/ptSrcFieldTD.m
@@ -3,7 +3,7 @@
% Author: Achilles Kappis
% e-mail: axilleaz@protonmail.com
%
-% Date: 13/09/2024 (DD/MM/YYYY)
+% Date: 14/09/2024 (DD/MM/YYYY)
%
% Copyright: MIT
% --------------------------------------------------
@@ -25,12 +25,23 @@
% 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, or a real
-% positive integer denoting the length of the
-% source signals. In case the length of the source
-% signals is given, the generated signals are
-% zero-mean, (approximately) uniformly distributed
-% and normalised to unity. [Default: 128].
+% 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"
@@ -72,11 +83,11 @@
% long source signals
%
% --------------------------------------------------
-function [rSig, rSigMean, rSigMtx, rSigMtxMean, Q] = ptSrcFieldTD(sPos, rPos, fs, Q, nTrials, c)
+function [rSig, rSigMean, rSigMtx, rSigMtxMean, Q] = ptSrcFieldTD(sPos, rPos, fs, Q, sigLen, nTrials, c)
% ====================================================
% Check for number of arguments
% ====================================================
- narginchk(3, 6);
+ narginchk(3, 7);
nargoutchk(0, 5);
% ====================================================
@@ -91,23 +102,36 @@ function [rSig, rSigMean, rSigMtx, rSigMtxMean, Q] = ptSrcFieldTD(sPos, rPos, fs
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 = 128;
+ Q = 1;
end
- if ~isscalar(Q)
+ 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 > 4 && ~isempty(nTrials) && isscalar(Q)
- validateattributes(nTrials, "numeric", {'scalar', 'positive', 'nonnan', 'nonempty', 'finite'}, mfilename, "Number of trials/sound field realisations", 5);
+ 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 > 5 && ~isempty(c)
- validateattributes(c, "numeric", {'scalar', 'nonempty', 'nonnan', 'finite', 'real', 'positive'}, mfilename, "The speed of sound", 6);
+ if nargin > 6 && ~isempty(c)
+ validateattributes(c, "numeric", {'scalar', 'nonempty', 'nonnan', 'finite', 'real', 'positive'}, mfilename, "The speed of sound", 7);
else
c = 343;
end
@@ -117,10 +141,17 @@ function [rSig, rSigMean, rSigMtx, rSigMtxMean, Q] = ptSrcFieldTD(sPos, rPos, fs
% Pre-process data
% ====================================================
% Generate source signals
- if isscalar(Q)
- Q = rand(Q, size(sPos, 1), nTrials);
- Q = Q - mean(Q);
- Q = Q./max(abs(Q));
+ 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
--
GitLab