FFT of a Simple Sinusoid

Our first example is an FFT of the simple sinusoid

$$x(n) = \cos(\omega_x n T)$$

where we choose $\omega_x=2\pi(f_s/4)$ (frequency $f_s/4$ Hz) and $T=1$ (sampling rate $f_s$ set to 1). Since we're using a Cooley-Tukey FFT, the signal length $N$ should be a power of $2$ for fastest results. Here is the Matlab code:

% Example 1: FFT of a DFT-sinusoid

% Parameters:
N = 64;              % Must be a power of two
T = 1;               % Set sampling rate to 1
A = 1;               % Sinusoidal amplitude
phi = 0;             % Sinusoidal phase
f = 0.25;            % Frequency (cycles/sample)
n = [0:N-1];         % Discrete time axis
x = A*cos(2*pi*n*f*T+phi); % Sampled sinusoid
X = fft(x);          % Spectrum

% Plot time data:
figure(1);
subplot(3,1,1);
plot(n,x,'*k');        
ni = [0:.1:N-1];     % Interpolated time axis
hold on; 
plot(ni,A*cos(2*pi*ni*f*T+phi),'-k'); grid off;
title('Sinusoid at 1/4 the Sampling Rate');
xlabel('Time (samples)'); 
ylabel('Amplitude');
text(-8,1,'a)');
hold off; 

% Plot spectral magnitude:
magX = abs(X);
fn = [0:1/N:1-1/N];  % Normalized frequency axis
subplot(3,1,2);
stem(fn,magX,'ok'); grid on;
xlabel('Normalized Frequency (cycles per sample))'); 
ylabel('Magnitude (Linear)');
text(-.11,40,'b)');

% Same thing on a dB scale:
spec = 20*log10(magX); % Spectral magnitude in dB
subplot(3,1,3);
plot(fn,spec,'--ok'); grid on;
axis([0 1 -350 50]);
xlabel('Normalized Frequency (cycles per sample))'); 
ylabel('Magnitude (dB)');
text(-.11,50,'c)');
cmd = ['print -deps ', '../eps/example1.eps'];
disp(cmd); eval(cmd);

Figure 8.1: Sampled sinusoid at frequency $f=f_s/4$. a) Time waveform. b) Magnitude spectrum. c) DB magnitude spectrum.

The results are shown in Fig.8.1. The time-domain signal is shown in the upper plot (Fig.8.1a), both in pseudo-continuous and sampled form. In the middle plot (Fig.8.1b), we see two peaks in the magnitude spectrum, each at magnitude $32$ on a linear scale, located at normalized frequencies $f= 0.25$ and $f= 0.75 = -0.25$. A spectral peak amplitude of $32 = (1/2) 64$ is what we expect, since

$$\hbox{\sc DFT}_k(\cos(\omega_x n)) \isdef \sum_{n=0}^{N-1} \frac{e^{j\omega_x n} + e^{-j\omega_x n}}{2} e^{-j\omega_k n},$$

and when $\omega_k=\pm\omega_x$, this reduces to

$$\sum_{n=0}^{N-1}\frac{e^{j 0 n}}{2} = \frac{N}{2}.$$

For $N=64$ and $\omega_x=2\pi f_s/4$, this happens at bin numbers $k = 0.25 N = 16$ and $k = 0.75N = 48$. However, recall that Matlab requires indexing from $1$, so that these peaks will really show up at indexes $17$ and $49$ in the magX array.

The spectrum should be exactly zero at the other bin numbers. How accurately this happens can be seen by looking on a dB scale, as shown in Fig.8.1c. We see that the spectral magnitude in the other bins is on the order of $300$ dB lower, which is close enough to zero for audio work $(\stackrel{\mbox{.\,.}}{\smile})$.