# Discrepancy between Gaussian FFT and its Fourier Transform

I am trying to do the FFT of a Gaussian signal and comparing it to the theoretical Fourier transform. For infinitely small time step $dt$ and infinitely long signal length $T$, the 2 should become identical right?

• Could this just be a rounding issue? I ran your code and it looks like the largest error was on the order of 10^-5 percent. – Atul Ingle Oct 10 '17 at 16:28
• Not really because usually round up errors are on the order of $10^{-16}$ and exhibit strong irregularities. – Ronan Tarik Drevon Oct 10 '17 at 16:37
• I'm not sure I agree that rounding errors are always on the order of 10^-16. – Atul Ingle Oct 10 '17 at 16:39
• Maybe but the curve is very smooth as well which is not typical of round of errors – Ronan Tarik Drevon Oct 10 '17 at 16:54
• can you please precisely state the error that you consider? And what was your expectation ? Since there are a number of error sources it's important to fix the one you want to discuss. – Fat32 Oct 10 '17 at 17:10

First consider the theoretical DTFT of $x[n]$ $$X(e^{j\omega}) = \sum_{n=-\infty}^{\infty} x[n] e^{-j \omega n} = \sum_{n=0}^{\infty} a^n e^{-j \omega n} = \frac{1}{1 - a e^{-j \omega} }$$

Then consider the N-point DFT of the truncated signal: $$X[k] = \sum_{n=0}^{N-1} x[n] e^{-j \frac{2\pi}{N} n k} = \sum_{n=0}^{N-1} a^n e^{-j \frac{2\pi}{N} n k} = \frac{1 - a^N}{1 - a e^{-j \omega} }$$

A simple inspection of those DTFT abd DFT yields the following conclusion: $$\boxed{ X[k] = (1-a^N) X(e^{j\omega_k}) }$$

Where the DTFT is evaluated at the DFT frequency of $w_k = \frac{2\pi}{N}k$.

And the error between $X[k]$ and $X(e^{j \omega_k})$ is: $$E[k] = X(e^{j \omega_k}) - X[k] = a^N X(e^{j \omega_k})$$

As can be seen, the error due to the finite length DFT truncation is multiplied by the constant $a^N$ and the frequency dependent variable $X(e^{j \omega_k})$ for each DFT frequency of $\omega_k$. It's obvious from this that the error plot will have a shape of the actual DTFT $X(w)$ scaled by $a^N$.

In order to make this error smaller than the floating point round-off errors (such as $10^{-16}$ for 64 bit IEEE floating point numbers) you have to make sure that for all frequencies $\omega_k$ for $0 \leq k \leq N-1$ you should satisfy $E[k] < 10^{-16}$ that is:$$E[k]=a^N X(e^{j \omega_k}) < 10^{-16}$$

Note that this error is frequency dependent due to $X(e^{j \omega_k})$ term: the higher the DTFT at that frequency, the higher will be the error and vice versa. So for a randomly fluctuating DTFT $X(e^{j \omega_k})$ the error will also randomly fluctuate.

In the analysis of the Gaussian signal $x[n] = e^{-\alpha n^2} = a^{n^2}$ where $0 < a = e^{-\alpha} < 1$ , an analytic expression for the error will be difficult to simplify as you have to evaluate the following sum:

$$X_g(e^{j\omega}) = \sum_{n=-\infty}^{\infty} a^{n^2} e^{-j \omega n} = \beta a^{\gamma w^2}$$ where $\beta$ and $\gamma$ should be determined.

And the N-point DFT of the truncated (and shifted by $N/2$) Gaussian signal: $$X_g[k] = \sum_{n=0}^{N-1} a^{(n-N/2)^2} e^{-j \frac{2\pi}{N} n k} = ?$$

Unfortunately since the series sum may not have a closed form expression (or it will be quite complicated) therefore a simple formulation of the error does not exist as were the case for the simpler signal $x[n] = a^n u[n]$.

EDIT: The following exposition shows the same results as attained above however from a different perspective which may be less intuitive to follow. I suggest the above perspective for a clear understanding and the below perspective a complementary case.

Consider the following infinite length signal $x[n] = a^n u[n]$, with $|a|<1$, where the discrete-time Fourier tranform (DTFT) $X(e^{j\omega})$ of $x[n]$ is:

$$X(e^{j\omega}) = \text{DTFT} \{x[n] \} = \sum_{n=0}^{\infty} a^n e^{-j \omega n} = \frac{1}{1-a e^{-j\omega} }$$

In order to practically compute the DTFT $X(e^{j\omega})$ using a digital computer, the relative transform, DFT (discrete Fourier transform) is used. For finite length signals it's well known that DFT $X[k]$ will be samples of DTFT $X(e^{j \omega})$. But for infinite length signals (such as $x[n]=a^n u[n]$) the computed N-point DFT (which should practically use a finite length truncated version of the infinite length signal) will be an approximation to the true theoretical DTFT. The question is ; what's the relation between the true DTFT spectrum $X(e^{j\omega})$ and the computed DFT samples $X[k]$ ?

Now it's obvious that for infinite length signals the values of DFT $X[k]$ will not be the samples of the function $X(e^{j\omega})$. For example given $x[n]=a^nu[n]$ and $X(e^{j \omega}) = \frac{1}{1-a e^{-j\omega}}$ then we see that $$X[k] \neq X(e^{j\frac{2 \pi}{N}k}) = \frac{1}{1-a e^{-j\frac{2 \pi}{N}k}}$$

