Select Git revision
obsFilt.m 7.05 KiB
%% Calculate optimal observation filters in the frequency domain
% --------------------------------------------------
% Author: Achilles Kappis
% e-mail: axilleaz@protonmail.com
%
% Date: 13/02/2025 (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.
%
% 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
% 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]
%
% snrVal [numeric] (Optional): This is the Root-Mean-Square (RMS)
% Signal-to-Noise Ratio (SNR) in the
% monitoring microphone signals. This must
% be either a scalar indicating the common
% SNR value of the microphone signals, or a
% vector of dimensions Kx1, containing the SNR
% values of each microphone. The values must
% be real or Inf for the noiseless case. See
% the Notes below for references on how this
% value is calculated. [Default: Inf].
%
% --------------------------------------------------
% 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. The matrix that is returned
% includes the regularisation term and the noise term.
%
% condNum [numeric]: The condition number of the regularised monitoring
% microphone cross spectrum matrix.
% % --------------------------------------------------
% Notes
%
% - The SNR at the monitoring microphones calculation is based on the
% paper: "Combining the remote microphone technique with head-tracking
% for local active sound control" by W. Jung, S. J. Elliott and J. Cheer.
%
% --------------------------------------------------
function [oOpt, Sme, Smm, condNum] = obsFilt(Pe, Pm, Svv, regFacs, snrVal)
% ====================================================
% Check for number of arguments
% ====================================================
narginchk(2, 5);
nargoutchk(0, 4);
% ====================================================
% 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(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
Svv = 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
if nargin > 4 && ~isempty(snrVal)
validateattributes(snrVal, "numeric", {'vector', 'real', 'nonnan', 'nonempty'}, mfilename, "Normalised Root-Mean-Square Singal-to-Noise Ratio of the monitoring microphones", 5);
% Make sure snrVal has correct dimensions
if ~isscalar(snrVal) && numel(snrVal) ~= size(Pm, 1)
error("SNR must be either a scalar or its length must match the number of monitoring microphones.");
end
if sum(snrVal == -Inf) > 0
error("SNR value cannot be equal to -Inf");
end
else
snrVal = Inf;
end
% ====================================================
% Calculate optimal filters
% ====================================================
% Calculate needed quantities
tmpVal = Svv * Pm';
Sme = Pe * tmpVal; % Virtual-Monitor mics cross-spectra Smm = Pm * tmpVal; % Monitor-Monitor mics cross-spectra
% Regularisation matrix
if isscalar(regFacs)
regMat = eye(size(Smm)) .* regFacs;
else
regMat = diag(regFacs);
end
% Signal-to-Noise Ratio at the monitoring microphones
if ~isinf(snrVal)
snrVal = 10^(-snrVal/10); % Convert deciBel to linear
snrMtx = (snrVal.^2) .* norm(Smm, 'fro')/(size(Pm, 1)^2); % Calculate the appropriate SNR amplitude values
else
snrMtx = 0;
end
snrMtx = snrMtx .* eye(size(Pm, 1));
Smm = Smm + snrMtx + regMat;
oOpt = Sme/Smm;
% Condition number of auto-correlation (power spectrum) matrix
if nargout > 3
condNum = cond(Smm);
end
end