Category: Data Acquisition

Circular Buffer in Python

A circular or ring buffer is a fixed-size data structure that is commonly used in real-time software applications to store a pre-defined number of values. The analogy of a ring with a fixed number of positions is quite useful to capture the FIFO (First-In-First-Out) nature of such data structure. Once the buffer is full, the first element that was written (“In” ) to the buffer is the first one to be overwritten (“Out”) by the next incoming element.

Circular buffers are particularly useful in situations where the data is continuously being sampled and calculations need to be done using a pre-defined sample size or continuous visualization is required.

The two images below represent a circular buffer with 10 positions. On the left, with four elements written to it being 5 the “first in”. On the right, the buffer is full where the value 3 was written. Note how the current write position pointer moves around as new values are added.

Additionally, a ring buffer has basically two states: not-yet- full (left image) and full (right image). The not-yet-full state occurs after the buffer is initialized and exists only until it becomes full. We will see later how this affects the Python code used to represent the buffer.

Once the buffer state changes to full, the FIFO nature of this type of data structure becomes evident. In the next two images, the next element “in” is the value 12, replacing the “first in” value 5, which therefore is the “first out”. The next value “in” is 2, replacing the “second in” value 10, and so on.

In the example figures above, we could imagine the buffer being used to display the average value of its elements (once the buffer reaches the full state). As time progresses, a new sample is written to the buffer, replacing the oldest sample. Each figure of the full buffer is a snapshot at a given time step of the sampling process. In this case the average values are 6.4, 7.1, and 6.3.

Python Implementation

The most elegant way to implement the circular buffer is by using a class. Also, a true circular buffer does not involve shifting its elements to implement the FIFO structure. The current position pointer is used to keep track of where the newest element is.

The circular buffer class should have the following attributes:

bufsize: the number of elements that the circular buffer can hold (integer)
data: the list of elements, i.e., the buffer itself
currpos: the current position pointer (integer)

The class should also have two methods:

add(): writes an element to the buffer
get(): returns a list containing the buffer elements

Let’s go over some of the features of the Python implementation below:

As mentioned earlier, the ring buffer has two states (not-yet-full and full). The first state only exists for a limited time until the buffer is filled up. With that in mind, it makes sense to define the __Full class within the RingBuffer class with the exact same methods. Once the buffer is full, the original class is permanently changed to the full class implementation, with its add() and get() methods now superseding the original ones.
The not-yet-full add() method just does a regular list append operation of new elements to the buffer list. It is also in this method that the class is changed to the __Full class once the buffer reaches its capacity.
The full add() method, on the other hand, writes the new element at the current (newest) element position and increments the current position by one unit (wrapping it around, or resetting it, once the buffer size value is reached).
More interesting however is the full get() method. It splits the data list at the current position into two lists which are then concatenated. The list going from currpos to the last element goes first, then comes the list from the first element to currpos (excluded). By returning this concatenated list, the method fully implements the circular buffer without shifting any of it elements! Inspect the data attribute as you add elements to a full buffer to see what the original list looks like.

class RingBuffer:
    """ Class that implements a not-yet-full buffer. """
    def __init__(self, bufsize):
        self.bufsize = bufsize
        self.data = []

    class __Full:
        """ Class that implements a full buffer. """
        def add(self, x):
            """ Add an element overwriting the oldest one. """
            self.data[self.currpos] = x
            self.currpos = (self.currpos+1) % self.bufsize
        def get(self):
            """ Return list of elements in correct order. """
            return self.data[self.currpos:]+self.data[:self.currpos]

    def add(self,x):
        """ Add an element at the end of the buffer"""
        self.data.append(x)
        if len(self.data) == self.bufsize:
            # Initializing current position attribute
            self.currpos = 0
            # Permanently change self's class from not-yet-full to full
            self.__class__ = self.__Full

    def get(self):
        """ Return a list of elements from the oldest to the newest. """
        return self.data


# Sample usage to recreate example figure values
import numpy as np
if __name__ == '__main__':

    # Creating ring buffer
    x = RingBuffer(10)
    # Adding first 4 elements
    x.add(5); x.add(10); x.add(4); x.add(7)
    # Displaying class info and buffer data
    print(x.__class__, x.get())

    # Creating fictitious sampling data list
    data = [1, 11, 6, 8, 9, 3, 12, 2]

    # Adding elements until buffer is full
    for value in data[:6]:
        x.add(value)
    # Displaying class info and buffer data
    print(x.__class__, x.get())

    # Adding data simulating a data acquisition scenario
    print('')
    print('Mean value = {:0.1f}   |  '.format(np.mean(x.get())), x.get())
    for value in data[6:]:
        x.add(value)
        print('Mean value = {:0.1f}   |  '.format(np.mean(x.get())), x.get())

The output shown next is produced if the code containing the class is executed. The values and the states of the ring buffer shown in the figures at the beginning of the post are recreated in an iterative process that simulates some data acquisition. Note how the ring buffer class changes once the buffer is full.

    <class '__main__.RingBuffer'> [5, 10, 4, 7]
    <class '__main__.RingBuffer.__Full'> [5, 10, 4, 7, 1, 11, 6, 8, 9, 3]

    Mean value = 6.4   |   [5, 10, 4, 7, 1, 11, 6, 8, 9, 3]
    Mean value = 7.1   |   [10, 4, 7, 1, 11, 6, 8, 9, 3, 12]
    Mean value = 6.3   |   [4, 7, 1, 11, 6, 8, 9, 3, 12, 2]

As a final remark, in the Pulse Rate Monitor post, the buffer that was used is not a true circular buffer, since the elements in the arrays are shifted by one position for every new value being sampled. Because the buffer is relatively small (200 elements) and the sampling period of the data is reasonably large (0.1s) this solution does not constitute a problem. You can check out the code that implements a true ring buffer for the monitor on my GitHub page.

Pulse Rate Monitor with Raspberry Pi

Pulse rate can be detected using a variety of physical principles. One way of doing it is by applying infrared light to one side of the tip of your finger and sensing the “output” on the other side with a phototransistor. The blood flow through the capillaries at the finger tip fluctuates with the same frequency (rate) of your pulse. That fluctuation can be detected after processing the signal captured by the sensor.

Keyestudio has a rudimentary pulse rate monitor (KS0015) that does the trick. In this post we will go over the signal processing and how to display the pulse rate in real time. Since the KS0015 outputs an analog signal, it has to be used with the MCP3008 analog-to-digital converter chip, as illustrated in the schematics below.

As noted in previous posts, adopting a reference voltage of 3.3V for the MCP3008 provides the best resolution for the analog-to-digital conversion. The sensor output therefore has to be no greater than 3.3V, which can be achieved by using 3.3 instead of 5V as the supply voltage for the KS0015.

The very first thing to do is to look at the signal coming from the sensor as the tip of my index finger is placed between the IR LED and the phototransistor. The plot on the left can be generated by running the Python code below.

By inspecting the program, the first thing to observe is the sampling period (tsample) of 0.02 s. Since we don’t know much about the signal, it’s always good to start with a smaller value. The second thing is the wait time of 5 seconds before collecting any data. After tinkering around for a bit, I realized that the sensor is very sensitive to finger movements. Therefore, waiting for the signal to “settle down” really helps.

# Importing modules and classes
import time
import numpy as np
import matplotlib.pyplot as plt
from gpiozero import MCP3008


# Defining function for plotting
def make_fig():
    # Creating Matplotlib figure
    fig = plt.figure(
        figsize=(8, 3),
        facecolor='#f8f8f8',
        tight_layout=True)
    # Adding and configuring axes
    ax = fig.add_subplot(xlim=(0, max(t)))
    ax.set_xlabel('Time (s)', fontsize=12)
    ax.grid(linestyle=':')
    # Returning axes handle
    return ax


# Creating object for the pulse rate sensor output
vch = MCP3008(channel=0, clock_pin=11, mosi_pin=10, miso_pin=9, select_pin=8)

# Assigning some parameters
tsample = 0.02  # Sampling period for code execution (s)
tstop = 30  # Total execution time (s)
vref = 3.3  # Reference voltage for MCP3008
# Preallocating output arrays for plotting
t = []  # Time (s)
v = []  # Sensor output voltage (V)

# Waiting for 5 seconds for signal stabilization
time.sleep(5)

# Initializing variables and starting main clock
tprev = 0
tcurr = 0
tstart = time.perf_counter()
# Running execution loop
print('Running code for', tstop, 'seconds ...')
while tcurr <= tstop:
    # Getting current time (s)
    tcurr = time.perf_counter() - tstart
    # Doing I/O and computations every `tsample` seconds
    if (np.floor(tcurr/tsample) - np.floor(tprev/tsample)) == 1:
        # Getting current sensor voltage
        valuecurr = vref * vch.value
        # Assigning current values to output arrays
        t.append(tcurr)
        v.append(valuecurr)
    # Updating previous time value
    tprev = tcurr

print('Done.')
# Releasing GPIO pins
vch.close()

# Plotting results
ax = make_fig()
ax.set_ylabel('Sensor Output (V)', fontsize=12)
ax.plot(t, v, linewidth=1.5, color='#1f77b4')
ax = make_fig()
ax.set_ylabel('Sampling Period (ms)', fontsize=12)
ax.plot(t[1::], 1000*np.diff(t), linewidth=1.5, color='#1f77b4')

Digital Filtering and Pulse Detection

In this stage of the signal processing, let’s put some of the content of previous posts to good use. More specifically digital band-pass filtering and event detection for the pulse signal.

First let’s remove the low and high frequency components of the signal that are outside the band of interest. I chose 0.5 and 5 Hz (30 and 300 beats per minute) as the cutoff frequencies. That seems to be a reasonable bandwidth to account for the filter attenuation at the cutoff frequencies, as well as for the human heart “range of operation”.

Then, let’s calculate the derivative of the signal so it’s easier to detect the pulse events. By normalizing the derivative based on its maximum value, the signal amplitude becomes more consistent and therefore a single value threshold can be used. In this case 25% of the maximum value.

The plot on the top shows the normalized derivative of the filtered signal. In red are the trigger events detected by the function find_cluster, which I placed inside the module utils.py . As usual, all the code for this post can be found on my GitHub page.

