- Email: [email protected]
Noise reduction in multichannel neural recordings using a new array wavelet denoising algorithm
Neurocomputing 38}40 (2001) 1687}1693
Noise reduction in multichannel neural recordings using a new array wavelet denoising algorithm Karim G. Oweiss *, David J. Anderson Electrical Engineering & Computer Science Department, University of Michigan, Ann Arbor, USA Biomedical Engineering Department, University of Michigan, Ann Arbor, USA Kresge Hearing Research Institute, University of Michigan, Ann Arbor, USA
Abstract We investigate a new technique for noise reduction in multichannel neural recordings based on the discrete wavelet transform. Starting with the denoising technique proposed by Donoho et al. (IEEE Trans. Inform. Theory 41 (1995) 613}627), we suggest a new thresholding method for the multiresolution decomposition of the multichannel data. The potential of this technique lies in the fact that thresholds at di!erent resolution levels of the wavelet transform are estimated spatially to account for signi"cant correlation of the wavelet coe$cients across channels. The method is applied to a simulated multichannel data as well as real silicon microprobe recordings obtained in our laboratory. Preliminary results show the ability of the technique to reduce both spatially correlated and uncorrelated noise components in the neural recordings. Results are compared to existing techniques and the overall performance is evaluated. 2001 Published by Elsevier Science B.V. Keywords: Multichannel recording; Silicon probes; Array processing; Wavelet denoising
1. Introduction With the recent advance in recording extracellular neural activity in the brain using multichannel microelectrode systems, the need for robust array signal processing algorithms became essential to take advantage of the high dimensionality in signal space. Extracellular neural recordings are considered event-related processes that can be decomposed into an invariant deterministic signal and an additive noise
* Corresponding author. E-mail addresses: [email protected] (K.G. Oweiss), [email protected] (D.J. Anderson). 0925-2312/01/$ - see front matter 2001 Published by Elsevier Science B.V. PII: S 0 9 2 5 - 2 3 1 2 ( 0 1 ) 0 0 5 3 3 - 1
1688
K.G. Oweiss, D.J. Anderson / Neurocomputing 38}40 (2001) 1687}1693
component. Noise reduction poses a main challenge in multichannel recordings since it obscures neural discharges from cells of interest. The noise sources in this environment can be divided into 2 main categories: 1. Thermal and electrical noise due to ampli"ers in the head stage of the associated circuitry which can be modeled as a spatio-temporally uncorrelated additive white Gaussian noise [1]. 2. High levels of background activity caused by neural cells far from the recording array. This type has 2 components, a spatially correlated noise process representing groups of cells close to the recording array but not enough for their "ring patterns to be detected and sorted by the recording array. The second component is a spatially uncorrelated noise process representing local but low level neural activity. Previous e!orts have been made to reduce background activity in neural recordings [2,6]. Of particular interest is the work done by Bierer et al. [2], to reduce the correlated part of the noise based on estimating the noise spatial covariance matrix and perform Wiener "ltering of the data. Gozani and Miller used software based linear and nonlinear "lters to optimally discriminate and classify neuronal action potential waveforms [6]. Almost all the e!orts to reduce noise in a similar recording environment dealt with reducing one type of noise but not the other(s). In some cases, reducing one type was achieved at the expense of boosting up the other component. This is particularly true in the case of reducing the correlated component. In this paper, we describe a new technique based on thresholding a discrete wavelet transform representation of the multichannel data. Preliminary results suggest that the technique has a strong capability to reduce both noise types simultaneously.
2. A brief background on denoising using the DWT The discrete wavelet transform (DWT) is a powerful signal processing tool to perform multi-resolution analysis. The decomposition consists of observing the signal at di!erent resolution levels and di!erent translations in time by bandpass "ltering using di!erent versions of a `mothera scaling and wavelet functions. Time localization is achieved by projecting the signal onto a subspace spanned by contracted versions of the wavelet function (high frequency), while frequency localization is achieved by implementing the same procedure using dilated (low frequency) versions of the mother scaling function. Thus, the strength of the obtained signal decomposition lies in using short high frequency basis functions and long low frequency ones to isolate di!erent characteristics of the signal such as discontinuities, etc. The mathematical formulation of the DWT will be omitted here and can be found elsewhere [8]. On the recent breakthroughs in noise reduction mechanisms, the work done by Donoho et al. [4,5], improves signi"cantly existing noise reduction techniques. The basic idea stems from the fact that the signal contributes few high amplitude coe$cients in some levels and some translations in time, while the noise contributes low
K.G. Oweiss, D.J. Anderson / Neurocomputing 38}40 (2001) 1687}1693
1689
amplitude coe$cients at all levels and all translations depending on the bandwidth and characteristics governing the noise process. By estimating the noise level at each resolution level and comparing it to the signal level, thresholding can be done by throwing away coe$cients attributed to noise while keeping coe$cients that contribute the signal. Di!erent methods have been proposed for threshold selection [4,5,7]. The main characteristic of all these existing methods is that they rely on threshold estimation using single channel data and the associated temporal characteristic of the noise process. As far as we know, none of them has been applied to the multichannel data case.
3. Proposed algorithm By observing the signi"cant spatial correlation between wavelet coe$cients at each decomposition level, we propose to perform a spatial weighted average of the decomposition structure. If x denotes the ith channel data vector of length N time G samples, then the weights are estimated at the jth level by performing an eigen analysis of the array wavelet covariance matrix, i.e. w "max (E[= =2]), (1) H H H H where = denotes the M;N/2H array wavelet transform matrix at the jth level H obtained by applying the N;N/2H orthonormal wavelet transformation matrix W to the snapshot of the M channels array data matrix X"[x x 2 x ]. This amounts + to selecting the eigen vector corresponding to the maximum eigen value obtained from the eigen decomposition of the array wavelet covariance matrix. Subsequently, thresholding takes place on the weighted average obtained and reconstruction using the inverse wavelet transformation matrix W\ follows next to obtain the denoised array data.
4. Simulation and experimental results 4.1. Simulation results The above-described algorithm was implemented in MATLAB. Three di!erent spike shapes were extracted from real data and the simulation was carried out by generating 3 spike trains from Poisson distributions with di!erent parameters. The spike trains were multiplied by a 4;3 mixing matrix that describes the amplitude distribution across 4 channels. This was varied arbitrarily to simulate the attenuation of multiple cell discharges recorded at variable distances from the electrode array. In the processing stage, we chose the Daubechies mother wavelet as an orthonormal basis with 6 decomposition levels. The Daubechies wavelet family have proven to maximize the smoothness of the mother wavelet function by maximizing the rate of decay of its Fourier transform [3].
1690
K.G. Oweiss, D.J. Anderson / Neurocomputing 38}40 (2001) 1687}1693
We report here the results obtained in the case of simulating signal#experimental noise (Fig. 1). The amplitude mixing matrix distribution of the signal and noise is also shown as a fraction of the full amplitude spikes. It is to be noted that spikes were scaled according to Fig. 1b while the noise across channels was used directly in the simulation without scaling.
4.2. Experimental multichannel recordings results Real multichannel data was acquired using Michigan type silicon substrate microprobes designed and fabricated in our laboratory [1]. We recorded data from Dorsal Cochlear Nucleus (DCN) of guinea pigs. A subset of the data collected using a 16channel microprobe is shown in Fig. 2 along with the output of the algorithm. Performance analysis was carried out and plotted to assess the SNR improvement in all cases (Fig. 2d).
Fig. 1. Signal#experimental noise simulation: (a) simulated tetrode data, (b) spike and noise amplitude distribution, (c) individual channel denoising, (d) array denoising.
K.G. Oweiss, D.J. Anderson / Neurocomputing 38}40 (2001) 1687}1693
1691
Fig. 2. Real tetrode recordings & performance analysis: (a) real tetrode recordings, (b) individual channel denoising, (c) array denoising, (d) SNR for: signal#uncorrelated noise, signal#real noise and real recordings.
5. Conclusion Some promising results have been presented to show the capability of the new array denoising algorithm in reducing noise sources in multichannel extracellular recordings. Spatial averaging is a simple and well-known technique to reduce cross correlation in coherent and correlated signal environments thereby permitting proper signal estimation and source separation.
1692
K.G. Oweiss, D.J. Anderson / Neurocomputing 38}40 (2001) 1687}1693
Wavelet coe$cients due to merely spatially uncorrelated noise will be small in magnitude and spatially uncorrelated as well [8]. Thus those will be "ltered out by the thresholding step. While this is the case in a single channel signal, the spatial averaging in the wavelet domain makes this "ltering operation more robust by reducing the variance in the amplitude of the noise coe$cients in the case where the uncorrelated noise variance is high on any one channel, thus improving the overall performance. In the case of the spatially correlated noise component, the wavelet domain structure is similar across channels. In otherwords, when two noise processes are spatially correlated, their wavelet decomposition exhibits high cross correlation but the coe$cient magnitudes will still be smaller than the signal coe$cients. This can be attributed to variations in distance between the recording array and populations of neurons from which this correlated background activity was observed. Thus, we expect the averaging technique to perform well in case where the SNR is high on any one channel. This was veri"ed experimentally but was not included here for the lack of space. It is thus expected that the increase in the number of channels in the recording device will directly lead to a more signi"cant improvement in the algorithm performance. Current work aims at assessing di!erent ways to re"ne and enhance the algorithm performance in di!erent cases. Exploring the correlation in the wavelet domain to assess di!erent spike shapes as a tool to sort di!erent neural cells is under investigation.
Acknowledgements This work was supported by NIH under contract number P 41 RR09754. Special thanks to Jamie Hetke, James Wiler and Steven Bierer for microprobe design and experimental setup.
References [1] D. J. Anderson, K. Naja", S. J. Tanghe, D. A. Evans, K. L. Levy, J. F. Hetke, X. Xue, J. J. Zappia, K. D. Wise, Batch-fabricated thin-"lm electrodes for stimulation of the central auditory system, IEEE Trans. BME (1989) 693}704. [2] S.M. Bierer, D.J. Anderson, Multi-channel spike detection and sorting using an array processing technique, NeuroComputing, Vol. 26}27, Elsevier, Amsterdam, 1999, pp. 947}956. [3] I. Daubechies, Ten Lectures on Wavelets, SIAM, Philadelphia, PA, 1992. [4] D.L. Donoho, De-noising by soft thresholding, IEEE Trans. Inform. Theory 41 (1995) 613}627. [5] D.L. Donoho, I. Johnstone, Ideal spatial adaptation via wavelet shrinkage, Biometrika 81, 425}455. [6] S.N. Gozani, J.P. Miller, Optimal discrimination and classi"cation of neuronal action potential waveforms from multiunit, multichannel recordings using software-based linear "lters, IEEE Trans. BME 41 (4) (1994). [7] I. Johnstone, B. Silverman Wavelet threshold estimators for data with correlated noise, J. Roy. Statist. Soc. (1996). [8] M. Vetterli, C. Herley, Wavelets and "lter banks: theory and design, IEEE Trans. Signal Process. 40 (9) (1992) 2207}2232.
K.G. Oweiss, D.J. Anderson / Neurocomputing 38}40 (2001) 1687}1693
1693
Karim G. Oweiss was born in Alexandria, Egypt in 1970. He obtained his B.Sc. and M.Sc. degrees from Alexandria University, Egypt in 1993 and 1995 respectively, all in Electrical Engineering. Since 1997 he has been with the Electrical Engineering & Computer Science Department at the University of Michigan, Ann Arbor as a research assistant where he is currently working towards his Ph.D. degree. His research interests are in the areas of array signal processing, multiresolution analysis and multichannel signal estimation with application to biomedical signal processing & neurophysiology.
David J. Anderson received the BSEE degree from Rensselaer Polytechnic Institute, Troy, NY, and the MS and Ph.D. degrees from the University of Wisconsin, Madison. After completing a postdoctoral traineeship with the Laboratory of Neurophysiology, University of Wisconsin Medical School, he joined the University of Michigan, Ann Arbor where he is now professor of Electrical Engineering and Computer Science, Biomedical Engineering and Otolaryngology. He is a fellow of the American Institute for Medical and Biological Engineering, a member of the IEEE, Neuroscience Society, Acoustical Society of America, the Association for Research in Otolaryngology, and the Barany Society.
Noise reduction in multichannel neural recordings using a new array wavelet denoising algorithm Karim G. Oweiss *, David J. Anderson Electrical Engineering & Computer Science Department, University of Michigan, Ann Arbor, USA Biomedical Engineering Department, University of Michigan, Ann Arbor, USA Kresge Hearing Research Institute, University of Michigan, Ann Arbor, USA
Abstract We investigate a new technique for noise reduction in multichannel neural recordings based on the discrete wavelet transform. Starting with the denoising technique proposed by Donoho et al. (IEEE Trans. Inform. Theory 41 (1995) 613}627), we suggest a new thresholding method for the multiresolution decomposition of the multichannel data. The potential of this technique lies in the fact that thresholds at di!erent resolution levels of the wavelet transform are estimated spatially to account for signi"cant correlation of the wavelet coe$cients across channels. The method is applied to a simulated multichannel data as well as real silicon microprobe recordings obtained in our laboratory. Preliminary results show the ability of the technique to reduce both spatially correlated and uncorrelated noise components in the neural recordings. Results are compared to existing techniques and the overall performance is evaluated. 2001 Published by Elsevier Science B.V. Keywords: Multichannel recording; Silicon probes; Array processing; Wavelet denoising
1. Introduction With the recent advance in recording extracellular neural activity in the brain using multichannel microelectrode systems, the need for robust array signal processing algorithms became essential to take advantage of the high dimensionality in signal space. Extracellular neural recordings are considered event-related processes that can be decomposed into an invariant deterministic signal and an additive noise
* Corresponding author. E-mail addresses: [email protected] (K.G. Oweiss), [email protected] (D.J. Anderson). 0925-2312/01/$ - see front matter 2001 Published by Elsevier Science B.V. PII: S 0 9 2 5 - 2 3 1 2 ( 0 1 ) 0 0 5 3 3 - 1
1688
K.G. Oweiss, D.J. Anderson / Neurocomputing 38}40 (2001) 1687}1693
component. Noise reduction poses a main challenge in multichannel recordings since it obscures neural discharges from cells of interest. The noise sources in this environment can be divided into 2 main categories: 1. Thermal and electrical noise due to ampli"ers in the head stage of the associated circuitry which can be modeled as a spatio-temporally uncorrelated additive white Gaussian noise [1]. 2. High levels of background activity caused by neural cells far from the recording array. This type has 2 components, a spatially correlated noise process representing groups of cells close to the recording array but not enough for their "ring patterns to be detected and sorted by the recording array. The second component is a spatially uncorrelated noise process representing local but low level neural activity. Previous e!orts have been made to reduce background activity in neural recordings [2,6]. Of particular interest is the work done by Bierer et al. [2], to reduce the correlated part of the noise based on estimating the noise spatial covariance matrix and perform Wiener "ltering of the data. Gozani and Miller used software based linear and nonlinear "lters to optimally discriminate and classify neuronal action potential waveforms [6]. Almost all the e!orts to reduce noise in a similar recording environment dealt with reducing one type of noise but not the other(s). In some cases, reducing one type was achieved at the expense of boosting up the other component. This is particularly true in the case of reducing the correlated component. In this paper, we describe a new technique based on thresholding a discrete wavelet transform representation of the multichannel data. Preliminary results suggest that the technique has a strong capability to reduce both noise types simultaneously.
2. A brief background on denoising using the DWT The discrete wavelet transform (DWT) is a powerful signal processing tool to perform multi-resolution analysis. The decomposition consists of observing the signal at di!erent resolution levels and di!erent translations in time by bandpass "ltering using di!erent versions of a `mothera scaling and wavelet functions. Time localization is achieved by projecting the signal onto a subspace spanned by contracted versions of the wavelet function (high frequency), while frequency localization is achieved by implementing the same procedure using dilated (low frequency) versions of the mother scaling function. Thus, the strength of the obtained signal decomposition lies in using short high frequency basis functions and long low frequency ones to isolate di!erent characteristics of the signal such as discontinuities, etc. The mathematical formulation of the DWT will be omitted here and can be found elsewhere [8]. On the recent breakthroughs in noise reduction mechanisms, the work done by Donoho et al. [4,5], improves signi"cantly existing noise reduction techniques. The basic idea stems from the fact that the signal contributes few high amplitude coe$cients in some levels and some translations in time, while the noise contributes low
K.G. Oweiss, D.J. Anderson / Neurocomputing 38}40 (2001) 1687}1693
1689
amplitude coe$cients at all levels and all translations depending on the bandwidth and characteristics governing the noise process. By estimating the noise level at each resolution level and comparing it to the signal level, thresholding can be done by throwing away coe$cients attributed to noise while keeping coe$cients that contribute the signal. Di!erent methods have been proposed for threshold selection [4,5,7]. The main characteristic of all these existing methods is that they rely on threshold estimation using single channel data and the associated temporal characteristic of the noise process. As far as we know, none of them has been applied to the multichannel data case.
3. Proposed algorithm By observing the signi"cant spatial correlation between wavelet coe$cients at each decomposition level, we propose to perform a spatial weighted average of the decomposition structure. If x denotes the ith channel data vector of length N time G samples, then the weights are estimated at the jth level by performing an eigen analysis of the array wavelet covariance matrix, i.e. w "max (E[= =2]), (1) H H H H where = denotes the M;N/2H array wavelet transform matrix at the jth level H obtained by applying the N;N/2H orthonormal wavelet transformation matrix W to the snapshot of the M channels array data matrix X"[x x 2 x ]. This amounts + to selecting the eigen vector corresponding to the maximum eigen value obtained from the eigen decomposition of the array wavelet covariance matrix. Subsequently, thresholding takes place on the weighted average obtained and reconstruction using the inverse wavelet transformation matrix W\ follows next to obtain the denoised array data.
4. Simulation and experimental results 4.1. Simulation results The above-described algorithm was implemented in MATLAB. Three di!erent spike shapes were extracted from real data and the simulation was carried out by generating 3 spike trains from Poisson distributions with di!erent parameters. The spike trains were multiplied by a 4;3 mixing matrix that describes the amplitude distribution across 4 channels. This was varied arbitrarily to simulate the attenuation of multiple cell discharges recorded at variable distances from the electrode array. In the processing stage, we chose the Daubechies mother wavelet as an orthonormal basis with 6 decomposition levels. The Daubechies wavelet family have proven to maximize the smoothness of the mother wavelet function by maximizing the rate of decay of its Fourier transform [3].
1690
K.G. Oweiss, D.J. Anderson / Neurocomputing 38}40 (2001) 1687}1693
We report here the results obtained in the case of simulating signal#experimental noise (Fig. 1). The amplitude mixing matrix distribution of the signal and noise is also shown as a fraction of the full amplitude spikes. It is to be noted that spikes were scaled according to Fig. 1b while the noise across channels was used directly in the simulation without scaling.
4.2. Experimental multichannel recordings results Real multichannel data was acquired using Michigan type silicon substrate microprobes designed and fabricated in our laboratory [1]. We recorded data from Dorsal Cochlear Nucleus (DCN) of guinea pigs. A subset of the data collected using a 16channel microprobe is shown in Fig. 2 along with the output of the algorithm. Performance analysis was carried out and plotted to assess the SNR improvement in all cases (Fig. 2d).
Fig. 1. Signal#experimental noise simulation: (a) simulated tetrode data, (b) spike and noise amplitude distribution, (c) individual channel denoising, (d) array denoising.
K.G. Oweiss, D.J. Anderson / Neurocomputing 38}40 (2001) 1687}1693
1691
Fig. 2. Real tetrode recordings & performance analysis: (a) real tetrode recordings, (b) individual channel denoising, (c) array denoising, (d) SNR for: signal#uncorrelated noise, signal#real noise and real recordings.
5. Conclusion Some promising results have been presented to show the capability of the new array denoising algorithm in reducing noise sources in multichannel extracellular recordings. Spatial averaging is a simple and well-known technique to reduce cross correlation in coherent and correlated signal environments thereby permitting proper signal estimation and source separation.
1692
K.G. Oweiss, D.J. Anderson / Neurocomputing 38}40 (2001) 1687}1693
Wavelet coe$cients due to merely spatially uncorrelated noise will be small in magnitude and spatially uncorrelated as well [8]. Thus those will be "ltered out by the thresholding step. While this is the case in a single channel signal, the spatial averaging in the wavelet domain makes this "ltering operation more robust by reducing the variance in the amplitude of the noise coe$cients in the case where the uncorrelated noise variance is high on any one channel, thus improving the overall performance. In the case of the spatially correlated noise component, the wavelet domain structure is similar across channels. In otherwords, when two noise processes are spatially correlated, their wavelet decomposition exhibits high cross correlation but the coe$cient magnitudes will still be smaller than the signal coe$cients. This can be attributed to variations in distance between the recording array and populations of neurons from which this correlated background activity was observed. Thus, we expect the averaging technique to perform well in case where the SNR is high on any one channel. This was veri"ed experimentally but was not included here for the lack of space. It is thus expected that the increase in the number of channels in the recording device will directly lead to a more signi"cant improvement in the algorithm performance. Current work aims at assessing di!erent ways to re"ne and enhance the algorithm performance in di!erent cases. Exploring the correlation in the wavelet domain to assess di!erent spike shapes as a tool to sort di!erent neural cells is under investigation.
Acknowledgements This work was supported by NIH under contract number P 41 RR09754. Special thanks to Jamie Hetke, James Wiler and Steven Bierer for microprobe design and experimental setup.
References [1] D. J. Anderson, K. Naja", S. J. Tanghe, D. A. Evans, K. L. Levy, J. F. Hetke, X. Xue, J. J. Zappia, K. D. Wise, Batch-fabricated thin-"lm electrodes for stimulation of the central auditory system, IEEE Trans. BME (1989) 693}704. [2] S.M. Bierer, D.J. Anderson, Multi-channel spike detection and sorting using an array processing technique, NeuroComputing, Vol. 26}27, Elsevier, Amsterdam, 1999, pp. 947}956. [3] I. Daubechies, Ten Lectures on Wavelets, SIAM, Philadelphia, PA, 1992. [4] D.L. Donoho, De-noising by soft thresholding, IEEE Trans. Inform. Theory 41 (1995) 613}627. [5] D.L. Donoho, I. Johnstone, Ideal spatial adaptation via wavelet shrinkage, Biometrika 81, 425}455. [6] S.N. Gozani, J.P. Miller, Optimal discrimination and classi"cation of neuronal action potential waveforms from multiunit, multichannel recordings using software-based linear "lters, IEEE Trans. BME 41 (4) (1994). [7] I. Johnstone, B. Silverman Wavelet threshold estimators for data with correlated noise, J. Roy. Statist. Soc. (1996). [8] M. Vetterli, C. Herley, Wavelets and "lter banks: theory and design, IEEE Trans. Signal Process. 40 (9) (1992) 2207}2232.
K.G. Oweiss, D.J. Anderson / Neurocomputing 38}40 (2001) 1687}1693
1693
Karim G. Oweiss was born in Alexandria, Egypt in 1970. He obtained his B.Sc. and M.Sc. degrees from Alexandria University, Egypt in 1993 and 1995 respectively, all in Electrical Engineering. Since 1997 he has been with the Electrical Engineering & Computer Science Department at the University of Michigan, Ann Arbor as a research assistant where he is currently working towards his Ph.D. degree. His research interests are in the areas of array signal processing, multiresolution analysis and multichannel signal estimation with application to biomedical signal processing & neurophysiology.
David J. Anderson received the BSEE degree from Rensselaer Polytechnic Institute, Troy, NY, and the MS and Ph.D. degrees from the University of Wisconsin, Madison. After completing a postdoctoral traineeship with the Laboratory of Neurophysiology, University of Wisconsin Medical School, he joined the University of Michigan, Ann Arbor where he is now professor of Electrical Engineering and Computer Science, Biomedical Engineering and Otolaryngology. He is a fellow of the American Institute for Medical and Biological Engineering, a member of the IEEE, Neuroscience Society, Acoustical Society of America, the Association for Research in Otolaryngology, and the Barany Society.