Freuency Modulation (FM)¶

Let's start out with a quote from John Chowning's Original Paper:

"[...] frequency modulation is shown to result in a surprising control of audio spectra. The technique provides a means of great simplicity to control the spectral components and their evolution in time. Such dynamic spectra are diverse in their subjective impressions and include sounds both known and unknown."

As said in the above quote, Frequency Modulation (FM) is very simple. You take two oscillators, plug one into the frequency control of the other and are more or less ready to play around, get rich spectra and interesting, sometimes unexpected behavior.

But how exactly this works, what exactly can be expected, what other things we can do with FM; This can get rather involved. Here, we first cover a practical use case and then present some of the more hideen underlying aspects. A quick overview of what is covered here:

  • Basic FM Usage
  • Frequencies of the resulting sidebands
  • Amplitudes of the resulting sidebands
  • Estimation of Bandwidth, or 'spectral richness' (Carson's Rule)
  • When do we get harmonic spectra out of FM? When inharmonic?
  • Chowning's rapid determination figure
  • Dynamic Spectra
  • About the $\approx$ sign in Eq. 1, naive ideas to implement FM and connections phase modulation.
  • FM in radio. How to demodulate FM.
  • References

Basic FM¶

Frequency modulation can be described by:

$$ FM(t) \approx A sin(\omega_c t + \beta sin(\omega_m t)) \tag{1} $$

Where:

  • $t$, t, is the time in seconds.
  • $A$, A is the output Amplitude.
  • $\omega_c$, ω_c is the carrier angular frequency in radians per second.
  • $\omega_m$, ω_m is the modulator angular frequency in radians per second.
  • $\beta$, β is the modulation index.
In [1]:
%pylab inline
π = pi
import scipy.signal as sig
from IPython.display import Audio
import matplotlib.pyplot as plt
plt.style.use('website')

Populating the interactive namespace from numpy and matplotlib
In [2]:
T = 5 # [s] Total Duration 
sr = 44100 # [Hz] sampling rate
N = int(T*sr) # [] NUmber of Samples.
ns = arange(N) # Sample index array
t = ns/sr # [s] time array
In [3]:
f_c = 100 # [Hz], carrier Frequency
f_m = 5 # [Hz], modulator Frequency
ω_m = f_m*2*π # [rad/s] Angular modulator frequency
ω_c = f_c*2*π # [rad/s] Angular carrier frequency
In [4]:
A = 1 # [] Output Amplitude.
β = 5 # [] modulation index.

Equation (1) again, for simplicity of comparison: $$ FM(t) \approx A sin(\omega_c t + \beta sin(\omega_m t)) \tag{1} $$

In [5]:
fm_sig = A*sin(ω_c*t + β*sin(ω_m*t))

Let's plot the result (note that the wave form get's stretched and sqeezed, indicating lower and higher apparent fundamental frequency)

In [6]:
plot(t, fm_sig)
xlim([0, 0.5]);
xlabel('$t$ [s]');

and listen to it:

In [7]:
Audio(fm_sig,rate = sr)
Out[7]:
Your browser does not support the audio element.

Using a low modulation frequnecy $f_m$ and a low modulation index $\beta$, we achived a vibrato sound. Increasing both, we can achieve typical FM sounds, rich spectra.

In [8]:
f_c = 100 # [Hz], carrier Frequency
f_m = 500 # [Hz], modulator Frequency
ω_m = f_m*2*π # [rad/s] Angular modulator frequency
ω_c = f_c*2*π # [rad/s] Angular carrier frequency
β = 1.5 # [] modulation index.
fm_sig = A*sin(ω_c*t + β*sin(ω_m*t))
In [9]:
f, psd = sig.welch(fm_sig, fs=sr, nperseg=2**13)
semilogx(f, psd)
title('Power spectral Density')
xlabel('$f$ [Hz]')
ylabel('Linear Amplitude');
In [10]:
Audio(fm_sig,rate = sr)
Out[10]:
Your browser does not support the audio element.

A closer Look¶

