Skip to content
Snippets Groups Projects
Select Git revision
  • aa60543427c9d43493b8810861c746ab8ea8958f
  • master default protected
2 results

README.md

Blame
  • 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