diff --git a/Active Noise Control/MATLAB/Functions/tonalAnc.m b/Active Noise Control/MATLAB/Functions/tonalAnc.m new file mode 100644 index 0000000000000000000000000000000000000000..af7b80b23d2e3b6ebb20d2017b68ed2e74f728e8 --- /dev/null +++ b/Active Noise Control/MATLAB/Functions/tonalAnc.m @@ -0,0 +1,328 @@ +%% 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 diff --git a/Active Noise Control/README.md b/Active Noise Control/README.md new file mode 100644 index 0000000000000000000000000000000000000000..7da78cd0bd18f13a76133789f005650934fcbc69 --- /dev/null +++ b/Active Noise Control/README.md @@ -0,0 +1,38 @@ +# 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