The plot on the bottom shows the identified trigger events (red dots) displayed on the filtered sensor output signal.

The Python program below is an extension of the one in the previous section, where the digital filter and the trigger detection features are now incorporated. The pulse rate is calculated using the median of the time difference between two consecutive trigger events. In this case, using the median instead of the mean is recommended since it makes the calculation more robust to outliers.

Running the code will generate the two graphs shown above and print the median pulse rate for the corresponding (30-second) data acquisition window. You will notice that the signal amplitude is fairly small (approx. 30 mV) and very sensitive to the finger tip location. It might take a few tries to find the correct position and pressure for the finger.

# Importing modules and classes
import time
import numpy as np
import matplotlib.pyplot as plt
from gpiozero import MCP3008
from utils import find_cluster


# Defining function for plotting
def make_fig():
    # Creating Matplotlib figure
    fig = plt.figure(
        figsize=(8, 3),
        facecolor='#f8f8f8',
        tight_layout=True)
    # Adding and configuring axes
    ax = fig.add_subplot(xlim=(0, max(t)))
    ax.set_xlabel('Time (s)', fontsize=12)
    ax.grid(linestyle=':')
    # Returning axes handle
    return ax


# Creating object for the pulse rate sensor output
vch = MCP3008(channel=0, clock_pin=11, mosi_pin=10, miso_pin=9, select_pin=8)

# Assigning some parameters
tsample = 0.1  # Sampling period for code execution (s)
tstop = 30  # Total execution time (s)
vref = 3.3  # Reference voltage for MCP3008
# Preallocating output arrays for plotting
t = np.array([])  # Time (s)
v = np.array([])  # Sensor output voltage (V)
vfilt = np.array([])  # Filtered sensor output voltage (V)

# Waiting for 5 seconds for signal stabilization
time.sleep(5)

# First order digital band-pass filter parameters
fc = np.array([0.5, 5])  # Filter cutoff frequencies (Hz)
tau = 1/(2*np.pi*fc)  # Filter time constants (s)
# Filter difference equation coefficients
a0 = tau[0]*tau[1]+(tau[0]+tau[1])*tsample+tsample**2
a1 = -(2*tau[0]*tau[1]+(tau[0]+tau[1])*tsample)
a2 = tau[0]*tau[1]
b0 = tau[0]*tsample
b1 = -tau[0]*tsample
# Assigning normalized coefficients
a = np.array([1, a1/a0, a2/a0])
b = np.array([b0/a0, b1/a0])
# Initializing filter values
x = [vref*vch.value] * len(b)  # x[n], x[n-1]
y = [0] * len(a)  # y[n], y[n-1], y[n-2]
time.sleep(tsample)

# Initializing variables and starting main clock
tprev = 0
tcurr = 0
tstart = time.perf_counter()
# Running execution loop
print('Running code for', tstop, 'seconds ...')
while tcurr <= tstop:
    # Getting current time (s)
    tcurr = time.perf_counter() - tstart
    # Doing I/O and computations every `tsample` seconds
    if (np.floor(tcurr/tsample) - np.floor(tprev/tsample)) == 1:
        # Getting current sensor voltage
        valuecurr = vref * vch.value
        # Assigning sensor voltage output to input signal array
        x[0] = valuecurr
        # Filtering signals
        y[0] = -np.sum(a[1::]*y[1::]) + np.sum(b*x)
        # Updating output arrays
        t = np.concatenate((t, [tcurr]))
        v = np.concatenate((v, [x[0]]))
        vfilt = np.concatenate((vfilt, [y[0]]))
        # Updating previous filter output values
        for i in range(len(a)-1, 0, -1):
            y[i] = y[i-1]
        # Updating previous filter input values
        for i in range(len(b)-1, 0, -1):
            x[i] = x[i-1]
    # Updating previous time value
    tprev = tcurr

print('Done.')
# Releasing pins
vch.close()

# Calculating and normalizing sensor signal derivative
dvdt = np.gradient(vfilt, t)
dvdt = dvdt/np.max(dvdt)
# Finding heart rate trigger event times
icl, ncl = find_cluster(dvdt>0.25, 1)
ttrigger = t[icl]
# Calculating heart rate (bpm)
bpm = 60/np.median(np.diff(ttrigger))
print("Heart rate = {:0.0f} bpm".format(bpm))

# Plotting results
ax = make_fig()
ax.set_ylabel('Normalized Derivative ( - )', fontsize=12)
ax.plot(t, dvdt, linewidth=1.5, color='#1f77b4', zorder=0)
for ik, nk in zip(icl, ncl):
    if ik+nk < len(t):
        ax.plot(t[ik:ik+nk], dvdt[ik:ik+nk], color='#aa0000')
ax = make_fig()
ax.set_ylabel('Filtered Output (V)', fontsize=12)
ax.plot(t, vfilt, linewidth=1.5, color='#1f77b4', zorder=0)
for ik, nk in zip(icl, ncl):
    ax.plot(t[ik:ik+nk], vfilt[ik:ik+nk], color='#aa0000')
    ax.scatter(t[ik+1], vfilt[ik+1], s=25, c='#aa0000')

Real-time Pulse Rate Detection

Now that we have the pulse rate detection and calculation working, we can modify the code from the previous section to add the real-time display of the pulse rate. For this step, we will use a 7-segment LED display with a TM1637 chip. In case you don’t have one, all the calls for the tm object in the code below should be removed and a simple

print("Heart rate = {:0.0f} bpm".format(bpm))

can be used to replace the code line

tm.number(int(bpm))

Also, notice how we are using a circular buffer (tbuffer) of 20 seconds that holds the signal for the pulse rate calculation. The buffer is updated every tsample seconds and the actual pulse rate calculation occurs every tdisp seconds. Through comparison with the previous program, you can see how the bpm calculation is now placed inside the execution loop, hence happening in real-time.

# Importing modules and classes
import time
import numpy as np
import matplotlib.pyplot as plt
import tm1637
from gpiozero import MCP3008
from utils import find_cluster


# Defining function for plotting
def make_fig():
    # Creating Matplotlib figure
    fig = plt.figure(
        figsize=(8, 3),
        facecolor='#f8f8f8',
        tight_layout=True)
    # Adding and configuring axes
    ax = fig.add_subplot(xlim=(0, max(t)))
    ax.set_xlabel('Time (s)', fontsize=12)
    ax.grid(linestyle=':')
    # Returning axes handle
    return ax


# Creating object for the pulse rate sensor output and LED display
vch = MCP3008(channel=0, clock_pin=11, mosi_pin=10, miso_pin=9, select_pin=8)
tm = tm1637.TM1637(clk=18, dio=17)

# Assigning some parameters
tsample = 0.1  # Sampling period for code execution (s)
tstop = 60  # Total execution time (s)
tbuffer = 20  # Data buffer length(s)
tdisp = 1  # Display update period (s)
vref = 3.3  # Reference voltage for MCP3008
# Preallocating circular buffer arrays for real-time processing
t = np.array([])  # Time (s)
v = np.array([])  # Sensor output voltage (V)
vfilt = np.array([])  # Filtered sensor output voltage (V)

# Waiting for 5 seconds for signal stabilization
time.sleep(5)

# First order digital band-pass filter parameters
fc = np.array([0.5, 5])  # Filter cutoff frequencies (Hz)
tau = 1/(2*np.pi*fc)  # Filter time constants (s)
# Filter difference equation coefficients
a0 = tau[0]*tau[1]+(tau[0]+tau[1])*tsample+tsample**2
a1 = -(2*tau[0]*tau[1]+(tau[0]+tau[1])*tsample)
a2 = tau[0]*tau[1]
b0 = tau[0]*tsample
b1 = -tau[0]*tsample
# Assigning normalized coefficients
a = np.array([1, a1/a0, a2/a0])
b = np.array([b0/a0, b1/a0])
# Initializing filter values
x = [vref*vch.value] * len(b)  # x[n], x[n-1]
y = [0] * len(a)  # y[n], y[n-1], y[n-2]
time.sleep(tsample)

# Initializing variables and starting main clock
tprev = 0
tcurr = 0
tstart = time.perf_counter()
# Running execution loop
tm.show('Hold')
print('Running code for', tstop, 'seconds ...')
print('Waiting', tbuffer, 'seconds for buffer fill.')
while tcurr <= tstop:
    # Getting current time (s)
    tcurr = time.perf_counter() - tstart
    # Doing I/O and computations every `tsample` seconds
    if (np.floor(tcurr/tsample) - np.floor(tprev/tsample)) == 1:
        # Getting current sensor voltage
        valuecurr = vref * vch.value
        # Assigning sensor voltage output to input signal array
        x[0] = valuecurr
        # Filtering signals
        y[0] = -np.sum(a[1::]*y[1::]) + np.sum(b*x)
        # Updating circular buffer arrays
        if len(t) == tbuffer/tsample:
            t = t[1::]
            v = v[1::]
            vfilt = vfilt[1::]
        t = np.concatenate((t, [tcurr]))
        v = np.concatenate((v, [x[0]]))
        vfilt = np.concatenate((vfilt, [y[0]]))
        # Updating previous filter output values
        for i in range(len(a)-1, 0, -1):
            y[i] = y[i-1]
        # Updating previous filter input values
        for i in range(len(b)-1, 0, -1):
            x[i] = x[i-1]
    # Processing signal in the buffer every `tdisp` seconds
    if ((np.floor(tcurr/tdisp) - np.floor(tprev/tdisp)) == 1) & (tcurr > tbuffer):
        # Calculating and normalizing sensor signal derivative
        dvdt = np.gradient(vfilt, t)
        dvdt = dvdt/np.max(dvdt)
        # Finding heart rate trigger event times
        icl, ncl = find_cluster(dvdt>0.25, 1)
        ttrigger = t[icl]
        # Calculating and displaying pulse rate (bpm)
        bpm = 60/np.median(np.diff(ttrigger))
        tm.number(int(bpm))
    # Updating previous time value
    tprev = tcurr

print('Done.')
tm.show('Done')
# Releasing pins
vch.close()

Joystick with Raspberry Pi

