Free Books

FFT of a Simple Sinusoid

Our first example is an FFT of the simple sinusoid

$\displaystyle 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:
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)');
hold off;

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

% Same thing on a dB scale:
spec = 20*log10(magX); % Spectral magnitude in dB
plot(fn,spec,'--ok'); grid on;
axis([0 1 -350 50]);
xlabel('Normalized Frequency (cycles per sample))');
ylabel('Magnitude (dB)');
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

$\displaystyle \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

$\displaystyle \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 array indexes in matlab start at $ 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})$.

Next Section:
FFT of a Not-So-Simple Sinusoid
Previous Section:
Why a DFT is usually called an FFT in practice