From: Michalis on
Hi, I'm trying to create an output signal using fixed volterra kernels and I would like to know if there is a faster (computationaly) way to do the following equation rather than using a double for loop:
y(n) = Σk1(i)x(n-i) + Σk(i,j)x(n-i)x(n-j)
the first part can be made with a single convolution but for the 2nd part I couldn't find a satisfying command that could do it and return me a vector instead of a matrix.
Thanks
From: Matt J on
"Michalis " <micpan84(a)gmail.com> wrote in message <hvlsu5$88r$1(a)fred.mathworks.com>...
> Hi, I'm trying to create an output signal using fixed volterra kernels and I would like to know if there is a faster (computationaly) way to do the following equation rather than using a double for loop:
> y(n) = &#931;k1(i)x(n-i) + &#931;k(i,j)x(n-i)x(n-j)
> the first part can be made with a single convolution but for the 2nd part I couldn't find a satisfying command that could do it and return me a vector instead of a matrix.
==================

Convolution with x can be represent by multiplication with a Toeplitz matrix X. You can then represent the above in matrix-vector form by

y=X*k1 + sum( (X*k).*X , 2)

The expressions X*k1 and X*k can be evaluated as single convolutions using conv().
The operation sum(z.*X,2) will require you to generate the Toeplitz matrix X.
You can use the toeplitz() function to do this. Or, if X is sparse, it is advisable to generate it using

http://www.mathworks.com/matlabcentral/fileexchange/26292-regular-control-point-interpolation-matrix-with-boundary-conditions
From: Bruno Luong on
"Michalis " <micpan84(a)gmail.com> wrote in message <hvlsu5$88r$1(a)fred.mathworks.com>...
> Hi, I'm trying to create an output signal using fixed volterra kernels and I would like to know if there is a faster (computationaly) way to do the following equation rather than using a double for loop:
> y(n) = &#931;k1(i)x(n-i) + &#931;k(i,j)x(n-i)x(n-j)

I wonder if what you have is Volterra integral equation, the equation I know is linear wrt x:
http://en.wikipedia.org/wiki/Volterra_integral_equation

Are you sure the equation is correctly written?

Bruno
From: Michalis on
"Bruno Luong" <b.luong(a)fogale.findmycountry> wrote in message <hvmucp$ga2$1(a)fred.mathworks.com>...
> "Michalis " <micpan84(a)gmail.com> wrote in message <hvlsu5$88r$1(a)fred.mathworks.com>...
> > Hi, I'm trying to create an output signal using fixed volterra kernels and I would like to know if there is a faster (computationaly) way to do the following equation rather than using a double for loop:
> > y(n) = &#931;k1(i)x(n-i) + &#931;k(i,j)x(n-i)x(n-j)
>
> I wonder if what you have is Volterra integral equation, the equation I know is linear wrt x:
> http://en.wikipedia.org/wiki/Volterra_integral_equation
>
> Are you sure the equation is correctly written?
>
> Bruno
There you can see only the zero'th and first order (or first and second if you prefer). I want to add an additional order. For the infinitive equation see below:
http://en.wikipedia.org/wiki/Volterra_series
From: Bruno Luong on
"Michalis " <micpan84(a)gmail.com> wrote in message <hvnavt$la8$1(a)fred.mathworks.com>...

> There you can see only the zero'th and first order (or first and second if you prefer). I want to add an additional order. For the infinitive equation see below:
> http://en.wikipedia.org/wiki/Volterra_series

Thanks, I learn new thing today, which is good.

Back to the question:

Here is two ways to carried out the sum:

n = 5;
K=rand(n);
x=rand(1,n);

% Engine 1
y1 = zeros(2*n-1,n);
for j=1:n
y1(:,j) = convn(K(:,j),x(:),'full');
end
y = zeros(2*n-1);
for i=1:2*n-1
y(i,:) = convn(x(:).',y1(i,:),'full');
end
y

% Engine 2, no loop
% Need FEX submission
% http://www.mathworks.com/matlabcentral/fileexchange/24504-fft-based-convolution
y1 = convnfft(K,repmat(x(:),1,n),'full',1);
y = convnfft(repmat(x(:).',2*n-1,1),y1,'full',2);
y

% Bruno