A two-axis joystick is an input device that can be used to simultaneously control two degrees of freedom of a system, such as roll and pitch on an aircraft, or the X and Y coordinates of a cartesian robot. The Keyestudio KS0008 joystick discussed in this post provides analog signals for the two axis (left-to-right and front-to-back), as well as a digital signal (downward). Like any analog input device used with the Raspberry Pi, the KS0008 requires the MCP3008 analog-to-digital converter chip.

The wiring schematic shows how to connect the joystick to the Pi with the aid of the MCP3008. The analog outputs X and Y are wired to the channels 0 and 1 of the chip. The digital output B can be connected directly to any GPIO pin. In this configuration the KS0008 uses the 3.3V supply voltage. The Python code below samples and displays the three signals in an execution loop.

As a side note, I used the suffix LR (left-to-right) for the joystick’s X axis, and FB (front-to-back) for the Y axis. Mostly to avoid confusion with the digital filter input (X) and output (Y) used later in this post.

# Importing modules and classes
import time
import numpy as np
from gpiozero import MCP3008, DigitalInputDevice

# Creating objects for the joystick outputs
joyLR = MCP3008(channel=0, clock_pin=11, mosi_pin=10, miso_pin=9, select_pin=8)
joyFB = MCP3008(channel=1, clock_pin=11, mosi_pin=10, miso_pin=9, select_pin=8)
joyB = DigitalInputDevice(18)
# Assigning some parameters
tsample = 0.02  # Sampling period for code execution (s)
tdisp = 0.5  # Output display period (s)
tstop = 30  # Total execution time (s)
vref = 3.3  # Reference voltage for MCP3008

# Initializing variables and starting main clock
tprev = 0
tcurr = 0
tstart = time.perf_counter()
# Running execution loop
print('Running code for', tstop, 'seconds ...')
while tcurr <= tstop:
    # Getting current time (s)
    tcurr = time.perf_counter() - tstart
    # Doing I/O and computations every `tsample` seconds
    if (np.floor(tcurr/tsample) - np.floor(tprev/tsample)) == 1:
        # Getting joy stick normalized voltage output
        valLRcurr = joyLR.value
        valFBcurr = joyFB.value
        # Calculating current time raw voltages
        vLRcurr = vref*valLRcurr
        vFBcurr = vref*valFBcurr
        # Getting the Z axis state
        Bcurr = joyB.value
        # Displaying output voltages every `tdisp` seconds
        if (np.floor(tcurr/tdisp) - np.floor(tprev/tdisp)) == 1:
            print("X = {:0.2f} V , Y = {:0.2f} V , B = {:d}".
            format(vLRcurr, vFBcurr, Bcurr))
    # Updating previous time value
    tprev = tcurr

print('Done.')
# Releasing pins
joyLR.close()
joyFB.close()
joyB.close()

The interactive window output in VS Code should look something like what’s shown next. Note that the center position of the joystick outputs approximately half of the 3.3V supply voltage, i.e., 1.65V. Also observe the change in the B digital value as the joystick is pressed.

Running code for 30 seconds ...
X = 1.64 V , Y = 1.64 V , B = 0
X = 1.64 V , Y = 1.64 V , B = 0
X = 1.64 V , Y = 1.64 V , B = 0
X = 1.64 V , Y = 1.64 V , B = 0
X = 3.30 V , Y = 3.30 V , B = 0
X = 3.30 V , Y = 3.30 V , B = 0
X = 3.30 V , Y = 3.08 V , B = 0
X = 1.65 V , Y = 1.65 V , B = 0
X = 1.64 V , Y = 1.65 V , B = 1
X = 1.64 V , Y = 1.64 V , B = 1
X = 1.64 V , Y = 1.64 V , B = 0
X = 1.00 V , Y = 1.64 V , B = 0
X = 0.01 V , Y = 1.64 V , B = 0
X = 0.01 V , Y = 1.64 V , B = 0
X = 0.00 V , Y = 1.64 V , B = 0
X = 0.00 V , Y = 1.64 V , B = 0
X = 1.64 V , Y = 1.64 V , B = 0
X = 1.63 V , Y = 1.35 V , B = 0
X = 1.64 V , Y = 0.00 V , B = 1
X = 1.64 V , Y = 0.00 V , B = 1
X = 1.64 V , Y = 0.00 V , B = 1
X = 1.63 V , Y = 0.00 V , B = 1
X = 1.64 V , Y = 0.00 V , B = 0
X = 1.64 V , Y = 1.64 V , B = 0
...
X = 1.48 V , Y = 0.34 V , B = 0
X = 1.22 V , Y = 0.00 V , B = 0
X = 1.23 V , Y = 0.00 V , B = 0
Done.

Joystick Output Filtering

Eventually, the joystick output signal can be a bit jittery, due to the fact that not all of us are born with perfect control over our thumb motion.

The plot was generated using the Python code below, where a digital low-pass filter with a cutoff frequency of 1Hz was applied to the joystick X and Y output signals. Of course, the cutoff frequency has to be adjusted based on the desired overall system response, which includes what the joystick is controlling.

The outputs were also normalized between -1 and 1. While any transfer function can be designed to map the original 0 to 3.3V range, -1 to 1 seems appropriate if a DC motor is being controlled between full speed reverse and full speed forward.

# Importing modules and classes
import time
import numpy as np
from utils import plot_line
from gpiozero import MCP3008

# Creating ADC channel objects for the joystick inputs
joyLR = MCP3008(channel=0, clock_pin=11, mosi_pin=10, miso_pin=9, select_pin=8)
joyFB = MCP3008(channel=1, clock_pin=11, mosi_pin=10, miso_pin=9, select_pin=8)
# Assigning some parameters
tsample = 0.02  # Sampling period for code execution (s)
tstop = 10  # Total execution time (s)
vref = 3.3  # Reference voltage for MCP3008
# Preallocating output arrays for plotting
t = []  # Time (s)
xLRn = []  # Joy stick X direction output (-1 to 1)
xFBn = []  # Joy stick Y direction output (-1 to 1)
yLRn = []  # Filtered X direction output (-1 to 1)
yFBn = []  # Filtered Y direction output (-1 to 1)

# First order digital low-pass filter parameters
fc = 1  # Filter cutoff frequency (Hz)
tau = 1/(2*np.pi*fc)  # Filter time constant (s)
# Filter difference equation coefficients
a1 = -tau/(tsample+tau)
b0 = tsample/(tsample+tau)
# Initializing filter values
xLR = [joyLR.value]  # x[n]
xFB = [joyFB.value]  # x[n]
yLR = xLR * 2  # y[n], y[n-1]
yFB = yLR * 2  # y[n], y[n-1]
time.sleep(tsample)

# Initializing variables and starting main clock
tprev = 0
tcurr = 0
tstart = time.perf_counter()
# Running execution loop
print('Running code for', tstop, 'seconds ...')
while tcurr <= tstop:
    # Getting current time (s)
    tcurr = time.perf_counter() - tstart
    # Doing I/O and computations every `tsample` seconds
    if (np.floor(tcurr/tsample) - np.floor(tprev/tsample)) == 1:
        # Getting joy stick normalized voltage output
        xLR[0] = joyLR.value
        xFB[0] = joyFB.value
        # Filtering signals
        yLR[0] = -a1*yLR[1] + b0*xLR[0]
        yFB[0] = -a1*yFB[1] + b0*xFB[0]
        # Updating output arrays with normalized output
        t.append(tcurr)
        xLRn.append(-1 + 2*xLR[0])
        xFBn.append(-1 + 2*xFB[0])
        yLRn.append(-1 + 2*yLR[0])
        yFBn.append(-1 + 2*yFB[0])
        # Updating previous filter output values
        yLR[1] = yLR[0]
        yFB[1] = yFB[0]
    # Updating previous time value
    tprev = tcurr

print('Done.')
# Releasing pins
joyLR.close()
joyFB.close()
# Plotting results
plot_line([t]*2, [xLRn, yLRn], yname='X Output', legend=['Raw', 'Filtered'])
plot_line([t]*2, [xFBn, yFBn], yname='Y Output', legend=['Raw', 'Filtered'])

A Note on Joystick Drift

Even though the KS0008 that I used didn’t have a noticeable drift, it is possible that the output values at the “center position” can drift over time. This means that, over time, your system may start moving around, even when it’s supposed to be at rest at the joystick’s “center position”.

One way to address the issue is to use a digital band-pass filter, which attenuates both the low and high frequency content of a signal. Since the lower cutoff frequency will remove the DC component of the signal (the value that is held constant), it has to be chosen carefully depending on the application. If you are controlling an RC car (where there are extended periods of wide-open-throttle conditions), the band-pass filter may start attenuating the WOT value, causing the car to slow down! On the other hand, if you are controlling a drone, where zero-position drifts are highly undesirable, the band-pass filter may be what you need, since there are no extended periods of moving forward-backward or side-to-side.

The following Python program implements a digital first-order band-pass filter for the joystick outputs, with low and high cutoff frequencies respectively of 0.005 and 2 Hz. The plots show the joystick’s X-Y motion traces for 30 seconds, as well as the “center-position” rest case. The latter helps illustrate how the DC component of the signal is almost fully removed.

# Importing modules and classes
import time
import numpy as np
from utils import plot_line
from gpiozero import MCP3008

# Creating ADC channel objects for the joystick inputs
joyLR = MCP3008(channel=0, clock_pin=11, mosi_pin=10, miso_pin=9, select_pin=8)
joyFB = MCP3008(channel=1, clock_pin=11, mosi_pin=10, miso_pin=9, select_pin=8)
# Assigning some parameters
tsample = 0.02  # Sampling period for code execution (s)
tstop = 30  # Total execution time (s)
vref = 3.3  # Reference voltage for MCP3008
# Preallocating output arrays for plotting
t = []  # Time (s)
xLRn = []  # Joy stick X direction output (-1 to 1)
xFBn = []  # Joy stick Y direction output (-1 to 1)
yLRn = []  # Filtered X direction output (-1 to 1)
yFBn = []  # Filtered Y direction output (-1 to 1)

