At what amplitude does spectral leakage start to be noticeable in the phase spectrum?

Question

I understand that discontinuities at the boundaries of a finite signal may introduce spectral leakage, that alter both the amplitude and phase spectra of the FFT, specially at low amplitudes. That is the explanation to several questions like (baffled by fft phase spectrum!"). In such questions, the examples are normally based on sine waves, where there is a dominant frequency, and the spectral leakage is noticeable in all the other frequencies.

Signals are often more complex, with a more continuous amplitude spectrum, but even in that case - specially if they are low pass-filtered, there is a point where there is a noticeable discontinuity in the amplitude and phase diagrams.

Example (in Octave):

N = 101;
x = exp(-linspace(-5,5,N).^2);
f = fftshift(fft(x));
subplot(211)
semilogy(abs(f))
subplot(212)
plot(unwrap(angle(f)))

Here the amplitude spectrum is bell-shaped as I expected at low frequencies. But when the amplitude goes below approximately 1e-11, the amplitude spectrum flattens.

At the frequencies where the amplitude spectrum is bell-shaped, the phase spectrum is a straight line with slope $\pi (N-1)/N$, which coincides with the effect of a time shift of $(N-1)/2$ from an even signal. At higher frequencies, it goes more or less linearly towards $\pi/2$ (although in this case the slope is not strictly constant).

Why there is such a harsh change in the spectra at that frequency? Is it a matter of numeric error in the computation of the FFT? The amplitude at the discontinuity is small (about 1e-11), but greater than machine error. And for other signals with narrower spectra, the amplitude at the discontinuity is higher. For instance, the same example but using a "wider" time signal:

x = exp(-0.5*linspace(-5,5*(N-1)/N).^2);

Answer 1

The issue is that you are not taking a large enough piece of your original x to get an accurate frequency domain representation.

Try setting

x = exp(-linspace(-20,20,N).^2);

If I do that, I get the bottom two plots below (the middle plots are for -10,10). The top plot is a replication of your values.

In the three pairs below, the left plot is the semilogy plot of x and the right plot is the semilogy plot of its FFT f.

---

Derivation of the FFT of a windowed Gaussian

To answer your question in the comments, we need to look at how a windowed version of: $$ x(t) = e^{-t^2/2} $$ can be analyzed in the frequency domain.

The windowed version is: $$ x_w(t) = w(t)e^{-t^2/2} $$ where $w(t) = \mbox{rect(t/T)}$ if we use a rectangular window of width T (span from $-T/2$ to $+T/2$).

The Fourier transform will be: $$ X_w(\omega) = \int_{-T/2}^{T/2} e^{-t^2/2} e^{j\omega t} dt\\ = 2 \int_{0}^{T/2} e^{-t^2/2} \cos(\omega t) dt\\ = \left . \sqrt{\frac{\pi}{2}} e^{-\omega^2/2} \left[ \mbox{erf}(\frac{t - j\omega}{\sqrt{2}}) + j \mbox{erfi}(\frac{ \omega - j t}{\sqrt{2}} ) \right] \right |^{T/2}_0\\ = \left . \sqrt{\frac{\pi}{2}} e^{-\omega^2/2} \left[ \mbox{erf}(\frac{t - j\omega}{\sqrt{2}}) + \mbox{erf}(\frac{ t + j\omega}{\sqrt{2}} ) \right] \right |^{T/2}_0 $$ where this integral is used and the identity $$ \mbox{erfi}(\frac{ \omega - j t}{\sqrt{2}} ) = \mbox{erf}(j \frac{ \omega - j t}{\sqrt{2}} ) / j = -j \mbox{erf}(\frac{ t + j\omega}{\sqrt{2}}) $$

Unfortunately, I'm not sure where that gets us. I'll leave it for now and see if I can doodle more later.

---

scilab code below (Sorry, I tend to use scilab rather than Octave).

// Q26236;
N = 101;
x1 = exp(-linspace(-5,5,N).^2);
f1 = fftshift(fft(x1));
x2 = exp(-linspace(-10,10,N).^2);
f2 = fftshift(fft(x2));
x3 = exp(-linspace(-20,20,N).^2);
f3 = fftshift(fft(x3));

clf
subplot(321)
plot2d("nl",abs(x1))
title('x1')
subplot(323)
plot2d("nl",abs(x2))
title('x2')
subplot(325)
plot2d("nl",abs(x3))
title('x3')

subplot(322)
plot2d("nl",abs(f1))
subplot(324)
plot2d("nl",abs(f2))
subplot(326)
plot2d("nl",abs(f3))

Answer 2

I think I have found an explanation of this phenomenon. Not a formal one, but perhaps someone else can make a smarter formulation out of it.

The situation is clearer with an even signal, like the following circular shift of the original example:

N = 101;
x = exp(-linspace(-5,5,N).^2);
x_even = [x(51:end), x(1:50)];
f_even = fftshift(fft(x_even));

The imaginary part of f_even should ideally be zero, although it is nonzero due to floating point errors. But the interesting thing is what happens to the real part, which is the one that determines the amplitude spectrum in this case. At high frequencies, exactly in the "flat" zone of the amplitude spectrum, real(f_even) oscillates between positive and negative values. See the following detail:

plot(real(f_even(1:34)))
hold on
plot(abs(f_even(1:34)),':')
legend({'Real','Abs'})

The source of such oscillation may be shown graphically by plotting the series that is summated in the computation of the DFT ($x_n \cdot e^{j2\pi\omega n/N}$):

% Only the real part, for e.g. w = 25
plot( x_even .* cos(2*pi*25*(0:100)/N), '.-' );

For each point in the "convoluted" series, at high frequencies there is another point in its vicinity of approximately opposite value, so the overall sum tends to zero as $\omega$ increases. But there is a finite number of samples, so that approximation is not smooth. The discontinuity occurs when the order of magnitude of the sum approaches the differences between consecutive samples in the original time series. (In this case the minimum nonzero value of diff(x_even) is 2.35e-11, which is about the amplitude where the discontinuity in the FFT occurs.) At higher frequencies, those minimum differences in the time series dominate the amplitude of the DFT, instead of the spectrum that could be expected from a continuous signal.

To complete the explanation, the phase spectrum is conditioned by the aggregation of this sampling effect and the numerical error seen in the imaginary part of f_even.