Introduction

The problem of detecting transient signals in a noisy environment has been studied for decades. In the Statistical Signal Detection Theory the presence of a useful signal in a background noise is normally cast as a hypothesis testing, where under the null hypothesis no signal is present. If the statistics of signal to be detected is not perfectly known, usually no uniformly most powerful (UMP) test exists, and the performance of a detector depends on signal representation [1]. In general, a signal representation can be model based and expansion based. When no appropriate model for the signal can be found, one usually resorts to a ``canonical set'' of basis function where the signal is projected, giving rise to expansion coefficients. We can think of these coefficients as of signal representation in a new coordinate system. Depending on the signal representation the detection problem can be formulated in various domains such as time domain, frequency domain, time-frequency domain, etc. In time-frequency domain, a signal is projected onto a basis of waveforms that are localized (subject to Heisenberg uncertainty principle) in both time and frequency, yielding a two-dimensional signal representation Tx(w, t) of a one-dimensional signal x(t). An example of this representation is a windowed Fourier Transform introduced by Gabor. A breakthrough in the theory of wavelets offered a powerful alternative to windowed Fourier Transform, where a one-dimensional signal x(t) is represented in time-scale domain by virtue of a wavelet transform Tx(a, b).

FAQ:

• What is wrong with the Fourier basis? Action potentials (APs) are short (transient) pulses, and belong to the family of time-limited signals. The Fourier Transform of a time-limited signal cannot be band-limited. Therefore, the representation of APs in the frequency domain will be spread and possibly highly overlapped with the spectrum of noise. Based on these observations, one can argue that the Fourier basis is suboptimal for representing APs.
• Why use the continuous wavelet transform (CWT)? Simply because we have a good idea about the duration of a typical AP. Namely, a large majority of APs last between 0.5 and 2.0 [ms]. Rather than using dyadic scales (typically associated with the discrete wavelet transform), where one hunts for signals at a wide range of scales, we can zoom in on the scales of interest--relevant scales. In addition, running the wavelet transform continuously, ensures that no AP will be missed due to its location. This property of the CWT is known as a translation-invariance.
• What types of wavelets? Naturally, the wavelets of compact support are preferred, because APs are signals of compact support. In particular, wavelets from biorthogonal family have bi-phasic shapes that are reminiscent of those of APs. Simply put, this allows a single AP to be represented with a few wavelet basis functions--sparse representation. Therefore, over the set of relevant scales, a bank of approximately matched filters is created, with the hope of getting a large wavelet coefficient once the wavelet function and AP are aligned.

This tutorial is an accompanying document to the computer code for detection of action potentials with the continuous wavelet transform. The details of the method can be found in [2], and the code is written in MATLAB^TM. The code consists of several functions and to run the main function Wavelet Toolbox needs to be installed. The main function is:

1. detect_spikes_wavelet.m

This is a stand-alone function for the detection of APs in noisy neural data. The auxiliary functions are:

1. get_score.m
2. lag_ts.m
3. lead_ts.m

and are provided for tutorial purposes only. Two synthetic data files are provided as well:

1. clean.mat
2. corrupted.mat

detect_spikes_wavelet.m returns the value of the vector TE of arrival times of detected spikes. The arrival times are given in samples and can be converted to true times if the sampling rate is known. This function has several input arguments and they will be explained in Section 2. Remark: if you were using the previous version of this function (previously posted at http://robotics.caltech.edu/~zoran/Research/detection.html), there are slight changes. Firstly, the set of relevant scales is now estimated more precisely, by means of a subroutine within detect_spikes_wavelet.m. Secondly, the number of relevant scales is now passed to the function as an input argument. Type help detect_spikes_wavelet in MATLAB^TMcommand prompt to learn more about this function.

get_score.m is an auxiliary function that is useful for comparing the true arrival times (if they are known) and the estimated arrival times given by TE. This function is useful since it returns the number of correctly detected APs, as well as the number of omissions and false alarms.

lag_ts.m and lead_ts.m are the auxiliary functions that are used by get_score.m.