# First order digital band-pass filter parameters
fc = np.array([0.005, 2])  # Filter cutoff frequencies (Hz)
tau = 1/(2*np.pi*fc)  # Filter time constants (s)
# Filter difference equation coefficients
a0 = tau[0]*tau[1]+(tau[0]+tau[1])*tsample+tsample**2
a1 = -(2*tau[0]*tau[1]+(tau[0]+tau[1])*tsample)
a2 = tau[0]*tau[1]
b0 = tau[0]*tsample
b1 = -tau[0]*tsample
# Assigning normalized coefficients
a = np.array([1, a1/a0, a2/a0])
b = np.array([b0/a0, b1/a0])
# Initializing filter values
xLR = [joyLR.value] * len(b)  # x[n], x[n-1]
xFB = [joyFB.value] * len(b)  # x[n], x[n-1]
yLR = [0] * len(a)  # y[n], y[n-1], y[n-2]
yFB = [0] * len(a)  # y[n], y[n-1], y[n-2]
time.sleep(tsample)

# Initializing variables and starting main clock
tprev = 0
tcurr = 0
tstart = time.perf_counter()
# Running execution loop
print('Running code for', tstop, 'seconds ...')
while tcurr <= tstop:
    # Getting current time (s)
    tcurr = time.perf_counter() - tstart
    # Doing I/O and computations every `tsample` seconds
    if (np.floor(tcurr/tsample) - np.floor(tprev/tsample)) == 1:
        # Getting joy stick normalized voltage output
        xLR[0] = joyLR.value
        xFB[0] = joyFB.value
        # Filtering signals
        yLR[0] = -np.sum(a[1::]*yLR[1::]) + np.sum(b*xLR)
        yFB[0] = -np.sum(a[1::]*yFB[1::]) + np.sum(b*xFB)
        # Updating output arrays with normalized output
        # (The filtered values have no DC component)
        t.append(tcurr)
        xLRn.append(-1 + 2*xLR[0])
        xFBn.append(-1 + 2*xFB[0])
        yLRn.append(2*yLR[0])
        yFBn.append(2*yFB[0])
        # Updating previous filter output values
        for i in range(len(a)-1, 0, -1):
            yLR[i] = yLR[i-1]
            yFB[i] = yFB[i-1]
        # Updating previous filter input values
        for i in range(len(b)-1, 0, -1):
            xLR[i] = xLR[i-1]
            xFB[i] = xFB[i-1]
    # Updating previous time value
    tprev = tcurr

print('Done.')
# Releasing pins
joyLR.close()
joyFB.close()
# Plotting results
plot_line([t]*2, [xLRn, yLRn], yname='X Output', legend=['Raw', 'Filtered'])
plot_line([t]*2, [xFBn, yFBn], yname='Y Output', legend=['Raw', 'Filtered'])

Digital Filtering

After a continuous signal goes through an analog-to-digital conversion, additional digital filtering can be applied to improve the signal quality. Whether the signal is being used in a real-time application or has been collected for a later analysis, implementing a digital filter via software is quite straightforward. We already saw them being casually used in some of my previous posts, such as MCP3008 with Raspberry Pi. Now, let’s go over these types of filters in more detail. For the examples, we will be using the signal processing functions from the the SciPy Python package.

One of the drawbacks of digital filtering (or analog filtering, for that matter) is the introduction of a phase delay in the filtered signal.

The plots on the left illustrate what an actual filtered signal and its phase lag would look like (top), in contrast with an ideal filtered signal with no lag (bottom). While the ideal situation cannot be achieved with real-time signal processing, it can be achieved when processing signals that are already sampled and stored as a discrete sequence of numerical values.

Since the sequence from start to finish is available, the idea is to filter it once (causing a phase shift) and then filter the resulting sequence backwards (causing a phase shift again). Since the filtering is done both forward and backwards, the phase shifts cancel each other, thus resulting in a zero-phase filtered signal!

To see how it’s done, check out the filtfilt function and the very good application examples on the SciPy documentation. You can also go for the code used to generate the plots above on my GitHub page. I should note that I use filtfilt all the time in my signal processing endeavors. There’s a lot of information that can be extracted from a signal after you clean it up a little bit. If the time alignment between events is critical, avoiding a phase shift is the way to go.

Before we get into the real-time application of digital filters, let’s talk briefly about how they’re implemented. I’ll focus on filters that can be described as the difference equation below:

whereis the input sequence (raw signal) andis the output sequence (filtered signal). Before you abandon this webpage, let me rewrite the equation so we can start bringing it to a more applicable form. Since the idea is to find the outputas a function of its previous values and the input, we can go with:

It’s common to normalize the coefficients so . Also, for simplicity, consider a first order filter which, in terms of difference equations, only depends on the values of the current and last time steps, or and. Thus:

So, to build a first-order low-pass digital filter, all we need to do is determine the coefficients , and. Luckily for us, the SciPy MATLAB-Style filter design functions return those coefficients, reducing our task to just the implementation of the filter using Python. Before we go over a couple of code examples, let’s examine a first-order filter that I use quite a bit.

Backward Difference Digital Filter

Going from the continuous domain to the discrete domain involves a transformation where we approximate a continuous variable by its discrete equivalent. I will start with a first-order low-pass filter, by analyzing its continuous behavior, more specifically the response to a unit step input.

The continuous system response to the constant input is given by , whereis the response time of the filter and is related to the filter cutoff frequencyby:

The response time of the filter is the time it takes for the output to reach approximately 63 % of the final value. In the case of the graph above, the response time is 0.1 seconds. If you’re wondering where that comes from, just calculate the system response using _.

The differential equation that corresponds to the first order continuous system in question is:

And here is where our discrete approximation ofusing backward difference (with sampling period) comes into play:

By also transformingandinto their discrete counterpartsand, we can arrive to the difference equation below, which represents the first-order low-pass digital filter, where we used backward difference to approximate. Remember that andare the discrete current and previous time steps.

Finally, solving forgives us an equation very similar to the one we saw at the end of the previous section. Note that the coefficients andare a function of the sampling rate and the filter response time (in this particular case).

What I like about this filter implementation is that it’s fairly straightforward to see that the outputis a weighed average of its previous valueand the current input value. Furthermore, the smaller the response time (), the faster the filter is (higher cutoff frequency) and the more it follows the input value. Conversely, the slower the filter (), the more the output takes into account the previous value.

The Python code below shows how to implement this filter by placing it inside an execution loop and running a step input excitation through it. It will produce the plot shown at the beginning of this section. I invite you to experiment with different response times (cutoff frequencies) and sampling periods.

import numpy as np
import matplotlib.pyplot as plt

# Creating time array for "continuous" signal
tstop = 1  # Signal duration (s)
Ts0 = 0.001  # "Continuous" time step (s)
Ts = 0.02  # Sampling period (s)
t = np.arange(0, tstop+Ts0, Ts0)

# First order continuous system response to unit step input
tau = 0.1  # Response time (s)
y = 1 - np.exp(-t/tau)  # y(t)

# Preallocating signal arrays for digital filter
tf = []
yf = []

# Initializing previous and current values
xcurr = 1  # x[n] (step input)
yfprev = 0  # y[n-1]
yfcurr = 0  # y[n]

# Executing DAQ loop
tprev = 0
tcurr = 0
while tcurr <= tstop:
    # Doing filter computations every `Ts` seconds
    if (np.floor(tcurr/Ts) - np.floor(tprev/Ts)) == 1:
        yfcurr = tau/(tau+Ts)*yfprev + Ts/(tau+Ts)*xcurr
        yfprev = yfcurr
    # Updating output arrays
    tf.append(tcurr)
    yf.append(yfcurr)
    # Updating previous and current "continuous" time steps
    tprev = tcurr
    tcurr += Ts0

# Creating Matplotlib figure
fig = plt.figure(
    figsize=(6.3, 2.8),
    facecolor='#f8f8f8',
    tight_layout=True)
# Adding and configuring axes
ax = fig.add_subplot(
    xlim=(0, max(t)),
    xlabel='Time (s)',
    ylabel='Output ( - )',
    )
ax.grid(linestyle=':')
# Plotting signals
ax.plot(t, y, linewidth=1.5, label='Continuous', color='#1f77b4')
ax.plot(tf, yf, linewidth=1.5, label='Discrete', color='#ff7f0e')
ax.legend(loc='lower right')

SciPy Butterworth Digital Filter

The non-functional code snippet below shows how to import the SciPy signal processing module an use it to create a first-order digital Butterworth filter. Note that the implementation of the filter in the DAQ execution loop is very similar to what was done in the example code above. In this particular case however, you should be able to immediately identify the difference equation for the filter

where the arrays containing the coefficients , and, are the output of the SciPy function signal.butter for the corresponding cutoff frequency and discrete sampling period.

# Importing SciPy module
from scipy import signal

#
# Doing other initialization
#

# Discrete signal parameters
tsample = 0.01  # Sampling period for code execution (s)
fs = 1/tsample  # Sampling frequency (Hz)

# Finding a, b coefficients for digital butterworth filter
fc = 10  # Low-pass cutoff frequency (Hz)
b, a = signal.butter(1, fc, fs=fs)

# Preallocating variables
xcurr = 0  # x[n]
ycurr = 0  # y[n]
xprev = 0  # x[n-1]
yprev = 0  # y[n-1]

tprev = 0
tcurr = 0
# Executing DAQ loop
while tcurr <= tstop:
    # Doing I/O tasks every `tsample` seconds
    if (np.floor(tcurr/tsample) - np.floor(tprev/tsample)) == 1:
        # Simulating DAQ device signal acquisition
        xcurr = mydaq.value
        # Filtering signal (digital butterworth filter)
        ycurr = -a[1]*yprev + b[0]*xcurr + b[1]*xprev
        yprev = ycurr
    # Updating previous values
    xprev = xcurr
    tprev = tcurr
    # Incrementing time step
    tcurr += Ts0

Finally, as I mentioned at the top, you can implement higher order or different types of filters. Take a look at the signal.iirnotch filter, which is a very useful design that removes specific frequencies from the original signal, such as the good old 60 Hz electric grid one.

Python DAC Class

Creating a class to represent a Raspberry Pi digital-to-analog converter is a good example on how to put together the concepts we have been exploring in some of the previous posts. More specifically, DAC with Raspberry Pi and Classes vs. Functions in Python. Besides a Raspberry Pi and a couple of resistors and capacitors, we need an additional DAQ device to calibrate our Pi DAC output.

