Continuous convolution and the (inverse) FFT

I'm a man deeply rooted in the continuous Fourier transform, and am generally altogether befuddled by the conventions of the discrete Fourier transform. Now, unfortunately, I'm faced with the following problem:

Calculate the continuous convolution of a function f(omega') and g(omega') at omega. Numerically, I have sampled these functions at negative&positive equidistantly spaced values omegan (from -Omega to Omega in N points) with associated values fn and gn. I want the convolution evaluated at omegan as well. Enter the FFT: from the convolution theorem, I know that I can jump back and forth, and the FFT should allow me to do so with a favorable computational scaling.

Unfortunately, I cannot make it work - perhaps, this is an issue of not understanding the (non)dimensionalized form assumed in the DFT, perhaps it is related to zeropadding and cyclic-vs-linear DFTs, or perhaps fftshift needs to be involved: regardless, my attempts do not resolve the issue.

Therefore: for the above problem, what must I do to find the convolution at omegan?

My naive attempt does not do the job (ignoring scaling factors):

fgn_conv = fft(ifft(fn).*ifft(gn))

Any assistance is much appreciated! Thanks!

[David Goodmanson]

Hello Thomas, I don't quite understand what you mean by getting the convolution 'evaluated at omegan' since the convolution should give you a vector of N values, with an independent variable that includes zero and has spacing omegan.

Do the functions fn and gn die out at each end of the array, or do they have significant nonzero values at either end?

[Thomas]

Hi David,

My functions fn and gn die out very rapidly (they are sums of delta-like functions) and "zeros" outside the sampled range omegan. In other words, I want to evaluate int{ f(omega')g(omega'-omega) domega'} (integrated from -inf to inf) at omega = omegan and I know that f(omega) and g(omega) are nonzero only at the value of omegan at zero everywhere else.

Does that make my issue clearer?

[David Goodmanson]

Hi Thomas, I confess I am not getting it yet. You refer to omegan as a sampled range, as if it would contain a set of points. But you also refer to 'the value of omegan' as if it is one quantity, and the original question says, 'sampled these functions at negative&positive equidistantly spaced values omegan (from -Omega to Omega in N points)' as if omegan is the spacing between points. Do your arrays fn and gn contain nonzero values only at isolated points?

[Thomas]

Sorry for the confusion: omegan is a set of points, i.e. an array of size N, that spans a range -|Omega| to Omega within which the functions f and g are nonzero (and outside which they vanish). The continuous functions f and g are sampled at omegan (an array) with associated value fn and gn.

Answer 1

[David Goodmanson]

Hi Thomas, I started out with the continuous transform as well, so going to the discrete form took some getting used to. The basic approach is simply to turn the sum into a Riemann integral by multiplying by the spacing of points of the independent variable, denoted below by delx.

The easy way to do convolution is with 'conv', the first method below. conv requires its own x array xconv since its output is approximately twice as long as the two inputs.

As you mentioned you can also do the fft-ifft approach but there is an extra multiplicative factor N involved, and a bunch of shifting (assuming you started with x = 0 at the center of the array). The fft does circular convolution. If you want it to look like regular convolution you just have to make sure that your function takes up a small enough portion of the x array that the translated function can't wrap around to the beginning with any significant overlap.

Here is some example code. The bookkeeping is easier when N, the number of array points, is odd rather than even. The result by the two conv methods (red and yellow lines) can't be distinguished on the plot

delx = .01;                        % x array spacing
x =     (-2000:2000)*delx;              
xconv = (-4000:4000)*delx;              
y = exp(-x+2).*(1+sign(x-2))/2;    % example: decaying exponential for x>=2 

% regular convolution
yyconv = conv(y,y)*delx;           % include delx factor for integration

sum(y)*delx                        % check, integrates to 1
sum(yyconv)*delx                   % check, also also integrates to 1

% convolution by fft
yshif = ifftshift(y);
N = length(A);
A = N*fft(ifft(yshif).*ifft(yshif));
% next take real part to avoid to avoid a nuisance tiny imaginary part; shift and scale the result
yyconv2 = real(fftshift(A))*delx;

% plot results
plot(x,y,xconv,yyconv,x,yyconv2)

[Thomas]

Dear David,

Thank you so much for the help and for your efforts; this has essentially solved my problem and simultaneously cleared up a lot of confusion for me. I really appreciate that you included a complete example - much appreciated! Thank you!

[Guillaume]

I've not understood the whole conversation (I'm not a mathematician, so Riemann integrals fly over my head), but this page may be relevant: linear and circular convolution

[David Goodmanson]

You're welcome, Thomas. In the answer I see that N = length(A) should be N = length(y) since A is not yet defined.

The informative link that Guillaume provided was a reminder that convolution by fft can also be done as

A =  ifft(fft(yshif).*fft(yshif));

(yshif same as before) in which case you don't need the factor of N.