Select Git revision
rectRoomPtSrc.m
-
Achilles Kappis authoredAchilles Kappis authored
rectRoomPtSrc.m 15.21 KiB
%% Calculate wavefield generated by point sourcesinside a rectangular room
% --------------------------------------------------
% Author: Achilles Kappis
% e-mail: axilleaz@protonmail.com
%
% Date: 29/05/2025 (DD/MM/YYYY)
%
% Copyright: MIT
% --------------------------------------------------
% Functionality: Calculate the wavefield generated by point sources at
% specified locations inside a rectangular room of specified
% dimensions and absorption.
% --------------------------------------------------
% Input
%
% sPos [numeric]: The position of the sources in Cartesian coordinates. The
% variable should be a matrix with dimensions 3xN where N
% is the number of sources and the rows represent the x, y
% and z coordinates respectively.
% IMPORTANT: The coordinates of the sources must NOT lie
% outside the room (set in "roomSize" argument).
%
% rPos [numeric]: The position of the receivers in Cartesian coordinates.
% The variable should be a matrix with dimensions 3xM where
% M is the number of receivers and the rows represent the
% x, y and z coordinates respectively.
% IMPORTANT: The coordinates of the receivers must NOT lie
% outside the room (set in "roomSize" argument).
%
% roomDims [numeric]: The dimensions of the room. This can either be a
% scalar, in which case the room will correspond to a
% cubic enclosure, or a vector with three elements
% corresponding to the x, y and z dimensions of the
% room. The values must be positive scalars.
%
% f [numeric]: The frequencies for which to calculate the wave field. It
% must be a positive scalar or vector of positive values.
%
% nModes [numeric] (Optional): The (approximate) number of modes to be used
% for the calculation of the sound field. This
% value can either be a positive scalar
% denoting the total number of modes that will
% be used or a 3xN matrix where the first
% dimension will contain the mode index along
% each dimension and N corresponds to the
% number of total modes. See notes for the
% mode indices used when a scalar is provided.
% [Default: 16^3 = 4096].
%
% dampRatio [numeric] (Optional): The damping ratio of the boundaries. This
% can be a scalar denoting that the damping
% ratio is identical for all frequencies or
% a vector with the number of elements
% matching the length of the "f" argument.
% The damping ratio at the natural
% frequencies of the room will be
% calculated through linear interpolation
% from the damping ratio values provided.
% [Default: 0.1].
%
% Q [numeric] (Optional): The source strengths. This argument can be a
% scalar indicating that all sources have the same
% strength per frequency, a vector with N elements,
% each corresponding to the source of each source
% per frequency, or a matrix with dimensions NxF,
% where each element will correspond to the
% strength of each source for each frequency.
% [Default: 1].
%
% c [numeric] (Optional): The speed of sound in m/s. [Default: 343].%
% rho [numeric (Optional): The density of the medium in kg/(m^3).
% [Default: 1.204].
%
% --------------------------------------------------
% Output
%
% modeIdx [numeric]: The mode of the indices in the order they are placed
% in the output arguments below (this helps to identify
% contribution of each mode). This is a 3xK matrix with
% the first dimension denoting the indices of each mode
% and K denoting the number of modes.
%
% waveField [numeric]: The complex pressures at the position of the
% receivers. The variable is an array of dimensions
% MxNxFxK where M is the number of receivers, N the
% number of sources and F the number of frequencies.
%
% srcWaveField [numeric]: The complex pressures at the positions of the
% receivers with the contribution of the individual
% modes being summed up. This should correspond to
% the contribution of each source to each receiving
% position per frequency. This is an MxNxF array.
%
% broadWaveField [numeric]: The complex pressures at the positions of the
% receivers per mode with the contribution per
% frequency being summed up. This corresponds to
% the broadband (in the frequencies provided)
% contribution of each mode at each position per
% source. This is an MxNxK array.
%
% srcBroadWaveField [numeric]: The complex pressures at the positions of
% the receivers with the contribution of each
% mode and frequency summed up. This is an MxN
% array and corresponds to the "total" (from a
% frequency perspective) transfer function
% from each source to each receiver
%
% sumWaveField [numeric]: The complex pressures at the positions of the
% receivers with the contribution of the sources
% being summed up. This corresponds to the sum of
% the contributions of each source to each mode per
% frequency and is an MxFxK array.
%
% sumSrcWaveField [numeric]: The complex pressures at the positions of the
% receivers with the contributions of the
% sources and the modes being summed up. This
% corresponds to the total complex pressure at
% each receiver position per frequency due to
% all the sources and modes and is an MxF array.
%
% sumBroadWaveField [numeric]: The complex pressures at the positions of
% the receivers with the contributions of the
% sources summed up for all frequencies. This
% correspond to the broadband (in the
% frequencies provided) contribution of each
% mode due to all sources and is an MxK array.
%
% totalWaveField [numeric]: This is the wavefield at each receiver position
% with the contribution of each mode, source and
% frequency summed up. In a sense this
% corresponds to the total complex pressure at
% each receiver position and is a vector with M
% elements.
%
% --------------------------------------------------
% Notes
%
% The actual number of modes calculated when a scalar is provided for
% "nModes" is equal to ceil(nthroot(nModes, 3))^3 - 1. This value will% result in a set of mode indices that will provide all available
% combinations of mode indices up to ceil(nthroot(nModes, 3)). As an
% example, if nModes = 5 we would get ceil(nthroot(5, 3))^3 - 1 == 7 and
% the modes used would be
%
% (1, 0, 0) | (0, 1, 0) | (0, 0, 1)
% (1, 1, 0) | (1, 0, 1) | (0, 1, 1)
% (1, 1, 1)
%
%
%
% Dependencies: - twoPtDist(): To calculate the distances between sources
% and receivers.
% --------------------------------------------------
function [modeIdx, waveField, srcWaveField, broadWaveField, srcBroadWaveField, sumWaveField, sumSrcWaveField, sumBroadWaveField, totalWaveField] = rectRoomPtSrc(sPos, rPos, roomDims, f, nModes, dampRatio, Q, c, rho)
% ====================================================
% Check for number of arguments
% ====================================================
narginchk(4, 9);
nargoutchk(0, 9);
% ====================================================
% Validate input arguments
% ====================================================
% Validate mandatory arguments
validateattributes(sPos, "numeric", {'2d', 'nrows', 3, 'finite', 'nonnan', 'nonempty', 'real'}, mfilename, "Source positions", 1);
validateattributes(rPos, "numeric", {'2d', 'nrows', 3, 'finite', 'nonnan', 'nonempty', 'real'}, mfilename, "Receiver positions", 2);
validateattributes(roomDims, "numeric", {'vector', 'finite', 'nonnan', 'nonempty', 'real', 'positive'}, mfilename, "Room dimensions", 3);
if ~isscalar(roomDims) && numel(roomDims) ~= 3
error("The dimensions argument must either be a scalar or a vector with three positive elements.");
end
validateattributes(f, "numeric", {'vector', 'finite', 'nonnan', 'nonempty', 'real', 'positive'}, mfilename, "Frequencies", 4);
% Validate optional arguments (and set their values)
if nargin > 8 && ~isempty(rho)
validateattributes(rho, "numeric", {'scalar', 'finite', 'positive', 'nonnan', 'nonempty', 'real'}, mfilename, "Medium density", 9);
else
rho = 1.204;
end
if nargin > 7 && ~isempty(c)
validateattributes(c, "numeric", {'scalar', 'finite', 'positive', 'nonnan', 'nonempty', 'real'}, mfilename, "Speed of propagation", 8);
else
c = 343;
end
if nargin > 6 && ~isempty(Q)
if isscalar(Q)
validateattributes(Q, "numeric", {'scalar', 'nonempty', 'real', 'finite', 'nonnan', 'nonnegative'}, mfilename, "Source strengths", 7);
elseif isvector(Q)
validateattributes(Q, "numeric", {'vector', 'numel', size(sPos, 2), 'nonempty', 'real', 'finite', 'nonnan', 'nonnegative'}, mfilename, "Source strengths", 7);
else
validateattributes(Q, "numeric", {'2d', 'size', [size(sPos, 2), length(f)] ,'nonempty', 'real', 'finite', 'nonnan', 'nonnegative'}, mfilename, "Source strengths", 7);
end
else
Q = 1;
end
if nargin > 5 && ~isempty(dampRatio)
if isscalar(dampRatio)
validateattributes(dampRatio, "numeric", {'scalar', 'finite', 'nonempty', '<=', 1, 'real', 'nonnan', 'positive'}, mfilename, "Damping ratio", 6);
else
validateattributes(dampRatio, "numeric", {'vector', 'numel', length(f), 'finite', 'nonempty', '<=', 1, 'real', 'nonnan', 'positive'}, mfilename, "Damping ratio", 6);
end
else
dampRatio = 0.1;
end
if nargin > 4 && ~isempty(nModes)
if ~isscalar(nModes)
validateattributes(nModes, "numeric", {'2d', 'nrows', 3, 'nonempty', 'integer', 'real', 'finite', 'nonnan', 'nonnegative'}, mfilename, "Modes", 5);
else
validateattributes(nModes, "numeric", {'scalar', 'nonempty', 'integer', 'real', 'finite', 'nonnan', 'nonnegative'}, mfilename, "Number of modes", 5);
end
else
nModes = 4096;
end
% Check validity of source and receiver positions
if sum(~isbetween(sPos(1, :), -roomDims(1)/2, roomDims(1)/2)) > 0 || ...
sum(~isbetween(sPos(2, :), -roomDims(2)/2, roomDims(2)/2)) > 0 || ...
sum(~isbetween(sPos(3, :), -roomDims(3)/2, roomDims(3)/2)) > 0
error("Sources outside the room where found.");
end
if sum(~isbetween(rPos(1, :), -roomDims(1)/2, roomDims(1)/2)) > 0 || ...
sum(~isbetween(rPos(2, :), -roomDims(2)/2, roomDims(2)/2)) > 0 || ...
sum(~isbetween(rPos(3, :), -roomDims(3)/2, roomDims(3)/2)) > 0
error("Receivers outside the room where found.");
end
% ====================================================
% Pre-process
% ====================================================
% Convert frequencies to circular frequencies
f = 2 * pi * f;
% Create room dimensions vector
if isscalar(roomDims)
roomDims = roomDims * ones(3, 1);
end
% Create matrix with mode indices
if isscalar(nModes)
modeIdx = ceil(nthroot(nModes, 3)) - 1;
modeIdx = combinations(0:modeIdx, 0:modeIdx, 0:modeIdx);
modeIdx = table2array(modeIdx);
modeIdx = modeIdx(2:end, :);
else
modeIdx = nModes;
end
modeIdx = modeIdx./roomDims; % Calculate the ratio of the index to the dimensions of the room (used in the calculations)
% Calculate normalisation factors
normFact = sqrt(prod((modeIdx > 0) + 1, 2));
% Create vector with damping ratios
if isscalar(dampRatio)
dampRatio = dampRatio * ones(length(f), 1);
end
% Create matrix with source strengths
if isscalar(Q)
Q = Q * ones(size(sPos, 2), length(f));
elseif isvector(Q)
Q = repmat(Q(:), 1, length(f));
end
% ====================================================
% Calculate sound field
% ====================================================
% Calculate eigenfrequencies
wEig = pi * c * sqrt(sum(modeIdx.^2, 2));
% Inter/extrapolate damping ratio at the eigenfrequencies
if length(f) > 1 dampRatio = interp1(f, dampRatio, wEig, "linear", "extrap");
end
% Calculate the complex resonance terms
A = repmat(f(:).', length(wEig), 1)./ ...
(2 * dampRatio .* wEig .* f(:).' - ...
1j * (repmat(wEig.^2, 1, length(f)) - repmat((f(:).').^2, length(wEig), 1)));
A = A * (rho * c^2)/prod(roomDims); % Multiply with the constant factor
% Calculate the eigenfunction values at the source and receiver positions
modeIdx = modeIdx * pi;
% Loop through frequencies
for fIdx = 1:length(f)
tmpA = A(:, fIdx);
% Loop through the receivers
for rIdx = size(rPos, 2):-1:1
rPsi = normFact .* prod(cos(modeIdx .* rPos(:, rIdx).'), 2);
% Loop through the sources
parfor sIdx = 1:size(sPos, 2)
waveField(rIdx, sIdx, fIdx, :) = rPsi .* tmpA .* normFact .* prod(cos(modeIdx .* sPos(:, sIdx).'), 2) * Q(sIdx, fIdx);
end
end
end
% ====================================================
% Return additional output arguments
% ====================================================
if nargout > 2
srcWaveField = squeeze(sum(waveField, 4));
end
if nargout > 3
broadWaveField = squeeze(sum(waveField, 3));
end
if nargout > 4
srcBroadWaveField = squeeze(sum(broadWaveField, 3));
end
if nargout > 5
sumWaveField = squeeze(sum(waveField, 2));
end
if nargout > 6
sumSrcWaveField = squeeze(sum(sumWaveField, 3));
end
if nargout > 7
sumBroadWaveField = squeeze(sum(sumWaveField, 2));
end
if nargout > 8
totalWaveField = squeeze(sum(sumBroadWaveField, 2));
end
end