Select Git revision
obsFiltTD.m
-
Achilles Kappis authoredFixed errors in header documentation, fixed calculations of mean O to be calculated by averaged Rme and Rmm and introduced separate (fractional) delays for each virtual microphone in either samples or seconds
Achilles Kappis authoredFixed errors in header documentation, fixed calculations of mean O to be calculated by averaged Rme and Rmm and introduced separate (fractional) delays for each virtual microphone in either samples or seconds
obsFiltTD.m 13.17 KiB
%% Calculate optimal observation filters in the time domain
% --------------------------------------------------
% Author: Achilles Kappis
% e-mail: axilleaz@protonmail.com
%
% Date: 13/09/2024 (DD/MM/YYYY)
%
% Copyright: MIT
% --------------------------------------------------
% Functionality: Calculate the optimal, in the least squares sense,
% observation filters for the Remote Microphone Technique in
% the time domain.
% --------------------------------------------------
% Input
%
% e [numeric]: The measurement(s) at the virtual microphone position(s).
% This must be an IxNxJ matrix with I the length of the
% measurements in samples, N the number of virtual microphones
% and J the number of trials/sound field realisations.
%
% m [numeric]: The measurement(s) at the monitoring microphone position(s).
% This must be an KxMxJ matrix with K representing the length
% of the measurements in samples, M the number of monitoring
% microphones and J the number of trials/sound field
% realisations. The number of trials must be the same as in
% the argument "e".
%
% beta [numeric] (Optional): The regularisation factor. This must be a real
% non-negative scalar. [Default: 0].
%
% filtLen [numeric] (Optional): The length of the observation filters in
% samples. This must be a positive scalar, at
% less than or equal to K, the number of
% samples in the monitoring microphone
% signals. [Default: size(m, 2)].
%
% delay [numeric] (Optional): The delay to be added to the virtual
% microphone signals, effectively implementing
% the delayed Remote Microphone Technique. This
% can be either a scalar, in which case the
% same delay is applied to all virtual
% microphone signals, or a vector with each
% element holding the delay for each virtual
% microphone signal respectively. The delays
% must be real non-negative values less than
% or equal to I, the length of the virtual
% microphone signals. If "fs" is provided,
% these values correspond to seconds otherwise
% they are in samples and must be integer.
% [Default: 0].
%
% fs [numeric] (Optional): The sampling frequency. If this value is given,
% the delay is treated as being in seconds,
% otherwise it is in samples. This argument must
% be a positive real integer.
%
% --------------------------------------------------
% Output
%
% O [numeric]: The observation filters. This is an NxfirLenxMxJ array.
%
% Rme [numeric]: The cross-correlation matrices between the monitoring and
% virtual microphone signals. This is an filtLenxMxNxJ
% array.
%
% Rmm [numeric]: The cross-correlation matrices of the monitoring
% microphone signals. This is an filtLenxfiltLenxMxMxJ
% array.
%
% Ovec [numeric]: This is a matrix with the observation filters% stacked/vectorised, so that they can be applied to a
% stacked/vectorised monitoring microphone measurements
% vector. The dimensions of this matrix are NxTotFiltLenxJ,
% where TotFiltLen is equal to M times filtLen. This way,
% stacking the monitoring microphone measured signals one
% can perform Ovec(:, :, jIdx) * dm (dm is the vectorised
% monitoring microphone signals vector) to calculate the
% estimated measurements at the N virtual microphone
% positions.
%
% RmeMtx [numeric]: This is the matrix with the cross-correlation vectors
% between the monitoring and virtual microphones
% stacked/vectorised. The dimensions of the matrix are
% NxTotFiltLenxJ.
%
% RmmMtx [numeric]: This is the matrix with the cross-correlation matrices
% of the monitoring microphone signals stacked together.
% The dimensions are TotFiltLenxTotFiltLenxJ.
%
% mMtx [numeric]: The monitoring microphone signals vectorised so that they
% can be used directly with either "Ovec" or "Oopt"
% arguments to perform the filtering/estimation.
%
% Omean [numeric]: The observation filters calculated with the correlation
% matrices RmeMtx and RmmMtx averaged over the
% trials/sound field realisations (these are the
% "RmeMtxMean" and "RmmMtxMean" output arguments). This is
% an NxfirLenxM array.
%
% RmeMean [numeric]: The cross-correlation matrices between the monitoring
% and virtual microphone signals averaged over the
% trials/sound field realisations. This is an
% NxMxfirLen array.
%
% RmmMean [numeric]: The cross-correlation matrices of the monitoring
% microphone signals averaged over the trials/sound
% field realisations. This is an MxMxfirLenxfirLen
% array.
%
% Oopt [numeric]: This is a matrix with the observation filters
% stacked/vectorised, so that they can be applied to a
% stacked/vectorised monitoring microphone measurements
% vector, calculated with the correlation matrices RmeMtx
% and RmmMtx averaged over all trials/sound field
% realisations (the "RmeMtxMean" and "RmmMtxMean" output
% arguments). The dimensions of this matrix are
% NxTotFiltLen.
%
% RmeMtxMean [numeric]: This is the matrix with the cross-correlation
% vectors between the monitoring and virtual
% microphones stacked/vectorised, averaged over the
% trials/sound field realisations. The dimensions of
% the matrix are NxTotFiltLen.
%
% RmmMtxMean [numeric]: This is the matrix with the cross-correlation
% matrices of the monitoring microphone signals
% stacked together and averaged over the trials/sound
% field realisations. The dimensions are
% TotFiltLenxTotFiltLen.
%
% --------------------------------------------------
% Notes
%
% - The calculations are based on the paper: "Modeling local active sound
% control with remote sensors in spatially random pressure fields" by
% Stephen J. Elliott and Jordan Cheer.
%
% --------------------------------------------------
function [O, Rme, Rmm, Ovec, RmeMtx, RmmMtx, mMtx, Omean, RmeMean, RmmMean, Oopt, RmeMtxMean, RmmMtxMean] = obsFiltTD(e, m, beta, filtLen, delay, fs)
% ==================================================== % Check for number of arguments
% ====================================================
narginchk(2, 6);
nargoutchk(0, 13);
% ====================================================
% Validate input arguments
% ====================================================
% Validate mandatory arguments
validateattributes(e, "numeric", {'3d', 'real', 'nonnan', 'nonempty', 'finite'}, mfilename, "Virtual microphone measurements", 1);
validateattributes(m, "numeric", {'3d', 'real', 'nonnan', 'nonempty', 'finite', 'size', [NaN, NaN, size(e, 3)]}, mfilename, "Monitoring microphone measurements", 2);
% Validate optional arguments
if nargin > 2 && ~isempty(beta)
validateattributes(beta, "numeric", {'scalar', 'nonnegative', 'nonnan', 'nonempty', 'finite', 'real'}, mfilename, "Regularisation factor", 3);
else
beta = 0;
end
if nargin > 3 && ~isempty(filtLen)
validateattributes(filtLen, "numeric", {'scalar', 'positive', 'nonnan', 'real', 'nonempty', 'finite', '<=', size(m, 1)}, mfilename, "Length of observation filters", 4);
else
filtLen = size(m, 2);
end
if nargin > 4 && ~isempty(delay)
validateattributes(delay, "numeric", {'vector', 'nonnegative', 'nonnan', 'finite', 'real', 'nonempty', '<=', size(e, 1)}, mfilename, "Delay to be added to virtual microphone signals to implement the delayed Remote Microphone Technique", 5);
if isscalar(delay)
delay = delay * ones(size(e, 2), 1);
elseif length(delay) ~= size(e, 2)
error("The 'delay' argument must be either a scalar, or a vector with number of elements equal to the number of virtual microphone signals.");
end
else
delay = zeros(size(e, 2), 1);
end
if nargin > 5 && ~isempty(fs)
validateattributes(fs, "numeric", {'scalar', 'integer', 'positive', 'nonnan', 'nonempty', 'finite', 'real'}, mfilename, "The sampling frequency", 6);
else
fs = [];
end
% ====================================================
% Pre-process data
% ====================================================
% Delay the virtual microphone signals
if isempty(fs)
for jIdx = size(e, 3):-1:1
e(:, :, jIdx) = delayseq(e(:, :, jIdx), delay);
end
else
for jIdx = size(e, 3):-1:1
e(:, :, jIdx) = delayseq(e(:, :, jIdx), delay, fs);
end
end
% ====================================================
% Calculate cross- and auto-correlation matrices
% ====================================================
% Go through the trials/sound field realisations
for jIdx = size(m, 3):-1:1
% Go through the monitoring microphones
for mIdx = size(m, 2):-1:1
% Calculate the cross-correlations between virtual and monitoring microphones
for eIdx = size(e, 2):-1:1
[corr, lags] = xcorr(m(:, mIdx, jIdx), e(:, eIdx, jIdx));
lIdx = find(lags == 0);
Rme(:, mIdx, eIdx, jIdx) = corr(lIdx:-1:lIdx-filtLen + 1);
end
% Go through the monitoring microphones to calculate the monitoring microphone correlation matrices
for mmIdx = mIdx:-1:1
% Auto-correlation matrices are Toeplitz symmetric
if mIdx == mmIdx
[corr, lags] = xcorr(m(:, mmIdx, jIdx), m(:, mmIdx, jIdx));
lIdx = find(lags == 0);
Rmm(:, :, mIdx, mmIdx, jIdx) = toeplitz(corr(lIdx:-1:lIdx - filtLen + 1));
else
[corr, lags] = xcorr(m(:, mIdx, jIdx), m(:, mmIdx, jIdx));
lIdx = find(lags == 0);
% Cross-correlation matrices
for iIdx = filtLen-1:-1:0
Rmm(:, iIdx + 1, mIdx, mmIdx, jIdx) = corr(iIdx + (lIdx:-1:lIdx - filtLen + 1));
end
Rmm(:, :, mmIdx, mIdx, jIdx) = squeeze(Rmm(:, :, mIdx, mmIdx, jIdx)).';
end
end
end
end
% ====================================================
% Post-process cross- and auto-correlation matrices
% ====================================================
% "Reshape" the data
RmeMtx = reshape(Rme, prod(size(Rme, [1, 2])), size(Rme, 3), size(Rme, 4));
RmmMtx = reshape(permute(Rmm, [1, 3, 2, 4, 5]), prod(size(Rmm, [1, 3])), prod(size(Rmm, [2, 4])), size(Rmm, 5));
% ====================================================
% Calculate observation filters
% ====================================================
for jIdx = size(RmmMtx, 3):-1:1
Ovec(:, :, jIdx) = RmeMtx(:, :, jIdx).'/(RmmMtx(:, :, jIdx) + beta * eye(size(RmmMtx, 1)));
end
% "Split" observation filter vector to observation filters per monitoring and virtual microphone
O = permute(reshape(permute(Ovec, [2, 1, 3]), filtLen, size(m, 2), size(e, 2), size(m, 3)), [3, 1, 2, 4]);
% ====================================================
% Provide additional output arguments
% ====================================================
% The monitoring microphone signals vectorised per trial/sound field realisation
if nargout > 6
mMtx = reshape(m, prod(size(m, [1, 2])), size(m, 3));
end
% Observation filter calculated with the mean Rme and Rmm over the trials/sound field realisation
if nargout > 7
% Average the RmeMtx and RmmMtx matrices
RmeMtxMean = mean(RmeMtx, 3);
RmmMtxMean = mean(RmmMtx, 3);
% Calculate the observation filter
Oopt = RmeMtxMean.'/(RmmMtxMean + beta * eye(size(RmmMtxMean)));
% Reshape
Omean = permute(reshape(Oopt.', filtLen, size(m, 2), size(e, 2)), [3, 1, 2]);
end
% Mean cross-correlations between monitoring and virtual microphones over trials/sound field realisations
if nargout > 8
RmeMean = mean(Rme, 4);
end
% Mean cross-correlations of monitoring microphones over trials/sound field realisations
if nargout > 9 RmmMean = mean(Rmm, 5);
end
end