The above is mostly what we need to know about FM in practice. But there are a couple of interesting questions left:

  • What are the frequencies of these sidebands
  • What are their Amplitudes?
  • Can we estimate the Bandwith?
  • When do we get harmonic spectra?
  • How to use FM to produce interesting sounds?
  • Eq (1) says FM is approximately this and that ($\approx$)???
  • How is FM Used in Radio?
  • What else?

What are the Frequencies of the Sidebands?¶

The frequencies are at: $$ f_c \pm f_m m $$ For any interger $m$. Let's test this.

In [11]:
f_c = 10_000 # [Hz], carrier Frequency
f_m = 1000 # [Hz], modulator Frequency
ω_m = f_m*2*π # [rad/s] Angular modulator frequency
ω_c = f_c*2*π # [rad/s] Angular carrier frequency
β = 5 # [] modulation index.
fm_sig = A*sin(ω_c*t + β*sin(ω_m*t))
In [12]:
f, psd = sig.welch(fm_sig, fs=sr, nperseg=2**10, scaling='spectrum')
plot(f, psd)
for m in arange(1, 10):
    axvline(f_c + f_m*m, ls='--', alpha=0.5)
    h1 = axvline(f_c - f_m*m, ls= '--', alpha = 0.5, label='$f_c \pm f_m m$')
h2 = axvline(f_c, label='$f_c$', ls= ':')
legend(handles=[h1, h2]);
xlabel('$f$ [Hz]')
ylabel('Linear Amplitude');

What are the Amplitudes?¶

The amplitude of the $m$th sideband (including the 0th sideband or carrier frequency) is determined by the Bessel function of order $m$ with argument $\beta$, the modulation index.

In [13]:
from scipy.special import jv #getting the bessel function 'jv'.
In [14]:
βs = linspace(0, 20, 100)
harms = arange(0, 4)
for h in harms:
    plot(βs, jv(h, βs), label=fr'$J_m(\beta), m={h}$')
grid(True)
legend()
xlabel(r'Modulation Index $\beta$');
title('Relative amplitude of sideband number $m$');
In [15]:
s = fft.fft(fm_sig, norm='forward')
faxis = fft.fftfreq(len(s), d=1/sr)
s = fft.fftshift(s)
faxis = fft.fftshift(faxis)

plot(faxis, abs(s))
plt.xlim([0, sr/2])

for m in arange(1, 10):
    axvline(f_c + f_m*m, ls='--', alpha=0.5)
    h1 = axvline(f_c - f_m*m, ls= '--', alpha = 0.5, label='$f_c \pm f_m m$')
h2 = axvline(f_c, label='$f_c$', ls= ':')

harms = arange(0, 10)
scaling = 0.5
for m in harms:
    scatter(f_c + f_m*m, abs(jv(m, β))*scaling, marker='x', c='r')
    h3 = scatter(f_c - f_m*m, abs(jv(m, β))*scaling, marker='x', c='r', label='Amp and freq prediction')
    
legend(handles=[h1, h2, h3])
xlabel('$f$ [Hz]')
ylabel('Linear Amplitude');

Can we estimate the Bandwith?¶

A Bandwidth estimation is given by the so called Carson's Rule: $$ BW \approx 2(d+f_m) \tag{2}$$ Where:

  • $d$, d is the peak frequency deviation given by $d=\beta f_m$
In [16]:
d = β*f_m
bw = 2*(d+f_m)
In [17]:
f, psd = sig.welch(fm_sig, fs=sr, nperseg=2**10, scaling='spectrum')
plot(f, psd)
axvline(f_c + bw/2, ls='--')
axvline(f_c - bw/2, ls='--')

ax = gca()
xstmark=f_c-bw/2
xenmark=f_c+bw/2
ystmark=0.07
an1=ax.annotate(' ',xy=(xstmark, ystmark), xycoords='data',
				xytext=(xenmark, ystmark-0.001),textcoords='data',
				arrowprops=dict(arrowstyle="<->",color='w'))
ax.annotate(r'$\approx BW$', xy=(((xstmark+xenmark)/2), ystmark+0.001), xycoords='data',fontsize=12.0,textcoords='data',ha='center');
xlabel('$f$ [Hz]')
ylabel('Linear Amplitude');

When do we get harmonic spectra?¶

