I am trying to learn how to implement the FFT as a way to approximate the continuous-time Fourier transform, and as a "nice easy example" I have chosen to test it with a simple Gaussian pulse in the time domain, given by

$$ x(t) = A \exp \bigg(-\frac{t^2}{2\sigma^2} \bigg), \tag{1} $$

and I know that the analytic continuous-time Fourier transform of this function is given by

$$ X(f) = A\sigma\sqrt{2\pi} \exp \Big( -2\pi^2\sigma^2f^2 \Big), \tag{2} $$

using the definition of the Fourier transform

$$ \mathcal{F}\{x(t)\} = X(f) = \int_{-\infty}^{+\infty} x(t) e^{-2\pi ift} \,\,\textrm{d}t. \tag{3} $$ (Since $X(f)$ is real, we have $|X(f)| = X(f)$, and $\textrm{arg}(X) = 0$. This is because the Gaussian in Eq. (1) is centered at $t = 0$. If the signal was shifted in time to $t=t_0$, this would introduce an addtional phase shift factor of $e^{-2\pi ift_0}$, and $\textrm{arg}(X)$ would be non-zero.)

I then use Matlab to sample the Gaussian input signal in Eq. (1) and try to numerically obtain the result of Eq. (2) by implementing the FFT algorithm.

s = 2;                                % Sigma of input Gaussian signal [s]
A = 3;                                % Amplitude of input Gaussian signal [V]

Fs = 2;                               % Sampling rate [Hz]
N = 50;                               % Number of samples to collect

Ts = 1/Fs;                            % Sampling interval [s]
T = N*Ts;                             % Record window length [s]

t = -(T/2):Ts:(T/2-dt);               % Generate centered time vector
x = A*exp(-t.^2/(2*s^2));             % Create sampled verison of Gaussian

NFFT = N;                             % This indicates no zero-padding will be implemented
df = Fs/NFFT;                         % Frequency resolution (bin separation) [Hz]
f = -(Fs/2):df:(Fs/2-df);             % Generate frequency vector
             
X = fftshift(fft(ifftshift(x)))*dt;   % Continuous Fourier transform approximated by FFT
X_abs = abs(X);                       % Complex magnitude of FT
X((X_abs<1e-8)) = 0;                  % kill values below threshold, so phase is well-behaved
X_phase = unwrap(angle(X));           % Phase of FT

figure; plot(t,x)
figure; plot(f,X_abs)
figure; plot(f,X_phase)

This code produces the following plots, which agree very nicely with the analytic expression above.

The problem is that this code works, but I am not entirely sure exactly why it works! What I believe I am actually doing is approximating the continuous Fourier transform of Eq. (3) by a Riemann sum

$$ X(f_k) = X\Big(\frac{kF_s}{K}\Big) \approx \sum_{i=0}^{N-1} x(nT_s) \exp\Big( -2\pi i \frac{kn}{K} \Big) \times T_s. $$

What confuses me is that I have had to use ifftshift(x) (not fftshift(x), although they are equivalent if the number of points is even) before passing my Gaussian to fft(). I was careful to define my time vector in the code in such a way that the zero of time falls at the index $N/2+1 = 26$. This is because I know then that after ifftshift() the zero point will then be moved to the first entry in the vector (since I am using an even number of samples), as seen in the figure below:

Can someone explain why we have to do this, and how to interpret the fact that points corresponding to negative times are now at the end of the vector? I know this is necessary because if I don't do it then I obtain a linear phase shift, consistent with the signal not being centered at $t=0$.

The most common justification I keep seeing is

"It's because the FFT 'assumes' the time zero point is the first value in your vector!"

Is there a nice way to see this in terms of the approximation to the continuous case, using the indexing and Riemann sum as above, without anthropomorphizing the FFT algorithm and resorting to the above comment?

Ideally, I would like to see a mathematical proof which loosely says:

"with a centered time signal, the output of the FFT will only approximate the continous Fourier transform if the two halves of the vector are swapped, and zero is moved to the front, otherwise the phase will be wrong".

