How can I plot a correct fft of cosine wave?

I'm starting with Matlab and our professor told us to do some exercises to get it going.

One of them it's to plot the Fourier Transform of a cosine (cos(2*pi*t)), but sampled with sample frequency of fs=12Hz.

So I did the following:

>> fs=12;

>> n=0:1/fs:11/fs;

>> x=cos(2*pi*n);

>> stem(n,x)

Which gives me the following:

Then I calculated the fft of x:

>> y=abs(fft(x));

>> k=0:fs:11*fs;

>> stem(k,y)

Which should give me the real part of the Fourier Transform of a cosine, but If I recall correctly, the FT of a cosine is two spikes, one at (wave frequency)/2*pi and another at -(wave frequency)/2*pi, but I got this:

Why I'm getting both spikes at positive side of the axis? One could argue the DFT is periodic, which is, but it should be periodic with 2*pi, no?

I think I'm really messing up DTFT and DFT of cosine, but I can't figure out what can I do to correct my plot.

Thanks in advance.

[Adam]

It isn't the easiest thing to find in the help (so much so that I gave up!), but the output of the fft is as follows:

0 positiveFrequencies nyquist negativeFrequencies

with there being n/2 - 1 of each of positive and negative frequencies. The negative frequencies can equally well be regarded as frequencies between nyquist and the sampling frequency though, which is where they sit in that output, it is just convention as to how you wish to display it, so long as you understand what order you start with.

Personally I always just take the first n/2 + 1 points when looking at the spectrum because I am not interested in the negative frequencies. Some people like fftshift to move 0 to the middle, but I have never used it as I am used to looking at a double-side spectrum this way now.

[Henrique Orlandini]

So my answer is correct, but I'm just showing the positive spike and second one is just the negative one, but shifted because of the periodicity of DFT ?

Answer 1

[Image Analyst]

You forgot to use fftshift(). Remember, with fft, the origin is at element 1 and then "wraps around" so to see the origin in the middle like you'd expect, and positive frequencies to the right and negative frequencies to the left, you need to call fftshift().

[Henrique Orlandini]

I did like this now:

>> n=-11/fs:1/fs:11/fs;

>> x=cos(2*pi*n);

>> stem(n,x)

>> y=fftshift(abs(fft(x)));

>> stem(k,y)

>> k=-11*fs:fs:11*fs;

>> stem(k,y)

And I got this:

Which is closer to what I would expect, except for those point between the spikes and the fact that those spikes are located at k= -24 and 24.

[Image Analyst]

It's because you don't have a signal infinite in the x direction. It's multiplied by a rect function, so your fft will be two delta functions convolved with sinc functions (which gives two sinc functions) instead of two delta functions like you might think.

Answer 2

