Understanding the FFT documentation

As I understand, for a signal x, the frequency spectrum X = fft(x) is two-sided - the first half of the array is frequencies 0 to FS/2 (half the sampling frequency, i.e: Nyquist limit), while the second half has negative frequencies from -FS/2 to the lowest non-zero frequency (-1/N for a signal with N points).

Meanwhile the documentation for FFT (found using "help fft" in MATLAB R2019a) says:

                     N
       X(k) =       sum  x(n)*exp(-j*2*pi*(k-1)*(n-1)/N), 1 <= k <= N.
                    n=1

This aligns with my belief that X(1) contains DC offset (0 frequency), however I reckon this would mean that X(N) contains frequency N-1. It would also mean that X(N/2) contains frequency (N/2 - 1), which actually aligns with my statement at the top of the question.

Can you help me with my confusion around the equation from the documentation?

Answer 1

[Star Strider]

Your understanding essentially equates to my understanding of fft. It produces a ‘two-sided’ Fourier transform with the first value being the value of the d-c offset (mean of the time-domain signal). The mid-frequency is the Nyquist frequency. The second half is the flipped (not transposed) complex-conjugate of the first half. So the frequency is actually defined for the first half from 0 to the Nyquist frequency, and for the second half, from the negative Nyquist frequency to 0.

The fftshift function creates an accurate representation of this, so the frequency vector of the result goes from the -Nyquist frequency to the +Nyquist frequency.

I doubt if I clarified anything however, I likely just re-stated it.

[Will Hardiman]

Thanks for your answer. That affirms that I've understood how the FFT function works, now I just need to understand the equation from the documentation!

[Star Strider]

It simply muktiplies the complex exponential (from the Euler identities) by the value of the time domain funciton.

[Will Hardiman]

I understand that part, however that would mean that X(N) contains frequency 2*Nyquist frequency, surely?

[Star Strider]

It computes a two-sided Fourier transform, with the frequencieds defined as I previously described. It is not 2*Nyquist frequency, instead the frequency vector for the first half and the complex-conjugate of the first half of the resuult are flipped and concatenated to the end of the first half.

[Will Hardiman]

The complex exponential has for exp(j * 2 * pi * t * f), where j is the imaginary unit, t is a time variable, and f is the frequency, or number of cycles in the range t = 0:1. The documentation I quote above has the following equation:

                     N
       X(k) =       sum  x(n)*exp(-j*2*pi*(k-1)*(n-1)/N), 1 <= k <= N.
                    n=1

when k = N, the exponential is

                    exp(-j * 2 * pi * (N-1) * (n-1)/N)

here the "time variable" is (n-1)/N, which is in the range 0:(N-1)/N. This means the frequency should be (N-1). This is obviously not the case with the implementation of FFT that we have in MATLAB. As you have said, it is a two-sided Fourier transform, and X(N) contains frequency 1 (i.e.: 1 cycle per N points).

What I am saying (and have been trying to say since the start) is that the documentation does not match the reality. This equation from the documentation is not the correct equation for FFT in MATLAB because it describes decomposition of a time-domain signal x into frequency components ranging from f = 0 to f = 2*Nyquist. The FFT function in MATLAB does not do what the documentation says it does.

[Star Strider]

Well, the upper limit of the frequency is in reality not ‘2*Nyquist’ (actually 0 to ) , since the frerquency folds on itself as I described. The fft function cares not one whit about the actual frequency. It simply does the Fourier transform and lets your handle the frequency vector appropriately with respect to the time-domain signal.

[Will Hardiman]

Oh, I see what I have done here - I didn't properly consider the effect of the discrete sampling, and how that causes the frequency to fold. Thank you for you patience and explanations!

[Star Strider]

As always, my pleasure!

It’s definitely different from the result the Symbolic Math Toolbox would produce for a Fourier transform:

syms t w 
FTsp(w) = int(1*exp(1j*w*t), t, -1, 1);
[FTspn(w),FTspd(w)]= numden(FTsp);
FTsplh(w) = diff(FTspn)/diff(FTspd);                            % L’Hospital
FTsp(w) = piecewise(w<0,FTsp(w), w==0,FTsplh(w), w>0,FTsp(w));

figure
fplot(abs(FTsp), [-10*pi 10*pi])
xlabel('\omega')
ylabel('|F(\omega)|')
title('Fourier Transform Of A Pulse')
grid

EDIT — (24 Oct 2020 at 16:33)

Corrected typographical error.