This would then justify the use of ifftshift().

EDIT 1 -------------------- Without shifting, it doesn't work -----------------------

At the request of DSP Rookie, here is an example to show that if the input is actually a shifted Gaussian signal (centered at $t=3$), then the expected linear phase is only obtained if ifftshift() is applied to the signal before the fft() command:

EDIT 2 ------------------------------- My own attempt ---------------------------------------

Here is my own attempt at what I am looking for, but I am struggling to finish it:

The goal is to calculate a numerical version of the Fourier transform: $$ X(f) = \int_{-\infty}^{+\infty} x(t) e^{-2\pi ift} \,\,\textrm{d}t $$ by using Matlab's fft() function.

Let's say I collect $N$ samples of $x(t)$ at specific times centered around zero-time with sampling period $T_s = 1/F_s$. The first sample is collected at time $t_1 = -T/2$, and the final sample collected at time $t_N=T/2 - T_s$, where $T = NT_s$. I can then define a sampled version of the signal which I will denote by $x[n]$ (I want to stick with Matlab's $n=1$ indexing convention for this exercise): $$ x[n] := x(t_n) = x\bigg( t=(n-1)\;T_s -\frac{T}{2}\bigg) \hspace{1cm} \textrm{for}\;\;\;\ n = 1,2,..., N $$ so that we have $$ \begin{align} x[1] &= x\bigg(t = -\frac{T}{2}\bigg)\\ x[2] &= x\bigg(t = -\frac{T}{2}+T_s\bigg)\\ &\vdots \\ x[N] &= x\bigg(t = \frac{T}{2}-T_s\bigg).\\ \end{align} $$ The sampled version of the signal would look like this:

I can put this sampled version into the CTFT formula (assuming I captured sufficiently the underlying signal) to obtain $$ \begin{align} X(f) &= \lim_{T_s\rightarrow0} \sum^{N}_{n=1} x[n] \exp \bigg( -2\pi i f \bigg[(n-1)\;T_s-\frac{T}{2}\bigg] \bigg)\cdot T_s\\ &\approx \sum^{N}_{n=1} x[n] \exp \bigg( -2\pi i f \bigg[(n-1)\;T_s-\frac{T}{2}\bigg] \bigg)\cdot T_s. \end{align} $$

Now in practice, the result is only going to give me values of $X(f)$ at certain frequencies $f_k$. Specifically, I am going to obtain $N$ values of the transform (neglecting zero-padding), centered around zero (because the CTFT that I am approximating runs from positive to negative), with the first frequency given by $f_1 = -F_s/2$ (the Nyquist). Each value in the vector will correspond to frequencies that increment by $\Delta f = F_s/N$, and so we can define a sampled verison of the transform by $$ X[k] := X(f_k) = X\bigg( f=(k-1)\;\Delta f -\frac{F_s}{2}\bigg) \hspace{1cm} \textrm{for}\;\;\;\ k = 1,2,..., N $$ so that we have $$ \begin{align} X[1] &= X\bigg(f = -\frac{F_s}{2}\bigg)\\ X[2] &= X\bigg(f = -\frac{F_s}{2}+\Delta f\bigg)\\ &\vdots \\ X[N] &= X\bigg(f = \frac{F_s}{2}-\Delta f\bigg).\\ \end{align} $$ Evaluating the approximation to the CTFT at these specific frequencies, I get $$ X[k] \approx \sum^{N}_{n=1} x[n] \exp \bigg( -2\pi i \bigg[(k-1)\;\Delta f-\frac{F_s}{2}\bigg] \bigg[(n-1)\;T_s-\frac{T}{2}\bigg] \bigg)\cdot T_s. \tag{4} $$ Now, this expression is clearly not that used by Matlab's fft() function, which is actually given here to be $$ Y[k] =\sum^{N}_{n=1} y[n] \exp \bigg( -2\pi i\; \frac{(n-1)(k-1)}{N} \bigg), \tag{5} $$ and so it is no surprise that simply passing the sampled signal straight to fft()*Ts without any modification does not give a way to approximate the CTFT of the orginal signal.

