winSincFracDel.m

%% Windowed-sinc fractional delay
% --------------------------------------------------
% Author: Achilles Kappis
% e-mail: axilleaz@protonmail.com
%
% Date: 17/04/2024 (DD/MM/YYYY)
%
% Copyright: MIT
% --------------------------------------------------
% Functionality: Calculate a fractional delay of a signal with
%                windowed-sinc FIR filter.
% --------------------------------------------------
% Input
% 
% delay [numeric]: The delay in samples. This must be a non-negative real
%                  scalar.
% 
% len [numeric] (Optional): This is the length of the filter. This must be
%                           a positive real integer. If this values is less
%                           than the "delay" value, the filter length
%                           equals the "delay" value rounded up.
%                           [Default: ceil(delay)].
% 
% win [numeric] (Optional): A window to apply to the FIR filter. This can
%                           be a vector in which case the same window will
%                           be applied to the filter for all signals, or it
%                           can be a matrix in which case each column will
%                           contain the window corresponding to the filter
%                           for each signal. The number of rows must be
%                           equal to "len" and the number of columns to the
%                           number of columns of "sig". If left empty, no
%                           window is applied. It has to be noted that the
%                           provided window has to be shifted by "delay"
%                           samples.
%                           [Default: []].
% 
% sig [numeric] (Optional): The signal(s) to be delayed in the time-domain.
%                           This can be a matrix with each column being a
%                           different signal. Only real-valued signals are
%                           accepted. [Default: []].
% 
% sigLen [string/char] (Optional): The length of the resulting signal. This
%                                  argument can be either "Full" or "Same"
%                                  (not case-sensitive). If "Same", the
%                                  delayed signals will have the same
%                                  length as the original. If "Full" is
%                                  provided, the full result of the
%                                  convolution is provided. If the number
%                                  of output arguments is less than 3, this
%                                  argument is ignored. [Default: "Full"].
% 
% --------------------------------------------------
% Output
% 
% delFilt [numeric]: The fractional delay windowed-sinc FIR filter(s). If
%                    multiple windows are provided then a matrix is
%                    returned with each column corresponding to the filter
%                    for each signal. Otherwise, this argument is a vector
%                    with the FIR filter.
% 
% intDel [numeric]: The causal latency of the filter. This value
%                   corresponds to half the length of the filter.
% 
% causalFilt [numeric]: The filter after left shifting "intDel" samples to
%                       the left.
% 
% delSig [numeric]: The delayed signal(s). The length of each signal is
%                   equal to that of the original. If the delay is
%                   longer than the original signal vector, a zeroed vector
%                   is returned (with possible artefacts from the filtering
%                   process). Only returned if a signal is provided,
%                   otherwise an empty array is returned.
% 
% causalSig [numeric]: The delayed signal(s) shifted to the left by the
%                      causal delay of the filter. 
% 
% --------------------------------------------------
% Notes
% 
% - The function uses time-domain convolution, which can be quite slow for
%   long signals. If the "delSig" is not asked for this calculation is
%   omitted so, when the signals are not needed (or want to perform
%   frequency-domain convolution manually), use only one output argument.
% --------------------------------------------------
function [delFilt, intDel, causalFilt, delSig, causalSig] = winSincFracDel(delay, len, win, sig, sigLen)
    % ====================================================
    % Check for number of arguments
    % ====================================================
    narginchk(1, 5);
    nargoutchk(0, 4);

    % ====================================================
    % Validate input arguments
    % ====================================================
    % Validate mandatory arguments
    validateattributes(delay, "numeric", {'scalar', 'real', 'nonnegative', 'nonnan', 'nonempty', 'finite'}, mfilename, "Delay in samples", 1);

    % Validate optional arguments
    if nargin > 1 && ~isempty(len)
        validateattributes(len, "numeric", {'scalar', 'integer', 'real', 'finite', 'nonnan', 'nonnegative'}, mfilename, "Filter length", 2);

        if len < delay
            warning("winSincFracDel(): The length of the filter cannot be less than the delay, setting equal to 'delay',");
            len = ceil(delay);
        end
    else
        len = ceil(delay);
    end

    if nargin > 2 && ~isempty(win)
        validateattributes(win, "numeric", {'nonnan', '2d', 'real', 'finite', 'nrows', len}, mfilename, "Window(s) to be applied", 3);
    else
        win = ones(len, 1);
    end

    if nargin > 3 && ~isempty(sig)
        validateattributes(sig, "numeric", {'2d', 'nonempty', 'real'}, mfilename, "Signal(s) to be delayed", 4);

        if size(win, 2) == 1
            win = repmat(win, 1, size(sig, 2));
        else
            if size(win, 2) ~= size(sig, 2)
                error("Number of windows must match the number of signals provided");
            end
        end
    else
        sig = zeros(len, 1);
    end

    if nargin > 4 && ~isempty(sigLen)
        validateattributes(sigLen, {'string', 'char'}, {'scalartext', 'nonempty'}, mfilename, "Length of delayed signal vectors/matrix", 5)
        validatestring(sigLen, ["Full", "Same"], mfilename, "Length of delayed signal vectors/matrix", 5);
    else
        sigLen = "Full";
    end

    % ====================================================
    % Generate FIR filter
    % ====================================================
    delFilt = -len + 1:len; delFilt = delFilt.'; % Filter sample indices
    
    intDel = -floor(delFilt(1)/2); % The causal delay of the filter
    delFilt = sinc(delFilt - delay); % Calculate the filter

    % Keep only the central part of the filter to match length asked for
    delFilt = delFilt(intDel + 1:intDel + len);

    % Replicate filter and apply window function
    delFilt = repmat(delFilt, 1, size(sig, 2)) .* win;
    

    % ====================================================
    % Filter the signals
    % ====================================================
    if nargout > 2
        for idx = size(sig, 2):-1:1
            delSig(:, idx) = conv(sig, delFilt, "full");
        end
    else
        delSig = [];
    end

    % Return the shifted filter
    if nargout > 3
        causalFilt = delFilt(intDel:end);
    end

    % Shift signal to the left by the amount of the "causal" delay of the filter
    if nargout > 3
        causalSig = delSig(intDel:end, :);
    end

    % Return the length of the signal asked for
    if strcmpi(sigLen, "Same")
        delSig = delSig(1:size(sig, 1), :);
        if nargout > 3
            causalSig = causalSig(1:size(sig, 1), :);
        end
    end
end