function [ys,times] = vcpd(u,model)
%VCPD Computes the output of a Volterra-CPD model (also known as 
%Volterra-PARAFAC model) for a given input.
%
%   YS = VCPD(U,MODEL,MEMORIES) returns an L x P array whose 
%   columns are the outputs of the homogeneous subsystems defined by the P
%   kernels whose CPD factors are stored in the cell array MODEL. Each
%   entry of MODEL stores a single matrix, which is assumed to represent
%   the single factor of a symmetric CPD of the corresponding symmetric
%   kernel. U is the input signal. MODEL{1} stores the first-order kernel
%   as a vector.
%
%  [YS,TIMES] = VCPD(U,MODEL,MEMORIES) returns also a vector TIMES
%  containing the time spent on the computation of each homogeneous
%  subsystem output.
%
%   See also VFDIAGS
%
%   CBVK Toolbox
%   JHMG, PMSB 2016

% Model order and memories
P = length(model);
N = max(cellfun(@(x)size(x,1),model));

% Input and output length
M = length(u);
L = M+N-1;

% Initialization
ys = zeros(L,P);
times = zeros(P,1);

inttimer = tic;

% 1st order
y1 = conv(u,model{1});
ys(:,1) = [y1; zeros(L-length(y1),1)];
times(1) = toc(inttimer);

% Higher orders
for p=2:P
    inttimer = tic;
    
    yp = conv2(u(:),model{p});
    yp = yp.^p;
    yp = sum(real(yp),2);
    ys(:,p) = [yp; zeros(L-length(yp))];
    
    times(p) = toc(inttimer);
end;    