We get harmonic spectra if there is a 'simple ratio' between the two frequencies $f_c$ and $f_m$. We define their ratio to be equals two integers $N_1$ and $N_2$: $$f_c/f_m = N_1/N_2 \tag{3}.$$ When $N_1$ and $N_2$ are integers, we get a harmonic spectrum. E.g. if $f_m = 2 f_c$, so the modulation frequncy is an octave above the carrier frequency, then we have $N_1 = 1$ and $N_2 = 2$. We can rearange to: $$ f_m = \frac{N_2}{N_1}f_c.$$ Setting $N_1$, we can call $N_2$ the ratio that is typically found in FM synthsizers.

Below are 2 examples, one with a ratio of $N_1/N_2 = 1/2$ and one with $N_1/N_2 = 1/\sqrt{2}$ since irrational numbers give inharmonic spectra.

In [18]:
f_c = 100 # [Hz], carrier Frequency
N1 = 1 # fixing to 1
N2 = 4 # ratio
f_m = (N2/N1)*f_c # [Hz], modulator Frequency
ω_m = f_m*2*π # [rad/s] Angular modulator frequency
ω_c = f_c*2*π # [rad/s] Angular carrier frequency
β = 2.5 # [] modulation index.
fm_sig = A*sin(ω_c*t + β*sin(ω_m*t))

f, psd = sig.welch(fm_sig, fs=sr, nperseg=2**12)
semilogx(f, psd)
title('Power spectral Density')
xlabel('$f$ [Hz]')
ylabel('Linear Amplitude');
for i in range(15):
    h2 = axvline(f_c + i*f_c, alpha = 0.5, ls=':', label='harmonic series')
legend(handles =[h2])
xlim([50, 3000])
Out[18]:
(50, 3000)
In [19]:
Audio(fm_sig, rate=sr)
Out[19]:
Your browser does not support the audio element.
In [20]:
f_c = 100 # [Hz], carrier Frequency
N1 = 1 # fixing to 1
N2 = sqrt(2) # ratio
f_m = (N2/N1)*f_c # [Hz], modulator Frequency
ω_m = f_m*2*π # [rad/s] Angular modulator frequency
ω_c = f_c*2*π # [rad/s] Angular carrier frequency
β = 2.5 # [] modulation index.
fm_sig = A*sin(ω_c*t + β*sin(ω_m*t))

f, psd = sig.welch(fm_sig, fs=sr, nperseg=2**12)
semilogx(f, psd)
title('Power spectral Density')
xlabel('$f$ [Hz]')
ylabel('Linear Amplitude');
for i in range(15):
    h2 = axvline(f_c + i*f_c, alpha = 0.5, ls=':', label='harmonic series')
legend(handles =[h2])
xlim([50, 3000])
Out[20]:
(50, 3000)
In [21]:
Audio(fm_sig, rate=sr)
Out[21]:
Your browser does not support the audio element.

Recreation of Chowning's 'rapid determination' Figure¶

Chowning in his paper gives a figure (Fig 8) that allows for estimation of bandwith. Let's recreate it.

In [22]:
nOrders = 900
nModInds = 900
βMax = 30
orderMax = 30
βs = linspace(0, βMax, nModInds)
orders = linspace(0, orderMax, nOrders)
bess = zeros([nOrders, nModInds])

for i,aβ in enumerate(βs):
    for j,o in enumerate(orders):
        bess[j, i] = abs(jv(o, aβ))
In [23]:
plt.imshow(bess, aspect='auto', extent = [0, βMax, 0, orderMax], origin='lower', cmap='jet')
xlabel(r'Modulation Index $\beta$')
ylabel('Number of side band')
grid(ls=':')
colorbar()
title('Amplitude of side frequencies');
In [24]:
X = βs
Y = orders
X, Y = np.meshgrid(X, Y)
Z = abs(jv(Y,X))

fig1 = plt.figure(figsize=[15,10], tight_layout=True)
ax=fig1.add_subplot(111,projection='3d')
ax.set_box_aspect((5, 5, 1))  # aspect ratio is 1:1:1 in data space
stride = 5
ax.plot_surface(X, Y, Z, vmin=Z.min() * 2, cmap='jet', rstride=stride, cstride=stride)
ax.view_init(elev=40., azim=30)

xlabel(r'Modulation Index $\beta$')
ylabel('Number of side band')
ax.set_zlabel('Relative Amplitude')
plt.show()

Dynamic Spectra¶

In [25]:
α = 1
αβ = 2
amp = exp(-t*α)
βt = exp(-t*αβ)

f_c = 100 # [Hz], carrier Frequency
N1 = 1 # fixing to 1
N2 = sqrt(2) # ratio
f_m = (N2/N1)*f_c # [Hz], modulator Frequency
ω_m = f_m*2*π # [rad/s] Angular modulator frequency
ω_c = f_c*2*π # [rad/s] Angular carrier frequency
β = 2.5 # [] modulation index.
fm_sig = amp*sin(ω_c*t + (β+βt)*sin(ω_m*t))

figure(figsize=[20,5])
subplot(121)
plot(t, amp, label='amp envelope')
plot(t, βt, label='mod envelope')
xlabel('$t$ [sek]')
legend()
subplot(122)
_ = specgram(fm_sig, Fs= sr, cmap='jet', vmin=-150)
xlabel('$t$ [sek]')
ylabel('$f$ [Hz]')
colorbar();
In [26]:
Audio(fm_sig, rate=sr)
Out[26]:
Your browser does not support the audio element.

About the '$\approx$' and Phase Modulation¶

Equation 1 says: $$ FM(t) \approx A sin(\omega_c t + \beta sin(\omega_m t)) \tag{1} $$

And thats more or less what the original paper by John Chowning says.

But one might say this looks like phase modulation because we know our oscillators as: $$ sin(\omega_0 t +\phi).$$ So what we added looks like it should be $\phi$.

Let's try what happens if we just really modulate the frequency. So let's see what happens if we just do this:

$$ FM_{maybe}(t) = sin \left( t 2 \pi (f_c + \beta f(t)) \right) $$ Where $f(t)$ is a modulating function. Looks reasonable, right?

In [27]:
T = 5 # [s] Total Duration 
sr = 44100 # [Hz] sampling rate
N = int(T*sr) # [] NUmber of Samples.
ns = arange(N) # Sample index array
t = ns/sr # [s] time array
In [28]:
f_c = 100 # [Hz], carrier Frequency
f_m = 5 # [Hz], modulator Frequency
β = 30 # [] modulation index.
In [29]:
f = sin(2*π*t*f_m)
fm_maybe = sin(t*2*π* (f_c + β*f))
In [30]:
plot(t, fm_maybe)
xlim([0, 0.5])
Out[30]:
(0.0, 0.5)

Looks great, right?¶

Let's look and listen

In [31]:
Audio(fm_maybe, rate=sr)
Out[31]:
Your browser does not support the audio element.
In [32]:
_=specgram(fm_maybe, Fs=sr, cmap='jet')
colorbar()
xlabel('$t$ [sek]')
ylabel('$f$ [Hz]')
Out[32]:
Text(0, 0.5, '$f$ [Hz]')

So this admittedly tempting version (that can also sometimes be found on the internet) does not quite work. If we want a combined signal, that is our frequency, where we can use a constant as a carrier frequency and add a signal to modulate, we need to think about instantaneous frequency. Basically we create this signal and intergrate it to produce a running phase.

Or as wikipedia puts it:

$$ {\begin{aligned}FM(t)&\ =\ A\,\sin \left(\,\int _{0}^{t}\left(\omega _{c}+B\,\sin(\omega _{m}\,\tau )\right)d\tau \right)\\&\ =\ A\,\sin \left(\omega _{c}\,t-{\frac {B}{\omega _{m}}}\left(\cos(\omega _{m}\,t)-1\right)\right)\\&\ =\ A\,\sin \left(\omega _{c}\,t+{\frac {B}{\omega _{m}}}\left(\sin(\omega _{m}\,t-\pi /2)+1\right)\right)\\\end{aligned}}$$

Here we can also see that our original Eq 1 is the same besides some phase offsets.