However, I do know that the command X[k] = fftshift(fft(ifftshift(x[n])))*Ts works, and produces the correct thing, because I have shown it in the graphs in the OP. Therefore, I must be able to show that

• $y[n]$ is an ifftshifted version of $x[n]$
• $X[k]$ is an fftshifted version of $Y[k]$.

In doing so, I should be able to express Eq. (4) in such a way that if ifftshift() is first applied to my signal $x[n]$, then Eq. (5) (the Matlab fft) is applied, and finally fftshift() is applied to the result, then I should obtain exactly the expression for $X[k]$, given by Eq. (4).

I can try to expand Eq. (4) as follows: $$ \begin{split} X[k] \approx &\sum^{N}_{n=1} x[n]\times ... & \\ &\exp \bigg( -2\pi i \bigg[(n-1)(k-1)\;\Delta f\; T_s - (k-1)\frac{T\Delta f}{2} -(n-1)\frac{T_sF_s}{2} + \frac{F_sT}{4}\bigg] \bigg)\cdot T_s \end{split} $$ and since $F_s = 1/T_s$, $F_sT=N$, and $\Delta f = F_s/N$, this simplifies to $$ \begin{split} X[k] &\approx \sum^{N}_{n=1} x[n] \exp \bigg( -2\pi i \bigg[\frac{(n-1)(k-1)}{N} - \frac{1}{2}(k-1) - \frac{1}{2}(n-1) + \frac{N}{4}\bigg] \bigg)\cdot T_s \\ &= \underbrace{\sum^{N}_{n=1} x[n] \exp \bigg( -2\pi i\; \frac{(n-1)(k-1)}{N} \bigg)}_{\text{This is the expression for Matlab's fft()}}\times... \\ &\hspace{3cm}\underbrace{\exp \bigg(-2\pi i\bigg[-\frac{(k-1)}{2} - \frac{(n-1)}{2} + \frac{N}{4}\bigg] \bigg)}_{\text{This must represent the ifftshift and fftshift somehow?}}\cdot T_s \hspace{1cm} (6) \end{split} $$

It is now in a form that I can recognise Matlab's fft() expression in there (although in the term I have indicated the fft() is applied to $x[n]$, not the shifted version), but this is where I am stuck - how can I show that this last equation is indeed calculated by the command

X[k] = fftshift(fft(ifftshift(x[n])))*Ts

I am not even sure how to mathematically represent these two operations (which essentially just swap the two halves of the vectors around). As a sanity check so far, I can also verify that Eq.(6) is correct up to this point, by manually evaluating the two things for a simple test vector:

x = [1 2 3 4 5 6];
N = length(x);

%%% Calculate fft using MATLAB's fft (with ifftshift / fftshift)
X_WithShifts = fftshift(fft(ifftshift(x)));

%%% Calculate fft using MATLAB's fft (without ifftshift / fftshift)
X_WithoutShifts = fft(x);

%%% Calculate X using manual formula, without using fft, ifftshift, or fftshift
X_Manual = zeros(1,N);
for k = 1:N
   for n = 1:N
       X_Manual(k) = X_Manual(k) + x(n)*exp(-2*pi*1i*(n-1)*(k-1)/N)*...
           exp( -2*pi*1i*( -(k-1)/2 - (n-1)/2 + N/4 ) );
   end    
end

which gives the same results for each of the shifted versions, and does not work correctly if you don't use the shifting functions:

X_WithShifts = 
      [ (3+0i) (-3-1.73i) (3+5.20i) (21+0i) (3-5.20i) (-3+1.73i) ]

X_Manual = 
      [ (3+0i) (-3-1.73i) (3+5.20i) (21+0i) (3-5.20i) (-3+1.73i) ]

X_WithoutShifts = 
      [ (21+0i) (-3+5.20i) (-3+1.73i) (-3+0i) (-3-1.73i) (-3-5.20i) ]