Prev: creating den num for fix point filter
Next: Help to calculate spectrogram ??
From: Bruno Luong on 21 Jun 2010 06:03
"Bruno Luong" <b.luong(a)fogale.findmycountry> wrote in message <hvnd1h$sq4$1(a)fred.mathworks.com>...
> "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);

Note that you can hack the code to perform REPMAT on fft(x), which save some unnecessary computation.

Bruno
From: Bruno Luong on 21 Jun 2010 06:51

> > % 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);

I overlook at a much simpler syntax:
y = convnfft(K,x(:)*x(:)','full')

However, I still recommend to separate the dimensions so that the fft can be performed only once.

Bruno
From: Michalis on 21 Jun 2010 06:57
I also learn a lot of things today thanks to you! Thanks a lot; I'll try to see how they work. I already have the construction with for loops (as I know doing better because I used to write at fortran and it is my first program on Matlab :-) ):
for i = 1:n
if i<kerlength
for k = 1:i
for j = 1:i
y2(i) = y2(i) + k2_true(k,j)*inp(i-k+1)*inp(i-j+1);
end
end
else
for k = 1:kerlength
for j = 1:kerlength
y2(i) = y2(i) + k2_true(k,j)*inp(i-k+1)*inp(i-j+1);
end
end
end
end
but I know that it's computationally big so I'll try your solutions.
Thanks once again
From: Bruno Luong on 21 Jun 2010 07:22
"Bruno Luong" <b.luong(a)fogale.findmycountry> wrote in message
> y = convnfft(K,x(:)*x(:)','full')

Or using stock function (slow for large kernel/data)

y = convn(K,x(:)*x(:)','full')

Bruno
From: Bruno Luong on 21 Jun 2010 07:27
And finally a must: David young's submission that should handle automatically the rank-1 property of x(:)*x(:)'

% http://www.mathworks.com/matlabcentral/fileexchange/22619-fast-2-d-convolution
y = convolve2(K,x(:)*x(:)','full')

Bruno
First  |  Prev  |  Next  |  Last
Pages: 1 2 3
Prev: creating den num for fix point filter
Next: Help to calculate spectrogram ??