10. Hands-On: Synthesizing AM and FM in Python // Analog Communications // bhaswanth

A few lines of NumPy reproduce the math we have been writing. Run this and you can see and hear what we have been talking about.

python

import numpy as np
import matplotlib.pyplot as plt
 
fs = 200_000           # 200 kHz sample rate (we will work below Nyquist)
T  = 1.0               # 1 second of signal
t  = np.arange(int(fs*T)) / fs
 
# Message: a 1 kHz tone
fm = 1000.0
m  = np.cos(2*np.pi*fm*t)
 
# Carrier at 20 kHz (so we can plot sensibly without ultra-high sample rate)
fc = 20_000.0
c  = np.cos(2*np.pi*fc*t)
 
# AM with modulation index 0.7
mu     = 0.7
s_am   = (1 + mu*m) * c
 
# DSB-SC
s_dsb  = m * c
 
# FM with peak deviation 5 kHz (so beta = df/fm = 5)
df     = 5000.0
phase  = 2*np.pi*fc*t + 2*np.pi*df*np.cumsum(m)/fs   # integrate m
s_fm   = np.cos(phase)
 
# Spectra
def spectrum(x, fs):
    X = np.fft.rfft(x * np.hanning(len(x)))
    f = np.fft.rfftfreq(len(x), 1/fs)
    return f, 20*np.log10(np.abs(X) + 1e-12)
 
for name, sig in [('AM', s_am), ('DSB-SC', s_dsb), ('FM', s_fm)]:
    f, S = spectrum(sig, fs)
    mask = (f > 10_000) & (f < 35_000)
    plt.figure()
    plt.plot(f[mask], S[mask])
    plt.title(f'{name} spectrum')
    plt.xlabel('Hz')
    plt.ylabel('dB')
 
plt.show()

What you will see:

AM: a tall carrier line at 20 kHz, two side spikes at 19 and 21 kHz (the sidebands at $f_c \pm f_m$ ).
DSB-SC: only the two sidebands, no carrier.
FM: many sidebands at $f_c \pm n f_m$ for $n = 0, 1, 2, \ldots$ , with amplitudes following the Bessel function pattern. With $\beta = 5$ , you will see roughly 6 to 7 significant sidebands on each side, covering ~30 kHz of total bandwidth, exactly Carson's rule.

Try varying $\mu$ in AM (push it past 1; the envelope distorts). Try varying $\beta$ in FM (small $\beta$ looks like AM; large $\beta$ shows the Bessel comb spreading). Try demodulating the FM by computing the analytic signal:

python

from scipy.signal import hilbert
analytic = hilbert(s_fm)
inst_phase = np.unwrap(np.angle(analytic))
inst_freq  = np.diff(inst_phase) / (2*np.pi) * fs
# inst_freq oscillates around fc with deviation ~ df * m(t)
# subtract fc to get the message
recovered_m = (inst_freq - fc) / df

That last block is, in essence, a complete software FM demodulator. The same code (with a finer carrier and 200 kHz sample rate appropriate to broadcast FM) runs inside every RTL-SDR-based receiver.