Merge branch 'BugsAndFeatures' into 'main'

Various bug fixes and updates

See merge request in-nova/in-nova!9

CHANGELOG.md

+### v0.3.2 ###
+**Signal Processing - Generic**\
+\+ Add option to use negative delay values with the `winSincFracDel()` function.
+
+**Signal Processing - Array Processing**\
+\* Fix error in `arrMan()` in the inner product of the source directions matrix and the sensor position vectors.\
+\* Fix error in `arrMan()` where the returned array manifold had incorrect dimensions. It was $N_{\textrm{d}} \times N_{\textrm{m}} \times N_{\textrm{k}}$, where $N_{\textrm{d}}$ is the number of source directions, $N_{\textrm{m}}$ the number of sensors and $N_{\textrm{k}}$ the number of frequencies. Now it is $N_{\textrm{m}} \times N_{\textrm{d}} \times N_{\textrm{k}}$.
+
+**Utilities - Geometries**\
+\* Fix bug in `virtMicGeo()` where translation and rotation were not applied to the geometry mesh coordinates array.
+
+**Virtual Sensing - Remote Microphone Technique**\
+\* Update `obsFiltEstTD()`, removing the observatoin filter replication when a single filter is provided for many trials/sound field realisations to reduce the memory footprint of the function.\
+\* Update `obsFiltTD()`, performing the delay operation for the virtual microphone signals only when needed and added check for $SNR = \infty$ to calculate the SNR matrix appropriately.
+\* Fix bug in `obsFiltTD()`, returning the condition number without taking into account the regularisation and SNR values.
+
+-------------------------------------------
+
+
 ### v0.3.1 ###
 **Virtual Sensing - Remote Microphone Technique**\
 \+ Add the option to add "noise" to the monitoring microphone signals in the time domain estimation of the optimal observation filters in `obsFiltTD()` to replicate the feature of `obsFilt()`. This change **breaks** backwards compatibility.
 
+-------------------------------------------  
+
 
 ### v0.3.0 ###
 **General**\
@@ -21,6 +42,8 @@
 - `ptSrcFieldSH()`
 - `ptSrcFieldTD()`
 
+-------------------------------------------
+
 
 ### v0.2.8 ###
 **Utilities - Generic**\
@@ -29,17 +52,23 @@
 **Utilities - Geometries**\
 \* Update the geometry generating functions `rcvGeo()`, `srcGeo()` and `virtMicGeo` to accommodate the changes in `rotMat3d()`.
 
+-------------------------------------------
+
 
 ### v0.2.7 ###
 **Signal Processing - Measurements**\
 \+ Introduction of the *Measurements* "topic" in the Signal Processing part of the codebase.\
 \+ Function to estimate the impulse responses from sweep measurements.
 
+-------------------------------------------
+
 
 ### v0.2.6 ###
 **Signal Processing - Generic**\
 \+ Add function to perform fractional octave smoothing of spectra.
 
+-------------------------------------------
+
 
 ### v0.2.5 ###
 **Utilities - Geometries**\
@@ -47,6 +76,8 @@
 \* Combine translation offsets to a single vector input argument for th `virtMicGeo()` MATLAB function.\
 \* Add the ability to rotate the receiver geometries in `rcvGeo()` MATLAB function.
 
+-------------------------------------------
+
 
 ### v0.2.4 ###
 **Virtual Sensing**\
@@ -60,6 +91,8 @@
 **Utilities - Geometries**\
 \* The virtual microphone geometry generation function has been heavily modified and its interface has changed in a **backward incompatible** way. It is now more generic and allows for "arbitrary" size of each Cartesian dimension of the geometries as well as different number of sensors along each dimension.
 
+-------------------------------------------
+
 
 ### v0.2.3 ###
 **Virtual Sensing**\
@@ -73,16 +106,22 @@
 **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.
 
+-------------------------------------------
+
 
 ### v0.2.0 ###
 **Utilities**\
@@ -132,6 +171,8 @@
 \+ Added function to calculate sound field generated by a point source in the Spherical Harmonics domain (truncated order).\
 \+ Added function to calculate sound field generated by a plane wave in the Spherical Harmonics domain (truncated order).
 
+-------------------------------------------
+
 
 ### v0.1.0 ###
-\* This is the initial version of the project.
+This is the initial version of the project.

README.md

