function y  = vfilterp(u, h, N, inds)
% Computes the output of a pth-order homogeneous Volterra filter 
% for a given input signal.
%
%   Y = VFILTERP(U, H, N, INDS) returns the output produced by the 
%   pth-order homogeneous Volterra filter having (vectorized and 
%   nonredundant) kernel H with memory N. The structure INDS stores
%   auxiliar indices used in the computation. It is given by the pth cell
%   array component of vfinds(Ns) when Ns(p) = N.
%
%   This function is called by vfilter in order to compute the partial
%   outputs of each Volterra filter subsystem.
%
%   See also VFILTER
%
%   CBVK Toolbox
%   JHMG, PMSB 2016

h = h(:);

L = length(h);

MM = length(u) + N - 1;
y = zeros(MM,1);

u = [u; zeros(MM-length(u),1)];

% Temporary vectors
xN = zeros(N,1);
xp = zeros(L,1);
for m=1:MM
    
    % Extracts column of input signal
    xN(1:min([m N])) = u(m:-1:max([1 m-N+1]));
    
    % Computes corresponding column of Khatri-Rao product
    xp(inds.inon1) = xp(inds.inonN);
    xp(inds.i1) = prod(xN(inds.Su(inds.i1,:)),2);
    
    % Adds to result
    y(m) = h'*xp;
end