In [33]:
ω_c = f_c*2*π # carrier angular freq
ω_m = f_m*2*π # modulator angular freq
f = β*sin(ω_m*t) #modulator signal
instantaneousFreq = ω_c + f # Fundamental and modulator combined. Instantaneous frequency.
ϕ = cumsum(instantaneousFreq)/sr # computing instantaneous phase by integration of instantaneous frequency
fm_sig = A*sin(ϕ)
In [34]:
_=specgram(fm_sig, Fs=sr, cmap='jet')
colorbar()
xlabel('$t$ [sek]')
ylabel('$f$ [Hz]')
Out[34]:
Text(0, 0.5, '$f$ [Hz]')
In [35]:
Audio(fm_sig, rate=sr)
Out[35]:
Your browser does not support the audio element.

Application in RF Transmission¶

Please note: the author of this is not an RF expert. Also, the following is highly simplified. Critique highly welcome.

FM for information transmission is often combined with so called I/Q-Modulation. IQ stands for 'In-phase / Quadrature'. Basically we generate our frequency modulated signal similar to the above 'instantaneous frequency'-approach. After that, we use an I/Q-Encoding scheme that actually shifts the spectrum to the desired frequency bands.

These two links helped me understand I/Q Modulation and demodulation and also present some of the more interesting ways to look at this operation:

  • https://pysdr.org/content/sampling.html
  • https://www.microwavejournal.com/ext/resources/pdf-downloads/IQTheory-of-Operation.pdf?1336590796

In short: its all a bit complicated.

Creating a Test Signal¶

In [36]:
import soundfile as sf
useSample = True #Use drum sample or stick to sine waves as test input.
In [37]:
f_c = 7000 # [Hz], carrier Frequency
f_m = 100 # [Hz], modulator Frequency
ω_c = f_c*2*π # [rad/s] Angular carrier frequency

if useSample:
    message, sr = sf.read('../data/drumLoop.aif')
    message = average(message, axis=1) # to mono
    N =len(message)
    T = N/sr
    ns = arange(N)
    t = ns/sr

else:
    T = 5 # [s] Total Duration 
    sr = 44100 # [Hz] sampling rate
    N = int(T*sr) # [] NUmber of Samples.
    ns = arange(N) # Sample index array
    t = ns/sr # [s] time array
    ω_m = f_m*2*π # [rad/s] Angular modulator frequency
    message = sin(ω_m*t)

Doing the actual Modulation¶

In [38]:
carrier = sin(ω_c*t)
int_x = cumsum(message)/sr;

df = 1

xi=cos(2*π*f_c*int_x*df)# The i component of carrier wave
xq=sin(2*π*f_c*int_x*df)# The q component of carrier wave.

iqEncoded = xi*cos(2*π*f_c*t) - xq*sin(2*π*f_c*t) #QI encoded single fm signal for transmission.

Plotting the spectrum of the Modulated signal¶

In [39]:
s_i = fft.fft(xi, norm='forward')


s = fft.fft(iqEncoded, norm='forward')
f = fft.fftfreq(len(s), d=1/sr)
f= fft.fftshift(f)
s = fft.fftshift(s)
s_i = fft.fftshift(s_i)

plot(f/1e3, 20*log10(abs(s)), label='IQ-Encoded FM Signal')
plot(f/1e3, 20*log10(abs(s_i)), label='I')
axvline(f_c/1e3, label='carrier Frequency $f_c$', ls='--', lw=0.5)
legend()
ylabel('Amplitude [dB]')
xlabel('$f$ [kHz]');

Adding Some noise and (more) inperfectness to the transmission¶

In [40]:
detune = 0.9 #detune local oscillator (LO) at receiver.
noiseAmount = 0.01 # add white noise in transmission.

f_d = f_c*detune
iqEncoded = iqEncoded+(random.rand(len(iqEncoded))-0.5)*noiseAmount

Starting the Demodulation (I/Q)¶

In [41]:
xiRec = iqEncoded*cos(2*π*f_d*t) #reconstructed I signal
xqRec = iqEncoded*-sin(2*π*f_d*t) #reconstructed Q signal

cutoff = 4000
ω = cutoff/(sr/2)
b,a = sig.butter(11, ω) # rather steep lowpass filter

