FFT-based convolution

Discrete convolution using FFT method

Description

As opposed to Matlab CONV, CONV2, and CONVN implemented as straight forward sliding sums, CONVNFFT uses Fourier transform (FT) convolution theorem, i.e. FT of the convolution is equal to the product of the FTs of the input functions.

In 1-D, the complexity is O((na+nb)*log(na+nb)), where na/nb are respectively the lengths of A and B.

Optional arguments to control the dimension(s) along which convolution is carried out.

Slightly less accurate than sliding sum convolution.

Good usage recommendation:
        In 1D, this function is faster than CONV for nA, nB > 1000.
        In 2D, this function is faster than CONV2 for nA, nB > 20.
        In 3D, this function is faster than CONVN for nA, nB > 5.

Acknowledgements

This file inspired Matching Pursuit For 1 D Signals and Conv2fft Reuse.

• MATLAB release: MATLAB 7.8 (R2009a)

[Luke Phai]

Thanks Bruno.

[Bruno Luong]

The GPU acceleration depends on user's hardware. It is impossible to give reliable number without what is your computer setup. A test I run long ago shows an acceleration between 3-5 times.

[Bruno Luong]

Sorry, Luc -> Luke

[Luke Phai]

Excellent work Bruno. Many thanks.
Quick question. How much faster it would be with the GPU jacket enabled?
I try GPUmat using fft2 to programme similar code but it turns out to be slower.
I am thinking to get MATLAB Parallel Computing Toolbox to run the GPU if it is a lot faster.
I hope it is.
Could you please give me some idea? Say
A = rand(1000,1000);
B = rand(1000,1000);
tic;C=convnfft(A, B, 'same', [1, 2],'false');toc
given me> 0.213153 seconds without GPU

tic;C=convnfft(A, B, 'same', [1, 2],'true');toc
time???
Thanks

[Nikolay S.]

Hello again. Apparently it's Matlab filter2 function fault: "Given a matrix X and a two-dimensional FIR filter h, filter2 rotates your filter matrix 180 degrees to create a convolution kernel."
Why the rotation is needed? hack knows...

[Nikolay S.]

Hi Bruno. Excellent contribution and elegant code.
A few small comments:
1) The speed up is smaller then the one you state- as conv is optimized both by both Мatlab and by CPU vendor (Intel in my case).
2) The convolution shift in your (and btw mine) is different form the one resulting from filter2 function. Try running the following code:

img=rgb2gray(imread('peppers.png'));
filt=zeros(50);
filt(1,1)=1;

convFFT=convnfft(img, filt, 'same', [1, 2]);
regConv=filter2(filt, img);

figure;
subplot(1, 2, 1);
imshow(uint8(convFFT));
title('FFT based convolution');
subplot(1, 2, 2);
imshow(uint8(regConv));
title('Matlab regular convolution');

[Jonathan]

 

[Edwin]

Moreover, when I checked the result form this code, there are some different between this and convn for two 3d matrix, the 2 input matrix are both positive.
the result from this code has negative values while convn does not.

[Bruno Luong]

Hi Michael,
May be the inplace times is no longer necessary for recent Matlab.

I remember implement that from a user request.

Thanks for the useful feedback.

[Michael Völker]

Hi Bruno, are you sure your inplaceprod() is (still) useful?

I'm pretty sure, MATLAB does A=A.*B in-place itself. I just compared my memory usage for very large A/B with both methods and there was no difference.
This post is from 2007: http://blogs.mathworks.com/loren/2007/03/22/in-place-operations-on-data/

In a heterogeneous environment, it is useful to avoid mex code for such small tasks. If you are a non-privileged user on a compute server it is really a mess when compiling fails due to compiler version, libraries or whatever.

A minor notice is that (i)fftn is faster than for-loops around 1D (i)fft calls. At least as long as the input and output are of the same size. So I got at least a little speed gain by replacing
for dim=dims
A = ifft(A,[],dim);
end
with
A = ifftn( A );
for the MATLAB ifft case.

Thank you for the code!

[Romesh]

I think this says it all...

>> tic;C = convn(Vs,Vs);toc;
Elapsed time is 473.103412 seconds.
>> tic;C2 = convnfft(Vs,Vs);toc;
Elapsed time is 1.351315 seconds.
>> max(max(max(abs(C-C2))))
ans =
5.2208e-15

Thanks so much for this!

[Arun Shamanna]

 

[Felipe G. Nievinski]

Well written (IMHO).

[Alexander]

Awesome function! My code runs 60x faster now (thanks to your GPU support).

[Eric]