winSincFracDel.m

%% Windowed-sinc fractional delay
% --------------------------------------------------
% Author: Achilles Kappis
% e-mail: axilleaz@protonmail.com
%
% Date: 09/03/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: []].
% 
% --------------------------------------------------
% 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.
% 
% delSig [numeric]: The delayed signal(s). The length 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 a
%                   zeroed-out vector is returned.
% --------------------------------------------------
% 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, delSig] = winSincFracDel(delay, len, win, sig)
    % ====================================================
    % Check for number of arguments
    % ====================================================
    narginchk(1, 4);
    nargoutchk(0, 2);

    % ====================================================
    % 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

    % ====================================================
    % Generate FIR filter
    % ====================================================
    delFilt = -floor((len - 1)/2):floor(len/2); delFilt = delFilt.'; % Filter sample indices
    delFilt = sinc(delFilt - mod(delay, floor(delay))); % Calculate the sinc for the fractional part of the delay
    delFilt = repmat(delFilt, 1, size(sig, 2)) .* win; % Replicate filter and apply window function


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