Appendix G: Measuring spikes


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));


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=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.


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.