2
$\begingroup$

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:

enter image description here

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?

$\endgroup$
1
  • 1
    $\begingroup$ Be careful comparing expressions in continuous-time and continuous-frequency to their counterparts in discrete-time and discrete-frequency. There are scaling issues both in the independent variable ($t$ and $n$ or $\omega$ and $k$) and in the dependent variable ($x(t)$ and $x[n]$ or $X(\omega)$ and $X[k]$). To make a good direct comparison, you might want to review the Riemann integral. That's a good way to compare the integrals and dirac impulses in the continuous-time/frequency domain to the summations and kronecker impulses in the discrete-time/frequency domain. $\endgroup$ Jun 6, 2023 at 17:59

1 Answer 1

4
$\begingroup$

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.

Both phases plotted.

If I plot the magnitudes:

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

Then I get the following plot.

Magnutude responses.


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.

$\endgroup$
14
  • $\begingroup$ I get confused by your use of the semicolon. ;. Especially the double ;;. $\endgroup$ Jun 6, 2023 at 18:21
  • 1
    $\begingroup$ @robertbristow-johnson That was in the original code. I didn't notice! Matlab just ignores it, in effect. $\endgroup$
    – Peter K.
    Jun 6, 2023 at 18:41
  • $\begingroup$ But if I make the change as you suggested, then amplitudes won't align. Is this just a matter of shifting the x axis around? It causes no end in confusion the way that MATLAB requires the fftshift(fft(x))... $\endgroup$
    – Darcy
    Jun 6, 2023 at 20:39
  • 3
    $\begingroup$ @Darcy Go back to basics: go to the matlab site I linked to and see how they've defined it. Matlab doesn't require anything of the sort (fftshift). You may want to use if for display of centered data. I rarely do. I just plot the positive frequencies if it's an issue. $\endgroup$
    – Peter K.
    Jun 6, 2023 at 20:54
  • 1
    $\begingroup$ You seem to have forgotten to mention fftshift(fft(x)) -> fft(x). Your edit makes it a standard FFT. $\endgroup$ Jun 7, 2023 at 4:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.