1
$\begingroup$

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)))

FFT amplitude and phase of exponential signal

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);
$\endgroup$
  • $\begingroup$ Why the downvote? $\endgroup$ – Jim Clay Oct 5 '15 at 15:31
  • 1
    $\begingroup$ Wondering that, too. Seems like a good question! (+1, BTW) $\endgroup$ – Peter K. Oct 5 '15 at 15:56
1
$\begingroup$

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.

enter image description here


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))
$\endgroup$
  • $\begingroup$ Well, the exponential function was just a toy example. This seemingly happens with any signal that has a similar shape (in my case, it is a physiological response to a impulse-like stimulus). You're right: longer signals, (with longer "silent" sections), have a significant spectrum at higher frequencies, and that effect vanishes. What makes me wonder is that -when it happens- it looks like a clear discontinuity in the spectrum; I would like to know if there is some way to predict where that discontinuity occurs. $\endgroup$ – Helios Oct 5 '15 at 21:49
  • $\begingroup$ @Helios: that's an interesting question. Let me doodle with it offline and see if I can get a satisfactory answer. $\endgroup$ – Peter K. Oct 5 '15 at 22:01
  • $\begingroup$ I'd say that the discontinuity seen in the FFT is caused by the fact that the signal is discrete in time (see my later answer). So the analysis of the continuous Fourier Transform may give interesting conclusions, but not the answer to that particular situation. $\endgroup$ – Helios Oct 6 '15 at 21:50
1
$\begingroup$

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'})

Detail of the real part of the FFT

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), '.-' );

Convoluted series

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.

$\endgroup$

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.