@@ -13,20 +13,20 @@ The current structure of the codebase is shown below. More information can be fo
 - [Utilities](https://gitlab.com/in-nova/in-nova/-/tree/main/Utilities?ref_type=heads)
   - Generic
   - Geometries
-- [Active Noise Control](https://gitlab.com/in-nova/in-nova/-/tree/main/Active%20Noise%20Control)
+- [Active Noise Control](https://gitlab.com/in-nova/in-nova/-/tree/main/Active%20Noise%20Control?ref_type=heads)
   - Optimal tonal control
     - Frequency domain
-- [Virtual Sensing](https://gitlab.com/in-nova/in-nova/-/tree/main/Virtual%20Sensing)
+- [Virtual Sensing](https://gitlab.com/in-nova/in-nova/-/tree/main/Virtual%20Sensing?ref_type=heads)
   - Remote Microphone Technique
     - Frequency domain
     - Time domain
-- [Optimisation](https://gitlab.com/in-nova/in-nova/-/tree/main/Optimisation/Memetic%20Algorithm)
+- [Optimisation](https://gitlab.com/in-nova/in-nova/-/tree/main/Optimisation/Memetic%20Algorithm?ref_type=heads)
   - Memetic Algorithm
-- [Signal Processing](https://gitlab.com/in-nova/in-nova/-/tree/main/Signal%20Processing)
+- [Signal Processing](https://gitlab.com/in-nova/in-nova/-/tree/main/Signal%20Processing?ref_type=heads)
   - Array Processing
   - Generic
   - Measurements
-- [Sound fields](https://gitlab.com/in-nova/in-nova/-/tree/main/Sound%20Fields)
+- [Sound fields](https://gitlab.com/in-nova/in-nova/-/tree/main/Sound%20Fields?ref_type=heads)
   - Generation
     - Frequency domain
     - Time domain
@@ -83,11 +83,11 @@ Part of the `ptsOnSphere()` MATLAB/Octave function is an adaptation of a functio
 
 
 ## License
-The project is licensed under the ***MIT License***, which is a rather unrestrictive license, effectively passing the code to the public scope. For more information see the [MIT License](https://mit-license.org/), or the [LICENSE](https://gitlab.com/in-nova/in-nova/-/blob/main/LICENSE) file of the project.
+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?ref_type=heads) file of the project.
 
 
 ## Versioning ##
-The project uses [semantic versioning](https://semver.org/) and the current version is **v0.3.1**.
+The project uses [semantic versioning](https://semver.org/) and the current version is **v0.3.2**. The project hasn't reached v1.0.0 yet so changes that break backwards compatibility may (and have been) introduced in any version update. For information about the changes in each version consult the [CHANGELOG](https://gitlab.com/in-nova/in-nova/-/blob/main/CHANGELOG.md?ref_type=heads).
 
 
 #### **Important**

Signal Processing/Array Processing/MATLAB/Functions/arrMan.m

@@ -3,7 +3,7 @@
 % Author: Achilles Kappis
 % e-mail: axilleaz@protonmail.com
 %
-% Date: 29/09/2024 (DD/MM/YYYY)
+% Date: 05/10/2024 (DD/MM/YYYY)
 %
 % Copyright: MIT
 % --------------------------------------------------
@@ -64,10 +64,10 @@ function am = arrMan(mPos, dirs, k)
     [x, y, z] = sph2cart(dirs(1, :), dirs(2, :), ones(size(dirs(1, :))));
 
     % Inner product of directions and sensor position
-    am = [x; y; z] * mPos;
+    am = [x; y; z].' * mPos;
 
     % Multiply with wavenumbers
-    am = am .* permute(k(:), [3, 2, 1]);
+    am = am.' .* permute(k(:), [3, 2, 1]);
 
     % Calculate the exponential
     am = exp(1i * am);

Signal Processing/Generic/MATLAB/Functions/winSincFracDel.m

@@ -3,7 +3,7 @@
 % Author: Achilles Kappis
 % e-mail: axilleaz@protonmail.com
 %
-% Date: 13/08/2024 (DD/MM/YYYY)
+% Date: 05/10/2024 (DD/MM/YYYY)
 %
 % Copyright: MIT
 % --------------------------------------------------
@@ -12,14 +12,13 @@
 % --------------------------------------------------
 % Input
 % 
-% del [numeric]: The delay in samples. This must be a non-negative real
-%                scalar.
+% del [numeric]: The delay in samples. This must be a real scalar.
 % 
 % len [numeric] (Optional): This is the length of the filter. This must be
-%                           a positive real integer. If this values is less
-%                           than the "delay" value, the filter length
-%                           equals the "delay" value rounded up.
-%                           [Default: ceil(delay)].
+%                           a positive real integer. If this value is less
+%                           than the absolute value of "del", the filter    
+%                           length equals the absolute value of the "del",
+%                           rounded up. [Default: ceil(abs(del))].
 % 
 % winFun [string/char, numeric] (Optional): Windowing functions to be
 %                                           applied to the truncated sinc
@@ -57,6 +56,8 @@
 % 
 % causDel [numeric]: The delay required to make the filter causal. This
 %                    value corresponds to half the length of the filter.
+%                    The causal part of the filter are at indices
+%                    (causDel:len).
 % 
 % 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
@@ -85,7 +86,8 @@ function [sincFilt, causDel, dSig] = winSincFracDel(del, len, winFun, sig, sigLe
     % Validate input arguments
     % ====================================================
     % Validate mandatory arguments
-    validateattributes(del, "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);
+    validateattributes(del, "numeric", {'scalar', 'real', 'nonnan', 'nonempty', 'finite'}, mfilename, "Delay in samples", 1);
 
     % Validate optional arguments
     if nargin > 1 && ~isempty(len)
@@ -93,10 +95,10 @@ function [sincFilt, causDel, dSig] = winSincFracDel(del, len, winFun, sig, sigLe
 
         if len < del
             warning("winSincFracDel(): The length of the filter cannot be less than the delay, setting equal to 'delay',");
-            len = ceil(del);
+            len = ceil(abs(del));
         end
     else
-        len = ceil(del);
+        len = ceil(abs(del));
     end
 
     if nargin > 2 && ~isempty(winFun)
@@ -145,7 +147,7 @@ function [sincFilt, causDel, dSig] = winSincFracDel(del, len, winFun, sig, sigLe
     % ====================================================
     idx = floor(-len/2 + 1):floor(len/2); % Filter sample indices
     
-    causDel = -idx(1); % The causal delay of the filter
+    causDel = -idx(1) + 1; % The causal delay of the filter
     sincFilt = sinc(idx - fracDel); % Calculate the filter
 
 
@@ -165,14 +167,20 @@ function [sincFilt, causDel, dSig] = winSincFracDel(del, len, winFun, sig, sigLe
         end
     else
         % Kaiser window
-        win = besseli(0, pi * winFun * sqrt(1 - (2 * (idx - fracDel)/len).^2))/besseli(0, pi * winFun);
+        win = besseli(0, pi * winFun * sqrt(1 - (2 * (idx - abs(fracDel))/len).^2))/besseli(0, pi * winFun);
     end
     
     % Apply window to the filter
     sincFilt = sincFilt .* win;
+    sincFilt = sincFilt(:);
 
     % Add the integral delay
-    sincFilt = [zeros(1, fix(del)), sincFilt(1:end - fix(del))];
+    % sincFilt = [zeros(1, fix(del)), sincFilt(1:end - fix(del))];
+    if del >= 0
+        sincFilt = [zeros(fix(del), 1); sincFilt(1:end - fix(del))];
+    else
+        sincFilt = [sincFilt(abs(floor(del)):end); zeros(abs(floor(del)), 1)];
+    end
 
 
     % ====================================================

Utilities/Geometries/MATLAB/Functions/virtMicGeo.m

@@ -3,7 +3,7 @@
 % Author: Achilles Kappis
 % e-mail: axilleaz@protonmail.com
 %
-% Date: 01/10/2024 (DD/MM/YYYY)
+% Date: 10/10/2024 (DD/MM/YYYY)
 %
 % Copyright: MIT
 % --------------------------------------------------
@@ -198,7 +198,7 @@ function [vPos, vPosMesh] = virtMicGeo(gType, geoDim, nSens, trans, rot, rotOrd)
 
     % For "Cube" geometry provide the coordinates in a "mesh format"
     if nargout > 1 && strcmpi(gType, "Cube")
-        vPosMesh = cat(3, x, y, z);
+        vPosMesh = reshape(vPos.', size(x, 1), size(x, 2), 3);
     elseif nargout > 1
         vPosMesh = [];
     end

Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/obsFiltEstTD.m

@@ -3,7 +3,7 @@
 % Author: Achilles Kappis
 % e-mail: axilleaz@protonmail.com
 %
-% Date: 17/09/2024 (DD/MM/YYYY)
+% Date: 13/10/2024 (DD/MM/YYYY)
 %
 % Copyright: MIT
 % --------------------------------------------------
@@ -87,14 +87,6 @@ function [estPerMic, est, err, estMean, errMean] = obsFiltEstTD(m, O, e)
     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
     % ====================================================
@@ -102,9 +94,16 @@ function [estPerMic, est, err, estMean, errMean] = obsFiltEstTD(m, O, e)
     for jIdx = size(m, 3):-1:1
         % Go through the virtual microphone positions
         for eIdx = size(O, 1):-1:1
+            % Pick the right set of filters (for this trial/"realisation")
+            if size(O, 4) ~= size(m, 3)
+                tmpO = squeeze(O(eIdx, :, :));
+            else
+                tmpO = squeeze(O(eIdx, :, :, jIdx));
+            end
+
             % 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));
+                estPerMic(:, eIdx, mIdx, jIdx) = filter(tmpO(:, mIdx), 1, m(:, mIdx, jIdx));
             end
         end
     end

Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/obsFiltTD.m

@@ -3,7 +3,7 @@
 % Author: Achilles Kappis
 % e-mail: axilleaz@protonmail.com
 %
-% Date: 05/10/2024 (DD/MM/YYYY)
+% Date: 14/10/2024 (DD/MM/YYYY)
 %
 % Copyright: MIT
 % --------------------------------------------------
@@ -222,6 +222,7 @@ function [O, Rme, Rmm, Ovec, RmeMtx, RmmMtx, condNum, mMtx, Omean, RmeMean, RmmM
     % Pre-process data
     % ====================================================
     % Delay the virtual microphone signals
+    if logical(delay)
         if isempty(fs)
             for jIdx = size(e, 3):-1:1
                 e(:, :, jIdx) = delayseq(e(:, :, jIdx), delay);
@@ -231,6 +232,7 @@ function [O, Rme, Rmm, Ovec, RmeMtx, RmmMtx, condNum, mMtx, Omean, RmeMean, RmmM
                 e(:, :, jIdx) = delayseq(e(:, :, jIdx), delay, fs);
             end
         end
+    end
 
     % ====================================================
     % Calculate cross- and auto-correlation matrices
@@ -279,13 +281,22 @@ function [O, Rme, Rmm, Ovec, RmeMtx, RmmMtx, condNum, mMtx, Omean, RmeMean, RmmM
     % Regularisation
     regMtx = beta * eye(size(RmmMtx, 1));
 
-
     for jIdx = size(RmmMtx, 3):-1:1
+        if ~isinf(snrVal)
             snrMtx = 10^(-snrVal/10); % Convert SNR to linear value
             snrMtx = (snrMtx^2) * norm(RmmMtx(:, :, jIdx), 'fro')/(size(m, 2)^2); % Calculate the appropriate SNR amplitude values
+        else
+            snrMtx = 0;
+        end
         snrMtx = snrMtx * eye(size(RmmMtx, 1)); % Make the matrix
+        denomMtx = RmmMtx(:, :, jIdx) + regMtx + snrMtx;
 
-        Ovec(:, :, jIdx) = RmeMtx(:, :, jIdx).'/(RmmMtx(:, :, jIdx) + regMtx + snrMtx);
+        Ovec(:, :, jIdx) = RmeMtx(:, :, jIdx).'/denomMtx;
+
+        % Condition number of RmmMtx for each trial/sound field realisation
+        if nargout > 6
+            condNum(jIdx) = cond(denomMtx);
+        end
     end
 
     % "Split" observation filter vector to observation filters per monitoring and virtual microphone
@@ -294,13 +305,6 @@ function [O, Rme, Rmm, Ovec, RmeMtx, RmmMtx, condNum, mMtx, Omean, RmeMean, RmmM
     % ====================================================
     % Provide additional output arguments
     % ====================================================
-    % Condition number of RmmMtx for each trial/sound field realisation
-    if nargout > 6
-        for jIdx = size(RmmMtx, 3):-1:1
-            condNum(jIdx) = cond(RmmMtx(:, :, jIdx));
-        end
-    end
-
     % The monitoring microphone signals vectorised per trial/sound field realisation
     if nargout > 7
         mMtx = reshape(m, prod(size(m, [1, 2])), size(m, 3));