Many real-time filters were built in Labview to extract the relevant features from the sensor data. Among the filtered parameters included: instantaneous values, zero-mean, rectification, beats (peak detection), beat timings (which give the information to compute the rate, or tempo), beat amplitudes, first derivatives, second derivatives, envelopes (smoothed instantaneous values), power spectral densities, spectral analyses and spectrograms (using both the Fast Fourier Transform and the Fast Hartley Transform), and noise analysis. Most of these were used in the final performance system; some were used to analyze the data in order to build appropriate mappings. For example, a number of the above-mentioned filters were written and used initially to determine if the frequency domain contained any useful information that differed from the time domain. After spending several weeks experimenting with various frequency-domain filters I determined that they did not provide any extra functionality, and so they were not included in the final performance system. Nonetheless, since they were built and used for testing purposes, they are included in the present discussion. The algorithms and details of all the filters are given below.
Choosing how to represent the instantaneous value of a signal is not always a straightforward matter, when various sampling rates and linear transforms are available to apply to the signal to make the data easier to apply. Sometimes it is advantageous to sample at a relatively low rate (relative to the frequency content of the signal or the distribution of frequencies with the most power), since it reduces the load on the processor and reduces multiplexing errors. However, if done in a way such that it does not satisfy the Nyquist sampling criteria it can introduce aliasing problems. The problem that was encountered with the ComputerBoards sampling hardware was that in order to time the sampling rate the data had to be read into a buffer; when the buffer is released and cleared, it often generated an instantaneous glitch. EMG sensors carry frequency information up through 500 Hz and therefore have to be sampled at 1 KHz in order to satisfy the Nyquist criteria. However, sometimes it seems to be practical for real-time use to sample the data at the 330 Hz board default in order to avoid the glitching problem. EMG signals that were needed for advanced filtering and feature recognition were sampled at 3 Khz, whereas some of the less crucial signals were sampled at 330 Hz. The challenge was to keep keep the sampling rates high and the window sizes small in order to ensure real-time responses.
In order to make the sensor data more usable for mappings, it was run through a zero-mean filter. This estimated the mean from a fixed number of samples and subtracted it from the signal. This yields a signal that has a baseline at zero. The mean was continually re-estimated at every new window of 150 samples.
After giving the data a zero-mean, full rectification was done by taking the absolute value of every sample. Half rectification (removing the negative values of the signal) was explored for EMG signals since it is computationally less expensive, but it was rejected because it left gaps in the signal that caused problems with the mappings.
Beat detection was done for both the right and left biceps signals. While it was suggested that I first low-pass filter the data before finding beats, I tested the idea by plotting the frequency densities in real-time; this experiment demonstrated that EMG peaks contained frequencies all the way through the spectrum. Therefore I decided that it was unnecessary to remove the higher frequencies and might unnecessarily slow down the processor. Instead, I used a peak detector with a fixed threshold of .1 volt (a value of 600 out of 65536) and a running window size of 150. The algorithm I used from Labview ("peak detector.vi") picked peaks by fitting a quadratic polynomial to sequential groups of data points. The number of data points used in the fit was specified by the width control (set at 50). For each peak the quadratic fit was tested against the threshold level (600); peaks with heights lower than the threshold were ignored. Peaks were detected only after about width/2 data points have been processed beyond the peak location; this delay had serious implications for realtime processing. If more than one peak was found in a 50-sample window then one beat would be sent. Beats were sent as a value of 1; otherwise the default value was 0. A huge pitfall for this beat detection algorithm became multiple triggers, which had to be filtered out. They occurred when a beat overlapped over adjacent sample buffers. Solving this problem was not as simple as increasing the buffer size, however, since that would slow down the overall response time of the system.
In order to gather a measure of the force that was used to generate each beat, I wrote a filter to detect beat intensity. It calculated amplitudes for each peak. When one or more peaks were detected in a window of 150 samples, this filter reviewed all the samples and gave the value of the greatest one. This value was frequently used to determine the volume of the notes for a given beat. It also had another property which was unexpectedly useful for visual debugging: if the beat detector misfired too often, the graphical output of this filter allowed me to quickly diagnose what noise source was generating the problem. If it was generating high amplitude values, then that usually meant that the sensor was lifting off the surface of the skin; if it generated low amplitude values, that indicated that there was some ambient electronic noise that was raising the general signal level.
It also became critical to develop a beat timer, which simply time-stamped the moment of a beat to the millisecond. The clock had no reference, so this time-stamp was just used to find the inter-onset interval in milliseconds between consecutive beats. Each successive time value was sent over to the next computer, where C++ code in Visual Developer Studio subtracted its predecessor from it to determine the millisecond interval between beats. This clock would regularly ‘turn over’ from its 32-bit maximum value to 0, which would generate a very large number in the second computer; this was easily dealt with by removing very high values.
During the development process I explored the possibilities for using first and second derivatives of the different signals. The Labview algorithm performs a discrete differentiation of the sampled signal as follows:
d/dt X[i] = (½*dt)(x[i+1] – x[i-1])
In the end, these measurements were not used in the real-time system because they did not make the signals more practical to use or reveal any extra information in the data.
One method for generating a smoother version of the EMG signal that closely follows its contour is to use a time-domain envelope follower. This is done by first convolving an impulse response with a decaying exponential. The result is then scaled and convolved with the EMG signal. I experimented with many aspects of this algorithm and was somewhat successful in removing unwanted spikes in the signal, but did require a lot of processing time and slowed down the application. The algorithm might also have been
improved upon by convolving the result with a low-pass filter. It also might have benefited from the method developed by Raul Fernandez for blood volume pressure (BVP) readings.
Then I developed a group of filters to perform and graph
spectral analyses of the signals. These allowed me to observe the frequency-domain
behavior of the EMG signals in real-time. The first computes a fast
fourier transform after the signal is zero-meaned (but not rectified,
since the loss of negative frequencies takes out half of the DC component).
If the input sequence is not a power of 2 (which this was not), it calls
an efficient RDFT (real discrete Fourier transform) routine. The output
sequence is complex and is returned in a complex array of real and imaginary
parts. Another filter is a utility for graphing the power spectrum
of an input signal. The formula Labview uses for this function is:
A third filter computes a spectrogram (sometimes called spectrum), a graphical representation of a signal in which each vertical slice represents the magnitude (intensity) of each of its component frequencies. In a file called spectrogram.vi, I wrote a real-time utility that computed the real FFT and plotted the magnitudes of the frequencies (organized in discrete bins and sorted according to frequency value) versus time. The greater the magnitude, the more the bin became red in color; the lower the magnitude, the closer to white the bin appeared. The full frequency spectrum was divided into either 256, 512, or 1024 bins, depending on how much delay was tolerable. I also needed to know how quickly the frequencies changed, so that I could determine what window size to use. The uncertainty principle is based on the fact that there is a tradeoff between frequency resolution and time resolution; in order to get great frequency resolution you need to use a large chunk of time-sample, but that slows you down. Conversely, to have fine knowledge about exactly what time a frequency changed, you can only have a few frequencies. I decided to make the compromise in the frequency domain, since my interest was to generate feedback in real-time. Since I picked small window sizes to optimize speed, I lost quite a bit of low frequency component in my signal.
A variant of my spectrogram utility was weirdspectrogram.vi, in which bins were sorted according to the strength of the frequencies within them; that is, magnitudes would be stacked in order of their value. Instead of giving an image of the variety in the behavior across the spectrum, it gave a sense of the overall magnitude across all the frequencies. This gave a smoother result that was easier to read and resembled a Power Spectral Density graph. In the end it turned out that the frequency domain correlated well with the amplitude domain; at a moment of high amplitude (i.e. at the moment of a beat), the frequencies would all increase markedly in magnitude. There were few instances where frequencies didn’t respond together; sometimes this provided a helpful indicator of noise: if I could see the magnitudes of the frequencies were elevated around 60 Hz, I knew that my sources of electrical noise were too high and needed to be removed.
In my last experiment with frequency domain graphing utilities and filters I explored the Fast Hartley Transform. The Discrete Fast Hartley Transform is similar to the Discrete Fast Fourier Transform except that the Hartley is its own inverse and therefore more symmetric. Hartley essentially removes the imaginary part of the calculation so that the calculations are all real, and not complex. Although there are straightforward ways to remove the imaginary part of the DFT, the Hartley is slightly less involved to implement. In hartleyvsfftspectrogram.vi, I compared the relative performance of the DFT and the DHT, but didn’t see significant enough differences to pursue it further.
Finally, I wrote a filter to perform noise analysis on the various signals. The file called probabilitydensity.vi fits the frequency distribution of the signal to a Gaussian curve; if they match closely, then the noise on the sensor is white and minimal. Otherwise, if the two distributions don’t fit, then there is a large source of either signal or noise. If the sensor is placed properly on the surface of the skin and the muscle is not contracting, then it is safe to say that the anomaly is from a source of noise and should be removed. This utility was useful for detecting and removing sources of noise that were not detectable from a visual analysis of the time-domain component of the signals.