function [ys, times] = vfilter(u, kernels, Ns, inds)
%VFILTER computes the output of a conventional Volterra filter for a given
%input signal.
%
%   YS = VFILTER(U, KRS, NS) returns an L x P array whose columns are the
%   outputs of the homogeneous subsystems defined by the P (vectorized) 
%   kernels contained in the cell array KERS. The memories of these kernels
%   are specified by the input vector NS. 
%
%   [YS, TIMES] = VFILTER(U, KRS, NS) returns also a vector TIMES
%   containing the time spent on the computation of each homogeneous
%   subsystem output.
%
%   YS = VFILTER(U, KRS, NS, INDS) allows providing a structure INDS which
%   contains auxiliar indices used in the output computation. This is
%   recomended for greater efficiency. INDS must be pre-computed with the
%   command: 
%      INDS = vfinds(NS)
%
%   See also VFILTERP
%
%   CBVK Toolbox
%   JHMG, PMSB 2016

if nargin < 4
    warning(['Indices for the computation of products of' ...
        ' input samples have not been provided. This slows ' ...
        'down vfilter. To suppress this warning, precompute these'...
        ' indices by using inds = vfinds(Ns).'])
    inds = vfinds(Ns);
end
    
% Model order and memories
P = length(kernels);

% Output signals and measured times
yps = cell(1,P);
times = zeros(1,P);

% Loops over pth-order subsystems
inttimer = tic;
yps{1} = conv(kernels{1},u);
times(1) = toc(inttimer);
for p=2:P
    inttimer = tic;
    yps{p} = vfilterp(u, kernels{p}, Ns(p), inds{p});
    times(p) = toc(inttimer);
end

% Arranges outputs in array with correct size
MM = max(cellfun(@length,yps));
ys = zeros(MM,P);
for p=1:P
    ys(:,p) = [yps{p}(:); zeros(MM-length(yps{p}(:)))];
end