Recently, I have posted a few blogs and projects which looked at signal processing or filters. I thought it would be a good idea to go into more detail about one of the commonly used elements of these blogs and projects which is the Fourier Transform.
We should be comfortable working in both the time and frequency domains, so hopefully this blog provides a good starting point for those who might not be.
Time or Frequency Domain
As engineers, we can analyze and manipulate signals in either time or frequency domains and knowing when to do either one is a core reason an engineer is required on a project.
The time domain enables engineers to analyze how a signal changes over time. Typically, in electronic systems, the signal in question is either a changing voltage, current, or frequency which has been output by a sensor or generated by another part of the system. Within this domain, we can measure a signal’s amplitude, frequency, and period in addition to more interesting parameters such as the rise and fall times of the signal. If we wish to observe a time domain signal in a laboratory environment, it’s common to use an oscilloscope or logic analyzer.
However, there are signal parameters that require analysis within the frequency domain to access the information contained within. These parameters are present within the frequency domain where we can identify the frequency components of the signal, their amplitudes, and the phase of each frequency. Working within the frequency domain also makes it much simpler to manipulate signals due to the ease with which convolution can be performed in the frequency domain. Similar to the time domain analysis, we must use a spectrum analyzer if we wish to observe a frequency domain signal within the laboratory environment.
For some applications, we will wish to work within the time domain like systems that monitor voltage or temperature, for example. While noise may be an issue, taking an average number of samples will be sufficient in many cases. For other applications, however, it is preferable to work within the frequency domain. For example, this could include signal processing applications which requires the filtering of one signal from another or to separate it from a noise source.
Of course, working within the time domain requires less processing on the quantized digital signal. Working within the frequency domain, on the other hand, first requires the application of a transform to the quantized data to convert from the time domain. Similarly, we need to perform the inverse transform again back to the time domain to output the post processed data from the frequency domain.
From the Time to Frequency Domain
Depending upon the type of signal (i.e., repetitive or non-repetitive, discrete or non-discrete), there are a number of methods that can be used to covert between time and frequency domains such as Fourier series, Fourier Transforms or Z Transforms. Within electronic signal processing and FPGA applications in particular, the most often used method is the Discrete Fourier Transform (DFT) which is a subset of the Fourier Transform. The DFT is used to analyze signals that are periodic and discrete meaning they consist of a number of N bit samples evenly spaced at a sampling frequency which, in many applications, are supplied by an ADC within the system.
Simply put, the DFT decomposes the input signal into two output signals which represent the sine and cosine components of that signal. Thus, for a time domain sequence of N samples, the DFT will return two groups of N/2+1 cosine and sine wave samples respectively referred to as the real and imaginary components as shown in figure one below. The real and imaginary sample width will also be N/2 for an input signal width of n bits.
The algorithm to calculate the DFT is pretty straight forward and we can use the equation shown below.
Where x[i] is the time domain signal, i ranges from 0 to N-1 and k ranges from 0 to N/2. The algorithm above is called the correlation method and it multiplies the input signal with the sine or cosine wave for that iteration to determine its amplitude.
At some point in our application, we’ll wish to transform back from the frequency domain into the time domain. To do this, we can use the synthesis equation which combines the real and imaginary waveforms to recreate a time domain signal.
ReX ̅ and ImX ̅ are, however, scaled versions of the cosine and sine waves so we need to scale them. ImX[k] or ReX[k] is divided by N/2 to determine the values for ReX ̅ and ImX ̅ in all cases except when ReX[0] and ReX[N/2], in this case, they are divided by N.
For obvious reasons this is called the Inverse Discrete Fourier Transform or IDFT for short.
We can use tools such as Python, MATLAB and even Excel to perform DFT calculations on captured data and many lab tools like oscilloscopes are capable of performing DFT on request.
Before we get any further, it’s worth pointing out that both the DFT and IDFT above are referred to as real DFT and real IDFT. This means that the input is a real number and not complex. It will soon become apparent why we need to know this.
We can demonstrate the DFT algorithm using Python in JupyterLab using the following code. We can also generate the plots below showing the real element using a signal of 10 MHz and sampling rate of 100 Msps.
def dft(signal):
"""
Compute the Discrete Fourier Transform (DFT) of the given signal.
"""
N = len(signal)
n = np.arange(N)
k = n.reshape((N, 1))
e = np.exp(-2j * np.pi * k * n / N)
return np.dot(e, signal)
# Compute the DFT
Y = dft(signal)
n = len(signal)
freq = np.arange(n) / (n*1/sampling_rate) # Frequency bins
# Plot the signal
plt.figure(figsize=(12, 6))
plt.subplot(2, 1, 1)
plt.plot(t, signal)
plt.title('Sample Signal')
plt.xlabel('Time')
plt.ylabel('Amplitude')
# Plot the magnitude spectrum
plt.subplot(2, 1, 2)
#plt.plot(freq, np.abs(Y))
plt.plot(freq[:n // 2], np.abs(Y)[:n // 2] * 1/n) # Plot only the positive frequencies
plt.title('Magnitude Spectrum')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude')
#plt.xlim([-2e7,2e7])
plt.show()
FPGA based Implementation
The implementation of the DFT and IDFT as described above is often implemented as a nested loop, each performing N calculations as shown in the Python example. As such, the time taken to implement the DFT calculation is as follows.
DFTtime=N*N*Kdft
Where Kdft is the processing time for each iteration to be performed, this can become quite time consuming to implement. It will be pretty fast on a high-performance desktop but much slower on a processor intended for the edge.
As a result, DFTs within FPGAs are normally implemented using an algorithm called the Fast Fourier Transform (FFT) to calculate the DFT. This algorithm has often been called the “the most important algorithm of our lifetime” because of its enabling impact on many industries. The FFT differs slightly from the DFT algorithms explained previously in that it calculates the complex DFT. This means that it expects real and imaginary time domain signals and produces results in the frequency domain which are n bits wide as opposed to n/2. To clarify further, when we wish to calculate a real DFT, we must first set the imaginary part to zero and move the time domain signal into the real part. We have two options if we wish to implement a FFT within our AMD FPGA. We can either write one from scratch using the HDL of our choice or use the FFT IP provided within the Vivado IP catalogue. Unless there are pressing reasons not to use the IP, I suggest using the IP core to allow for reduced development time.
The basic approach of the FFT is that it decomposes the time domain signal into a number of single point time domain signals. This is often called bit reversal because the samples are reordered. If a bit reversal algorithm is not used as a short cut, the number of stages that it takes to create these single point time domain signals is calculated by Log2 N where N is the number of bits.
These single point time domain signals are then used to calculate the frequency spectra for each point. This is pretty straight forward since the frequency spectra is equal to the single point time domain.
It’s in the recombination of these single frequency points that the FFT algorithm gets complicated and we must recombine these spectra points one stage at a time. This is the opposite of the time domain decomposition and will again take Log2 N stages to recreate the spectra. This is where the famous FFT butterfly comes into play.
When compared with the DFT execution time, the FFT takes
FFTtime=Kfft*N Log2 N
This results in significant improvements in execution time to calculate a DFT.
When implementing an FFT within our FPGA, we must also take into account the size of the FFT because this will determine the noise floor below which we can’t see signals of potential interest. The FFT size will also determine the spacing of the frequency bins. The equation below is used for this.
FFTNoise Floor (dB)=6.02n+1.77+10log10(FFTsize/2)
Where n is the number of quantized bits within the time domain and FFT size is the FFT size. For FPGA-based implementations, this is normally a power of 2 (e.g., 256, 512, 1024 etc.). The frequency bins will be evenly spaced at
Bin Width=(Fs/2)/FFTsize
To illustrate a very simple example, a FS of 100 MHz with a FFT size of 128 would have a frequency resolution of 0.39 Hz. This means frequencies with 0.39 Hz of each other cannot be distinguished.
Again, we can look at the FFT in JupyterLabs and observe the output for the same configuration as for the DFT.
# Compute the FFT
n = len(signal) # Length of the signal
freq = np.fft.fftfreq(n, d=1/sampling_rate) # Frequency bins
Y = np.fft.fft(signal)/n # FFT and normalization
# Plot the signal
plt.figure(figsize=(12, 6))
plt.subplot(2, 1, 1)
plt.plot(t, signal)
plt.title('Sample Signal')
plt.xlabel('Time')
plt.ylabel('Amplitude')
# Plot the magnitude spectrum
plt.subplot(2, 1, 2)
plt.plot(freq, np.abs(Y))
plt.title('Magnitude Spectrum')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude')
plt.xlim([-2e7,2e7])
plt.show()
Of course, now that we understand a little more about FFTs and DFT we might want to implement them in our FPGA. If you have a PYNQ board, you can find detailed step-by-step instructions on how to do so for MPSoC or Zynq 7000 boards below.
This should help you understand how to go about implementing an FFT in a FPGA and give you good starting point to working with them further.
Workshops and Webinars
If you enjoyed the blog why not take a look at the free webinars, workshops and training courses we have created over the years. Highlights include
Professional PYNQ Learn how to use PYNQ in your developments
Introduction to Vivado learn how to use AMD Vivado
Ultra96, MiniZed & ZU1 three day course looking at HW, SW and Petalinux
Arty Z7-20 Class looking at HW, SW and Petalinux
Mastering MicroBlaze learn how to create MicroBlaze solutions
HLS Hero Workshop learn how to create High Level Synthesis based solutions
Perfecting Petalinux learn how to create and work with PetaLinux OS
Embedded System Book
Do you want to know more about designing embedded systems from scratch? Check out our book on creating embedded systems. This book will walk you through all the stages of requirements, architecture, component selection, schematics, layout, and FPGA / software design. We designed and manufactured the board at the heart of the book! The schematics and layout are available in Altium here Learn more about the board (see previous blogs on Bring up, DDR validation, USB, Sensors) and view the schematics here.
Sponsored by AMD
OOPS - After my first post I couldn't see a way to view the rest of the post, which had the signal function. I assumed that the post length was limited and reposted the signal function and plot format statement. Now I know to refresh following a new post.
Nice article, Adam. My last DSP activity was my DSP class in college. So I decided to play with your code a bit.
I discovered that the code doesn't run as published. It needs the following:
A signal needs to be defined
Library declarations for matplotlib, numpy
Here is the missing code, which will allow both examples to run. As I'm a python hack and haven't done much python coding, the code is likely not optimal. But it works :)
Not sure how to insert code, but trying ``` markdown. Don't add these backticks if if they show up in this post.
Library declarations
```
import math
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
```
Signal generator