3
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Atul Ingle Oct 10 '17 at 16:28
  • $\begingroup$ Not really because usually round up errors are on the order of $10^{-16}$ and exhibit strong irregularities. $\endgroup$ – Ronan Tarik Drevon Oct 10 '17 at 16:37
  • $\begingroup$ I'm not sure I agree that rounding errors are always on the order of 10^-16. $\endgroup$ – Atul Ingle Oct 10 '17 at 16:39
  • $\begingroup$ Maybe but the curve is very smooth as well which is not typical of round of errors $\endgroup$ – Ronan Tarik Drevon Oct 10 '17 at 16:54
  • $\begingroup$ 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. $\endgroup$ – Fat32 Oct 10 '17 at 17:10
3
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Ronan Tarik Drevon Oct 10 '17 at 23:57
  • $\begingroup$ 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. $\endgroup$ – Fat32 Oct 11 '17 at 8:18
  • $\begingroup$ you can upvote the answer too if you found it useful. $\endgroup$ – Fat32 Oct 11 '17 at 8:21
  • $\begingroup$ 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. $\endgroup$ – Ronan Tarik Drevon Oct 11 '17 at 10:38
  • $\begingroup$ I'll put a few extra lines and discuss the error pattern in the answer asap. $\endgroup$ – Fat32 Oct 11 '17 at 12:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.