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

Placed optimal tonal control into Active Noise Control section and added README [WIP]

parent e51dc9c5
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Active Noise Control/MATLAB/Functions/tonalAnc.m 0 → 100644
+ 328
− 0
View file @ 0da6d9a9
%% Tonal control in the frequency domain
% --------------------------------------------------
% Author: Achilles Kappis
% e-mail: axilleaz@protonmail.com
%
% Date: 22/02/2024 (DD/MM/YYYY)
%
% Copyright: MIT
% --------------------------------------------------
% Functionality: Perform tonal active noise control in the frequency
% domain.
% --------------------------------------------------
% Input
%
% Pe [numeric]: The NxM matrix holding the tranfer functions from the
% primary sources to the error/virtual microphone positions.
% N is the number of error microphones and M the number of
% primary sources.
%
% Pm [numeric]: The NxM matrix holding the transfer functions from the
% primary sources to the monitoring microphones. N is the
% number of monitoring microphones and M the number of
% primary sources, which must match the number of columns
% of Pe.
%
% Px [numeric]: The NxM matrix holding the transfer functions from the
% primary sources to the reference microphone positions. N is
% the number of reference microphones and M the number of
% primary sources, which must match the number of columns of
% Pe.
%
% Ge [numeric]: The NxM matrix holding the transfer functions from the
% secondary sources to the error/virtual microphone
% positions. N is the number of error/virtual microphones
% which must match the number of rows of "Pe" and M the
% number of secondary sources.
%
% Gm [numeric]: The NxM matrix holding the transfer functions from the
% secondary sources to the monitoring microphone positions.
% N is the number of monitoring microphones which must match
% the number of rows of "Pm" and N the number of secondary
% sources, which must match the number of columns of "Ge".
%
% GeEst [numeric] (Optional): The NxM matrix holding the estimated transfer
% functions from the secondary sources to the
% error/virtual microphone positions. N is the
% number of error/virtual microphones which
% must match the number of rows in Pm and M the
% number of secondary sources, which must match
% the number of columns of Ge. [Default: Ge].
%
% GmEst [numeric] (Optional): The NxM matrix holding the estimated transfer
% functions from the secondary sources to the
% monitoring microphone positions. N is the
% number of monitoring microphones which must
% match the number of rows in Pm and M the
% number of secondary sources, which must match
% the number of columns in Ge. [Default: Gm].
%
% v [numeric] (Optional): The primary source strengths or primary source
% cross power spectral density matrix. If the
% argument is a vector, its number of elements
% must match the number of columns in Pe and
% corresponds to the source strengths. If it is a
% matrix it must be square with each dimension
% matching the number of columns of Pe and will be
% treated as the cross spectral density matrix of
% the primary sources. [Default: eye(size(Pe, 2))].
%
% O [numeric] (Optional): The observation filter in case the remote
% microphone technique is to be used. This argument
% must be a matrix of dimensions NxM, with N being
% the number of error/virtual microphones matching
% the number of rows in Pe and M the number of
% monitoring microphones matching the number of
% rows of Pm. [Default: []].
%
% gRegFac [numeric] (Optional): A scalar regularisation factor used for
% the inversion of the plant responses
% matrix (G' * G). [Default: 0].
%
% xRegFac [numeric] (Optional): A scalar regularisation factor used for
% the inversion of the reference signal
% cross spectral density matrix (Sxx).
% [Default: 0].
%
% filterType [char/string] (Optional): Whether the optimal filter, in the
% least-squares sense, will be used or
% the FxLMS implementation. Possible
% values are "Optimal"/"opt" and
% "FxLMS". The values are NOT case
% sensitive. [Default: "Optimal"].
%
% Peval [numeric] (Optional): The transfer function matrix between the
% primary sources and the positions of sound
% field evaluation. The dimensions must be NxM
% where N is the number of evaluation positions
% and M the number of primary sources, which
% must match the columns of "Pe". If left empty
% or not provided, the values in "Pe" will be
% used. [Default: []].
%
% Geval [numeric] (Optional): The transfer function matrix between the
% secondary sources and the positions of sound
% field evaluation (plant responses). The
% dimensions must be NxM where N is the number
% of evaluation positions, which must be equal
% to the positions provided in "Peval", and M
% is the number of secondary sources which must
% be equal to the number of columns of "Ge". If
% left empty or not provided "Ge" will be used.
% [Default: []].
%
% --------------------------------------------------
% Output
%
% W [numeric]: The NxM filter matrix. N is the number of secondary
% sources matching the number of columns in Ge and M the
% number of reference microphones, matching the number of
% rows of Px.
%
% micErr [numeric]: This is a vector with number of elements equal to the
% number of error/virtual microphones (equal to the
% number of rows of Pe), holding the control error at
% each error/virtual microphone.
%
% normSqMicErr [numeric]: The normalised squared errors at the
% error/virtual microphones. This vector holds the
% squares of "micErr" normalised to the power
% spectral density of the true sound field at the
% position of each error/virtual microphone.
%
% evalErr [numeric]: The vector holding the error at the evaluation
% positions. This is empty if "Peval" is not provided.
%
% normSqEvalErr [numeric]: The vector holding the normalised squared error
% at the evaluation positions. This is empty if
% "Peval" is not provided.
%
% gCondNum [numeric]: The inversion condition numbers of the plant
% response matrices to be inverted
% [G' * G + regFac * eye()].
%
% --------------------------------------------------
% Notes
%
% - If O is not provided or given as an empty array, the remote
% microphone technique will not be used. In this case Pe, Ge and GeEst
% are ignored in the control phase and only monitoring microphones are
% assumed. Thus, only Pm, Gm and GmEst are used for control. However,
% Pe and Ge are still used for the estimation at the error/virtual
% microphone positions. This allows the evaluation of control without
% virtual sensing at remote locations.
%
%
% - The implementation is based on:
% 1) Robust performance of virtual sensing methods for active noise
% control by J. Zhang, S. J. Elliott and J. Cheer.
% 2) Modeling local active sound control with remote sensors in
% spatially random pressure fields by S.J. Elliott and J. Cheer.
%
% --------------------------------------------------
function [W, micErr, normSqMicErr, evalErr, normSqEvalErr, gCondNum] = tonalAnc(Pe, Pm, Px, Ge, Gm, GeEst, GmEst, v, O, gRegFac, xRegFac, filterType, Peval, Geval)
% ====================================================
% Check for number of arguments
% ====================================================
narginchk(5, 14);
nargoutchk(0, 6);
% ====================================================
% Validate input arguments
% ====================================================
% Validate mandatory arguments
validateattributes(Pe, "numeric", {'nonempty', 'finite', 'nonnan'}, mfilename, "Tranfer function matrix from primary sources to error/virtual the microphones.", 1);
validateattributes(Pm, "numeric", {'nonempty', 'finite', 'nonnan', 'ncols', size(Pe, 2)}, mfilename, "Tranfer function matrix from primary sources to the monitoring microphones.", 2);
validateattributes(Px, "numeric", {'nonempty', 'finite', 'nonnan', 'ncols', size(Pe, 2)}, mfilename, "Tranfer function matrix from primary sources to the reference microphones.", 3);
validateattributes(Ge, "numeric", {'nonempty', 'finite', 'nonnan', 'nrows', size(Pe, 1)}, mfilename, "Plant response matrix from the secondary sources to the error/virtual microphones", 4);
validateattributes(Gm, "numeric", {'nonempty', 'finite', 'nonnan', 'nrows', size(Pm, 1)}, mfilename, "Plant response matrix from the secondary sources to the monitoring microphones", 5);
% Validate optional arguments
if nargin > 5 && ~isempty(GeEst)
validateattributes(GeEst, "numeric", {'finite', 'nonnan', 'nrows', size(Pe, 1), 'ncols', size(Ge, 2)}, mfilename, "Estimated plant response matrix from the secondary sources to the error/virtual microphones", 6);
else
GeEst = Ge;
end
if nargin > 6 && ~isempty(GmEst)
validateattributes(GmEst, "numeric", {'finite', 'nonnan', 'nrows', size(Pm, 1), 'ncols', size(Ge, 2)}, mfilename, "Estimated plant response matrix from the secondary sources to the monitoring microphones", 7);
else
GmEst = Gm;
end
if nargin > 7 && ~isempty(v)
if isvector(v)
validateattributes(v, "numeric", {'nonnan', 'finite', 'real', 'numel', size(Pe, 2)}, mfilename, "Primary source strengths.", 8);
Svv = diag(v.^2);
else
validateattributes(v, "numeric", {'2d', 'square', 'finite', 'nonnan', 'real', 'ncols', size(Pe, 2)}, mfilename, "Primary source cross spectral density matrix.", 8);
Svv = v;
end
else
Svv = eye(size(Pe, 2));
end
if nargin > 8 && ~isempty(O)
validateattributes(O, "numeric", {'2d', 'finite', 'nonnan', 'nrows', size(Pe, 1), 'ncols', size(Pm, 1)}, mfilename, "Observation filter.", 9);
else
O = [];
end
if nargin > 9 && ~isempty(gRegFac)
validateattributes(gRegFac, "numeric", {'scalar', 'nonnegative', 'finite', 'nonnan', 'real'}, mfilename, "Plant response regularisation factor.", 10);
else
gRegFac = 0;
end
if nargin > 10 && ~isempty(xRegFac)
validateattributes(gRegFac, "numeric", {'scalar', 'nonnegative', 'finite', 'nonnan', 'real'}, mfilename, "Reference signal cross spectral density matrix regularisation factor.", 11);
else
xRegFac = 0;
end
if nargin > 11 && ~isempty(filterType)
validateattributes(filterType, {'char', 'string'}, {'scalartext'}, mfilename, "Filter type", 12);
validatestring(filterType, ["Optimal", "Opt", "FxLMS"], mfilename, "Filter type", 12);
else
filterType = "optimal";
end
if nargin > 12 && ~isempty(Peval)
validateattributes(Peval, {'numeric'}, {'2d', 'ncols', size(Pe, 2)}, mfilename, "Primary path transfer functions to evaluation positions", 13);
else
Peval = [];
end
if nargin > 13 && ~isempty(Geval)
if isempty(Peval)
error("'Peval' cannot be empty if 'Geval' is provided.");
end
validateattributes(Geval, {'numeric'}, {'2d', 'ncols', size(Ge, 2)}, mfilename, "Secondary path transfer functions to evaluation positions (plant responses)", 14);
else
Geval = [];
end
% ====================================================
% Calculated needed quantities
% ====================================================
% Calculate cross spectral densities
Sxm = Pm * Svv * Px';
Sxx = Px * Svv * Px';
% Calculate plant response matrices
if isempty(O)
G = GmEst;
else
G = GeEst + O * (Gm - GmEst);
end
% ====================================================
% Calculate filters
% ====================================================
if isempty(O)
% No remote microphone technique
invG = (G' * G);
if gRegFac ~= 0
invG = invG + gRegFac * eye(size(invG));
end
% Optimal solution and FxLMS are the same
W = -invG\G';
W = W * Sxm/(Sxx + xRegFac * eye(size(Sxx)));
else
% With remote microphone technique
if strcmpi(filterType, "FxLMS")
invG = GeEst' * G;
if gRegFac ~= 0
invG = invG + gRegFac * eye(size(invG));
end
W = -invG\GeEst' * O * Sxm/(Sxx + xRegFac * eye(size(Sxx)));
else
invG = G' * G;
if gRegFac ~= 0
invG = invG + gRegFac * eye(size(invG));
end
W = -invG\G' * O * Sxm/(Sxx + xRegFac * eye(size(Sxx)));
end
end
% ====================================================
% Calculate errors
% ====================================================
% Control error
if nargout > 1
micErr = Pe + Ge * W * Px;
end
% Normalised squared control error
if nargout > 2
normSqMicErr = diag(micErr * Svv * micErr')./diag(Pe * Svv * Pe');
end
% Error at evaluation positions
if nargout > 3 && ~isempty(Peval)
evalErr = Peval + Geval * W * Px;
else
evalErr = [];
end
% Normalised squared error at evaluation positions
if nargout > 4 && ~isempty(evalErr)
normSqEvalErr = diag(evalErr * Svv * evalErr')./diag(Peval * Svv * Peval');
else
normSqEvalErr = [];
end
% Condition number of inverse of (G' * G + regFac * eye())
if nargout > 5
gCondNum = cond(invG);
end
end
\ No newline at end of file
Active Noise Control/README.md 0 → 100644
+ 38
− 0
View file @ 0da6d9a9
# Active Noise Control
The *Remote Microphone Technique* (RMT) was introduced by Roure and Albarrazin \[[1](#references)\] to control a sound field at locations remote from sensors. The sound field is estimated at a remote position based on the measurements at the physical sensors through an observation filter, constituting the transfer function from the physical sensors to the remote location termed *virtual microphone*. An optimal, in the least-squares sense, observation filter is estimated in a preliminary *identification phase* and consequently used during the control phase.
More information on the formulation used in the development of this project can be found in the [*Documentation*]() file.
## Dependencies
The implementations are dependent on functions from the [Utilities](https://gitlab.com/in-nova/in-nova/-/tree/main/Utilities/MATLAB?ref_type=heads). The necessary files are listed below:
- ptsOnSphere.m
- twoPtDist.m
- rotMat3d.m
## Examples
A condensed example of a very simple estimation case is shown below and more examples can be found in the [*Documentation*](https://gitlab.com/in-nova/in-nova/-/tree/main/Documentation/latex?ref_type=heads), along with a complete description of all the functions.
```matlab
% Place primary sources, physical and virtual microphones
mPos = rcvGeo("UCA", 4, 0.5); % Omni mic positions
vPos = srcGeo("Diffuse", 0, 3, 10); % Primary source positions
ePos = virtMicGeo("Grid", 2, 0, 0, 0, 11^2, "xy"); % Virtual/error mic positions on a grid
% Calculate transfer functions
Pm = ptSrcField(vPos, mPos, 7.5e2); % Primary-to-monitoring transfer function (f = 750 Hz)
Pe = ptSrcField(vPos, ePos, 7.5e2); % Primary-to-virtual transfer function (f = 750 Hz)
% Perform estimation and acquire the normalised square estimation error
[~, ~, ~, ~, ~, ~, normSqrErr] = obsFilt(Pe, Pm, [], 1e2); % Use regularisation factor 100
```
#### Acknowledgments
None.
## References
[1] A. Roure, A. Albarrazin, *"The remote microphone technique for active noise control"*, Inter-Noise and Noise-Con Congress and Conference Proceedings, Active99, Fort Lauderdale FL, pp. 1233-1244(12).
\ 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