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}
$$

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.

```lang-matlab
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. 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.)

[![enter image description here][1]][1]
[![enter image description here][2]][2]

The problem is that this code works, but I am not entirely sure exactly *why* it works!

What confuses me is that I have had to use ```ifftshift(x)``` before passing my Gaussian to ```fft()```, and ```fftshift(x)``` afterwards. 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:

[![enter image description here][4]][4]

Can someone explain why we have to do this. The most common justification I keep seeing is

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

or

**"It's just because of the periodicity of the DFT!"**

Is there a nice way to see this *explicitly* in terms of the approximation to the continuous case, without anthropomorphizing the FFT algorithm and resorting to the above comments?  

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

*"if your time-domain signal is centered, the output of the FFT will only approximate the continous-time Fourier transform if the two halves of the vector are swapped, and you first move your zero to the front, and then move it back to the middle again afterwards. Otherwise the phase obtained will be wrong".* 

A proof of this nature would then justify the use of ```ifftshift()``` and ```fftshift()``` for me. 

**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, and ```fftshift()``` applied afterwards:
[![enter image description here][5]][5]

**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:
[![enter image description here][6]][6]


I can write the CTFT formula as a Riemann sum, and insert this sampled signal into it (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][3] 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, and so this command must be equivalent to 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 have a feeling that there might be a way to use the shift theorem to show that if you circularly shift the vectors through ifftshift and fftshift then you can cancel the extra phase factors, and the fft can be applied.

As a quick sanity check on Eq. (6) so far, I can verify it by manually evaluating the two things for a simple test vector of six numbers:
```lang-matlab
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 Eq. (6) 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) ] % WRONG
```
  [1]: https://i.sstatic.net/VfkhS.png
  [2]: https://i.sstatic.net/pYCLB.png
  [3]: https://www.mathworks.com/help/matlab/ref/fft.html
  [4]: https://i.sstatic.net/q4fdg.png
  [5]: https://i.sstatic.net/7bbp7.png
  [6]: https://i.sstatic.net/UksYv.png