Skip to content
Snippets Groups Projects
Commit c14444d8 authored by Achilles Kappis's avatar Achilles Kappis
Browse files

Add function to calculate the optimal observation filters in the time domain

parent 0b321aef
No related branches found
No related tags found
1 merge request!3Update to v0.2.3
Hide whitespace changes
Inline Side-by-side
Virtual Sensing/Remote Microphone Technique/MATLAB/Functions/obsFiltTD.m 0 → 100644
+ 247
− 0
View file @ c14444d8
%% Calculate optimal observation filters in the time domain
% --------------------------------------------------
% Author: Achilles Kappis
% e-mail: axilleaz@protonmail.com
%
% Date: 12/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
% must be a real non-negative scalar less than
% or equal to I, the length of the virtual
% microphone signals. [Default: 0].
%
% --------------------------------------------------
% Output
%
% O [numeric]: The observation filters. This is an firLenxNxMxJ 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 averaged over the trials/sound
% field realisations. This is an NxMxfirLen 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, averaged over the trials/sound field
% realisations. The dimensions of this matrix are
% NxTotFiltLen. These filters are denoted "optimal" since
% they are averaged over many realisations (if J > 1).
%
% 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)
% ====================================================
% Check for number of arguments
% ====================================================
narginchk(2, 5);
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'}, mfilename, "Regularisation factor", 3);
else
beta = 0;
end
if nargin > 3 && ~isempty(filtLen)
validateattributes(filtLen, "numeric", {'scalar', 'positive', 'nonnan', '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", {'scalar', 'nonnegative', 'nonnan', 'finite', 'nonempty', '<=', size(e, 1)}, mfilename, "Delay to be added to virtual microphone signals to implement the delayed Remote Microphone Technique", 5);
else
delay = 0;
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(delay + 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
% Mean observation over trials/sound field realisations
if nargout > 7
Omean = mean(O, 4);
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
% Average of optimal observation filters over trials/sound field realisations
if nargout > 10
Oopt = mean(Ovec, 3);
end
% Mean cross-correlation matrix between monitoring and virtual microphones over trials/sound field realisations
if nargout > 11
RmeMtxMean = mean(RmeMtx, 3);
end
% Mean cross-correlation matrix of monitoring microphones over trials/sound field realisations
if nargout > 12
RmmMtxMean = mean(RmmMtx, 3);
end
end
\ No newline at end of file
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment