JohnAlan (Jack) Keegan and I have been doing some interesting work on real-time electroencephalogram (EEG) analysis recently and we had a really productive practical session in the lab today. EEG is a recording of electrical activity in the brain, usually measured using electrodes on the scalp. Jack is doing some fascinating work on Visual Evoked Potentials (VEPs), which are EEG patterns that occur in response to visual stimuli. They are recorded from the visual cortex (at back of the head). Jack is developing new types of visual stimulus which we hope will elicit VEPs more easily and reliably.
The system Jack uses to record EEG is the BioSemi ActiveTwo (shown above), which is a high-end biopotential acquisition system used by many researchers all over the world. It’s expensive and a little bit cumbersome to use, but it does a really great job of recording biopotentials and does it reliably. For a long time, we have been using the software that comes with the ActiveTwo (called ActiView) to record data during our experiments and then analysing the resulting data offline. For example, in Jack’s experiments, he would wire up a subject, display several visual stimuli and record the resulting EEG to files. Only once the experiment was complete would he be able to analyse the data and find out if any of the stimuli worked as expected. ActiView displays the time-domain signals in real time, but it doesn’t perform much real-time analysis.
So… we’ve had this long-standing intention to implement real-time data analysis and display (of frequency spectra and the like), so that it’s possible to manipulate stimuli in real-time and see what effect the changes are having on the neurological response. A year or more ago, Jack developed some software to access data directly from the BioSemi libraries. It was a bit of a struggle, but it was finally more or less working. However, it had worked out to be quite fiddly, and I think we both just kind of ran out of steam on that approach.
What we did today was to access the data in a completely different way by connecting to ActiView via TCP/IP (which is something it allows). It turned out to be incredibly easy to connect to ActiView via TCP/IP using Python and parse the data packets it streams over the network. As it happens, we were actually running ActiView and the Python program on the same PC, so we connected via the localhost loopback interface (which is always IP address 127.0.0.1). All the Python program does is connect to ActiView, which needs to be running and recording data. The Python program parses incoming data packets, storing only the data from channel 1. Once it has a complete window of data (we chose a window length of 512 arbitrarily) it performs an FFT and plots the spectrum of the signal (actually just the real part of the DFT). It repeats this 50 times, updating the same plot each time.
We didn’t have the BioSemi connected to a subject while testing the program today, so the data we recorded is just noise, but we’re optimistic that the program is working correctly and we’re looking forward to trying it on a real EEG signal!
This is a snapshot of the DFT (real part only) that was displayed by the Python program using the matplotlib library. Because the data we recorded was just noise, no particularly interesting features are evident in the spectrum. However, the symmetry we expect to see in the DFT of a real-valued time-domain signal is clearly present.
ActiView needs to be running to control the recording and stream the data. During recording, it provides a scrolling real-time display of the data. In this case, there was no subject connected to the equipment, so it just recorded low amplitude noise:
This is the TCP Server configuration tab in ActiView. In this screenshot, the Python program is already connected via TCP, so the corresponding green light is illuminated. Once you select the number of channels to stream over the network, the information labels at the bottom of the screen display the packet size that will be used and how many samples per channel will be included in each packet. Jack’s Python code is configured to receive 384-byte packets, containing 16 samples per channel for 8 channels. Each sample is 3 bytes long.
Jack runs the Python program from within iPython, as shown below:
This is the complete source code for the current program.
# # test_plot.py - Written by Jack Keegan # Last updated 16-1-2014 # # This short Python program receives data from the # BioSemi ActiveTwo acquisition system via TCP/IP. # # Each packet received contains 16 3-byte samples # for each of 8 channels. The 3 bytes in each sample # arrive in reverse order (least significant byte first) # # Samples for all 8 channels are interleaved in the packet. # For example, the first 24 bytes in the packet store # the first 3-byte sample for all 8 channels. Only channel # 1 is used here - all other channels are discarded. # # The total packet size is 8 x 16 x 3 = 384. # (That's channels x samples x bytes-per-sample) # # 512 samples are accumulated from 32 packets. # A DFT is calculated using numpy's fft function. # the first DFT sample is set to 0 because the DC # component will otherwise dominates the plot. # The real part of the DFT (all 512 samples) is plotted. # That process is repeated 50 times - the same # matplotlib window is updated each time. # import numpy # Used to calculate DFT import matplotlib.pyplot as plt # Used to plot DFT import socket # used for TCP/IP communication # TCP/IP setup TCP_IP = '127.0.0.1' # ActiView is running on the same PC TCP_PORT = 778 # This is the port ActiView listens on BUFFER_SIZE = 384 # Data packet size (depends on ActiView settings) # Open socket s = socket.socket(socket.AF_INET, socket.SOCK_STREAM) s.connect((TCP_IP, TCP_PORT)) # Create a 512-sample signal_buffer (arange fills the array with # values, but these will be overwritten; We're just using arange # to give us an array of the right size and type). signal_buffer = numpy.arange(512) # Calculate spectrum 50 times for i in range(50): # Parse incoming frame data print("Parsing data") # Data buffer index (counts from 0 up to 512) buffer_idx = 0 # collect 32 packets to fill the window for n in range(32): # Read the next packet from the network data = s.recv(BUFFER_SIZE) # Extract 16 channel 1 samples from the packet for m in range(16): offset = m * 3 * 8 # The 3 bytes of each sample arrive in reverse order sample = (ord(data[offset+2]) << 16) sample += (ord(data[offset+1]) << 8) sample += ord(data[offset]) # Store sample to signal buffer signal_buffer[buffer_idx] = sample buffer_idx += 1 # Calculate DFT ("sp" stands for spectrum) sp = numpy.fft.fft(signal_buffer) sp = 0 # eliminate DC component # Plot spectrum print("Plotting data") plt.plot(sp.real) plt.hold(False) plt.show() # Close socket s.close()