ptSrcFieldTD.m

%% Signal generated by point sources in the time-domain
% --------------------------------------------------
% Author: Achilles Kappis
% e-mail: axilleaz@protonmail.com
%
% Date: 15/09/2024 (DD/MM/YYYY)
%
% Copyright: MIT
% --------------------------------------------------
% Functionality: Calculate the signals generated by ideal point sources in
%                the time-domain.
% --------------------------------------------------
% Input
% 
% sPos [numeric]: The Cartesian coordinates of the source positions. This
%                   must be an Nx3 matrix where N is the number of sources.
% 
% rPos [numeric]: The Cartesian coordinates of the receiver positions. This
%                 must be an Mx3 matrix where M is the number of receivers.
% 
% fs [numeric]: The sampling frequency. This is required in order to
%               calculate the time-of-flight in samples and must be a
%               positive real scalar.
% 
% Q [numeric] (Optional): The source signals. This can be either an IxNxJ
%                         array, where I is the length of the source
%                         signals in samples and J is the number of
%                         trials/sound field realisations, a vector of
%                         positive real values holding the source strengths
%                         (this is just an amplitude scaling factor), or a
%                         real positive integer denoting the common
%                         strength of all the source signals. [Default: 1].
% 
% sigLen [numeric] (Optional): The length of the source signals. This must
%                              be real positive integer. The generated
%                              signals are zero-mean, (approximately)
%                              uniformly distributed and normalised to
%                              unity. If "Q" is an array, "sigLen" is used
%                              to set the length of the signals. If
%                              "sigLen" is smaller than I, the signals are
%                              truncated and if it is larger the signals
%                              are padded with zeros. Leave empty to have
%                              the signals unchanged.
%                              [Default: 128 or if Q is matrix size(Q, 1)].
% 
% nTrials [numeric] (Optional): This is the number of trials/sound field
%                               realisations that will be generated. If "Q"
%                               is provided, this argument is ignored. It
%                               must be a positive real scalar.
%                               [Default: 1]
% 
% c [numeric] (Optional): The speed of sound. [Default: 343].
% 
% --------------------------------------------------
% Output
% 
% rSig [numeric]: The signals at the receiver positions for all
%                 trials/sournd field realisations. This is an IxMxJ array.
% 
% rSigMean [numeric]: The signals at the receiver positions averaged over
%                     all trials/sound field realisations. This is an IxM
%                     matrix.
% 
% rSigMtx [numeric]: The signals at each receiver position due to each
%                    source for each trial/sound field realisation. This is
%                    an IxMxNxJ array.
% 
% rSigMtxMean [numeric]: The signals at each receiver position due to each
%                        source averaged over the trials/sound field
%                        realisations. This is an IxMxN array.
% 
% Q [numeric]: The source signals. If Q is provided, the same variable is
%              returned here. If Q is generated internally, the signals are
%              returned.
% 
% --------------------------------------------------
% Notes
% 
% Dependencies: - twoPtDist(): To calculate the distances between sources
%                              and receivers.
% 
% - It is the responsibility of the user to pick (or provide) adequately
%   long source signals
% 
% --------------------------------------------------
function [rSig, rSigMean, rSigMtx, rSigMtxMean, Q] = ptSrcFieldTD(sPos, rPos, fs, Q, sigLen, nTrials, c)
    % ====================================================
    % Check for number of arguments
    % ====================================================
    narginchk(3, 7);
    nargoutchk(0, 5);

    % ====================================================
    % Validate input arguments
    % ====================================================
    % Validate mandatory arguments
    validateattributes(sPos, "numeric", {'2d', 'real', 'nonnan', 'nonempty', 'finite', 'ncols', 3}, mfilename, "Cartesian coordinates of the source positions", 1);
    validateattributes(rPos, "numeric", {'2d', 'real', 'nonnan', 'nonempty', 'finite', 'ncols', 3}, mfilename, "Cartesian coordinates of the receiver positions", 2);
    validateattributes(fs, "numeric", {'scalar', 'real', 'nonnan', 'nonempty', 'finite', 'positive'}, mfilename, "The sampling frequency", 3);

    % Validate optional arguments
    if nargin > 3 && ~isempty(Q)
        if isscalar(Q)
            validateattributes(Q, "numeric", {'scalar', 'real', 'nonnan', 'nonempty', 'finite', 'positive', 'integer'}, mfilename, "Length of source signals in samples", 4);
        elseif isvector(Q) && size(Q, 2) == size(sPos, 1)
            validateattributes(Q, "numeric", {'vector', 'real', 'nonnan', 'nonempty', 'finite', 'ncols', size(sPos, 1)}, mfilename, "Source signals", 4);
        else
            validateattributes(Q, "numeric", {'3d', 'real', 'nonnan', 'nonempty', 'finite', 'ncols', size(sPos, 1)}, mfilename, "Source signals", 4);
        end
    else
        Q = 1;
    end

    if nargin > 4 && ~isempty(sigLen)
        validateattributes(sigLen, "numeric", {'scalar', 'nonempty', 'nonnan', 'finite', 'real', 'positive', 'integer'}, mfilename, "The length of the source signals", 5)
    else
        if ~isvector(Q)
            sigLen = size(Q, 1);
        else
            sigLen = 128;
        end
    end


    if ~ismatrix(Q)
        nTrials = size(Q, 3);
    elseif nargin > 5 && ~isempty(nTrials)
        validateattributes(nTrials, "numeric", {'scalar', 'positive', 'nonnan', 'nonempty', 'finite', 'integer'}, mfilename, "Number of trials/sound field realisations", 6); 
    else
        nTrials = 1;
    end

    if nargin > 6 && ~isempty(c)
        validateattributes(c, "numeric", {'scalar', 'nonempty', 'nonnan', 'finite', 'real', 'positive'}, mfilename, "The speed of sound", 7);
    else
        c = 343;
    end
    

    % ====================================================
    % Pre-process data
    % ====================================================
    % Generate source signals
    if isvector(Q)
        tmp = rand(sigLen, size(sPos, 1), nTrials);
        tmp = tmp - mean(tmp);
        tmp = tmp./max(abs(tmp));
        Q = Q(:).' .* tmp;
    elseif sigLen ~= size(Q, 1)
        if sigLen > size(Q, 1)
            Q = cat(1, Q, zeros(sigLen - size(Q, 1), size(Q, 2), size(Q, 3)));
        else
            Q = Q(1:sigLen, :, :);
        end
    end

    % Calculate source-receiver distances
    dist = twoPtDist(sPos, rPos); % Distances
    del = dist/c; % Delays in seconds

    % ====================================================
    % Calculate signals
    % ====================================================
    % Go through the trials/sound field realisations
    for jIdx = nTrials:-1:1
        % Go through the sources (calculate for all receivers in one go)
        for sIdx = size(dist, 1):-1:1
            rSigMtx(:, :, sIdx, jIdx) = (1./dist(sIdx, :)) .* delayseq(Q(:, sIdx, jIdx), del(sIdx, :), fs);
        end
    end

    % Sum source signals at each receiver position
    rSig = reshape(sum(rSigMtx, 3), size(rSigMtx, 1), size(rSigMtx, 2), size(rSigMtx, 4));

    % ====================================================
    % Calculate output arguments
    % ====================================================
    % Mean receiver signal over all trials/sound field realisations
    if nargout > 1
        rSigMean = squeeze(mean(rSig, 3));
    end

    % Mean receiver signal due to each source over all trials/sound field realisations
    if nargout > 3
        rSigMtxMean = squeeze(mean(rSigMtx, 4));
    end
end