Why is there a difference in the phase when using fft versus a manual DFT on a cosine signal?

Question

Apologies if this is a duplicate, but I can't find a good answer anywhere.

If I have a time-domain cosine signal of the form:

$$x(t) = A\mathrm{cos}(\omega_0 t)$$

Then this will result in an Fourier transform with spectrum:

$$X(\omega)=\frac{A}{2}\delta(\omega-\omega_0) + \frac{A}{2}\delta(\omega+\omega_0) $$

I can generate these results in MATLAB using either fft or doing a summation of DFT. For example:

f = 1; %cosine frequency
A = 20; %cosine amplitude

w = 2*pi*f;
fs = 16; %sample rate in Hz
N = 2; %number of seconds
t = 0:1/fs:(fs*N)-1; %time vector

L = length(t); %length of time
df = fs/L; %frequency spacing
fAxis = (0:df:(fs-df)) - (fs-mod(L,2)*df)/2; %frequency vector


%time domain signal
x = A*cos(w*t);

%FFT of time domain signal
Xfft = fftshift(fft(x)/L);

%DFT of time domain signal
n=(-floor(L/2):1:floor(L/2))';
k=n';
Xdft = (1/L)*sum(x.*exp(-1i*2*pi*n.*k/L),2);

If I take the difference between the modulus of each spectrum, the difference is on the order of 10^-15. So the amplitudes are identical as expected.

However, I was surprised to find that the phases are quite different:

plot(fAxis,angle(Xdft)*180/pi,'--r');; hold on
plot(fAxis,angle(Xfft)*180/pi,'-k','LineWidth',2); hold on
grid on;
xlabel('Frequency (Hz)');
ylabel('Phase Angle (deg)')
legend('DFT','FFT')
set(gca,'FontSize',14)

Result:

The phase of the DFT jumps around between -180 and 180 degrees, which makes sense since the phase should be either -180 or 180 (it also has some zeros as well?). However, I'm confused by the fft phase angle. It looks like it is a constantly increasing function which phase wraps.

I understand that the real and imaginary parts are $\approx 0$ at all values not equal to $\omega_0$. So is this just a case of numerical imprecision at calculating the angle for small numbers? If this were the case, then at 1 Hz (i.e. where the delta spike is) I would expect the fft phase angle to match the DFT phase angle. But they don't there either. Is there something about the fft algorithm that causes the phase to be different when compared to DFT? Or am I making some error somewhere in my code/reasoning?

Answer 1

Check the Matlab definition of their FFT.

If I change the code

%DFT of time domain signal
n=(-floor(L/2):1:floor(L/2))';
k=n';
Xdft = (1/L)*sum(x.*exp(-1i*2*pi*n.*k/L),2);

to

%DFT of time domain signal
n = 1:L;
k=n';
Xdft = (1/L)*sum(x.*exp(-1i*2*pi*(n-1).*(k-1)/L),2);

and then plot the two

plot(fAxis,angle(Xfft)*180/pi,'-k','LineWidth',20); hold on
plot(fAxis,angle(Xdft)*180/pi,'r+');; hold on
grid on;
xlabel('Frequency (Hz)');
ylabel('Phase Angle (deg)')
legend('DFT','FFT')
set(gca,'FontSize',14)

Then I get the two overlapping precisely.

Note: I do not use fftshift in my plots.

If I plot the magnitudes:

plot(abs(Xfft));
hold on;
plot(abs(Xdft),'r+')

Then I get the following plot.

---

Why?

So let's look at the differences between the OP code and my code.

The first difference is between:

Xfft = fftshift(fft(x)/L);

and

Xfft = fft(x)/L;

This means that my labels (which I didn't change) are incorrect. So instead of the frequency range being:

fAxis = (0:df:(fs-df)) - (fs-mod(L,2)*df)/2; %frequency vector

it's really

fAxis = (0:df:(fs-df)); %frequency vector

And this leads to an error in my plots (which I am not currently going to correct, so you can see the error).

The other difference is in the direct calculation of the DFT

n=(-floor(L/2):1:floor(L/2))';
k=n';
Xdft = (1/L)*sum(x.*exp(-1i*2*pi*n.*k/L),2);

which is

\begin{align*} X_{\tt dft}[k] &= \frac{1}{L} \sum_{n=-\lfloor L/2 \rfloor}^{\lfloor L/2 \rfloor} x[n + \lfloor L/2 \rfloor] e^{-\jmath 2\pi n k /L}\\ &= \frac{1}{L} \sum_{n'=0}^{L-1} x[n'] e^{-\jmath 2\pi (n' - \lfloor L/2 \rfloor) k /L}\\ &= e^{\jmath 2\pi \lfloor L/2 \rfloor k/L} \frac{1}{L} \sum_{n'=0}^{L-1} x[n'] e^{-\jmath 2\pi n'k /L}\\ &= e^{\jmath \pi k} \frac{1}{L} \sum_{n'=0}^{L-1} x[n'] e^{-\jmath 2\pi n'k /L} \end{align*}

(assuming $L$ is even) versus

n = 1:L;
k=n';
Xdft = (1/L)*sum(x.*exp(-1i*2*pi*(n-1).*(k-1)/L),2);

which is

$$ X_{\tt dft}[k] = \frac{1}{L} \sum_{n=0}^{L-1} x[n] e^{-\jmath 2\pi n k /L} $$

where, in both cases, I'm assuming $x[n]$ is non-zero for $0\le n \le L-1$.

Reproducing with FFT

k = n' means the input is also fftshifted. Changing

Xfft = fftshift(fft(x)/L);

to

Xfft = fftshift(fft(ifftshift(x))/L);

yields

The few mismatches are likely explained by float limitations. Interestingly, something else happens if we do fftshift(x) instead, but that's for another day.