function y = convsi(hd, u, spiopts)
%CONVSI Encapsulates the convolution routine 
% with the possibility of employing a sparse-interpolated diagonal hd.
%
% Y = CONVSI(HD, U, []) computes the convolution between hd and u.
% 
% Y = CONVSI(HD, U, SPIOPTS) computes the convolution between h and u,
% where h(1+(m-1)*L) = hd(m), with L being defined by spiopts.L
%
%  Input arguments:
%
%  hd : kernel diagonal
%  u  : input signal
%  spiopts : structure storing sparse interpolation options:
%  spiopts.L : interpolation factor
%
%  Outputs:
%   y : result of convolution taking into account sparse interpolation
%   options
%
%   CBVK Toolbox
%   JHMG, PMSB 2016

% If no options are declared, assumes no sparse interpolation is applied
if isempty(spiopts)
    y = conv(u,hd);
else
    
    y = zeros(length(u)+spiopts.L*length(hd)-1,1);
    lr = rem(length(y),spiopts.L);
    
    for l=1:spiopts.L
        if l > lr
            y(l:spiopts.L:end) = conv(hd,u(l:spiopts.L:end));
        else
            y(l:spiopts.L:end) = [conv(hd,u(l:spiopts.L:end)); 0];
        end
    end
    
end