First lets consider the following, assume you have artificially obtained $N$ uniform samples from the function $X(e^{j\omega})=\frac{1}{1-a e^{-j\frac{2 \pi}{N}k}}$ as
$$X_f[k] = \frac{1}{1-a e^{-j\frac{2 \pi}{N}k}}$$ where $X_f[k]$ stands for the DFT samples of a finite length signal $\tilde{y}[n]$ whose N-point DFT would be $X_f[k]$. Therefore we can assert that the inverse DFT of the $X_f[k]$ would yield the N-point discrete-time sequence $\tilde{y}[n]$ as:

$$\tilde{y}[n] = IDFT\{ X_f[k] \} = \frac{1}{N} \sum_{k=0}^{N-1} \frac{1}{1-a e^{-j\frac{2 \pi}{N}k}} e^{ j \frac{2\pi}{N} k n}$$

It can be shown from the periodicity assumption of the IDFT that the result is also equivalent to: $$\tilde{y}[n] = IDFT\{ X_f[k] \} = \sum_{r=-\infty}^{\infty} x[n-rN]$$

Where the signal $x[n]$ is nothing but the theoretical infinite length inverse of the DTFT $X(e^{j\omega}) = \frac{1}{1-ae^{-j\omega}}$ which is as we know $$x[n] = IDTFT \{\frac{1}{1-ae^{-j\omega}} \} = \frac{1}{2\pi} \int_{-\pi}^{\pi} \frac{1}{1-a e^{-j\omega}} e^{j\omega n} d\omega = a^n u[n]$$

Therefore we have the following:

$$\tilde{y}[n] = \sum_{r=-\infty}^{\infty} x[n-rN] = \sum_{r=-\infty}^{\infty} a^{n-rN}u[n-rN]$$

It will be interesting to consider the value of $\tilde{y}[n]$ for the range $0 \leq n \leq N-1$ for which we have:

$$\tilde{y}[n] = a^n \sum_{r=-\infty}^{\infty} a^{-rN}u[n-rN] = a^n \sum_{r=-\infty}^{0} a^{-rN} = a^n \sum_{r=0}^{\infty} (a^N)^r = a^n \frac{1}{1-a^N}$$

Which shows the relation between finite length signal $\tilde{y }[n]$ and the infinite length signal $x[n] = a^n u[n]$ as: $$\tilde{y}[n] = \frac{x[n]}{1-a^N}$$ where $x[n] = a^n u[n]$.

One can see that the modification involved on the computed value of finite length approximation to the infinite length signal $x[n]$ is a constant weight of $K = \frac{1}{1-a^N}$ for the signal $x[n]= a^n u[n]$

Now using the linearity properties of DFT one can deduce the following: Given the pair of N-point sequence $\tilde{y}[n]$ and its N-point DFT $X_f[k]$ which is the uniform samples of the true theoretical DTFT $X(e^{j\omega})$ associated with infinite length signal $x[n]=a^n u[n]$ we have the following

$$\tilde{y}[n] \leftrightarrow X_f[k]$$ $$\tilde{y}[n] = K x[n] \leftrightarrow X_f[k] \Rightarrow x[n] \leftrightarrow \frac{1}{K} X_f[k] = (1-a^N) X_f[k]$$

This last statment is concluded as follows: given an N-point finite length truncation $x_N[n]$ of the infinite length signal $x[n] = a^n u[n]$ for the range $0 \leq n< N$, then the N-point DFT of the N-point sequence $x_N[n]$ will be $(1-a^N)X[k]$ where $X[k]$ are the samples of the true DTFT spectrum $X(e^{j\omega}) = \frac{1}{1 -a e^{-j \omega}}$ associated with the infinite length signal $x[n]=a^n u[n]$. So the error associated with the DFT computation is the fixed weight $1-a^N$ per sample of the DFT.

You should perform a similar analysis for the Gaussian signal $x[n] = \beta e^{-\alpha n^2}$ to find the error associated per computed DFT sample.

• That's looks rigorous indeed but could you just say that the computed DFT being the finite sum of geometrical series is then $\Sigma_1^N (a\Omega)^n=\frac{1- (a\Omega)^N}{1-a\Omega}$ where $\Omega=e^{-j\omega}$? So you confirm that the difference DTFT-DFT should rapidly get below round-off errors. – Ronan Tarik Drevon Oct 10 '17 at 23:57
• Th edifference $X(e^{jw}) - X[k] = X(e^{jw}) - (1-a^N)X(e^{jw}) = a^N X(e^{jw})$ So the difference (error) depends on 1-) the factor $a^N$, 2-) the value of $X(e^{jw})$ at the frequency of investigation for the signal of $x[n] = a^n u[n]$. So for a roundoff error you look for $a^N X(e^{jw}) = 10^{-16}$ at a given frequency. – Fat32 Oct 11 '17 at 8:18
• you can upvote the answer too if you found it useful. – Fat32 Oct 11 '17 at 8:21
• That was useful indeed, although I realised that line "df = 1/((N-1)dt) " in the code should be replaced by "df = 1/(N*dt)" so the error decrease uniformly with increasing sampling length N. – Ronan Tarik Drevon Oct 11 '17 at 10:38
• I'll put a few extra lines and discuss the error pattern in the answer asap. – Fat32 Oct 11 '17 at 12:04