In the post DAC with Raspberry Pi, we used an MCP3008 to measure (and calibrate) the DAC output. This time around I will follow my recommendation of not using the same DAQ to generate and measure the output, using instead a LabJack U3. While a LabJack is not a calibration DAQ device, it is far more accurate than a Pi, therefore moving us in the right direction when it comes to choosing a calibration device.

The setup is quite simple: On the Pi side, a GPIO pin used as PWM input to the the low-pass filter and a GND pin. On the U3 side, the RC filter output (which is our DAC output) is connected to the analog input port AIN0, while the GND port is connected to the same GND as the Pi.

The U3 is connected to a separate computer via USB and will be running a variation of the code that can be found on my GitHub page. The Python package that I made containing the LabJack U3 class can be installed from the PyPI (Python Package Index) page.

As before, the cascaded second order filter is made of two first order low-pass RC filters, where R = 1 kΩ and C = 10 μF.

Of course, if you don’t have a LabJack, a digital multimeter should do the trick and allow you to measure the DAC output in order to do its calibration. I chose to use a LabJack to gain some insight on the DAC output signal by using its 50 kHz streaming capability (and because an oscilloscope is not among my few possessions).

The next step is to layout what our DAC class should look like. You may want to put down a “skeleton” of the class, as shown below, with a brief description of what the attributes are and what the methods do. Using the reserved word pass allows for placeholders for the methods to be easily created and filled out later.

class DAC:

    def __init__(self):
        # Class constructor
        self._vref = 3.3  # Reference voltage output
        self._slope = 1  # Output transfer function slope
        self._offset = 0  # Output transfer function intercept

    def __del__(self):
        # Class destructor
        # - makes sure GPIO pins are released
        pass

    def reset_calibration(self):
        # Resets the DAC slope and offset values
        pass

    def set_calibratio(self):
        # Sets the DAC slope and offset with calibration data
        pass

    def get_calibration(self):
        # Retrieves current slope and offset values
        pass

    def set_output(self):
        # Sets the DAC output voltage to the desired value
        pass

One of the nice things about using classes is that it’s easy to add (or remove) attributes and methods as you code away. And even better, if you’re using the interactive session of VS Code, saved updates to your methods are immediately reflected in any instance of the class that is present in the session’s “workspace” (to use a MATLAB term, for those familiar with it). In other words, if you were to call that method again using dot notation, it would behave based on the latest modifications that you saved.

In the case of our example, the attributes need to hold the values of a reference voltage (Pi’s 3.3 V if you’re not using any OpAmp with your filter) and the calibration parameters which are used so the DAC can output the correct set voltage value.

from gpiozero import PWMOutputDevice


class DAC:

    def __init__(self, dacpin=12):

        # Checking for valid hardware PWM pins
        if dacpin not in [12, 13, 18, 19]:
            raise Exception('Valid GPIO pin is: 12, 13, 18, or 19')
        # Assigning attributes
        self._vref = 3.3  # Reference voltage output
        self._slope = 1  # Output transfer function slope
        self._offset = 0  # Output transfer function intercept
        # Creating PWM pin object
        self._dac = PWMOutputDevice(dacpin, frequency=700)

    def __del__(self):

        # Releasing GPIO pin
        self._dac.close()

    def reset_calibration(self):
        self._slope = 1
        self._offset = 0

    def set_calibration(self, slope, offset):
        self._slope = slope
        self._offset = offset

    def get_calibration(self):
        print('Slope = {:0.4f} , Offset = {:0.4f}'.format(
            self._slope, self._offset))

    def set_output(self, value):

        # Limiting output
        output = self._slope*value/self._vref + self._offset
        if output > 1:
            output = 1
        if output < 0:
            output = 0
        # Applying output to GPIO pin
        self._dac.value = output

As far as the methods are concerned, a constructor and a destructor, being the latter a good idea so the GPIO pin can be released once you’re done using the class instance. Also, a group of methods to deal with the DAC calibration and a single method to set the DAC output value. Notice that the attribute names start with an underscore. That’s a loose Python convention to define private properties (or methods, if I were to use the underscore as the first character of the method name), unlike other programming languages which are more rigorous about it by making you explicitly define what’s private. The general idea of private properties (or methods) is that they’re not supposed to be accessible directly by the user of the program, therefore, being accessed internally by the code.

DAC Calibration

First, we create an instance of the DAC class on GPIO pin 18: mydac = DAC(dacpin=18)

Then, the output calibration is done with two output settings (0.5 and 3.0 V) using the method: mydac.set_output(0.5) and mydac.set_output(3.0). Each time the output is read using the streaming feature of the LabJack U3, as shown on the left. The high frequency noise in the signal is quite apparent and is a most likely due to the limitations of the Raspberry Pi hardware. For the two desired outputs of 0.5 and 3.0 V, the corresponding mean actual values are 0.614 and 2.854 V. Because the system response is linear, two data points are sufficient.

The slope and offset are calculated as slope = (3.0 – 0.5) / (2.854 – 0.614) = 1.116, and offset = 0.5 – (3.0 – 0.5) / (2.854 – 0.614) x 0.614 = -0.185. The calibrated values are applied to the DAC instance using the method mydac.set_calibration(1.116, -0.185). If a new calibration has to be done, the default values of 1 and 0 can be restored by using the method mydac.reset_calibration(). As mentioned earlier, the private attributes _slope and _offset are not dealt with directly by the user but instead through different methods. Of course, in this simple example, they could be changed directly by doing mydac._slope = 1.116 and mydac._offset = -0.185. In more complex situations there might be a few good reasons to avoid direct access to an instance’s attributes.

With the calibrated DAC, the new outputs to the same settings as before (0.5 and 3.0 V) give the results on the left, where the actual mean values are 0.493 and 3.006 V, respectively.

DAC Test

Let’s put the DAC device to the test! The Python code below can be used to generate a random sequence of steps every 0.5 seconds between 0 and 3 V (every 0.5 V). The data is collected using the LabJack U3 streaming feature and plotted in the next figure. Not too bad for our poor man’s DAC.

import time
import numpy as np
from gpiozero_extended import DAC

# Assigning parameter values
tstep = 0.5  # Interval between step changes (s)
tstop = 10  # Total execution time (s)

# Creating DAC object on GPIO pin 18
mydac = DAC(dacpin=18)
mydac.set_calibration(1.116, -0.185)

# Initializing timers and starting main clock
tprev = 0
tcurr = 0
tstart = time.perf_counter()

# Running execution loop
print('Running code for', tstop, 'seconds ...')
while tcurr <= tstop:
    # Updating analog output every `tstep` seconds with
    # random voltages between 0 and 3 V every 0.5 V
    if (np.floor(tcurr/tstep) - np.floor(tprev/tstep)) == 1:
        mydac.set_output(np.round(6*np.random.rand())/2)
    # Updating previous time and getting new current time (s)
    tprev = tcurr
    tcurr = time.perf_counter() - tstart

print('Done.')
# Releasing pin
del mydac

DAC with Raspberry Pi

Let’s put some of the concepts from the Digital-to-Analog Conversion post to work with a Raspberry Pi. To get the most out of it, you will need a few resistors and capacitors, as well as an MCP3008. The former will be used to build the low-pass analog filters, while the latter will be used to measure the DAC output. As a general guideline, you don’t want to generate and measure a signal on the same DAQ device, since all sorts of extraneous loops and noise could potentially show up. But because we’re not doing any instrumentation-critical work, we won’t be following that guideline.

The figure shows the layout of the components on a breadboard. The MCP3008 wires will use the same pins that were used in the MCP3008 with Raspberry Pi.

The PWM wire should be connected to GPIO12, which is one of the hardware-implemented PWMs on the Pi (GPIO13, 18, and 19 also have hardware PWM). As a side note, the software-implemented PWM is good for 100 to 200 Hz. The PWM period starts to fall apart above 500 Hz. The hardware PWM can go up to 8 kHz! However, that would be an overkill for our application and we will settle down for a lower value, as shown later on.

It’s also possible to identify on the board two cascaded low-pass RC filters, whose output is connected to pin 1 on the MCP3008 (channel 0 on the GPIO Zero MCP3008 class).

Let’s aim for a step response time (τ) of 0.01 s for our filters. We can then use the formula τ = RC to determine the RC product and choose suitable R and C values. In our case, we should select R = 1 kΩ and C = 10 μF. When putting everything together, make sure the negative side of the capacitors are connected to ground!

The following Python code (mostly derived from the example in Prescribed PWM duty cycle) will run a sequence of PWM duty cycle steps that can be used to better understand the DAC output. Check the list of comments below before digging into the code.

Even though the hardware PWM can go up to 8 kHz, it seemed reasonable to use 700 Hz for it. Change this number to lower values (such as 50 Hz) to see the impact on the signal ripple.
Both the PWM values as well as the ADC measurements are normalized between 0 and 1 in GPIO Zero. Therefore, they are scaled throughout the code to V_ref = 3.3 V.
The PWM duty cycle output steps (that correspond to a DAC output) span the entire range. However, the 0 and 1 values are omitted. That’s because they fully bypass the PWM, applying respectively a constant 0 and V_ref to the output.
Even though the sampling rate of the hardware-based ADC can go up to 6.25 kSamples/s, the signals are sampled at a rate of 500 Samples/s. That’s sufficient for visualization and provides plenty of margin in the execution loop.

import time
import numpy as np
from utils import plot_line
from gpiozero import PWMOutputDevice, MCP3008

# Assigning parameter values
pwmfreq = 700 # PWM frequency (Hz)
pwmvalue = [0.05, 0.2, 0.4, 0.6, 0.8, 0.95, 0.9, 0.7, 0.5, 0.3, 0.1]
tsample = 0.002 # Sampling period for data acquisition (s)
tstep = 0.5 # Interval between PWM output step changes (s)
tstop = tstep * len(pwmvalue) # Total execution time (s)
vref = 3.3  # Reference voltage for the PWM and MCP3008

# Preallocating output arrays for plotting
t = []  # Time (s)
dacvalue = []  # Desired DAC output
adcvalue = []  # Measured DAC output

# Creating DAC PWM and ADC channel (using hardware implementations)
dac0 = PWMOutputDevice(12, frequency=pwmfreq)
adc0 = MCP3008(channel=0, clock_pin=11, mosi_pin=10, miso_pin=9, select_pin=8)

# Initializing timers and starting main clock
tprev = 0
tcurr = 0
tstart = time.perf_counter()

# Initializing other variables used in the loop
count = 1  # PWM step counter
dacvaluecurr = pwmvalue[0]  # Initial DAC output value
dac0.value = dacvaluecurr  # Sending initial value to GPIO pin

# Running execution loop
print('Running code for', tstop, 'seconds ...')
while tcurr <= tstop:
    # Updating PWM output every `tstep` seconds
    # with values from prescribed sequence
    if (np.floor(tcurr/tstep) - np.floor(tprev/tstep)) == 1:
        dacvaluecurr = pwmvalue[count]
        dac0.value = dacvaluecurr
        count += 1
    # Acquiring digital data every `tsample` seconds
    # and appending values to output arrays
    # (values are converted from normalized to 0 to Vref)
    if (np.floor(tcurr/tsample) - np.floor(tprev/tsample)) == 1:
        t.append(tcurr)
        dacvalue.append(vref*dacvaluecurr)
        adcvalue.append(vref*adc0.value)
    # Updating previous time and getting new current time (s)
    tprev = tcurr
    tcurr = time.perf_counter() - tstart

print('Done.')
# Releasing pins
dac0.value = 0
dac0.close()
adc0.close()

# Plotting results
plot_line([t]*2, [dacvalue, adcvalue], yname='DAC Output (V)')
plot_line(t[1::], 1000*np.diff(t), yname='Sampling Period (ms)')

The plots show the output generated using plot_line from the utils.py module. Make sure you have a copy of it that can be imported. Observing the DAC response, it’s possible to notice that for lower PWM duty cycles the stead-state value is higher than the set point. Conversely, a higher duty cycle produces a steady-state value lower than the set point. Later on, we’ll see that it is possible to calibrate the DAC output so a one-to-one relationship can be obtained. It’s also a good practice to plot the actual sampling period that the execution loop is running on. Just to make sure it can keep up.

Output Ripple and Step Response

The next figures show the step response of our Pi DAC for a PWM duty cycle change from 25 to 75 % (0.825 to 2,457 V), respectively for PWM frequencies of 50 and 700 Hz.

For the 50 Hz PWM, the ripple in the output is quite significant and can be easily identified. We can also see that there’s no offset between the desired and the actual steady-state DAC output signals. In the other hand, for the 700 Hz PWM, the ripple in the output is much smaller and isn’t related to the PWM frequency. Unlike the lower frequency case, there’s now an offset between desired and actual outputs.

That offset is mostly caused by the transitions from low-to-high (and vice-versa) in the PWM signal. Because they’re not instantaneous (think of them as very steep ramps instead), the PWM signal going into the low-pass filter will assume values other than low (0V) and high (V_ref). These intermediate values become more relevant the higher the frequency is. Interestingly enough, the offset will be zero at 50% duty cycle, regardless of f_PWM.

The response time is not affected by the PWM frequency. The theoretical value of 0.02 s (for the combined cascaded low-pass filters used in this example) wasn’t observed on the actual signal. The higher response time of 0.05 s is due in part to the PWM not being a perfectly “square” wave form. One should also speculate other delays in the Raspberry Pi hardware itself, as well as the effect of input/output impedances.

DAC Calibration

The relationship between the set point and the actual analog output values is quite linear and can be determined with the help of a modified version of the Python code above, which also performs a linear regression on the collected data. The code stores (and averages) analog output data for the steady-state portion of every step of the input sequence. The data points are then used in a linear regression to determine the transfer function between desired and actual analog output values. The graphs below show the linear regression and the calibrated output. Check out the code on my GitHub repository for more details.

Final Remarks

While most instrumentation and DAQs these days don’t require analog inputs to do what they’re supposed to do, some of the older pieces of equipment out there still do. Op Amps can be used to expand the output voltage range of our simple Pi DAC to a more common 0 to 5 V.

In this post, I create a DAC class that can be used like any other GPIO Zero class (such as PWMOutputDevice or MCP3008) making it much easier to work with the code. The exercise of building the class also comes in handy as we explore more efficient ways to write DAQ software in Python.

Digital-to-Analog Conversion

As the counterpart of analog-to-digital conversion, DAC will take a digital signal and convert it to an analog one. Paraphrasing what I mentioned in previous posts, ADC and DAC walk hand-in-hand bridging the gap between the continuous and discrete domains that modern machines and devices have to deal with.

Even though there are several types of DACs, I will talk about the PWM-based one, where a reference input voltage is switched (in the form of a digital train of pulses) into a low-pass analog filter, producing an analog output. Before we get there, let’s first talk about Pulse Width Modulation, a very clever way to go from the discrete to the continuous domain, invented in the mid-seventies. A PWM signal has two main characteristics: its frequency and its duty cycle.

The figure illustrates the frequency and the concept of duty cycle for a 10 Hz PWM signal. In this case, the period of the PWM is 1/f_PWM = 0.1 s. The chosen duty cycle is 25 %, which corresponds to the time the output is “on” (0.025 s). During the remainder of the duration of the PWM period (0.075 s) the output is “off”. As a matter of fact, it’s the modulation of the duty cycle that will control the output level of the DAC system.

Since it’s either fully “on” or fully “off”, the PWM signal is fundamentally a digital one. It’s the frequency and duty cycle of that “on-off” train of pulses, coupled with the response of the sub-system it’s interacting with, that will result in an analog behavior of the output of the overall system, in our case the digital-to-analog converter.

As a side note, the Python code used to generate the graphs in this post can be found on my GitHub page. Make sure to check it out and experiment with the DAC parameters.

The simplest DAC implementation is shown on the left. It consists of a low-pass first-order passive RC filter connected to the PWM output signal V_ref. Typically, for a TTL circuit, the voltage can be either 0 (low) or 3.3V (high).

The cutoff frequency (rad/s) of the filter is f_c = 1/(RC). Along with the PWM frequency, it can be used to obtain the desired behavior of this simple DAC.

Output Ripple

The PWM signal can be thought as a sequence of low-to-high and high-to-low step inputs into the RC filter. Even with a constantly switching input, this first order system will eventually produce a near constant output. The figures below show the effect of the RC filter cutoff frequency f_c as well as the PWM frequency f_PWM on the amplitude of the output ripple.

As it can be observed in the plots above, the output voltage will converge to (duty cycle) x V_ref , which for the duty cycle of 25 % (used in the example) produces an output of 0.825 V when V_ref = 3.3 V. Also, reducing the filter cutoff frequency or increasing the PWM frequency are both effective in reducing the output signal ripple. While having a slow RC filter can greatly reduce the ripple, that will negatively affect the DAC response time when it’s subject to a change in the desired output value.

DAC Step Response

In the case of our DAC, its response time is essentially the response time of the RC filter, i.e. 1/f_c = RC (with f_c in rad/s). Leaving the details for a future post on analog and digital filtering, for the first order system under consideration, the response time can be defined as the time it takes for the output to reach about 63 % of its final (steady-state) value, after a step input is applied.

The next figure shows how an appropriate combination of a higher PWM frequency and a higher filter cutoff frequency can produce a faster response time, while maintaining roughly the same amplitude of the DAC output ripple. For our example, the response time based on the RC value goes from 0.16 s down to 0.04 s.

Cascaded Filters

The output ripple can be further improved by cascading two first-order RC filters, as shown below. The newly formed second-order filter will have twice the attenuation (dB/decade) above the cutoff frequency. It’s an easy-to-do improvement which is often adopted.

Practical Considerations

Just like its ADC counterpart, DACs also have an inherent quantization in the signal. The digital PWM, generated from the microprocessor clock, has a finite resolution (duty cycle levels) given by the number of bits of the PWM. Typical resolutions range from 8 to 12 bits.

Actual PWM frequencies for DAC devices are much higher than the ones used in the examples above, which were chosen mainly to illustrate the concepts. For applications involving testing of physical systems or automation, frequencies between 500 and 2000 Hz are reasonable. In audio applications, several times the highest frequency of the audible range (20 kHz) seems to be the case. That is, 100 to 200 kHz. However, just as we observed when using higher order filtering, there’s wiggle room to reduce these numbers while still having a good compromise between response time and ripple.

It’s important to make a distinction at this point: if the PWM is used to drive a DC motor directly (or with some power amplification), the motor itself will behave approximately as a first order system. Hence, there’s no need for the RC filter. In this type of application, PWM frequencies of 100 to 200 Hz are just fine.

Amplification of the DAC output beyond the typical 3.3V V_ref, can be achieved by using an Op Amp (Operational Amplifier). For a simple non-inverting amplifier, the output voltage is given by:

V₂ = V₁(1+R₁/R₂)

By choosing R₁ = R₂ = 10 kOhm, we can double the output to accommodate a more typical 0 to 5V requirement.

In my next post, we will explore the implementation of a DAC using a Raspberry Pi and a few resistors and capacitors. The Pi has hardware enabled PWM on two of its GPIO pins, which can reach frequencies in the kHz range! That should be an interesting project.

MCP3008 with Raspberry Pi

Since the Raspberry Pi is fundamentally a digital device, any I/O that is done through it’s GPIO pins will happen through high (one) and low (zero) states. When it comes to input signals that are analog, most likely from a transducer, they need to be converted to the digital domain so the Raspberry Pi can understand them. The MCP3008 is a 10-bit 8-channel analog-to-digital converter chip that has a very straightforward API implementation in GPIO Zero. The MCP3008 device uses an SPI (Serial Peripheral Interface) communication protocol that is fully taken care by the API.

SPI requires four pins for the communication:

A “SCLK” pin which provides the clock (pulse train) for the communication.
A “MOSI” pin (Master Out, Slave In) which the Pi uses to send information to the device.
A “MISO” pin (Master In, Slave Out) which the Pi uses to receive information from the device.
A “CEX” pin which is used to start and end the communication with the device, usually for one data sample at a time. Because the Raspberry Pi can have more than one device sharing the same SCLK, MOSI, and MISO pins, the CE0 or CE1 pin will indicate which device to use.

Additionally, the MCP3008 requires two separate voltage inputs ranging from 2.7 to 5.5V. The supply voltage at 5V allows for a higher data sampling rate of 200 kSamples/s (which should be plenty for the Raspberry Pi). A reference voltage of 3.3V improves the absolute resolution to about 3 mV instead of 5 mV since, for a 10-bit device, 1LSB = V_REF/1024. However, the choice of V_REF also has to satisfy the analog input range requirement of the transducer. If the input is higher that 5V, a resistive voltage divider can be used. Finally, in the case of a Raspberry Pi, both the digital and analog grounds on the device can share the same GND pin.

Important! You must enable the SPI for the MCP3008 communication to work. To do so, go to the Raspberry Pi Configuration menu and select Enable for the SPI in the Interfaces tab.

Hardware vs. Software SPI implementation

The Raspberry Pi has both a hardware and a software implementation of the SPI protocol. It is worth investigating what kind of sampling rates can be achieved with those two types.

The figure on the left shows a comparison between software and hardware implementation results of the sampling period per channel as a function of the number of analog inputs. The code used to collect the data, as well as GPIO pin configurations, can be found on my GitHub page.

For each point on the graph, data was sampled as fast as possible, i.e., with no time step clock, where the only I/O occurring in the execution loop was reading the MCP3008 output. Computational time storing the data in an array was negligible. Not surprisingly, the sampling period increases linearly with the number of input channels and, on my Raspberry Pi 4B, the hardware implementation is in average 3.6 times faster than the software one.

The hardware sampling rate is approximately 6.25 kSamples/s and the software one is approximately 1.72 kSamples/s. Note that it’s common to use Samples (or kSamples) per second for a DAQ device, where the number usually gives the maximum throughput of the device. For example, 5 channels with the hardware implementation still gives 6.25 kSamples/s. However the rate is 1.25 kSamples/s per channel.

Still on the 5 channels example, you can very likely get away without an anti-aliasing filter. Noise above 625 Hz (half of the 1250 Hz sampling rate) should have negligible amplitude and shouldn’t show up as an alias in the actual signal. That doesn’t mean that the analog signal being collected with the MCP3008 doesn’t have (relatively) low frequency noise, which in turn can be filtered using a digital low-pass filter in the DAQ software.

Example with a Potentiometer

A potentiometer is basically a variable resistive voltage divider that will output a voltage between a high and a low voltage value, in our case V_REF and GND. Because potentiometers convert a linear or angular position input to a voltage output, they do fall into the transducer category.

In this example, the analog output of the potentiometer will be sampled by the DAQ software. To add a bit of visual excitement, the sampled voltage will be used to change the PWM output that controls the intensity of the LED.

The wires on the right hand side of the figure are all connected to the appropriate pins on the Raspberry Pi. Since hardware SPI will be used:

SCLK = GPIO11
MISO = GPIO9
MOSI = GPIO10
CEO = GPIO8

PWM is connected to GPIO16. Note that the LED requires a series resistor, where R can be anywhere between 300 to 1000 Ohm.

Also, pay attention to the semi-circular notch on the MCP3008 chip, so it can be inserted with the right orientation into the breadboard.

The Python code for the example is shown below, where there are a few things worth mentioning:

The sampling period is 0.02 s, which gives a sampling rate of 50 Samples/s. That is far below the limit for the hardware SPI implementation (6.25 kSamples/s) and even the software one (1.72 kSamples/s).
The execution loop concept with a time step clock is used to ensure a fairly accurate sampling period. Both I/O and computation tasks are performed within a time step with plenty of time to spare.
Due to the noise in the analog signal, most likely due to the questionable quality of the potentiometer, a digital low-pass filter is implemented in the code. The cutoff frequency of 2 Hz can be modified for some extra exploration of its effect.
For the plotting of the sampled signal to work, a module named utils.py has to be accessible for importing as discussed in my previous post A useful Python plotting module.

import time
import numpy as np
from utils import plot_line
from gpiozero import PWMOutputDevice, MCP3008

# Creating LED PWM object
led = PWMOutputDevice(16)
# Creating ADC channel object
pot = MCP3008(channel=0, clock_pin=11, mosi_pin=10, miso_pin=9, select_pin=8)
# Assining some parameters
tsample = 0.02  # Sampling period for code execution (s)
tstop = 10  # Total execution time (s)
vref = 3.3  # Reference voltage for MCP3008
# Preallocating output arrays for plotting
t = []  # Time (s)
v = []  # Potentiometer voltage output value (V)
vfilt = []  # Fitlered voltage output value (V)
# First order digital low-pass filter parameters
fc = 2  # Filter cutoff frequency (Hz)
wc = 2*np.pi*fc  # Cutoff frequency (rad/s)
tau = 1/wc  # Filter time constant (s)
c0 = tsample/(tsample+tau)  # Digital filter coefficient
c1 = tau/(tsample+tau)  # Digital filter coefficient
# Initializing filter previous value
valueprev = pot.value
time.sleep(tsample)
# Initializing variables and starting main clock
tprev = 0
tcurr = 0
tstart = time.perf_counter()

# Execution loop
print('Running code for', tstop, 'seconds ...')
while tcurr <= tstop:
    # Getting current time (s)
    tcurr = time.perf_counter() - tstart
    # Doing I/O and computations every `tsample` seconds
    if (np.floor(tcurr/tsample) - np.floor(tprev/tsample)) == 1:
        # Getting potentiometer normalized voltage output
        valuecurr = pot.value
        # Filtering value
        valuefilt = c0*valuecurr + c1*valueprev
        # Calculating current raw and filtered voltage
        vcurr = vref*valuecurr
        vcurrfilt = vref*valuefilt
        # Updating LED PWM output
        led.value = valuefilt
        # Updating output arrays
        t.append(tcurr)
        v.append(vcurr)
        vfilt.append(vcurrfilt)
        # Updating previous filtered value
        valueprev = valuefilt
    # Updating previous time value
    tprev = tcurr

print('Done.')
# Releasing pins
pot.close()
led.close()
# Plotting results
plot_line([t]*2, [v, vfilt], yname='Pot Output (V)', legend=['Raw', 'Filtered'])
plot_line([t[1::]], [1000*np.diff(t)], yname='Sampling Period (ms)')

The graphs below show the results of one run on my Raspberry Pi 4B. Observe how the low-pass filter removes the noise from the raw signal as I twirl around the potentiometer knob. That of course comes with a price: the filtered signal now lags the original signal. In DAQ applications where the signal doesn’t need to be filtered in real-time, there are “zero-lag” filtering techniques that can be applied after the signal was acquired. Also note the consistent sampling time step of 20 ms.

A note on Quantization

What would happen if I never touched the potentiometer during the program execution, leaving its output at a constant value? The result is shown in the next graph, where the quantization effect is quite clear on the (very low amplitude) signal noise riding on top of the output voltage. The discrete nature of the output can now be observed, and the values are exactly V_REF/2¹⁰ = 3.3/1024 = 0.0032 V apart!

Analog-to-Digital Conversion

Also known as ADC, is the process of converting a physical quantity, usually represented by a voltage at this point, to a set of bits that’s more suitable for a computer to digest. Bearing in mind that the world we live in is analog in nature and the computers we use to interact with it are digital, going from the continuous domain to the discrete domain is an integral part of most machines and devices out there.

Back in the day, when micro-processors and computers weren’t around or were still in their infancy, many systems would function entirely in the analog domain. Take for instance an old record player. From the needle running in the record tracks (a continuous profile of peaks and valleys) to the sound coming out of the speakers, the entire sequence of changes in physical quantities occurred in the analog domain. Nowadays, still on the same example, music is stored as digital content and only converted back to an analog signal when it’s time to make it to the speakers, so we can listen to it with our good old analog ears. This last process of going back from the discrete to the analog domain is handled through a digital-to-analog conversion (DAC) and will be a topic of another post. The entire process can be illustrated with the diagram below.

As I alluded to at the very top, a sensor (or more specifically a transducer) will convert the measured physical quantity to an electrical signal, and ultimately a voltage. The next step in the process is the signal conditioning where some sort of filtering will be applied to the measured voltage. Not only to avoid aliasing, which we will discuss later on, but some times to accomplish other goals, like for example remove a DC component or a very slow drift in the signal.

The next block is the analog-to-digital conversion (ADC), which has mainly two tasks to it: sampling and quantization. Once they are performed, the continuous signal will be broken down into packets of bits that the computer can now handle. We will see that the signal is not only discrete in time but also each sample can only have a finite number of values.

At the bottom half of the diagram is the flip side of the coin: the digital-to-analog (DAC) conversion, an eventual power amplification and, generally speaking, the electro-mechanical actuators.

You should also note that if we are talking about data acquisition only, the focus is the top half of the diagram. If automation is the case, then the entire diagram is in scope.

Sampling

A physical system changes its states continuously over time and it would be impractical to collect continuous data to be used by an inherently discrete system such as the computer. Using the music example, and going back several decades, data would be collected and recorded on magnetic tapes so it could be later edited, mixed, and reproduced. With the computer entering the landscape, those analog systems had their days numbered.

The sampling process occurs in an integrated circuit (IC), where an electronic switch will bring a sample of the signal into the AD converter at a fixed sampling period.

The figure on the right shows 2 seconds of an arbitrary continuous signal sampled at 0.1 seconds (or 10 Hz). The first thing to notice is that we can now store 20 data points instead of an impractical infinite number of continuous points.

There are other considerations that dictate the sampling period (or frequency) than just the amount of data that can be manipulated. The main one is the frequency content of the signal. In other words, how “fast” does the signal change over time. For good signal reconstruction in the time domain, a sampling rate of 10 to 20 times the frequency content of the signal should be used. Most importantly, it must never be below 2 times the highest frequency component of the measured signal. Otherwise, aliasing will potentially occur to the sampled data.

As a quick note, all the plots used in this post were generated using Python and the Matplotlib package. You can get the code on my GitHub page.

Nyquist Theorem and Aliasing

According to the Nyquist sampling theorem, the sampling rate should be at least two times the maximum frequency component of the measured signal. The figure below shows a 20 Hz sinusoidal signal sampled at 21 Hz. Any sampling rate below 40 Hz will result in an alias of the original signal. In this case a new 1 Hz sinusoidal wave will show up.

The frequency of the aliased signals (yes, there’s more than one alias) is not exactly intuitive, as it appears to be in the previous figure. The same original 20 Hz signal sampled at 7 Hz would also produce a 1 Hz aliased signal! Check out the Python code on my GitHub page for a couple of formulas and to explore different sampling rates.

In practical terms, it is hard to guarantee that the signal being sampled won’t contain stray frequencies above the Nyquist frequency. The best way to avoid aliasing is to use an analog low-pass filter before the signal is sampled. Even a simple passive RC filter should suffice, provided a reasonable frequency transition band to accommodate for the poor filter attenuation right above its cut-off frequency.

Back to the music example, audio is usually sampled at 44.1 kHz. Since the human ear can respond to frequencies up to about 20 kHz, sound should be sampled at least at 40 kHz. 44.1 kHz gives a Nyquist frequency of 22.05 kHz and therefore a transition band of 2.05 kHz for a filter with a cut-off frequency of 20 kHz. However, it’s not uncommon to record audio at 48 kHz (or even 96 kHz). The former gives a Nyquist frequency of 24 kHz and therefore a wider 4 kHz transition band for a low-pass filter to attenuate the signal above it’s cut-off frequency of 20 kHz.

For the classically inclined, how does 48 kHz compare with the frequency contents of the music that’s being recorded? A soprano can hit C6 (1046.5 Hz), a violin A7 (3520 Hz), and a piano C8 (4186.01 Hz). Those frequencies are well bellow the Nyquist frequency of 24 kHz and even 10 times (or more) lower than the 48 kHz sampling frequency, allowing also for good signal reconstruction in the time domain.

Quantization

The second task in the AD processing is the quantization of the sampled values. Even though your computer can deal with floating numbers (they’re still a bunch of bits in disguise), the ADC device still ships out a packet of bits for every analog sample. The analog value will be discretized based on a reference voltage and the number of bits used for the quantization.

The ideal transfer function (TF) between the analog input voltage and the digital output code is a straight line where for every possible input there’s a unique output. In other words, the ideal ADC has an infinite resolution. Because the output of the conversion is a digital value, the transfer function of a perfect ADC is in fact a staircase function where each step is 1 LSB (Least Significant Bit).

As shown in the figure, let’s consider the case with a reference voltage V_REF = 2 V and a resolution of 3 bits. The step size (1 LSB) is 0.25 V, or in the general case V_REF / 2ⁿ, where n is the number of bits of the ADC.

Any analog input within a step will be coded to the same digital output. For example, 0.625 to 0.875 V will be coded as 011 (or 3 in decimal representation). Using the reference voltage of 2 V, that would give us 3V_REF / 8 = 0.75 V. Therefore, the quantization error goes from +0.125 V to -0.125 V. In other words +/- 0.5 LSB. This would be the case for the perfect ADC shown in the figure, more specifically a single-ended mode configuration with adjusted quantization. A good source explaining other configurations and also sources of errors that occur on an actual ADC can be found here. The full scale (FS in the figure) is V_REF – 1LSB, which in the case of this example would be 1.75 V.

If a 4-bit converter is used, the quantization error for the same reference voltage drops to +/- 0.0625 V (still +/- 0.5 LSB) and the FS goes up to 1.875 V.

As we move up in bit resolution, the ADC output starts to approximate the ideal transfer function. In the case of a 10-bit ADC, the quantization error becomes +/- 0.5V_REF/1024. Which in our case is approximately +/- 0.001 V, and in the general case +/- 0.05 % of the reference voltage. For a high accuracy 16-bit ADC, the error is approx. +/- 0.0008 % of the reference voltage.

Note that there are additional errors introduced in the quantization process due to offsets, non-linearities, and just noise in general in the electronic circuitry. A 12-bit converter will typically have +/- 0.13 % (instead of +/- 0.012 %) and a 16-bit a typical +/- 0.01%.

Practical Considerations

So, how high should you go with the sampling rate and number of quantization bits? How about the anti-aliasing filter? Keep in mind that faster, higher-resolution hardware comes with a steeper cost. And high sampling rates bring the burden of extra storage space and processing power.

If you are working with audio, definitely 48 kHz and 16-bit ADC for sound recording. 96 kHz (or even 172 kHz) and 24-bit if you’re a pro. And most definitely an anti-alias filter. By the way, the fact that 96 kHz is two times the 48 kHz base rate makes it very easy to down sample a recording for later reproduction: just throw away every other data point.

If you’re toying around with a Raspberry Pi, the sampling rate is the one that suits your needs. For slow changing temperatures 1 Hz is plenty. Pressures maybe 10 to 20 Hz. A 10-bit ADC MCP3008 chip is inexpensive and has plenty of resolution. An anti-aliasing filter is probably overkill. The higher frequency content of noise is also usually low amplitude and therefore shouldn’t affect you signal very much. However, it’s always a good idea to sample your data first at a higher frequency and inspect it as a function of time. Also, be careful with some typical electrical noise like the 60 Hz one caused by the AC outlet.

In a lab or instrumentation setting, 12 to 16-bit and an anti-aliasing filter, such as a simple passive RC one, should be considered. Again, try different sampling rates and filters. As I mentioned before, it’s always important to look at your data as a function of time first with as high as a sampling rate that can be run with the DAQ device at hand. Once you understand what the signal is all about, then a more conscious decision about the hardware requirements can be made.

Finally, just to have some fun with Python, below is the arbitrary signal of this post with two different sampling rates and resolutions. Observe the effects of both the sampling and the quantization on the original signal.

DAQ, SBC, etc.

DAQs, SBCs, Embedded Systems, and computers. How do they all fit together? Let’s first spell out the acronyms: DAQ stands (surprisingly) for Data Acquisition in general. It includes the software and hardware side of things. For this discussion, let’s stick to the hardware side. SBC stands for Single Board Computer. Something like your PC, but where everything that’s inside the box is on one single integrated circuit board. Furthermore, the SBCs we are interested in also have GPIO (General Purpose Input and Output), which we’ll get to later.

When we talk about data acquisition (in the broader sense) and automation, we are looking to interface the real world with a digital system, through sensors and actuators. Information will be collected or sensed from the real world, some sort of analysis or decisions will be made, and if we are doing automation, an action will be taken back into the real world.

As illustrated below, there are mostly three ways to put all the pieces together. Starting from the left, the most flexible, and usually most expensive arrangement, is a computer running on top of a DAQ device. In this type of application, it’s important to be able reconfigure your system to whatever your data acquisition needs are. A lab is usually where you will find this type of system, where the accuracy of what you’re measuring also plays a big role. A short list of DAQ devices should include LabJack, National Instruments and Yokogawa.

In the middle are the SBCs. Usually far less powerful than your computer, but still a fully functional computer. The SBCs of interest to us are the ones that also have GPIO pins. These General Purpose Input and Output pins enable a connection between the sensors/actuators and the software in the computer interacting with them. They are very cost-friendly and are used in education or in DIY projects, where accuracy and reliability are less critical. Also, the SBCs will mostly handle only digital signals. Any conversion from the real world analog signal to something the SBC can understand has to happen in the sensor or a separate ADC (Analog-to-Digital Conversion) component. Most DAQs, on the other hand, have the capability to handle both digital and analog signals, doing the needed conversions and providing the computer with a signal it can understand. A popular example of an SBC is the Raspberry Pi.

On the right hand side of the figure are the embedded systems, which by the way, are all around us. From the very sophisticated automotive ones (engine control units, transmission controllers, etc.) to the less glamorous ones, like the one in your dishwasher, these are dedicated computers, or more specifically micro-processors that only do one group of tasks and are designed to work with a specific set of sensors and actuators. They are more cost effective due to mass production and some of them are extremely fast, performing real-time calculations so all the actuation decisions can be made on-the-go. Depending on the application, reliability is quite critical. You don’t want the micro-processors in that air plane having to be rebooted amid flight. A good example of a very affordable embedded system is the Arduino.

A note on speed

When it comes to collecting data, it’s not uncommon to have an execution loop in the code that controls the acquisition. That is especially the case with automation. The tasks within an execution loop can be primarily grouped into either computations or I/O. The first group contains all activities where some processing of the data is done. The second group is where the inputs from and outputs to the real world are occurring, through a DAQ device or built-in I/O.

How fast an execution loop can run depends on how fast the tasks within each of those two groups can run. To some extent, it will come down to the hardware that’s being used for either the computation or the I/O. Let’s not forget that the number of data channels that are being handled is also an important factor. The better (and usually more expensive) the hardware, the faster it can run. Generally speaking, a computer with a DAQ device should have no problem with the computation part but may be held back by the DAQ device I/O.

On the other end of the spectrum, an embedded system, due to cost constraints, doesn’t have the computing power of a regular PC and may face some limits on how many computations can be performed. That is an important aspect, especially because these types of systems are used in real-time applications. The time it takes to get the sensor input data, do computations, and set the output actuator data (for one execution step) has to fit within the time frame of what’s happening with the real world that the embedded system is interacting with.

In terms of numbers (depending on complexity of calculations, hardware, sensors, and number of channels), a computer with a DAQ device should be able to run execution loops in the 10 to 100 Hz range. SBCs, which have the built-in I/O, can bring that number to 200 Hz or more. Embedded systems are all over the execution frequency range. Remember: they have to do the real-time task they’re designed for in the most cost-effective and reliable way possible. Real-time requirements for an autonomous vehicle are not the same as the ones for the thermostat controlling the temperature in your home.

Python Implementation

Share this:

Digital Filtering and Pulse Detection

Real-time Pulse Rate Detection

Share this:

Joystick Output Filtering

A Note on Joystick Drift

Share this:

Backward Difference Digital Filter

SciPy Butterworth Digital Filter

Share this:

DAC Calibration

DAC Test

Share this:

Output Ripple and Step Response

DAC Calibration

Final Remarks

Share this:

Output Ripple

DAC Step Response

Cascaded Filters

Practical Considerations

Share this:

Hardware vs. Software SPI implementation

Example with a Potentiometer

A note on Quantization

Share this:

Sampling

Nyquist Theorem and Aliasing

Quantization

Practical Considerations

Share this:

A note on speed

Share this: