Merge branch 'ObsFiltUpdate' into 'main'

Update the MATLAB implementation of the optimal frequency domain observation filter

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

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.1**.
+The project uses [semantic versioning](https://semver.org/) and the current version is **v0.2.2**.
 
 
 #### **Important**

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

@@ -3,7 +3,7 @@
 % Author: Achilles Kappis
 % e-mail: axilleaz@protonmail.com
 %
-% Date: 01/09/2024 (DD/MM/YYYY)
+% Date: 04/09/2024 (DD/MM/YYYY)
 %
 % Copyright: MIT
 % --------------------------------------------------
@@ -22,10 +22,10 @@
 %               where K is the number of monitoring microphones and M is
 %               the number of sources.
 %
-% srcCsd [numeric] (Optional): The source cross spectral density matrix.
-%                              This must be a square (MxM) symmetric matrix
-%                              with the cross power spectral density of the
-%                              sources. [Default: eye(M)]
+% Svv [numeric] (Optional): The source cross spectral density matrix. This
+%                           must be a square (MxM) symmetric matrix with
+%                           the cross power spectral density of the
+%                           sources. [Default: eye(M)]
 % 
 % regFacs [numeric] (Optional): The regularisation factors used for the
 %                               inversion of the cross-spectra of the
@@ -75,7 +75,7 @@
 %   for local active sound control" by W. Jung, S. J. Elliott and J. Cheer.
 % 
 % --------------------------------------------------
-function [oOpt, Sme, Smm, condNum] = obsFilt(Pe, Pm, srcCsd, regFacs, snrVal)
+function [oOpt, Sme, Smm, condNum] = obsFilt(Pe, Pm, Svv, regFacs, snrVal)
     % ====================================================
     % Check for number of arguments
     % ====================================================
@@ -90,10 +90,17 @@ function [oOpt, Sme, Smm, condNum] = obsFilt(Pe, Pm, srcCsd, regFacs, snrVal)
     validateattributes(Pm, "numeric", {'nonnan', 'nonempty', 'finite'}, mfilename, "Monitoring microphone pressure", 2);
 
     % Validate optional arguments
-    if nargin > 2 && ~isempty(srcCsd)
-        validateattributes(srcCsd, "numeric", {'2d', 'nonnegative', 'nonnan', 'finite'}, mfilename, "Source cross spectral density matrix", 3)
+    if nargin > 2 && ~isempty(Svv)
+        validateattributes(Svv, "numeric", {'2d', 'nonnan', 'finite'}, mfilename, "Source cross spectral density matrix", 3)
+
+        % Check for correct dimensions
+        if diff(size(Svv))
+            error("The source power spectral density matrix must be a square matrix");
+        elseif size(Svv, 1) ~= size(Pe, 1)
+            error("The number of rows of the source power spectral density matrix must be equal to the number of sources");
+        end
     else
-        srcCsd = eye(size(Pe, 2));
+        Svv = eye(size(Pe, 2));
     end
 
     if nargin > 3 && ~isempty(regFacs)
@@ -127,7 +134,7 @@ function [oOpt, Sme, Smm, condNum] = obsFilt(Pe, Pm, srcCsd, regFacs, snrVal)
     % Calculate optimal filters
     % ====================================================
     % Calculate needed quantities
-    tmpVal = srcCsd * Pm';
+    tmpVal = Svv * Pm';
     
     Sme = Pe * tmpVal; % Virtual-Monitor mics cross-spectra
     Smm = Pm * tmpVal; % Monitor-Monitor mics cross-spectra

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

@@ -3,7 +3,7 @@
 % Author: Achilles Kappis
 % e-mail: axilleaz@protonmail.com
 %
-% Date: 23/08/2024 (DD/MM/YYYY)
+% Date: 04/09/2024 (DD/MM/YYYY)
 %
 % Copyright: MIT
 % --------------------------------------------------
@@ -28,10 +28,10 @@
 %                          matrix must be NxM, where N is the number of
 %                          virtual microphones and M the number of sources.
 %
-% srcCsd [numeric] (Optional): The source cross spectral density matrix.
-%                              This must be a square (MxM) symmetric matrix
-%                              with the cross power spectral density of the
-%                              sources. [Default: eye(M)]
+% Svv [numeric] (Optional): The source cross spectral density matrix. This
+%                           must be a square (MxM) symmetric matrix with
+%                           the cross power spectral density of the
+%                           sources. [Default: eye(M)]
 % 
 % --------------------------------------------------
 % Output
@@ -60,7 +60,7 @@
 % Notes
 % 
 % --------------------------------------------------
-function [est, err, errSqr, normErrSqr, See] = obsFiltEst(Pm, O, Pe, srcCsd)
+function [est, err, errSqr, normErrSqr, See] = obsFiltEst(Pm, O, Pe, Svv)
     % ====================================================
     % Check for number of arguments
     % ====================================================
@@ -82,10 +82,17 @@ function [est, err, errSqr, normErrSqr, See] = obsFiltEst(Pm, O, Pe, srcCsd)
         Pe = NaN;
     end
 
-    if nargin > 3 && ~isempty(srcCsd)
-        validateattributes(srcCsd, "numeric", {'2d', 'nonnegative', 'nonnan', 'finite'}, mfilename, "Source cross spectral density matrix", 3)
+    if nargin > 3 && ~isempty(Svv)
+        validateattributes(Svv, "numeric", {'2d', 'nonnan', 'finite'}, mfilename, "Source cross spectral density matrix", 3)
+
+        % Check for correct dimensions
+        if diff(size(Svv))
+            error("The source power spectral density matrix must be a square matrix");
+        elseif size(Svv, 1) ~= size(Pe, 1)
+            error("The number of rows of the source power spectral density matrix must be equal to the number of sources");
+        end
     else
-        srcCsd = eye(size(Pe, 2));
+        Svv = eye(size(Pe, 2));
     end
 
 
@@ -103,7 +110,7 @@ function [est, err, errSqr, normErrSqr, See] = obsFiltEst(Pm, O, Pe, srcCsd)
     % Sum of squared estimation errors
     if nargout > 2
         for eIdx = size(Pe, 1):-1:1
-            errSqr(eIdx) = err(eIdx, :) * srcCsd * err(eIdx, :)';
+            errSqr(eIdx) = err(eIdx, :) * Svv * err(eIdx, :)';
         end
 
         errSqr = errSqr.';
@@ -112,7 +119,7 @@ function [est, err, errSqr, normErrSqr, See] = obsFiltEst(Pm, O, Pe, srcCsd)
     % Normalised squared errors
     if nargout > 3
         for eIdx = length(errSqr):-1:1
-            See(eIdx) = Pe(eIdx, :) * srcCsd * Pe(eIdx, :)';
+            See(eIdx) = Pe(eIdx, :) * Svv * Pe(eIdx, :)';
             normErrSqr(eIdx) = errSqr(eIdx)/See(eIdx);
         end