Add windowed sinc fractional delay function in Generic Signal Processing

Signal Processing/Generic/MATLAB/Functions/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
\ No newline at end of file