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 (assume the signal is in units of Volts), 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. $$

There are two things that confuse me:

1. The sum in the definition of the FFT in Matlab runs from $1$ to $N$ (or $0$ to $N-1$, but let's ignore the Matlab indexing details for now!), but the Fourier transform is defined from $-\infty$ to $+\infty$. Does anyone have a nice intuitive way to see why the shift in the limits works? Is it simply because what we are really looking at here is the area under a curve, and it doesn't matter if you change where this area is positioned?

2. The second thing which 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 -----------------

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: