[Introduction] [Signal arithmetic] [Signals and noise] [Smoothing] [Differentiation] [Peak Sharpening] [Harmonic analysis] [Fourier convolution] [Fourier deconvolution] [Fourier filter] [Wavelets] [Peak area measurement] [Linear Least Squares] [Multicomponent Spectroscopy] [Iterative Curve Fitting] [Hyperlinear quantitative absorption spectrophotometry] [Appendix and Case Studies] [Peak Finding and Measurement] [iPeak] [iSignal] [Peak Fitters] [iFilter] [iPower] [List of downloadable software] [Interactive tools]

*S**pikes*, narrow pulses with a width
of only one or a few points, are sometimes encountered
in signals as a result of an electronic "glitch" or
stray pickup from nearby equipment, and they can easily
be eliminated by the use of a "median"
filter. But it is possible that in some
experiments the spikes *themselves* might be the
important part of the signal and that it is required to
count or measure them. This situation was recently
encountered in a research application, and it opens up
some interesting twists on the usual procedures.

As a demonstration, the Matlab/Octave script SpikeDemo1.m creates a waveform (top panel of figure below) in which a series of spikes are randomly distributed in time, contaminated by two types of noise: white noise and a large-amplitude oscillatory interference simulated by a swept-frequency sine wave. The objective is to count the spikes and locate their position on the x (time) axis. Direct application of findpeaks or iPeak to the raw signal does not work well.

A single-point spike, called a *delta function* in mathematics,
has a power spectrum that is
*flat*; that is, it has *equal power at
all frequencies*, just like white noise. But
the oscillatory interference in this case is located in
a *specific range
of frequencies*, which opens some
interesting possibilities. One approach would be to use
a Fourier filter, for example
a notch or band-reject filter, to remove the
troublesome oscillations selectively. But if the
objective of the measurement is only to count the spikes
and measure their times, a simpler approach would be to
(1) compute the second derivative (which
greatly amplifies the spikes relative to the
oscillations), (2) smooth
the result (to limit the white noise amplification
caused by differentiation), then (3) take the absolute
value (to yield positive-pointing peaks). This can
be done in a *single line* of nested
Matlab/Octave code:

y1=abs(fastsmooth((deriv2(y)).^2,3,2));

The result, shown the lower panel of the figure on the left above, is an almost complete extraction of the spikes, which can then be counted with findpeaksG.m or peakstats.m or iPeak.m:

P=ipeak([x;y1],0,0.1,2e-005,1,3,3,0.2,0);

The second example, SpikeDemo2.m, is similar except that in
this case the oscillatory interference is caused by *two* fixed-frequency
sine waves at a higher frequency, which completely
obscure the spikes in the raw signal (top panel of the
left figure below). In the power
spectrum (bottom panel, in red), the oscillatory
interference shows as two sharp peaks that dominate the
spectrum and reach to y=10^{6}, whereas the spikes show as
the much lower broad flat plateau at about y=10. In this
case, use can be made of an interesting property of
sliding-average smooths, such as the boxcar, triangular,
and Gaussian smooths;
their frequency responses exhibit a series of deep cusps at frequencies
that are inversely proportional to their filter widths.
So this opens up the possibility of suppressing specific
frequencies of oscillatory interference by adjusting the
filter widths until the cusps occur at or near the
frequency of the oscillations. Since the signal in this
cases are spikes that have a flat power spectrum, they
are simply smoothed by this operation, which will reduce
their heights and increase their widths, but will have
little or no effect on their number or x-axis positions.
In this case a 9 point pseudo-Gaussian smooth puts the
first (lowest frequency) cusp right in between the two oscillatory frequencies.

In the figure on
the right, you can see the effect of applying this filter; the
spikes, which were *not even visible* in the original
signal, are now cleanly extracted (upper panel), and you can
see in the power spectrum (right lower panel, in red) that the
two sharp peaks of oscillatory interference is *reduced by
about a factor of about*
*1,000,000**!* This operation
can be performed by a single command-line function,
adjusting the smooth width (the second input argument,
here a 9) by trial and error to minimize the oscillatory interference:

y1=fastsmooth(y,9,3);

(If the interference varies
substantially in frequency across the signal, you could
use a segmented smooth rather
than the standard fastsmooth). Alternatively, the
segmented
Fourier spectrum could be used to visualize
this signal, and a Fourier filter in "notch"
mode could be employed to specifically eliminate
the interfering frequencies.

The extracted peaks can then be counted with any of the peak finding functions, such as:

P=findpeaksG(x,y1,2e-005,0.01,2,5,3);

or

P=findpeaksplot(x,y1,2e-005,0.01,2,5,3);

or

PS=peakstats(x,y1,2e-005,0.01,2,5,3,1);

The simple script "iSignalDeltaTest" demonstrates the power spectrum of the smoothing and differentiation functions of iSignal by a applying them to a delta function. Change the smooth type, smooth width, and derivative order and other functions in order to see how the power spectrum changes.

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.