obsFilt.m

%% Calculate optimal observation filters in the frequency domain
% --------------------------------------------------
% Author: Achilles Kappis
% e-mail: axilleaz@protonmail.com
%
% Date: 10/03/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), a scalar
%                               which will result in the same
%                               regularisation factor for all microphones,
%                               or a square matrix with dimensions KxK.
%                               [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, a
%                              vector of dimensions Kx1, containing the SNR
%                              values of each microphone or a square matrix
%                              with dimensions KxK. 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", {'2d', 'nonnan', 'nonempty', 'finite'}, mfilename, "Regularisation factors", 4);

        % Make sure regFacs has correct dimensions
        if ~isscalar(regFacs)
            if isvector(regFacs) && numel(regFacs) ~= size(Pm, 1)
                error("If the regularisation factors are given as a vector the number of elements must match the number of monitoring microphones.");
            elseif ismatrix(regFacs) && ~isequal(size(regFacs), ones(1, 2) .* size(Pm, 1))
                error("If the regularisation factors are given as a matrix, it must be a square matrix with each dimension equal to the number of monitoring microphones.");
            end
        end
    else
        regFacs = 0;
    end

    if nargin > 4 && ~isempty(snrVal)
        validateattributes(snrVal, "numeric", {'2d', 'real', 'nonnan', 'nonempty'}, mfilename, "Normalised Root-Mean-Square Singal-to-Noise Ratio of the monitoring microphones", 5);
        
        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

        % Make sure snrVal has correct dimensions
        if ~isscalar(snrVal)
            if isvector(snrVal) && numel(snrVal) ~= size(Pm, 1)
                error("If the SNR values are given as a vector the number of elements must match the number of monitoring microphones.");
            elseif ismatrix(snrVal) && ~isequal(size(snrVal), ones(1, 2) .* size(Pm, 1))
                error("If the SNR values are given as a matrix, it must be a square matrix with each dimension equal to the number of monitoring microphones.");
            end
        end

        if nnz(snrVal == -Inf)
            error("SNR value cannot be equal -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;
    elseif isvector(regFacs)
        regMat = diag(regFacs);
    else
        regMat = 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

    if isscalar(snrMtx)
        snrMtx = snrMtx .* eye(size(Pm, 1));
    elseif isvector(snrMtx)
        snrMtx = diag(snrMtx);
    end

    Smm = Smm + snrMtx + regMat;
    oOpt = Sme/Smm;

    % Condition number of auto-correlation (power spectrum) matrix
    if nargout > 3
        condNum = cond(Smm);
    end
end