diff --git a/Signal Processing/Generic/MATLAB/Functions/winSincFracDel.m b/Signal Processing/Generic/MATLAB/Functions/winSincFracDel.m index c2e5f3b29052731b46157ca2bb9a57b3d78f074c..3455345b0c309c557a8faf4eb5ea057f04228f51 100644 --- a/Signal Processing/Generic/MATLAB/Functions/winSincFracDel.m +++ b/Signal Processing/Generic/MATLAB/Functions/winSincFracDel.m @@ -12,8 +12,8 @@ % -------------------------------------------------- % Input % -% delay [numeric]: The delay in samples. This must be a non-negative real -% scalar. +% del [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 @@ -21,18 +21,16 @@ % 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: []]. +% winFun [string/char, numeric] (Optional): Windowing functions to be +% applied to the truncated sinc +% function. The available options +% are "Rectangular", "Hann" and +% "Hamming" and are not +% case-sensitive. Additionally, a +% non-negative scalar can be +% given to be used as the +% parameter for a Kaiser window. +% [Default: "Rectangular"]. % % sig [numeric] (Optional): The signal(s) to be delayed in the time-domain. % This can be a matrix with each column being a @@ -52,74 +50,75 @@ % -------------------------------------------------- % 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. +% sincFilt [numeric]: The fractional delay windowed-sinc FIR filter(s).This +% argument is a vector with the FIR filter. This is the +% causal version of the filter. If the non-causal +% version is required, shift to the left by "causDel". % -% intDel [numeric]: The causal latency of the filter. This value -% corresponds to half the length of the filter. +% causDel [numeric]: The delay required to make the filter causal. 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. +% dSig [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. % % -------------------------------------------------- % 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 +% long signals. If the "dSig" is not required 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) +function [sincFilt, causDel, dSig] = winSincFracDel(del, len, winFun, sig, sigLen) % ==================================================== % Check for number of arguments % ==================================================== narginchk(1, 5); - nargoutchk(0, 5); + nargoutchk(0, 3); % ==================================================== % Validate input arguments % ==================================================== % Validate mandatory arguments - validateattributes(delay, "numeric", {'scalar', 'real', 'nonnegative', 'nonnan', 'nonempty', 'finite'}, mfilename, "Delay in samples", 1); + validateattributes(del, "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 + if len < del warning("winSincFracDel(): The length of the filter cannot be less than the delay, setting equal to 'delay',"); - len = ceil(delay); + len = ceil(del); end else - len = ceil(delay); + len = ceil(del); end - if nargin > 2 && ~isempty(win) - validateattributes(win, "numeric", {'nonnan', '2d', 'real', 'finite', 'nrows', len}, mfilename, "Window(s) to be applied", 3); + if nargin > 2 && ~isempty(winFun) + validateattributes(winFun, {'char', 'string', 'numeric'}, {'nonempty'}, mfilename, "Window to be applied to the truncated sinc function (or the alpha parameter for a Kaiser window)", 3); + + if ischar(winFun) || isStringScalar(winFun) + validatestring(winFun, ["Rectangular", "Hann", "Hamming"], mfilename, "Window to be applied to the truncated sinc function", 3); + winFun = string(winFun); + elseif isscalar(winFun) + validateattributes(winFun, {'numeric'}, {'nonempty', 'nonnan', 'finite', 'nonnegative', 'real', 'scalar'}, mfilename, "Alpha parameter for a Kaiser window to be applied to the truncated sinc function", 3); + end else - win = ones(len, 1); + winFun = "rectangular"; 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)); + if size(winFun, 2) == 1 + winFun = repmat(winFun, 1, size(sig, 2)); else - if size(win, 2) ~= size(sig, 2) + if size(winFun, 2) ~= size(sig, 2) error("Number of windows must match the number of signals provided"); end end @@ -135,46 +134,50 @@ function [delFilt, intDel, causalFilt, delSig, causalSig] = winSincFracDel(delay end % ==================================================== - % Generate FIR filter + % Generate FIR filter for fractional delay % ==================================================== - delFilt = -len + 1:len; delFilt = delFilt.'; % Filter sample indices + idx = (-len + 1:len).'; % Filter sample indices - intDel = -floor(delFilt(1)/2); % The causal delay of the filter - delFilt = sinc(delFilt - delay); % Calculate the filter + causDel = -floor(idx(1)/2); % The causal delay of the filter + sincFilt = sinc(idx - del); % 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 + % Apply windowing % ==================================================== - % Return the shifted filter - if nargout > 2 - causalFilt = delFilt(intDel:end); - end + fracDel = abs(rem(del, 1)); % Total filter delay - if nargout > 3 - for idx = size(sig, 2):-1:1 - delSig(:, idx) = conv(sig, delFilt, "full"); + % Generate window + if isstring(winFun) + switch lower(winFun) + case "hann" + win = sin(pi * (idx - fracDel)/len).^2; + case "hamming" + alpha = 25/46; + win = alpha - (1 - alpha) * cos(2 * pi * (idx - fracDel)/len); + otherwise % Rectangular + win = ones(length(idx), 1); end + else + % Kaiser window + win = besseli(0, pi * winFun * sqrt(1 - (2 * (idx - del)/len).^2))/besseli(0, pi * winFun); end - - % Shift signal to the left by the amount of the "causal" delay of the filter - if nargout > 4 - causalSig = delSig(intDel:end, :); + + % Shift window function to centre on sinc's peak + if ~isscalar(winFun) + win = delayseq(win, floor(del) + causDel); end - % Return the length of the signal asked for - if strcmpi(sigLen, "Same") - if nargout > 3 - delSig = delSig(1:size(sig, 1), :); - end - if nargout > 4 - causalSig = causalSig(1:size(sig, 1), :); + % Apply window function and keep only the causal part up to the length of the filter + sincFilt = sincFilt(causDel + 1:causDel + len) .* win(causDel + 1:causDel + len); + + + % ==================================================== + % Filter the signals + % ==================================================== + if nargout > 2 + for idx = size(sig, 2):-1:1 + dSig(:, idx) = conv(sig, sincFilt, lower(sigLen)); end end end \ No newline at end of file