Convoluton of two functions $f_1$ and $f_2$, with $f_1=f_2= e^{-\pi x^2}$ using MatLab does not match WolframAlpha result?

Question

I have performed a convolution using the $\texttt{fft}$ function in MatLab. However the result I get does not match the same convolution if I perform it using the Fourier transform in WolframAlpha. The functions $f_1$ and $f_2$ that I am convolving are (for simplicity) given by $$ f_1 = f_2 = e^{-\pi x^2}. $$

With WolframAlpha I use the convolution theorem so I have $$ F[f_1\ast f_2] = F[f_1]\cdot F[f_2] = (F[f_1]^2), $$ and therefore $$ f_1\ast f_2 = F^{-1}(F[f_1]^2). $$ WolframAlpha gives an explicit answer of $$ f_1\ast f_2 = \frac{e^{-\pi t^2/2}}{2 \sqrt{\pi}}, $$ and the plot looks like

However this does not match the image I get from my MatLab code (which can be found below).

Here is my MatLab code

L = 2;                              % length
M = 200;                            % number of samples
dx = L/M;                           % sample interval

x = -L/2:dx:L/2-dx;             % x coordinates
f_1 = exp(-pi*(x.^2));          % Gaussian 1
f_2 = exp(-pi*(x.^2));          % Gaussian 2 (same as 1)

figure(1)
plot(x,f_1,x,f_2,'--'); title('functions');
xlabel('x (m)');

%% Step 2. Perform convolution
F_1 = fft(f_1);                     % transform f_1
F_2 = fft(f_2);                     % transform f_2
F_0 = F_1.*F_2;                     % multiply pointwise
f_0 = ifft(F_0)*dx;                 % inverse transform and scale
f = fftshift(f_0);                  % center result

figure(1)
hold on
plot(x,f); title('Convolution');
xlabel('x')

So why do the images not match? The formulas and code seem very straightforward, and I've rechecked them several times, so I can't see why the images are inconsistent. Can anyone see what's going wrong?

Answer 1

Short answer. Both computations are correct, and they differ simply because the result from WolframAlpha is not the convolution of your functions.

---

Longer answer. Fourier transform is not a single transform. It is actually a family of transforms involving some parameters. For instance, Mathematica implements the following formulation

\begin{align*} \mathcal{F}[f](\omega) &= \sqrt{\frac{|b|}{(2\pi)^{1-a}}} \int_{-\infty}^{\infty} f(t) e^{ib\omega t} \, \mathrm{d}t, \\ \mathcal{F}^{-1}[F](t) &= \sqrt{\frac{|b|}{(2\pi)^{1+a}}} \int_{-\infty}^{\infty} F(\omega) e^{-ib\omega t} \, \mathrm{d}\omega \end{align*}

with the default parameter $(a, b) = (0, 1)$. Under this formulation, the convolution is Fourier-transformed to

$$ \mathcal{F}[f * g] = \sqrt{\frac{(2\pi)^{1-a}}{|b|}} \mathcal{F}[f] \mathcal{F}[g]. $$

Notice that we have extra constant factor. That said, in order for your formulation to give a correct answer, you should choose a correct parameter so that the above constant factor reduces to $1$. (And unfortunately, the default choice does not work for you.)

One convenient choice is $(a, b) = (0, -2\pi)$, which corresponds to the following unitary Fourier transform

$$ \mathcal{F}[f](\omega) = \int_{-\infty}^{\infty} f(t) e^{-2\pi i \omega t} \, \mathrm{d}t, \qquad \mathcal{F}^{-1}[F](t) = \int_{-\infty}^{\infty} F(\omega) e^{2\pi i \omega t} \, \mathrm{d}\omega. $$

Under this choice, you should get the correct answer as follows:

Answer 2

The fundamental problem with MatLab and Mathematica and similar software is that it's not open source and therefore what it does "under the hood" is not known. In usual practice this does not matter since they will give the expected answers in all but "edge cases". You have run into an edge case. It could be that one or both of the results you get are not what you expected and the best you can do is to experiment and see what is actually returned as a result and adapt to it. It could also be a case of a software bug but that would require more experimentation. Bugs can happen in all software. I suggest studying the documentation very carefully to try to determine what the software is supposed to return as a result.

That is the situation in general. For Fourier transforms there is the added problem that conventions can vary and so you have to be extra careful to find out exactly what conventions the system you are using has in place. That is why you need to read the documentation and experiment as you did. As mentioned in another answer, Mathematica, as usual, allows you to specify "options" as given in the documentation, and in particular $\texttt{FourierParameters}$, so it is important to get that right.