How to plot the periodic digital spectrum?

Question

In a previous answer, Dan Boschen wrote (emphasis mine):

First let me explain the "unfolded" digital spectrum: If you allow the frequency axis of the sampled signal to extend to $\pm \infty$, instead of limiting to the unique digital frequency range of $\pm > F_s/2$ (where $F_s$ is the sampling rate), you will see replicas of the original spectrum that is centered about 0 (DC) to also be similarly centered around every multiple of $F_s$. This is why we only need to show the spectrum from $\pm F_s/2$ (or even $0$ to $F_s/2$ for real signals) since this replicates everywhere else. However, I find this visualization helps immensely in understanding many concepts in multi-rate signal processing as well as bridging analog and digital systems.

Indeed this was an(other) eye-opener to me. Unfortunately, all the GUI tools I have at hand at this time only plot digital spectrum in the $0$, $F_s/2$ range--and I always thought it was the full spectrum. But as mentioned by Matt in a comment below "the spectrum is periodic".

Do you know how I can plot several periods of the unfolded digital spectrum using Matlab/Octave or Python Scipy?

Answer 1

As a simple example, consider the discrete-time signal

$$x[n]=a^nu[n],\qquad |a|<1\tag{1}$$

where $u[n]$ is the unit step function. The discrete-time Fourier transform (DTFT) of $x[n]$ is given by

$$X(e^{j\omega})=\frac{1}{1-ae^{-j\omega t}}\tag{2}$$

Note that $\omega$ is a normalized angular frequency:

$$\omega=2\pi f/f_s\tag{3}$$

where $f_s$ is the sampling frequency.

Also note from $(2)$ that $X(e^{j\omega})$ is a $2\pi$-periodic function in $\omega$, i.e., it is periodic with period $f_s$ in the frequency variable $f$. This periodicity of the spectrum is a general property of discrete-time signals.

Now you can plot the spectrum of $x[n]$ (i.e., $X(e^{j\omega})$) as a function of the frequency $f$, and depending on the chosen range of $f$ and the chosen sampling frequency, you will plot a certain number of periods. This is shown in the following Octave/Matlab script:

a = 0.9;
fs = 1;
fmin = 0;
fmax = 2.5;
f = linspace(fmin,fmax,1000);
X = 1 ./ (1 - a * exp(-1i*2*pi*f/fs) );

plot(f,abs(X))
xlabel('f')
title('|X(f)|')
grid on