Spikes, 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));
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=106, 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.
