FM Radio & Demodulation¶

FM Radio is a technology that can be understood much more easily in the frequency domain, using techniques we have covered in this class. Using a Software Defined Radio (SDR), I will be analyzing a local FM station, KAMU 90.9MHz.

This station uses Wideband FM, meaning audio is carried with a 20 kHz bandwidth and maximum 75 kHz deviation from the carrier frequency (90.9 MHz). This means we should expect to see bands ranging from 90.825 MHz - 90.975 MHz.

import IPython.display as ipd
ipd.Image(filename="fmradio.png")

The above image is a spectrum analysis of the radio station when no audio (silence) is being broadcast. It is important to note that this frequency domain plot does not directly corresond to the audio signal since some demodulation has to be done. Nevertheless, the various subcarrier frequencies are abundantly clear, and as predicted, they do indeed range from 90.825 MHz - 90.975 MHz.

Additionally, the "blocks" on either side of the spectrum are related to HD radio. The specific details are out of the scope of this project, but their presence is interesting to note.

What would this spectrum look like in the time domain?¶

To see this, I have recorded the raw signal (before demodulation) centered at 90.9 MHz.

import matplotlib.pyplot as plt
import numpy as np
import librosa
from scipy.fftpack import fft, fftfreq, fftshift, ifft, rfft
from scipy.signal import decimate

x, sr = librosa.load('audio2.wav', sr=None)
samples = x.shape[0]
length = samples / sr
xcopy = x

time = np.linspace(0., length, samples)
plt.figure(figsize=(14,7))
plt.plot(time, xcopy)
plt.xlabel("Time (sec)")
plt.ylabel("Amplitude")
plt.title("Raw FM Signal - Time Domain")
plt.show()

The above plot shows the time domain signal of the raw, modulated FM signal. It looks relatively normal, but notice that the amplitude stays about the same at all times. This is one of the fundamental ideas of FM radio.

To get a better idea of what is happening in this signal, we can zoom in on a smaller chunk:¶

plt.figure(figsize=(14,7))
plt.plot(time[10000:20000], xcopy[10000:20000])
plt.xlabel("Time (sec)")
plt.ylabel("Amplitude")
plt.title("Raw FM Signal - Time Domain")
plt.show()

From this portion of the signal, you can clearly see the frequency modulation. The wavelength decreases around t=0.006 seconds.

Next, let's take a look at the frequency domain:¶

N = samples
T = 1.0 / sr
Xjw = fftshift(fft(x))
xf = fftshift(fftfreq(N,T)) / 1000.0
plt.figure(figsize=(14,7))
plt.plot(xf, np.abs(Xjw))
plt.xlabel("Frequency (KHz) from baseband (90.9 MHz)")
plt.ylabel("Magnitude")
plt.xlim([-200,200])
plt.ylim([0,25000])
plt.title("Raw FM Signal - Frequency Domain")
plt.show()

As expected, we can see the "data" of the signal is centered at zero and spans about +/- 75 kHz. This looks very similar to the screenshot of the spectrum from my SDR software (as expected).

Demodulation¶

Now, I will attempt to demodulate the recorded FM signal. I will first shorten the length of the recording by some, to reduce the processing load, then decimate by a factor of 8 to capture the "data" of the signal. The decimation function also applies an anti-aliasing filter.

from math import ceil
h = ceil(len(x)/2)
xcopy = x
xcopy = xcopy[:h]
xcopy = decimate(xcopy,8,None,'fir')
sr=ceil(sr/8)

The below function demodulates the FM signal. It is based on the equation:

$f(n)=\frac{f_s}{2\pi}[\Delta\theta_{rad}(n)]$ Hz

where $f_s$ is the sampling frequency.

def fm_demod(xa, df=1.0, fc=0.0):
    ''' Perform FM demodulation of complex carrier.

    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(xa))
    rx = xa*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(df*np.unwrap(phi)/(2*np.pi))

    return y

Now, we will call the demodulation function and plot the result in the frequency domain.

xnew = fm_demod(xcopy,10000,125000)
N = xnew.shape[0]
T = 1.0 / sr
Xjw = fftshift(fft(xnew))
xf = fftshift(fftfreq(N,T)) / 1000.0
plt.figure(figsize=(14,7))
plt.plot(xf, np.abs(Xjw))
plt.xlabel("Frequency (KHz)")
plt.ylabel("Magnitude")
plt.title("Demodulated FM Signal - Frequency Domain")
plt.show()

Clearly this demodulated signal is vastly different from the original signal in the frequency domain. We have a collection of impulses with a small layer of noise across all frequencies.

Another interesting thing to note is that this signal is perfectly symmetric (and therefore real). When the demodulation is done, any phase variation (imaginary components) are eliminated, and the result is a real, symmetric signal in the frequency domain.

ipd.Image("fmspectrum.png")

Now let's see how the time domain compares to the raw signal. I will again decimate the signal to capture only the mono audio portion of the signal (see the below diagram which shows the frequency make-up of the demodulated signal).

xfinal = decimate(xnew,10,None,'fir')
sr=ceil(sr/10)
samples = xfinal.shape[0]
length = samples / sr
time = np.linspace(0., length, samples)
plt.figure(figsize=(14,7))
plt.plot(time,xfinal)
plt.title("Demodulated FM Signal - Time Domain")
plt.xlabel("Time (s)")
plt.ylabel("Amplitude")
plt.show()

Clearly we have made some change to the signal: the amplitude now varies like we would expect from a normal audio signal. The result is not perfect, but with more time the parameters of the demodulation function could be changed to produce a cleaner signal. There is also a limit to what can be done due to the bandwidth limitations of the SDR dongle I used as well as the processing power needed to demodulate a signal with millions of samples.

Listening to the below player (turn down the volume first), you can hear the original signal faintly: a classical music broadcast. There is quite a bit of noise, but the music is indeed present. With more time and better hardware, the noise could be reduced even further.

ipd.Audio(xfinal,rate=sr)

Finally, as a comparison, I have included a pre-demodulated signal from the same station and its frequency and time domain plots below. This signal was demodulated directly in the SDR software I used and clearly has much higher fidelity.

The frequency axis uses a log scale to more easily display the present frequencies. Another interesting note, the type of audio in this signal is human voice, specifically that of a female. This is clearly reflected in the frequency plot, with most of the amplitude present around 200-500 Hz.

y, sr_y = librosa.load('audio1.wav', sr=None)
samples2 = y.shape[0]
length2 = samples / sr_y
time2 = np.linspace(0., length2, samples2)
fig, (ax1,ax2) = plt.subplots(1,2,figsize=(14,7))
ax1.plot(time2,y)
ax1.set_xlabel("Time (s)")
ax1.set_ylabel("Amplitude")
ax1.set_title("Time Domain")
Yjw = fftshift(rfft(y))
yf = fftshift(fftfreq(samples2,1/sr_y))
ax2.plot(yf,Yjw)
ax2.set_xlabel("Frequency (Hz)")
ax2.set_ylabel("Amplitude")
ax2.set_title("Frequency Domain")
ax2.set_xscale("log")
fig.suptitle("Sample Pre-Demodulated Radio Signal")
plt.show()
ipd.Audio(y,rate=sr_y)

References¶

https://en.wikipedia.org/wiki/Frequency_modulation
https://www.embedded.com/dsp-tricks-frequency-demodulation-algorithms/
http://www.aaronscher.com/wireless_com_SDR/RTL_SDR_AM_spectrum_demod.html