Ever seen those beautiful images Real-Time Spectrum Analyzers produce that show the FFT as an intensity graded persistence display? Well, I wanted to create a plotting tool that would allow you to create high-resolution static imagery for signals.
There is already an amazing tool aptly named Fosphor for GNURadio that is capable of real-time plotting of this intensity graded display utilizing the GPU processing but I could not create a high-resolution static output using it. Below is an example screen capture of the Fosphor display from the Osmocom website.
So I decided to create one in Python that could be customized to create static high-resolution plots that could be used in presentations or poster sessions. The added benefit being able to create artwork based on these beautiful colored images!
Background
Intensity-graded spectral plots are more commonly seen in waterfall plots. These display a 2D image that shows the frequency domain information relative to the time axis, with frequency magnitude displayed as a color z-axis. There are already many spectrogram functions in python that allow you to display these plots. The python module Matplotlib.pyplot provides the specgram() method which takes a signal as an input and plots the spectrogram.
Successive FFTs are computed each of which corresponds to one unit of the time axis. These are then stacked up to show how the frequency spectrum changes over time. The higher the magnitude of the FFT bin corresponds to a brighter or hotter color.
Digital Phosphor FFT
A lesser-known spectral intensity graded display are those used by really fast Real-Time Spectrum Analyzers to show an FFT plot that changes with time but which overlays multiple FFT computations overtop of each other. Here each FFT is mapped onto a grid, each grid hit from the FFT is summed with the previous hits.
The image is then colored based on this grid of hit points where brighter or hotter colors correspond to more FFT points that fall in the same grid location. This is then combined with a decaying factor to show how persistent a signal is in your band of interest.
Python Function
We only need a couple of python libraries to create these plots, make sure you have numpy, matplotlib, and pyrtlsdr if you have an RTLSDR to capture real RF samples. You can install these with pip.
pip install numpy
pip install matplotlib
pip install pyrtlsdrBasics
In order to compute an FFT in python, we can utilize the wonderful Numpy library, but it’s always a good idea to learn the basics of frequency-domain processing. For that Dr. Marc Lichtman has created PySDR, a wonderful resource for using python processing with software-defined radios. For this project, the section on Frequency Domain Python Processing was particularly useful.
I used the powerful matplotlib python library for plotting the heatmap data into a color-graded representation. A general resource can be found by Oreilly at Python Data Science Plotting
I will be using the imshow() function of matplotlib to plot the heatmap data. One of the most important parameters to configure when plotting heat maps is the color map. The color map assigns the different magnitudes in an array to various colors in the plot. matplotlib has a variety of options to choose from which can be viewed at Matplotlib color maps.
Collect or Generate Data
To build the function I needed to have some RF sample data to work with. I used an RTLSDR to collect a set of samples at the 915MHz center frequency as that band had some transient signals that would be a good test of the temporal nature of the intensity graded FFT processing. Thankfully the is already a python library, pyrtlsdr, that allows me to easily capture some samples with a few lines of python code.
from rtlsdr import RtlSdr
sdr = RtlSdr()
# configure device
sdr.sample_rate = 2.048e6 # Hz
sdr.center_freq = 915e6 # Hz
sdr.gain = 'auto'
samples = sdr.read_samples(256*10240)
sdr.close()To cover my bases I also wanted to synthetically generate some data to test with so those without an RTLSDR to capture data could still run the function to test their environment. The generated samples were created based on PySDR QPSK signal that has been pulse shaped to reduce the occupied bandwidth.
Generating this signal is a bit more involved than the RTLSDR capture as we need to generate random bits, modulate them into QPSK, design and apply the filter, and lastly apply some noise so we can see some variation in the spectrum on subsequent FFT calculations.
#Test QPSK Signal
num_symbols = 256*10240
sps = 2
x_int = np.random.randint(0, 4, num_symbols) # 0 to 3
x_int = np.repeat(x_int, sps, axis=0)
x_degrees = x_int*360/4.0 + 45 # 45, 135, 225, 315 degrees
x_radians = x_degrees*np.pi/180.0 # sin() and cos() takes in radians
x_symbols = np.cos(x_radians) + 1j*np.sin(x_radians) # this produces our QPSK complex symbols
# Create our raised-cosine filter
num_taps = 101
beta = 0.35
Ts = sps # Assume sample rate is 1 Hz, so sample period is 1, so *symbol* period is 8
t = np.arange(-51, 52) # remember it's not inclusive of final number
h = np.sinc(t/Ts) * np.cos(np.pi*beta*t/Ts) / (1 - (2*beta*t/Ts)**2)
# Filter our signal, in order to apply the pulse shaping
x_shaped = np.convolve(x_symbols, h)
n = (np.random.randn(len(x_shaped)) + 1j*np.random.randn(len(x_shaped)))/np.sqrt(2) # AWGN with unity power
noise_power = 0.001
r = x_shaped + n * np.sqrt(noise_power)
samples = rStatic Function Output
The static intensity graded FFT function takes a NumPy array of samples and splits them in time based on the length of your desired FFT. These FFTs are then processed into a hit map grid one on top of the other to build up the pixel intensity. So the final plot is a superposition of all the FFTs layered on top of one another so you are able to see the spectrum usage for a given length of sample time.
#Function
def fft_intensity_plot(samples: np.ndarray, fft_len: int = 256, fft_div: int = 2, mag_steps: int = 100, cmap: str = 'plasma'):
num_ffts = math.floor(len(samples)/fft_len)
fft_array = []
for i in range(num_ffts):
temp = np.fft.fftshift(np.fft.fft(samples[i*fft_len:(i+1)*fft_len]))
temp_mag = 20.0 * np.log10(np.abs(temp))
fft_array.append(temp_mag)
max_mag = np.amax(fft_array)
min_mag = np.abs(np.amin(fft_array))
norm_fft_array = fft_array
for i in range(num_ffts):
norm_fft_array[i] = (fft_array[i]+(min_mag))/(max_mag+(min_mag))
mag_step = 1/mag_steps
hitmap_array = np.random.random((mag_steps+1,int(fft_len/fft_div)))*np.exp(-10)
for i in range(num_ffts):
for m in range(fft_len):
hit_mag = int(norm_fft_array[i][m]/mag_step)
hitmap_array[hit_mag][int(m/fft_div)] = hitmap_array[hit_mag][int(m/fft_div)] + 1
hitmap_array_db = 20.0 * np.log10(hitmap_array+1)
figure, axes = plt.subplots()
axes.imshow(hitmap_array_db, origin='lower', cmap=cmap, interpolation='bilinear')
return(figure)Animated Function
After creating the static function I decided to experiment with the animation capabilities in matplotlib to see if I could create something similar to the GNURadio Fosphor module but without GNURadio complexity. This would also allow me to create video files that could be embedded into presentations.
The animation function takes in the same NumPy array of samples for a given time and split them into segments of FFT length. Now the difference between this function and the static function is that the FFTs get played on top of each other as frames of the video. The whole hit map array has a decay function applied to it so the transient frequency domain points appear and disappear based on how persistent they are in the spectrum.
The animation function calls a separate function that performs the hit map calculation and delay so that the frames can play in sequence.
#Animate function that interates over fft data
def animate(i, im, norm_fft_array, mag_steps, fft_len, fft_div):
mag_step = 1/mag_steps
if i == 0:
hitmap_array = im.get_array()*np.exp(-10)
else:
hitmap_array = im.get_array()*np.exp(-0.04)
for m in range(fft_len):
hit_mag = int(norm_fft_array[i][m]/mag_step)
hitmap_array[hit_mag][int(m/fft_div)] = hitmap_array[hit_mag][int(m/fft_div)] + .5
#hitmap_array_db = 10.0 * np.log10(hitmap_array)
im.set_array(hitmap_array)
return [im]
#Function
def fft_intensity_animate(samples: np.ndarray, fft_len: int = 256, fft_div: int = 2, mag_steps: int = 100):
num_ffts = math.floor(len(samples)/fft_len)
fft_array = []
for i in range(num_ffts):
temp = np.fft.fftshift(np.fft.fft(samples[i*fft_len:(i+1)*fft_len]))
temp_mag = 20.0 * np.log10(np.abs(temp))
fft_array.append(temp_mag)
max_mag = np.amax(fft_array)
min_mag = np.abs(np.amin(fft_array))
norm_fft_array = fft_array
for i in range(num_ffts):
norm_fft_array[i] = (fft_array[i]+min_mag)/(max_mag+min_mag)
#animation place holder
fig = plt.figure()
a = np.random.random((mag_steps+1,int(fft_len/fft_div)))
im = plt.imshow(a, origin='lower', cmap='plasma', interpolation='bilinear')
#compute animation
anim = FuncAnimation(fig, animate, frames=1000, fargs = (im,norm_fft_array,mag_steps,fft_len,fft_div), interval=1, blit=True)
return(anim)Performance
These are not optimized functions so the plotting can take a significant amount of time to produce a graph or animation. If we want a real-time solution we need to use another framework to plot the data in a GUI-like environment. Coming Soon!
Jupyter Notebooks
I created Jupyter Notebooks that demonstrate these functions in a GitHub repo. For each of the static and animated functions, I suppled a notebook that operated on a generated QPSK signal and one that uses an RTLSD to capture samples. You can clone the repo at the link below.
I also linked Github Gists of the notebooks below so you can get a quick idea of the function structure.
Static RTLSDR Notebook
Static QPSK Notebook
Animated RTLSDR Notebook
Animated QPSK Notebook
Intensity Graded FFT Desk Mat for your Home Office
I took some of the cool plots that my code created and placed them on a desk mat material using the RedBubble service. Check them out if you would like one for your own work-from-home setup! Ordering one would also help support this blog and allow me to create more open-source content.
1 thought on “Python Intensity Graded FFT Plots”
Comments are closed.