Appendix Y: Real-time signal
processing
All of the signal processing techniques
covered so far make the assumption that you have acquired and
stored the data in computer memory before beginning processing. In
some cases, however, it is necessary to do the signal processing
in "real time", that is, point-by-point as the data are acquired
from the sensor or instrument. That requires some modification of
the software, but the main ideas still apply. In this section we
will look at ways to perform real-time data plotting, smoothing,
differentiation, peak detection, harmonic analysis (frequency
spectra), and Fourier filtering. Because the details of the data
acquisition itself varies with each individual experimenter and
instrumental setup, these demonstration scripts will simulate real-time
data so that you can run them immediately on
your computer without additional hardware. I'll do this in either
of two ways:
(a) by using mouse-clicks to generate
each data point, using Matlab's "ginput" function, or
(b) by pre-calculating simulated data and then
accessing it point-by-point in a loop.
The first method is illustrated by the
simple script realtime.m.
When you run this script, it displays a graphical coordinate
system; position your mouse pointer along the y (vertical)
axis and click to enter data points as you move the mouse pointer
up and down. The "ginput" function waits for each click of the
mouse button, then the program records the y coordinate
position and counts the number of clicks. Data points are assigned
to the vector y (line 17), plotted on the graph as black
points (line 18), and print out in the command window (line 19). The script realtimeplotautoscale.m is an
expanded version that changes the graph scale as the data come in.
If the number of data points exceeds 20 ('maxdisplay'), the x
axis maximum is re-scaled to twice that (line 32). If the
data amplitude equals or exceeds ('maxy'), the y axis is
re-scaled to 1.1 times the data amplitude (line 36).
The script realtimeplotautoscale2.m
uses the second method to simulate real-time
data, using pre-calculated data (loaded from
disk in line 13) that is accessed point-by-point in lines 25 and
26 (animation shown on the right).
The script realtimeplotdatedtime.m demonstrates
one way to use Matlab's 'clock' function to record the data and
time of each data point that is acquired by clicking. You could
also have the computer control the time of data acquisition by
reading the clock in a loop until the desired time and date
arrives, then take a data point. Of course, a Windows machine is
not ideal for high-speed, precisely-timed data acquisition,
because there are typically so many interrupts and other processes
going on in the background, but it's adequate for low-speed
applications. For higher speeds, specialized
hardware and software is
available.
Smoothing. The script RealTimeSmoothTest.m
demonstrates real-time smoothing,
plotting the raw unsmoothed data as a black line and the smoothed
data in red. In this case the script pre-calculates
simulated data in line 28 and then accesses the data
point-by-point in the processing 'for' loop (lines 30-51).
The total number of data points is controlled by 'maxx' in line 17
(initially set to 1000) and the smooth width (in points) is
controlled by 'SmoothWidth' in line 20. (To do this with real time
data from your sensor, comment out line 29 and replace line 32
with the code that acquires one data point from your sensor).
As you can see in the
screen animation on the left, the smoothed
data (in red) is delayed compared to the raw data, because a
smoothed data point can not be computed until a number of data
points equal to the smooth width has been acquired - 21 points, or
about 0.5 seconds in this example. (However, knowing the smooth
width, you can correct the recorded y-axis positions of signal
features, such as maxima, minima, peaks, or inflection points).
This particular example implements a sliding
average smooth, but other smooth shapes can be implemented
simply by uncommenting line 24 (rectangular smooth), 25 (triangular smooth),
or 26 (Gaussian smooth), which requires
that the functions 'triangle'
and 'gaussian' be
in the Matlab/Octave path.
A practical application of a sliding average smooth like this is
in a control system where a noisy signal turns on a valve, switch
or alarm signal whenever the signal exceeds a certain value; an
example is shown in the figure on the right, where the
threshold value is 0.5 and the signal is so noisy that smoothing
is required to prevent the signal from prematurely triggering the
control. Too much smoothing, however, will cause an unacceptable
delay in operation.
On a standard desktop PC (Intel Core i5 3 Ghz) running
Windows 10 home, the smooth operation adds about 2 microseconds per data point to
the data acquisition time (without plotting, PlottingOn=0) and 20 milliseconds
per point (50 Hz) with point-by-point plotting (PlottingOn=1). With plotting off, the script acquires, smooths, and stores the smoothed data in the
variable "sy" in real time, then plots the data only after data
acquisition is complete, which is much faster than plotting in
real time.
Differentiation. The script RealTimeSmoothFirstDerivative.m
demonstrates real-time smoothed differentiation, using a simple
adjacent-difference algorithm (line 47) and plotting the raw data
as a black line and the first derivative data in red. The
demonstration script RealTimeSmoothSecondDerivative.m computes the smoothed second derivative
by using a central difference algorithm (line 47). Both of
these scripts pre-calculate the simulated data in line 28 and
then accesses the data point-by-point in the processing
loop (lines 31-52). In both cases the
maximum number of points is set in
line 17 and the smooth width is set in line
20. Again,
the derivatives are delayed
compared to the original
signal. Any derivative order
can be calculated this way
using the derivative
coefficients in the
Matlab/Octave derivative
functions listed on Differentiation.html#Matlab.
Peak detection. The little script realtimepeak.m
demonstrates simple real-time peak
detection
based on
derivative zero-crossing, using mouse clicks to simulate data. Each
time your mouse clicks form a peak (that is, go up and then down
again), the program will register and label the peak on the
graph and print out its n and y values.
Peak detected at n=13 and y=7.836
Peak detected at n=26 and y=1.707
In this case, a peak
is defined as any data point that has lower amplitude points
adjacent to it on both sides, which is determined by the nested
'for' loops in lines 31-36. Of course, a peak can not be
registered until the point following the peak is recorded,
because there is no way to predict ahead of time whether that
point will be lower or higher than the previous point.
If the data are noisy, it's better to smooth the data stream
before detecting the peaks, which is exactly what RealTimeSmoothedPeakDetection.m
does, which reduces the chance of false peaks due
to random noise but has
the disadvantage of delaying the peak detection further. Even
better, the script RealTimeSmoothedPeakDetectionGauss.m
uses the technique described on PeakFindingandMeasurement.htm#findpeaks;
it locates the positive peaks in a noisy data set that rise
above a set amplitude threshold ("AmpThreshold" in line 55),
performs a least-squares curve-fit of a Gaussian
function to the top part of the raw data peak (in line
58), identifies each peak (line 59), computes the position,
height, and width (FWHM) of each peak from that least-squares
fit, and prints out each peak found in the command window. The peak parameters are measured on the raw
data, so they are not distorted by smoothing. If you
look at the animation on the right, you can see that the "peak"
label appears next to each detected peak just a fraction of a
second after the top of the peak, but the actual peak times
listed are not delayed. In this example, the actual peak times
are x=500, 1000, 1100, 1200, 1400. (Also note that the first
visible peak, at x=300, is not detected because it falls below
the amplitude threshold, which is 0.1
in this case).
Peak detected at x=500.1705, y=0.42004, width=
61.7559
Peak detected at x=1000.0749, y=0.18477, width= 61.8195
Peak detected at x=1100.033, y=1.2817, width= 60.1692
Peak detected at x=1199.8493, y=0.36407, width= 63.8316
Peak detected at x=1400.1473, y=0.26134, width= 58.9345
The script
additionally numbers
the peaks and saves the
peak parameters of all the
peaks in a matrix called
PeakTable, which you can
interrogate as each peak is
encountered if you are looking
for particular peak
patterns. See PeakFindingandMeasurement.htm#UsingP
for some ideas on using
Matlab/Octave notation and
functions to do this.
Peak sharpening. The script RealTimePeakSharpening.m
demonstrates real-time peak
sharpening using the second derivative technique. It uses pre-calculated simulated
data in line 30 and then accesses the data
point-by-point in the processing loop (lines 33-55).
In both cases the maximum number of points is set in
line 17, the smooth width is set in line
20, and the weighting factor (K1) is set in line 21. In this
example on the left, the smooth width
is 101 points, which accounts for the
the delay in the sharpened peak
compared to the original (about 1.12
seconds without real-time plotting and
3.8 second with real-time
plotting).
Real-Time Frequency
Spectrum. The script RealTimeFrequencySpectrumWindow.m
computes and plots the Fourier
frequency spectrum of a signal. Like
the scripts above, it loads the
simulated real-time data from a
".mat file" and then accesses that
data point-by-point in the
processing loop. A critical variable
in this case is "WindowWidth", the
number of data points taken to
compute each frequency spectrum. The
larger this number, the fewer the
number of spectra that will be
generated, but the higher will be
the frequency resolution. On a
standard desktop PC (Intel Core i5 3
Ghz running Windows 10 home), this
script generates about 50 spectra
per second with an average data rate
(points per seconds) of about 50,000
Hz. Smaller spectra (i.e. lower
values of WindowWidth) generate
proportionally lower average data
rates (because the signal stream is
interrupted more often to calculate
and graph a spectrum).
If the
data stream is an audio
signal, it's also possible
to play the sound through
the computer's sound system
synchronized with the
display of the frequency
spectra; to do this, set PlaySound=1. Each
segment of the signal is
played, just before the
spectrum of that segment
is displayed, as shown on
the right. The sound
reproduction will not be
not perfect, because of
the slight delay while the
computer computes and
displays the spectrum
before going on to the
next segment. In this
demonstration script, the
data file is in fact an
audio recording of an
8-second excerpt of the
'Hallelujah Chorus' from
Handel's Messiah with a
sampling rate of 8192 Hz,
which is included in the
Matlab distribution. The
figure above shows one of
the 70 spectra generated
with a WindowWidth of
1024. You can adjust the
argument of the 'pause'
function for your computer
to minimize this problem
and to make the sound play
at the correct pitch.
Real-Time Fourier Filter.
The script RealTimeFourierFilter.m
is a demonstration of a real-time
Fourier filter. It pre-computes a
simulated signal starting in line 38,
then access the data point-by-point
(line 56, 57), and divides up the data
stream into segments to compute
each filtered
section.
"WindowWidth" (line 55) is the number
of data points in each segment. The larger this number, the
fewer the number of segments that will
be generated, but the higher will be
the frequency resolution within each
segment. On a standard desktop PC
(Intel Core i5 3 Ghz running Windows
10 home), with a window width of 1000
points, this script generates about 35
filtered segments per second with an
average data rate (points per second)
of about 34,000 Hz. Smaller segments
(i.e. lower values of WindowWidth)
generate proportionally lower average
data rate (because the signal stream
is interrupted more often to calculate
and graph the filtered spectrum). The
result of applying the filter to each
segment is displayed in real time
during the data acquisition, and then
at the end the script compares the
entire raw signal input to the
reassembled filtered output, shown the
the figure on the left.
In this demonstration, a bandpass
filter is used to detect a 500 Hz sine
wave (the
frequency is controlled by 'f' in line
28) that occurs in the middle third of
a very noisy signal (line 32), from
about 0.7 sec to 1.3 sec. The 500
Hz sine wave
is so weak it
can not be
seen at all in
the raw signal
(upper panel
of the figure
on the left).
The
filter center frequency
(CenterFrequency) and width
(FilterWidth) are set in lines 46 and
47.
Real-time classical least squares.
The classical
least squares technique can be
applied in real time to 2D
chromatography with array detectors
that can acquire a complete spectrum
multiple times per second over the
entire chromatogram. This is explored
in "Spectroscopy
and chromatography combined:
time-resolved Classical Least
Squares".
To apply any of these examples
to real-time data from your sensor or
instrument, you need only the main
processing 'for' loop, replacing the
first lines after the 'for' statement
with a call to a function that
acquires a single point of raw data
and assigns it to y(n). If you don't
need the data plotted out
point-by-point in real time, you can
speed things up greatly by removing
the "drawnow" statement at the end of
the 'for' loop or by removing all the
plotting code.
In the examples here, the output of the
processing operation is used to plot or to print out the processed
data point-by-point, but of course it could also be used as the
input to another processing function or to a digital-to-analog converter or simply to
trigger an alarm if certain specified results are obtained (e.g.
if the signal exceeds a certain value for a specified length of
time, or if a peak is detected at a specified position or height,
etc).
This page is part of "A Pragmatic Introduction to Signal
Processing", created and maintained by Prof. Tom O'Haver ,
Department of Chemistry and Biochemistry, The University of Maryland
at College Park. Comments, suggestions and questions should be
directed to Prof. O'Haver at toh@umd.edu. Updated July, 2022.