# Demonstrating the Periodic Spectrum of a Sampled Signal Using the DFT

One of the basic DSP principles states that a sampled time signal has a periodic spectrum with period equal to the sample rate. The derivation of can be found in textbooks [1,2]. You can also demonstrate this principle numerically using the Discrete Fourier Transform (DFT).

The DFT of the sampled signal x(n) is defined as:

$$X(k)=\sum_{n=0}^{N-1}x(n)e^{-j2\pi kn/N} \qquad (1)$$

Where

X(k) = discrete frequency spectrum of time sequence x(n)

n = time index

k = frequency index

N = number of samples of x(n)

The time and frequency variables are related to n and k as follows:

t = nT_{s }(2)

f= k/(NT_{s}) = kf_{s}/N (3)

where T_{s} is the sample time interval in seconds and f_{s} = 1/T_{s} is the sampling frequency in Hz.

This article is available in PDF format for easy printing.

While n has a range of 0 to N-1, the range of k depends on the frequency range over which we want to compute X(k). For example, if we let k = 0 to N-1, Equation 3 yields a frequency range of f = 0 to f_{s}(N-1)/N, which is the usual range used for the DFT. For our demonstration, we’ll evaluate X(k) over the wider range of k= -2N to 2N-1, which gives a frequency range of f = -2fs to f_{s}(2N-1)/N.

The Appendix lists Matlab code that uses Equation 1 to compute X(k) for an example real-valued time sequence of length N= 32. Running the code generates Figure 1, which shows the time sequence, the magnitude of the DFT, and the dB-magnitude of the DFT. As advertised, the spectrum is periodic, with period f_{s}. Looking at Equation 1, we see that the value of the complex exponent repeats every time k crosses a multiple of +/-N, which coincides with frequencies that are multiples of +/-f_{s}.

Now that we have demonstrated the periodicity of X(k), it’s obvious why the DFT is normally evaluated over only N values of k: all of the information about the spectrum is contained in that range. However, we have a choice over the particular N values of k that we use, as shown in Figure 2. The top plot shows our demonstration X(k). The middle plot shows just the samples of X(k) for the usual range of k = 0 to N-1, while the bottom plot shows just the samples for k = -N/2 to N/2-1, an equally valid range.

Let’s look a little closer at the DFT evaluated over k = -N/2 to N/2-1. Figure 3 shows the real part, imaginary part, and magnitude of the DFT. The top two plots illustrate another property of the DFT: for a real time sequence, the DFT has a real part that is an even function and an imaginary part that is an odd function. This property also holds for the DFT evaluated over k = 0 to N-1; but in that case, the even and odd properties are defined with respect to f_{s}/2 Hz, instead of 0 Hz.

Finally, we should note that Equation 1, while useful for our demonstrations, is not the most efficient way to find the DFT. For efficient computation, we would use the Fast Fourier Transform (FFT) algorithm [3].

Figure 1. Periodic Spectrum of a sampled time signal.

Top: sampled time signal Middle: Spectrum magnitude Bottom: Spectrum dB-magnitude

Figure 2. DFT magnitude |X(k)| for different frequency ranges.

Top: k = -2N to 2N-1 results in f = -2f_{s} to f_{s}(2N-1)/N

Middle: Conventional DFT. k = 0 to N-1 results in f = 0 to f_{s}(N-1)/N

Bottom: DFT centered at f = 0. k= -N/2 to N/2-1 results in f = -f_{s}/2 to f_{s}(N/2-1)/N

Figure 3. DFT evaluated over k = -N/2 to N/2-1 or f = -f_{s}/2 to f_{s}(N/2-1)/N

Top: Real part, showing even symmetry.

Middle: Imaginary part, showing odd symmetry.

Bottom: Magnitude.

## Appendix Matlab Code to Evaluate DFT over several periods of its spectrum

For this example, the time signal is a Hann, or Hanning, pulse [4, 5] with the formula

x(n) = K*(1 – cos(2πn/P)),

where P = 9 and n= 0: P. This pulse is then padded with leading and following zeros to length N= 32.

% extended_dft.m 3/4/19 Neil Robertson % compute dft using its definition % k = -2*N:2*N-1f = -2fs:2fs % fs= 100; % Hz sample frequency (arbitrary value) N= 32; % number of time samples % sampled time signal of length N x= [0 0 0 0 7 24 43 55 55 43 24 7 0 0 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0]/258; % Compute DFT M= 4*N; % number of frequency samples for k= -M/2:M/2-1 % frequency index f(k+ M/2+ 1)= k*fs/N; % frequency vector. index range is 1:M sum_n= 0; for n= 0:N-1 % time index sum_n= sum_n + x(n+1)*exp(-j*2*pi*k*n/N); % DFT sum for X(k) end X(k+ M/2+ 1)= sum_n; % DFT vector. index range is 1:M end Xreal= real(X); Ximag= imag(X); Xmag= sqrt(Xreal.^2 + Ximag.^2); % DFT magnitude vector XdB= 20*log10(Xmag); % DFT dB-magnitude vector % % % % plot x, Xmag, and XdB subplot(311),stem(0:N-1,x),grid axis([0 N-1 0 .25]),xlabel('n') subplot(312),stem(f,Xmag,'markersize',4),grid axis([-2*fs 2*fs 0 1.1]) xticks([-2*fs -3*fs/2 -fs -fs/2 0 fs/2 fs 3*fs/2 2*fs]) xticklabels({'-2fs','-3fs/2','-fs','-fs/2','0','fs/2','fs','3fs/2','2fs'}) subplot(313),plot(f,XdB,'.-','markersize',7),grid axis([-2*fs 2*fs -50 5]) xticks([-2*fs -3*fs/2 -fs -fs/2 0 fs/2 fs 3*fs/2 2*fs]) xticklabels({'-2fs','-3fs/2','-fs','-fs/2','0','fs/2','fs','3fs/2','2fs'}) ylabel('dB')

