# Discrete Fourier Transform of the Gaussian

Cross-posted from here

I encountered the following question in a Digital Image Processing examination:

Find the 2D DFT of $$\frac{1}{2 \pi \sigma^2} e^{-\frac{(x - x_0)^2 + (y - y_0)^2}{2 \sigma^2}}$$ where $$x_0, y_0$$ are integers and $$x_0, y_0 \in \left[ -\frac{N}{2}, \frac{N}{2} \right]$$. We are considering an image where pixel intensities can be real but the pixels are discrete and in the range $$[0, N-1] \times [0, N-1]$$

# My approach

I used the formula which was given in the question paper: $$F(u, v) = \frac{1}{N} \sum_{x = 0}^{N-1} \sum_{y = 0}^{N-1} f(x, y)e^{-\frac{2 \pi j}{N}(ux + vy)}$$ For simplicity, I tried to first solve for the case where $$x_0, y_0 \geq 0$$ \begin{align} \therefore F(u, v) &= \frac{1}{N} \sum_{x = 0}^{N-1} \sum_{y = 0}^{N-1} \frac{1}{2 \pi \sigma^2} e^{-\frac{(x - x_0)^2 + (y - y_0)^2}{2 \sigma^2}}e^{-\frac{2 \pi j}{N}(ux + vy)} \\ \therefore F(u, v) &= \frac{1}{2 \pi \sigma^2 N} \sum_{x = 0}^{N-1} \sum_{y = 0}^{N-1} e^{-\frac{(x - x_0)^2 + (y - y_0)^2}{2 \sigma^2}}e^{-\frac{2 \pi j}{N}(ux + vy)} \\ \therefore F(u, v) &= \frac{1}{2 \pi \sigma^2 N} \left( \sum_{x = 0}^{N-1} e^{-\frac{(x - x_0)^2}{2 \sigma^2} - \frac{2 \pi j u x}{N}} \right) \left( \sum_{y = 0}^{N-1} e^{-\frac{(y - y_0)^2}{2 \sigma^2} - \frac{2 \pi j v y}{N}} \right) \\ \therefore F(u, v) &= \frac{1}{2 \pi \sigma^2 N} g(u) g(v) \end{align} Hee, we define $$g(\cdot)$$ as $$g(t) = \sum_{x = 0}^{N-1} e^{-\frac{(x - x_0)^2}{2 \sigma^2} - \frac{2 \pi j t x}{N}}$$ Now, I tried to evaluate $$g$$ using some of the standard methods used to evaluate summations in DFT. I tried:

• Trying to write the summation as a GP. However, due to the square term in the exponent, this is clearly not possible.
• Trying to telescope the series
• Trying to differentiate and see if the derivative is easier to evaluate

# The given solution

While the professor hasn't given a solution, he said that the DFT of the Gaussian is the Gaussian with the variance as the multiplicative inverse of the original Gaussian. While I know that this property is true for the Fourier Transform, I could not find any references online or in the reference texts provided that claim the same.

Can someone please help me understand what is the closed form solution and if it is indeed the Gaussian, how can I derive it? Thanks!

I'm not sure that

the DFT of the Gaussian is the Gaussian

is correct.

As Cris says:

The Gaussian is almost, but not quite exactly, band-limited, so sampling it will not produce an exact Gaussian in the frequency domain.

Looking at the paper

Brian Conolly and I. J. Good , A Table of Discrete Fourier Transform Pairs SIAM Journal on Applied Mathematics Vol. 32, No. 4 (Jun., 1977), pp. 810-822 (13 pages) Published By: Society for Industrial and Applied Mathematics

the closest they come is F15

where $$r$$ is the time index and $$s$$ is the frequency index and $$t$$ is the duration.

The notes section suggests that the proof is not trivial:

where reference 13 is:

I.J. Good, Analogues of Poisson's summation formula, Amer. Math. Monthly, 69 (1962), pp. 259-266.

As Cris notes in the comments, this might not make sense... so perhaps the next one is closer?

• Note that $\omega^{r^2}$ is not even related to the Gaussian. Is $\omega$ $^{-1} = e$? Then the transform of the Gaussian is proportional to $e^{s^2}$ (no minus sign in the exponent)? I’m not sure this makes sense? Jan 15 at 13:23
• @CrisLuengo Good point! I've added the next one, F16, which might be closer... but still not quite there. Thoughts? Part of the point is that, for the DFT, the Gaussian is not it's own transform.
– Peter K.
Jan 15 at 13:51
• I don’t know if an exact close form solution is possible. An infinite sum of Gaussians, shifted by $kN$ would transform into an infinite sum of Gaussians shifted by $kN$ ($N$ the number of samples). The point in my answer is that, if the Gaussians are small enough compared to $N$ both in the spatial and the frequency domain, you can ignore the overlap of the other Gaussians. But I don’t know how to transform that into an exact transform for a single Gaussian. Jan 15 at 14:58
• We need to start with a Gaussian multiplied by a square window, which transforms (continuous domain) to a Gaussian convolved with a sinc. Sampling the time domain leads to an infinite sum of that Gaussian convolved with a sinc, each copy shifted by $kN$. Jan 15 at 15:01
• Right, and that's the takeaway here: there is no closed form solution as there is in the continuous time case.
– Peter K.
Jan 15 at 15:14

The (continuous-domain) FT of a Gaussian is a Gaussian, as OP knows. To get the DFT pair from a FT pair we need to sample and crop the time domain, which is equivalent to replicating and sampling the frequency domain.

Knowing that the Gaussian goes to almost zero quite quickly, though it never actually becomes zero, we can sample the time domain such that the repetitions in the frequency domain have little overlap, and we can crop such that we remove only the very thin tails.

The error we’re making in the frequency domain is then very limited. But we are making an error, and so the correspondence is not exact.

Thus: given the full Gaussian bell is sampled (at least out to 3 sigma on either side), and is sampled sufficiently densely (at least with a spacing of sigma), then the DFT of that sampled Gaussian is approximately the sampled FT of that Gaussian.

• How does that provide any help in finding the closed-form solution for the Fourier transform of the Gaussian? Am I missing something? Jan 14 at 0:15
• A closed form for the continuous Fourier Transform exists. Not sure for the DFT. Jan 14 at 3:21
• @ZaellixA What I’m saying is that, starting from the continuous-domain pair, you can sample and crop the time domain, which replicates and samples the frequency domain. Ignoring the tails of the Gaussian in both domains, this gives you the approximate discrete transform pair. Jan 14 at 5:30
• Why not just sample the Gaussian in the Frequency Domain. Usually we want to apply a Gaussian Filter. When doing it spatially we just sample the Gaussian spatially. For this who want to apply Gaussian Filter using frequency domain, then it would be great to be able to calculate in a closed form (Programming wise) the best approximation of the ideal Gaussian Filter.
– Royi
Jan 14 at 10:01
• Thanks for the clarification @CrisLuengo. It wasn’t clear to me from the answer. Jan 14 at 11:55