[John D'Errico]

What did you plot? Hint: look at the code you wrote:

k=0:fs:11*fs;
stem(k,y)

What is plotted as the x axis? The numbers 0:fs:11*fs.

Did MATLAB plot what you told it to plot? It did give you two spikes as you expected.

Answer 3

[Paul]

The problem statement is to plot the Fourier Transform of a cosine (cos(2*pi*t)), but sampled with sample frequency of fs=12Hz.

The discrete-time signal is therefore:

x[n] = cos(2*pi*n/Fs).

Because we are operating in discrete-time, the "Fourier Transform" in the problem statement (probably?) implies the Discrete-Time Fourier Transform (DTFT). The DTFT of x[n] is a periodic summation of a pair of (scaled) Dirac delta functions, which can't be computed using numerical methods.

Because x[n] is periodic (not every discrete-time sinusoid is periodic) with period N = 12, we can take the Discrete Fourier Transform (DFT) (as computed by fft) of one period of x[n] to compute complex, Discrete Fourier Series (DFS) coefficients of x[n].

N = 12;
n = 0:N-1;
Fs = 12; Ts = 1/Fs;
x = cos(2*pi*n*Ts); % x[n] = cos(2*pi*n*Ts)
X = fft(x);
max(abs(imag(X))) % imaginary part is roundoff error, discard

ans = 2.0500e-15

X = real(X);

If we interpret the DFT as a function of k, we have

figure
k = 0:N-1;
stem(k,X),xlabel('k')

Alternatively, we can interpret the DFT as a funtion of discrete frequency f_k = k/N*Fs, in which case we have

f_k = k/N*Fs;
figure
stem(f_k,X),xlabel('f_k (Hz)')

As shown, the frequencies with non-zero elements are at 1 and 11 Hz. In particular, the second non-zero element of X[k] in the "upper half" of the DFT corresponds to a positive frequency. And that's o.k. because the original signal can be reconstructed by summing a complex sinusoid at 1 Hz and another at 11 Hz, each scaled by the corresponding X[k] and divided by N per the convention of the DFT. Showing this explicitly, we have

sympref('FloatingPointOutput',false);
xr(sym('n')) = (6*exp(1j*2*sym(pi)/N*1*sym('n')) + 6*exp(1j*2*sym(pi)*11/N*sym('n')))/N

xr(n) = 

which can be shown to be exactly equal to x[n]. Negative frequency is not needed.

assume(sym('n'),'integer')
simplify(xr - cos(2*sym(pi)/Fs*sym('n')))  % xr[n] - x[n]

ans(n) = 

0

So where does this notion of the upper half of the output of fft correspond to negative frequencies come from? Well, both of the complex sinuosids are 2*pi periodic, so we can add or subtract integer multiples of 2*i*pi to either and obtain the same reconstruction. Subtracting 2*i*pi*n from the second we have

xr(sym('n')) = (6*exp(1j*2*sym(pi)/N*1*sym('n')) + 6*exp(1j*2*sym(pi)*11/N*sym('n') - 2i*sym(pi)*sym('n')))/N

xr(n) = 

which is

xr = simplify(xr)

xr(n) = 

as expected. Hence, it's true that we can interpret X[11] as corresponding to f_11 = -1 Hz, but that's not really necessary when working in discrete-time.

If we choose to interpret the upper half the of the DFS coefficients as corresponding to negative frequencies, then we can adjust the DFT using fftshift and recompute the f_k accordingly. Because N is even, we have

f_kshift = ((-N/2):(N/2-1))/N*Fs;
figure
stem(f_kshift,fftshift(X)),xlabel('f_k (Hz)')

●

and we see the non-zero elements at +-1 Hz, which might be more visually appealing, but is not mathematically necessary in discrete-time.

The alternative approach for dealing with x[n] is to multiply it with a window function w[n] and then take the transform (DTFT or DFT) of the product. If we have w[n] as

w[n] = 1, 0 < n < N-1, 0 otherwise,

then we have the values of x[n; 0 < n < N-1] and again take the DFT, which would yield the same result as above. However, in this case the underlying signal x[n]*w[n] is not periodic, so we can't interpret the DFT as DFS coefficients of x[n]*w[n] (but we could interpret the DFT as DFS coefficients of the N-periodic extension of x[n]*w[n], which is just x[n]). Instead, the approropiate interpretation is that the X[k] are frequency domain samples of the DTFT of x[n]*w[n]. Because x[n]*w[n] is causal and finite duration, we evaluate its DTFT using freqz. The DTFT is 2*pi periodic, so compute and display it for a few periods, and overlay the DFT samples.

[h,f] = freqz(x,1,linspace(-Fs,2*Fs,1000),Fs);
figure
plot(subplot(211),f,real(h)),xlabel('f (Hz)'),grid
hold on
stem(f_k,real(X))
plot(subplot(212),f,imag(h)),xlabel('f (Hz)'),grid
hold on
stem(f_k,imag(X))

Again, because of the inherent periodicity of the DTFT it's no problem to interpet the discrete frequencies f_k of the DFT as being all positive. The only thing that matters is that the DFT samples cover one period of the DTFT.

Now, suppose we want to use discrete-time processing of samples of a continuous-time signal x(t) to estimate the Continuous-Time Fourier Transform (CTFT) of the signal. The CTFT is (in general) not periodic. However, the DTFT of the sampled signal is a (appropriately scaled) periodic extension of the CTFT of x(t), so we can estimate the portion of the CTFT from -Fs/2 to Fs/2 by computing the "central" period of the DTFT of x[n] = x(n*Ts) from -pi to pi, or use the DFT to get frequency domain samples of that DTFT and then fftshift to get the central period.

Let x(t) = cos(2*pi*t)*r(t), where r(t) is 1 for 0 < t < 1, and 0 otherwise.

syms t f real
syms x(t)
x(t) = cos(2*sym(pi)*t)*rectangularPulse(0,1,t);

Its CTFT is

XX(f) = simplify(fourier(x(t),t,2*sym(pi)*f))

XX(f) = 

Compare to the fftshifted DFT

figure
fplot(subplot(211),real(XX(f)),[-Fs/2,Fs/2]),grid,xlabel('f (Hz)')
hold on
stem(f_kshift,fftshift(real(X))*Ts)
fplot(subplot(212),imag(XX(f)),[-Fs/2,Fs/2]),grid,xlabel('f (Hz)')
hold on
stem(f_kshift,fftshift(imag(X))*Ts)

So in this case of using discrete-time processing to estimate the transform of a continuous-time signal, we interpret the upper portion of the even-length DFT as corresponding to the portion of the CTFT of x(t) from -Fs/2 to -Fs/N.