## References

- Oppenheim, Alan V. and Shafer, Ronald W.,
__Discrete-Time Signal Processing__, Prentice Hall, 1989, section 3.2. - Rice, Michael,
__Digital Communications, a Discrete-Time Approach__, Pearson Prentice Hall, 2009, section 2.6.1. - Lyons, Richard G. ,
__Understanding Digital Signal Processing, 2__., Prentice Hall, 2004, Chapter 4.^{nd}Ed - Mathworks website, https://www.mathworks.com/help/signal/ref/hann.html
- Oppenheim and Shafer, Op. Cit., p. 447.

March 2019 Neil Robertson

**Previous post by Neil Robertson:**

Compute the Frequency Response of a Multistage Decimator

## Comments:

- Comments
- Write a Comment Select to add a comment

Hi Neil,

That is neat and useful. By the way I happened to compare DFT equation Vs iDFT. Apart from scaling,they look exactly same except for j/-j operator. Or am I wrong? So I tried both equations on a single tone @ + 0.1 Fs and can see DFT gives single line @ +0.1Fs as expected while iDFT gives same but @ -0.1Fs.

I wonder if the sign of (j) operator is actually chosen on arbitrary basis between DFT and iDFT. Moreover, I find it a bit hard to imagine how come one way the equation converts to frequency domain while the other way it converts to time domain.

Thanks

Kaz

Hi Kaz,

I'm glad you liked the article.

The continuous fourier transform is the integral of f(t)exp(-jwt)dt.

The inverse is the integral of 1/(2pi) * F(w)exp(jwt)dt.

If I ever knew how this arises, I have forgotten. Here is an article on computing the IDFT, by Rick Lyons:

Thanks Neil,

swapping Re/Im at input/output to fft, or inverting (Im) in effect forces +f of input to -f then back to +f at output... and this agrees with my test above. In practice we do rewiring of Re/Im as it costs nothing. Other three methods are costly but interesting.

Thus it confirms my conclusion that DFT and iDFT are almost same apart from transpose of Re/im input/outputs or cos/sin terms.

I can imagine DFT correlates set of cos/sin and picks up maximum as new samples indexed from zero to N-1. Hence it measures frequency content.

iDFT does same but using transposed cos/sin. How come this means un-corelating and converting from frequency domain to time domain. I can't "imagine".

Regards

Kaz

The choice of sign for which is forward and which is backward is by convention. So is the normalization factor.

The reason that the convention of using a negative sign is a little preferable is how the arg of the bin value relates to the phase of the signal. This is clearest in the complex case:

$$ x[n] = e^{ \omega n + \phi } $$

Suppose the DFT frame is a whole number of cycles, and $ k = \omega \cdot \frac{N}{2\pi} $, then

$$ X[k] = (value) \cdot e^{\phi} $$

By using the negative sign in the forward DFT the signs on $\phi$ align, that make sense.

What I have an issue with is the normalization factor convention. I can understand a factor of "1" in library routines to save a multiply, but when you are doing the math, it is less sensible.

The "true" normalization factor should be $ \frac{1}{\sqrt{N}} $, but that is a pain to deal with in practicality. It is the "true" one because it applies to the forward and inverse DFT making the operations "multiplications of unitary matrices" in Linear Algebra terminology. The common convention of using "1" for the forward and "1/N" for the inverse is backward in my opinion. It should be the reverse.

If you use the reverse, as I advocate, then the $(value)$ in the above equation becomes "1" for complex signals and "1/2" for real ones,

**making the interpretation of bin values independent of the sample count**. This is much more sensible and why I stubbornly use a "1/N" in all my articles even though it is immaterial for frequency calculations since those are amplitude independent.

In addition, in the cases when your DFT frame covering a whole number of cycles of a periodic signal, then the coefficients for the Fourier Series for the continuous function can be read straight from the DFT bin values. This should be the clincher.

Ced

Thanks Ced,

For FPGA/ASIC platforms as opposed to math or soft platforms we got commercial fft/ifft cores that scale by sqrt(N) for either direction. Still not practical due to bit growth limitation. I normally scale output back by 1/sqrt(N) using precomputed single constant plus one multiplier.

I am still not sure how to mentally explain the ifft equation...

Kaz

Conceptually (and literally), the bin values are the coefficients to the Sine and Cosine functions that reconstruct your function.

The inverse DFT is that reconstruction.

To post reply to a comment, click on the 'reply' button attached to each comment. To post a new comment (not a reply to a comment) check out the 'Write a Comment' tab at the top of the comments.

Registering will allow you to participate to the forums on ALL the related sites and give you access to all pdf downloads.