xiRec = sig.lfilter(b,a,xiRec) #reconstructed I signal, filtered
xqRec = sig.lfilter(b,a,xqRec) #reconstructed Q signal, filtered
In [42]:
plot(t, xq, label= 'Q')
plot(t, xqRec, label= 'Q recovered')
plot(t, xi, label= 'I')
plot(t, xiRec, label= 'I recovered')
xlim([0,0.01])
xlabel('$t$')
legend();
In [43]:
f, psd = sig.welch(xiRec, fs=sr,nperseg=2**14)
semilogx(f, 10*log10(psd), label='reconstructed I')
f, psd = sig.welch(xi, fs=sr, nperseg=2**14)
semilogx(f, 10*log10(psd), label= 'original I')
legend()
xlabel('$f$ [Hz]')
ylabel('Amplitude [dB]')
ylim([-200,None]);

FM-Demodulation¶

In [44]:
def fm_demod(x, df=1.0, fc=0.0):
    ''' Perform FM demodulation of complex carrier. There are other ways to do this (e.g. using a hilbert filter), but this works more or less.
    from: https://stackoverflow.com/questions/60193112/python-fm-demod-implementation

    Args:
        x (array):  FM modulated complex carrier.
        df (float): Normalized frequency deviation [Hz/V].
        fc (float): Normalized carrier frequency.

    Returns:
        Array of real modulating signal.
    '''

    # Remove carrier.
    n = np.arange(len(x))
    rx = x*np.exp(-1j*2*np.pi*fc*n)

    # Extract phase of carrier.
    phi = np.arctan2(np.imag(rx), np.real(rx))

    # Calculate frequency from phase.
    y = np.diff(np.unwrap(phi)/(2*np.pi*df))

    return y
In [45]:
complexFm = xiRec + 1j*xqRec
dem = fm_demod(complexFm, 1, fc=f_d/sr)*10
b,a = sig.butter(1, 20/(sr/2), btype='high') #dc offset removal
final = sig.lfilter(b,a, dem)

Comparing the Original to the Demodulated Signal¶

In [46]:
subplot(121)
plot(t, message, label='original')
plot(t[:-1], final, label='recovered', ls='-', c='r', alpha=0.5)
# xlim([0, 0.3])
xlabel('$t$ [sek]')
legend()

subplot(122)
f,psd = sig.welch(message, nperseg=2**12)
semilogx(f, 10*log10(psd), label='original')
f,psd = sig.welch(final, nperseg=2**12)
semilogx(f, 10*log10(psd),label='recovered', ls='-', c='r', alpha=0.5)
xlabel('$f$ [Hz]')
legend();
In [47]:
Audio(dem, rate=sr)
Out[47]:
Your browser does not support the audio element.

What Else?¶

  • What about using other waveforms than sine waves?
  • What about Multiple-Carrier FM (MC FM) to create formants?
  • What about Multiple-Modulator FM (MM FM)?
  • What about MM FM in parallel vs in Series?
  • What about Feedback FM?

Please refer to Roads, C. (1996). The computer music tutorial. MIT press..

References¶

  • 'The Synthesis of Complex Audio Spectra by Means of Frequency Modulation', ', John Chowning 1973, https://web.eecs.umich.edu/~fessler/course/100/misc/chowning-73-tso.pdf
  • Roads, C. (1996). The computer music tutorial. MIT press.
  • https://en.wikipedia.org/wiki/Frequency_modulation_synthesis#CITEREFChowning1973
  • https://en.wikipedia.org/wiki/Carson_bandwidth_rule#:~:text=In%20telecommunication%2C%20Carson's%20bandwidth%20rule,rather%20than%20a%20single%20frequency.
  • https://stackoverflow.com/questions/60193112/python-fm-demod-implementation
  • https://freqdemod.readthedocs.io/en/latest/intro.html#scientific-overview
  • https://github.com/shiva34/Digital-Modulation-Schemes/blob/main/FM/FM.m
  • https://pysdr.org/content/sampling.html
  • https://en.wikipedia.org/wiki/In-phase_and_quadrature_components
  • https://www.allaboutcircuits.com/technical-articles/should-i-and-q-combining-and-separation-be-done-digitally-or-analogly/
  • https://www.microwavejournal.com/ext/resources/pdf-downloads/IQTheory-of-Operation.pdf?1336590796