Select Git revision
-
Ben Anderson authoredBen Anderson authored
obsFilt.m 6.13 KiB
%% Calculate optimal observation filters
% --------------------------------------------------
% Author: Achilles Kappis
% e-mail: axilleaz@protonmail.com
%
% Date: 13/02/2024 (DD/MM/YYYY)
%
% Copyright: MIT
% --------------------------------------------------
% Functionality: Calculate the optimal, in the least squares sense,
% observation filters for the Remote Microphone Technique.
% --------------------------------------------------
% Input
%
% Pe [numeric]: The transfer function from sources to the virtual
% microphones. The dimensions of the matrix must be NxM,
% where N is the number of virtual microphones and M the
% number of sources.
%
% Pm [numeric]: The transfer function from sources to the monitoring
% microphones. The dimensions of the matrix must be KxM,
% 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)]
%
% regFacs [numeric] (Optional): The regularisation factors used for the
% inversion of the cross-spectra of the
% monitoring microphones. Can be a vector
% with number of elements equal to K (number
% of monitoring microphones) or a scalar
% which will result in the same
% regularisation factor for all microphones.
% [Default: 0]
%
% --------------------------------------------------
% Output
%
% oOpt [numeric]: The optimal observation filters for each virtual
% microphone position. The dimensions of the matrix are
% NxM, where N is the number of virtual positions and M the
% number of monitoring microphones.
%
% Sme[numeric]: The cross-spectra between monitoring and virtual
% microphones. The matrix has dimensions NxK where N is the
% number of virtual microphones and M the number of
% monitoring microphones.
%
% Smm[numeric]: The cross-spectra between monitoring microphones. The
% matrix has dimensions KxK, where K is the number of
% monitoring microphones.
%
% est [numeric]: The estimated complex pressure(s) at the virtual
% microphone position(s).
%
% err [numeric]: The error of the estimated complex pressure(s).
%
% errSqr [numeric]: The sum of the squared errors of the estimated
% complex pressure(s).
%
% normErrSqr [numeric]: The sum of the squared errors of the estimated
% complex pressure(s) normalised to the squared
% pressure(s) of the virtual microphones.
%
% See [numeric]: The normalisation factor (cross spectral density of the
% error microphones) for each virtual microphone position.
% "errSqr"./"See" gives "normErrSqr".
%
% condNum [numeric]: The condition number of the regularised monitoring
% microphone cross spectrum matrix.
%
% --------------------------------------------------
% Notes
%
% --------------------------------------------------
function [oOpt, Sme, Smm, est, err, errSqr, normErrSqr, See, condNum] = obsFilt(Pe, Pm, srcCsd, regFacs)
% ====================================================
% Check for number of arguments
% ====================================================
narginchk(2, 4);
nargoutchk(0, 9);
% ====================================================
% Validate input arguments
% ====================================================
% Validate mandatory arguments
validateattributes(Pe, "numeric", {'nonnan', 'nonempty', 'finite'}, mfilename, "Virtual microphone pressure", 1);
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)
else
srcCsd = eye(size(Pe, 2));
end
if nargin > 3 && ~isempty(regFacs)
validateattributes(regFacs, "numeric", {'vector', 'nonnegative', 'nonnan', 'nonempty', 'finite'}, mfilename, "Regularisation factors", 4);
% Make sure regFacs has correct dimensions
if ~isscalar(regFacs) && numel(regFacs) ~= size(Pm, 1)
error("Regularisation factors must be either a scalar or its length must match the number of monitoring microphones.");
end
else
regFacs = 0;
end
% ====================================================
% Calculate optimal filters
% ====================================================
% Calculate needed quantities
tmpVal = srcCsd * Pm';
Sme = Pe * tmpVal; % Virtual-Monitor mics cross-spectra
Smm = Pm * tmpVal; % Monitor-Monitor mics cross-spectra
if isscalar(regFacs)
regMat = diag(ones(size(Smm, 1), 1)) .* regFacs; % Regularisation matrix
else
regMat = diag(regFacs);
end
% Calcualte observation filters
invQty = Smm + regMat;
oOpt = Sme/invQty;
% Calculate pressure estimates
if nargout > 3
est = oOpt * Pm;
else
return;
end
% Calculate errors
if nargout > 4
err = Pe - est; % Estimation error
else
return;
end
if nargout > 5
errSqr = diag(err * srcCsd * err'); % Sum of squared estimation errors
else
return;
end
% Normalised squared errors
if nargout > 6
See = diag(Pe * srcCsd * Pe'); % Error microphone cross spectral density
normErrSqr = errSqr./See; % Normalised error
else
return
end
% Condition number of auto-correlation (power spectrum) matrix
if nargout > 8
condNum = cond(invQty);
end
end