This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TNSRE.2016.2525829, IEEE Transactions on Neural Systems and Rehabilitation Engineering
1
Automatic parametrization of somatosensory evoked potentials with chirp modeling Eero Väyrynen, Kai Noponen, Ashwati Vipin, Member, IEEE, Thow Xin Yuan, Hasan Al-Nashash, Senior Member, IEEE, Jukka Kortelainen, Member, IEEE, and Angelo All Abstract—In this paper, an approach using polynomial phase chirp signals to model somatosensory evoked potentials (SEPs) is proposed. SEP waveforms are assumed as impulses undergoing group velocity dispersion while propagating along a multipath neural connection. Mathematical analysis of pulse dispersion resulting in chirp signals is performed. An automatic parameterization of SEPs is proposed using chirp models. A Particle Swarm Optimization algorithm is used to optimize the model parameters. Features describing the latencies and amplitudes of SEPs are automatically derived. Rat model is then used to evaluate the automatic parameterization of SEPs in two experimental cases, i.e. anesthesia level and spinal cord injury (SCI). Experimental results show that chirp based model parameters and the derived SEP features are significant in describing both anesthesia level and SCI changes. The proposed automatic optimization based approach for extracting chirp parameters offers potential for detailed SEP analysis in future studies. The method implementation in Matlab technical computing language is provided online. Index Terms—Anesthesia, biological system modeling, parameter estimation, particle swarm optimization, spinal cord injury
I. INTRODUCTION
S
OMATOSENSORY evoked potentials (SEPs) represent a widely used neurophysiological method when examining the functional integrity of ascending sensory pathways. Today, they are clinically used for various diagnostic purposes such as the assessment of brain injury after cardiac arrest and intraoperative monitoring during high risk surgeries [1-3]. For studies involving laboratory animals, SEPs provide a convenient minimally-invasive method to quantitatively assess the functionality of various portions of the somatosensory Manuscript received July 21, 2015. This work was supported in part by grants 2013-MSCRFII-0109-00 from the Maryland Stem Cell Research Fund as well as R-175-000-121-133, R-175-000-121-733, R-175-000-122-112, and R-711-201-026-133 from the National University of Singapore. E. Väyrynen, K. Noponen, and J. Kortelainen are with the Physiological Signal Analysis Team, Center for Machine Vision and Signal Analysis, University of Oulu, Finland. A. Vipin and X. Thow are with the Singapore Institute for Neurotechnology (SINAPSE), National University of Singapore, Singapore. H. Al-Nashash is with the Department of Electrical Engineering, American University of Sharjah, United Arab Emirates. A. All is with the Department of Orthopedic Surgery and Biomedical Engineering, National University of Singapore, Singapore and the Department of Neurology and Biomedical Engineering, Johns Hopkins University, USA. Address correspondence to Jukka Kortelainen (e-mail: [email protected]).
pathways. However, the analysis of SEPs still relies mainly on the determination of simple parameters such as amplitudes and latencies from the averaged signal waveform [4] even though several more sophisticated measures for the signal quantification have been proposed [5,6]. Amplitude and/or frequency modulated non-stationary signals are common in many engineering problems (e.g. radar, communications, and biomedical imaging) as well as in natural signals (e.g. echolocation with natural sonar). The parameter estimation of chirps and in general polynomialphase signals has therefore been widely studied [7,8]. The estimation is generally based on models aiming for continuous chirp signal analysis using methods such as Maximum Likelihood (ML) estimators [6], polynomial phase transform, [9], variations of the higher order ambiguity functions [10], polynomial Wigner-Ville distributions [11], and cubic phase functions [12]. Furthermore, hybrid methods [13] to increase performance and phase unwrapping approaches [14] have recently been studied. These methods yield theoretically very effective parameter estimates for continuous chirp signals with Signal-to-Noise Ratio (SNR) over 0 dB, i.e. approach the Cramer-Rao lower bound for Mean Squared Error (MSE). For SEP measurements, however, the SNR is usually very low, e.g. -15 – 0 dB. In order to counteract the poor SNR, multiple trial average measurements are typically needed to analyze the signals. A growing interest has been for single trial measures of SEPs [15] to enable faster evaluation, e.g. important in spinal surgery, and also access the variability of parameters lost in the averaging processes with the expense of SNR. To work in the low SNR method using spatial measures such as slope analysis and morphological features [4,16], spectral domain features such as spectral coherence [5], wavelet and/or independent component analysis [15], and matching pursuit time-frequency analysis [17] have been proposed. In this paper, an analysis approach that incorporates a chirp based modeling in the normal SEP analysis structure is considered. In order to achieve robust operation, the aim is not to perform perfect reconstructions of the SEP signals with chirps but rather to find and extract best fitting polynomialphase waveform components, and associated parameters, of the SEPs that are consistent with a pulse propagation and dispersion model (i.e. imposes specific time-frequency transition characteristics). First, we discuss the mathematical implications of pulse propagation in dispersive channels
1534-4320 (c) 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TNSRE.2016.2525829, IEEE Transactions on Neural Systems and Rehabilitation Engineering
2 (Section II). We assume that the aggregate cortical signal in the sensor, caused by impulse-like stimuli to the peripheral nerve, can be explained by pulses propagating in a multipath dispersive channel with amplitude dampening resulting in chirp signals. Second, a chirp model for SEP is constructed and a Particle Swarm Optimization (PSO) model fitting technique is used in order to provide automatic measures for the chirp model parameters and also to estimate a set of more traditional amplitude/latency parameters (Section IV; A, B, and C). Third, the resulting parametrization is then evaluated with experimental rat model data in which SEPs are affected by two different factors, anesthetic level and spinal cord injury (SCI) (Section IV; D). Finally, statistical testing is used to identify significant parameters and features in both of the experimental cases (Section IV; E). II. PULSE DISPERSION AND CHIRP SIGNALS Empirically, the general functional form of the SEPs resembles decaying sine waves with varying frequency. Generally, they all seem to be rather closely modelled by dampened chirp signals which are a certain class of frequency modulated sinusoids with amplitude envelope decay 𝑥(𝑡) = 𝐴(𝑡)𝑒 𝑗φ(𝑡) ,
(1)
where 𝐴(𝑡) is the envelope function and 𝜑(𝑡) is the phase function of the model. While the individual neural transmissions are not chirp signals (rather, the neural transmissions more closely resemble to short pulses), the aggregated electric fields in the EEG sensors can be explained by a multipath propagation of pulses undergoing dispersion and amplitude dampening. A mathematical evaluation is provided next to show exactly the connection of dispersed pulses and chirp signals. A. Pulse Dispersion Pulse propagation in dispersive media is a process in which different spectral components of the pulse travel at different velocities. This may happen, for example, due to multipath propagation through channels having different passbands. In the frequency domain, this dispersion is seen as phase shifts of the different spectral components of the pulse. Let us denote the dispersive phase shift at the angular frequency 𝜔 by 𝜙(𝜔). In the frequency domain, a phase shift can be modeled as a linear time-invariant (LTI) system that has the frequency response of 𝐺(𝜔) = 𝑒 −𝑗𝜙(𝜔) .
(2)
Since |𝐺(𝜔)| = 1, the amplitude spectrum will stay unmodified as the signal passes such a phase modifying system. 1) Approximating the dispersion system If the Fourier transform of any signal 𝑥 is 𝑋, and it has the inverse Fourier transform
𝑥(𝑡) =
1 ∞ ∫ 𝑋(𝜔)𝑒 𝑗𝜔𝑡 𝑑𝜔, 2𝜋 −∞
(3)
then the output of the system after phase modification is 𝑦(𝑡) =
1 ∞ ∫ 𝑋(𝜔)𝑒 −𝑗𝜙(𝜔) 𝑒 𝑗𝜔𝑡 𝑑𝜔 . 2𝜋 −∞
(4)
Assuming 𝜙 is well-behaved; it can be approximated by a polynomial function. For example, any continuous 𝜙 can be uniformly approximated as closely as desired by polynomials over a closed interval [𝜔1 , 𝜔2 ], which is assured by the Weierstrass Approximation Theorem [18,19]. What is more, if 𝜙 is smooth enough, it can be expanded e.g. as a truncated Taylor series around a point 𝜔0 [18]. More generally, 𝜙 may be expanded e.g. in terms of Chebyshev polynomials that provide near optimal construction for a given degree [19]. This process is similar to Fourier analysis, but with a polynomial basis instead of the usual trigonometric one. In the following, we consider the approximation 𝜙(𝜔) ≈ 𝑎 + 𝑏ω + 𝑐𝜔2 ,
(5)
where 𝑎, 𝑏, 𝑐 ∈ ℝ are constants specific to 𝜙. In other words, we assume that the phase delay function may be adequately approximated by a second degree polynomial. 2) Impulse response of dispersion system Let us consider the impulse response of the dispersion system. Hence, let the input 𝑥(𝑡) = 𝛿(𝑡) for which 𝑋(𝜔) = 1. The response is now 𝑦(𝑡) = ≈
1 ∞ ∫ 1𝑒 −𝑗𝜙(𝜔) 𝑒 𝑗𝜔𝑡 𝑑𝜔 2𝜋 −∞ 2 1 ∞ 𝑒 −𝑗𝑎 ∫−∞ 𝑒 𝑗𝜔(𝑡−𝑏)−𝑗𝑐𝜔 𝑑𝜔 , 2𝜋
(6)
where the polynomial approximation was used. To evaluate the integral, we follow [20] and consider the Fourier integral ∞
1
∫ 𝑒 𝑗𝜔𝑡−(𝛼+𝑗𝛽)𝜔 2𝜋 −∞
2 /2
𝑑𝜔 =
1 √2𝜋(𝛼+𝑗𝛽)
𝑒
𝑡2 2(𝛼+𝑗𝛽)
−
,
(7)
where 𝛼, 𝛽 ∈ ℝ, and 𝛼 ≥ 0 showing that the transform of a real or complex Gaussian is also a Gaussian; real or complex correspondingly. Continuing with the approximation, we do a change of variable 𝑡 ↦ 𝑡 − 𝑏, set 𝛼 = 0, 𝛽 = 2𝑐, and choose from the polar representation 𝛼 + 𝑗𝛽 = 𝑟𝑒 𝑗𝜃 the square root 1
𝜃
𝑟 2𝑒 𝑗2 . Now since 𝛼 = 0, 𝛽 = 2𝑐, we have 𝜋 2 2 𝑟 = √0 + (2𝑐) = 2|𝑐|, and 𝜃 = ± depending on the sign of 𝑐. Thus √𝛼 + 𝑗𝛽 = √2|𝑐|𝑒 𝑦(𝑡) ≈
1 2√𝜋|𝑐|
𝜃
±𝑗 𝜃 4
𝑒 −𝑗(𝑎−𝑠𝑖𝑔𝑛(𝑐)4 ) 𝑒 𝑗
2
, and
(𝑡−𝑏)2 4𝑐
,
(8)
which is a linear chirp signal with a delay of 𝑏, phase shift of 𝜃 1 𝑎 − 𝑠𝑖𝑔𝑛(𝑐) , and amplitude scaling of . 4
2√𝜋|𝑐|
1534-4320 (c) 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TNSRE.2016.2525829, IEEE Transactions on Neural Systems and Rehabilitation Engineering
3 3) Gaussian pulse response of the dispersion system The Fourier integral (7) gives us also the Fourier transform of Gaussian pulses 𝑥(𝑡) =
1 √2𝜋𝜏
𝑒
𝑡2 2𝜏
−
(9)
as 𝜔2
𝑋(𝜔) = 𝑒 −𝜏 2 ,
1
∞
∫ 𝑒 2𝜋 −∞
≈ 𝑒 −𝑗𝑎
1
−𝜏
𝜔2 2
∞
𝑒 −𝑗𝜙(𝜔) 𝑒 𝑗𝜔𝑡 𝑑𝜔
∫ 𝑒 2𝜋 −∞
𝜏 2
𝑗𝜔(𝑡−𝑏)−𝑗(𝑐+ )𝜔2
𝑑𝜔 ,
(11)
which, based on (7-8), we can see to actually be a linear chirp with a modified instantaneous frequency slope. 4) Amplitude modulation Exponentially decreasing amplitude can be expressed as ℎ(𝑡) = 𝑒 −𝜆|𝑡| .
(12)
It has the Fourier transform 𝐻(𝜔) =
2𝜆
(13)
𝜆2 +𝜔2
that decays rather rapidly with increasing frequency (low-pass type characteristics). A chirp signal 𝑓(𝑡) that undergoes decaying amplitude may then be expressed as the product 𝑥(𝑡) = 𝑓(𝑡)ℎ(𝑡).
(14)
Since a multiplication in the time domain is a convolution in the frequency domain, we see that ∞
𝑋(𝜔) = 𝐹(𝜔) ∗ 𝑋(𝜔) = ∫−∞ 𝐹(𝜓)𝑋(𝜔 − 𝜓)𝑑𝜓.
(15)
This causes a blurring of spectral peaks. We may also combine the causality of the chirp (𝑓(𝑡) = 0 when 𝑡 ≤ 0) by using a unit step function 0 𝑢(𝑡) = { 1
,𝑡 < 0 . ,𝑡 ≥ 0
.
(18)
𝑦(𝑡) =
∞ 1 ∞ ∫ ∫ 𝐻(𝜓)𝑒 −𝑗𝜙(𝜔−𝜓) 𝑋(𝜔 2𝜋 −∞ −∞
− 𝜓)𝑒 𝑗𝜔𝑡 𝑑𝜓 𝑑𝜔 ,(19)
where 𝐻(𝜔) is the spectral smearing function, for which ∡𝐻(𝜔) is constant and 𝐻(𝜔) decreases with increasing |𝜔|, 𝜙(𝜔) is the phase dispersion function (|𝑒 −𝑗𝜙(𝜔) | = 1, ∡𝑒 −𝑗𝜙(𝜔) = −𝜙(𝜔)), and X(𝜔) is the Fourier transform of the input signal. For a single frequency signal we can now analyze the response. A sinusoid 𝑓(𝑡) = sin 𝜔0 𝑡 has the Fourier transform 𝜋 𝐹(𝜔) = [𝛿(𝜔 − 𝜔0 ) − 𝛿(𝜔 + 𝜔0 )]. So the convolution with 𝑗
the spectral smearing basically gives just the spectral smearing function which is then phase distorted (or if dispersion is first, the sinusoids phase is just altered before the smearing). An important finding is that, considering narrower band input signals, we can interpret the spectral smearing here as introduction of new spectral components. I.e. the original input frequency will be accompanied by others near it with dispersed phases. Furthermore, the approximations for dispersion and amplitude modulation reveals that a chirp model accounting for both linear and exponential amplitude as well as frequency components can be used to accommodate both effects in a system model. Chirp models with the desired properties are explored next. C. Chirp signals To accommodate the system model we next consider some suitable linear and nonlinear chirp signals and envelope decay functions. 1) Linear chirp A so-called linear chirp is described with the quadratic phase function that can be produced by integration: 𝜙𝑙𝑖𝑛 (𝑡) = 𝜙0 + 2𝜋 ∫[𝑓0 + ∫ 𝑘 𝑑𝑡]𝑑𝑡 𝑘
= 𝜙0 + 2𝜋(𝑓0 𝑡 + 𝑡 2 ), 2
(20)
where 𝜙0 is the initial phase, 𝑓0 is the initial frequency, (16)
So the output we are considering is actually of the form 𝑢(𝑡)𝑓(𝑡). Hence, we may combine the unit step with the amplitude decay part as the multiplicative factor ℎ(𝑡) = 𝑢(𝑡)𝑒 −𝜆𝑡 ,
1 𝜆+𝑗𝜔
B. System model (based on the aforementioned) To summarize the findings above, if the system consists of phase distortion and spectral smearing (convolution with a low-pass like in the frequency domain), it can exhibit decaying chirp signal responses e.g. to impulses and Gaussian pulses. The system output can be modeled as
(10)
which can easily be seen by setting 𝛼 = 𝜏 and 𝛽 = 0. So the output for a Gaussian pulse input would be 𝑦(𝑡) =
𝐻(𝜔) =
(17)
where 𝑢(𝑡) is the unit step function. The Fourier transform of ℎ is now also a low-pass type function
𝑓 −𝑓
𝑘 = 1 0 defines the chirp rate, and 𝑇 is the chirp duration. 𝑇 The complex representation of a linear chirp is thus 𝑘 2 )]
𝑓𝑙𝑖𝑛 (𝑡) = 𝑒 𝑗[𝜙0 +2𝜋(𝑓0𝑡+2𝑡
.
(21)
In addition to this fundamental frequency, the frequency modulation gives also rise to additional frequency components. The utilization of complex valued signal is useful for computations, but it should be noted that the actual observed chirp is the real part of 𝑓, i.e.
1534-4320 (c) 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TNSRE.2016.2525829, IEEE Transactions on Neural Systems and Rehabilitation Engineering
4
𝑘
ℛℯ 𝑓𝑙𝑖𝑛 (𝑡) = cos [𝜙0 + 2𝜋 (𝑓0 𝑡 + 𝑡 2 )].
(22)
2
2) Polynomial nonlinear chirp A polynomial chirp is given by expansion of the linear chirp 𝑐𝑛 𝑛 phase function integral of the form ∑𝑁 𝑡 where 𝑐𝑛 are the 𝑛=0 𝑛!
polynomial coefficients. The first coefficient 𝑐0 corresponds to the initial phase 𝜙0 , 𝑐1 is the initial frequency 𝑓0 , and 𝑐2 = 𝑘 is the linear chirp rate. The linear chirp is thus the special case when 𝑁 = 3. Now with the chirp length 𝑇 and the linear part separated we get the complex representation 𝑓𝑝𝑜𝑙𝑦 (𝑡) = 𝑒
𝑗[𝜙0 +2𝜋(𝑓0 𝑡+
𝑐𝑛 𝑘 2 𝑛 𝑡 +∑𝑁 𝑛=3𝑛!∙𝑇𝑛−1 𝑡 )] 2𝑇
,
ℛℯ 𝑓𝑝𝑜𝑙𝑦 (𝑡) = cos [𝜙0 + 2𝜋 (𝑓0 𝑡 +
∑𝑁 𝑛=3
+ 𝑐𝑛
𝑛!∙𝑇 𝑛−1
𝑡 𝑛 )].
(24)
3) Exponential nonlinear chirp An exponential chirp of the form 𝑘 𝑡 can be expressed in an exponential form by noting that 𝑘 = 𝑒 ln 𝑘 so that 𝑘 𝑡 = 𝑡 (𝑒 ln 𝑘 ) = 𝑒 𝑡 ln 𝑘 . The exponential form can be written as a Maclaurian series 𝑒 𝑥 = ∑∞ 𝑛=0
𝑥𝑛 𝑛!
corresponding to a special (ln 𝑘)𝑛 𝑡 𝑛
case of infinite polynomial series 𝑒 𝑡 ln 𝑘 = ∑∞ where 𝑛=0 𝑛! the polynomial coefficients have a uniform structure 𝑐𝑛 = (ln 𝑘)𝑛 . Using an exponential instantaneous frequency function of 𝑓 the form 𝑓0 𝑘𝑡/𝑇 , where 𝑘 = 1 , we get, by integration and 𝑓0
substitution 𝑢 = 𝜏/𝑇 𝑡
𝜏
𝑡
∫0 𝑓0 𝑘 𝑇 𝑑𝜏 = 𝑓0 𝑇 ∫0 𝑘 𝑢 𝑑𝑢 = 𝑓0 𝑇
𝑘 𝑡/𝑇 −1 ln 𝑘
.
(25)
Now, the complex representation of the exponential chirp becomes 𝑓𝑒𝑥𝑝 (𝑡) = 𝑒
𝑘𝑡/𝑇 −1 )] ln 𝑘
𝑗[𝜙0 +2𝜋𝑓0 𝑇(
(26)
and the corresponding time-domain signal is then ℛℯ 𝑓𝑒𝑥𝑝 (𝑡) = cos [𝜙0 + 2𝜋𝑓0 𝑇 (
𝑘 𝑡/𝑇 −1 ln 𝑘
)].
(𝑎−𝑏) 𝑇
𝑡.
(28)
An exponential decay defined as 𝐴𝑒𝑥𝑝 (𝑡) = 𝐶𝑒 −𝜆𝑡 is a monotonically decreasing function, when 𝐶 is a positive constant and 𝜆 > 0 is the decay rate. Because 𝑒 0 = 1 and we can subtract 𝑒 −𝜆𝑇 to set the function to 0 at time 𝑇, the function can be expressed on the amplitude interval [𝑎, 𝑏] and (𝑎−𝑏) time interval [0, 𝑇] with a scale term 𝐶 = −𝜆∙0 −𝜆∙𝑇 = 𝑒
(𝑎−𝑏) 1−𝑒
−𝑒
−𝜆𝑇 . Now, with the scale term 𝐶, the end amplitude 𝑏, and
subtracting 𝑒 −𝜆𝑇 , the exponential dampening can be expressed as 𝐴𝑒𝑥𝑝 (𝑡) = 𝐶(𝑒 −𝜆𝑡 − 𝑒 −𝜆𝑇 ) + 𝑏 to get an exponential decay
(23)
where 𝑐𝑛 [𝑛 ≥ 3] are the coefficients for the higher nonlinear polynomials. The time-domain signal is then
𝑘 2 𝑡 2
𝐴𝑙𝑖𝑛 (𝑡) = 𝑎 −
(27)
4) Linear and exponential decay The observed chirps show a decay of signal. A linear decay 𝐴𝑙𝑖𝑛 (𝑡) = 𝐶 − 𝑟𝑡, where 𝐶 is a constant and 𝑟 is the rate of decay, can be defined in an amplitude interval [𝑎, 𝑏] at time interval [0, 𝑇] as
𝐴𝑒𝑥𝑝 (𝑡) = (𝑎 − 𝑏) [
𝑒 −𝜆𝑡 −𝑒 −𝜆𝑇 1−𝑒 −𝜆𝑇
] + 𝑏.
(29)
III. DATA A. Data collection For the SEP model creation and testing, experimental data were collected using a rat model with two varying factors, anesthetic level and SCI, both known to affect the signal [21,22] (Fig. 1). The experimental protocol was approved by the Institutional Animal Care and Use Committee (IACUC) at the National University of Singapore. The data were collected from 21 adult female SpragueDawley rats (200-220 g). Four screw electrodes were implanted to the cranium of rats above the somatosensory cortex of different limbs (Fig. 2). A reference electrode was positioned in the right parietal area close to lambda. Carboxylate dental cement (3M) was used in fixing the electrodes to the cranium. The animals were anesthetized with ketamine (75mg/kg) and xylazine (10 mg/kg) cocktail during the electrode placement. Baseline SEP recoding was carried out at least seven days after the electrode implantation. During the recording, rats were anesthetized with a mixture of isoflurane and 100% oxygen delivered at a rate of 1.3 L/min. A rodent-size anesthesia mask connected to a diaphragm with a C-Pram circuit designed to deliver and evacuate the gas was used for maintaining the anesthesia. To assess the SEPs at different anesthetic levels, the recording was carried out as the isoflurane dosage was increased in a step-like manner from 1.5% to 2.5% using 0.5% increments. The dosage was kept fixed for 5 min at each step before recording the SEPs to guarantee the equilibrium. A neuro-electrophysiology monitoring setup (Tucker-Davis Technology) was used for multilimb acquisition of SEP [23]. During the recording, four pairs of stainless-steel stimulating needle electrodes (Safelead F-E3-48, Grass Technologies) were placed in proximity to the tibial nerve in both hindlimbs and to the median nerve in both forelimbs. The needle electrodes were then connected to a stimulus generator (Digitimer Ltd.) and the cranial screw electrodes to an amplifier. The stimulus generator was used to generate 3.5 mA stimulating pulses with a pulse width of 200 µs. Pulses were generated at a rate of 1 Hz sequentially to each of the four limbs approximately 150 times per limb. During the
1534-4320 (c) 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TNSRE.2016.2525829, IEEE Transactions on Neural Systems and Rehabilitation Engineering
5 Furthermore, when assessing the effects of SCI, the signal analysis was restricted to the data of severe injury group.
Fig. 1. The effect of anesthesia and spinal cord injury on SEP. The figure shows the SEP of a single rat before (left column) and after (right column) severe spinal cord injury at three different anesthetic levels (different rows anesthetic effect increasing downwards).
B. Preprocessing In the preprocessing procedure applied, individual single sweep signals are averaged to provide a mean SEP signal. A Band Pass (BP) filter is applied to remove low frequency drifting of the mean signals and noise above and below the frequency band of interest. The corner frequencies of 𝑓𝑙 = 10 Hz and 𝑓ℎ = 280Hz were chosen to remove much of the noises, EEG channels from uncorrelated brain activities, and also all signal components not fitting in the analysis window (100 ms). A zero-phase IIR filtering method [26] with reversed and flipped signal adjustment conserving the DC value is used to preserve the absolute amplitude and phase information of the SEP pulses. Maximum absolute value is searched within a preset minimum and maximum delay envelope to identify most likely first SEP peak direction in order to limit the optimization search space and to increase the robustness of the model fitting. IV. METHODS
Fig. 2. Electrode locations for somatosensory evoked potential (SEP) measurement. The forelimb recording sites (circles) are located 0.2 mm posterior and 3.8 mm lateral to the bregma. The hindlimb recording sites (stars) are located 2.5 mm posterior and 2.8 mm lateral to the bregma. Reference electrode (square) is located 3 mm lateral to lambda.
stimulation, SEPs were recorded simultaneously from all four cortical screw electrodes via the amplifier and data acquisition setup. In the analysis, only the recordings from the corresponding electrode of the stimulated limb were used. Furthermore, for the current study, the signal analysis was restricted to the hindlimb data. After the baseline SEP recording, the animals were divided into three SCI groups (n = 5, 6, and 5 for mild, moderate, and severe injury, respectively) and one control group (n = 5). For the rats in the SCI groups, a controlled spinal cord contusion was performed using the NYU-MASCIS (New York University - Multicenter Animal Spinal Cord Injury Study) weight-drop impactor [24]. In this procedure, a 10 g rod with a flat circular impact surface was dropped from a certain height (6.25, 12.5, and 25 mm for mild, moderate, and severe injury groups, respectively) on the midline of the spinal cord exposed by removing the lamina at T8. A more detailed description of the controlled spinal cord contusion procedure can be found in e.g. [25]. For the rats in the control group, only the laminectomy without contusion was performed. After the operation, SEPs were recorded following the above described protocol weekly for 4 weeks. In this paper, only the data from baseline recordings and 28 days after the operation are used.
An automatic parameterization of SEPs is proposed using a chirp model where the resulting SEP waveforms are assumed as impulses undergoing group velocity dispersion while propagating along a multipath neural connection. The dispersion is assumed to cause changes to the signal phase spectra that manifests in the observations as chirp-like result waveforms. The chirp waveforms are assumed as non-periodic narrowband signals with strictly monotonic decay in time and frequency. The method for the automatic parametrization of SEPs with chirp modeling was implemented in Matlab™ technical computing language and is available online1. A. Proposed Chirp Models Chirp waveforms are modeled as monotonically strictly decreasing sinusoids using linear (22), 3rd or 4th degree polynomial (24) or exponential (27) instantaneous phase functions. A phase model with both linear and exponential component was also considered. The amplitude dampening weight function 𝐴(𝑡) based on the linear and/or exponential dampening (29), a unit step function (16), and added with an onset delay of 𝜏 is 𝐴(𝑡) = {
0 (𝑎 − 𝑏) [
,𝑡 < 𝜏 𝑒 −𝜆𝑤 (𝑡−𝜏) −𝑒 −𝜆𝑤 𝑇 1−𝑒 −𝜆𝑤 𝑇
]+𝑏
,𝑡 ≥ 𝜏
,
(30)
where 𝑎 = initial weight, 𝑏 = end weight, 𝑡 = time, 𝜏 = onset delay, 𝜆𝑤 = amplitude nonlinearity coefficient (exponential weight), 𝑇 = chirp duration. An illustrative fitted example dampening function (30) with parameters indicated is presented in Fig. 3. A chirp phase model using a linear (22) and an exponential chirp (27) component weighted by 𝜆𝑓 , and with an onset delay 1
Available: http://www.ee.oulu.fi/~jukortel/SEP.html
1534-4320 (c) 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TNSRE.2016.2525829, IEEE Transactions on Neural Systems and Rehabilitation Engineering
6 B. Model fitting and optimization Parameters (𝜏, 𝑎, 𝑏, 𝜆𝑤 , 𝑓0 , 𝑓1 , and 𝜆𝑓 ) are fitted for the preprocessed mean SEP signals using a PSO method [27]. The PSO Toolbox [28] was used to implement a Trelea type 1 PSO. Random initial seed positions for 150 particles with a 10000 epoch maximum iteration time was used for robust finding of the model parameters. The initial peak direction was fixed to positive side to narrow the search space into half. Search ranges for the parameters are: 𝜏 = 0.005 – 0.02 s, 𝑎 = ½·max|𝑥| – 1½·max|𝑥|, 𝑏 = 0 – 0.3·max|𝑥|, 𝜆𝑤 = 0.0001 – 5, 𝑓0 = 30 – 120 Hz, 𝑓1 = 0.01 – 15 Hz, 𝜆𝑓 = 0 – 1. Chirp duration 𝑇 = 80 ms was used. The goodness of fit (objective function to be minimized) used is the Mean Squared Error (MSE) that is linearly weighted (weight decreases towards the ends of the recordings) to emphasize fitting of the beginning parts of the chirps. The linear weight 𝑤𝑙𝑖𝑛 is defined as Fig. 3. Chirp amplitude model. The dashed line indicates a real fitted amplitude function A example.
𝑤𝑎
, 𝑡(𝑖) < 𝜏
𝑤𝑙𝑖𝑛 (𝑖) = {𝑤𝑎 − (𝑤𝑎 − 𝑤𝑏 ) 𝑤𝑏
𝑡(𝑖)−𝜏 𝑇
, 𝜏 ≤ 𝑡(𝑖) < 𝜏 + 𝑇,(33) , 𝑡(𝑖) > 𝜏 + 𝑇
where 𝑖 = sample index, 𝑤𝑎 = linear weight start, 𝑤𝑏 = linear weight stop. Linear weights 𝑤𝑎 = 1 and 𝑤𝑏 = 0.1 was used. The linearly weighted 𝑀𝑆𝐸𝑤 𝑙𝑖𝑛 is then defined as: 2
1
𝑀𝑆𝐸𝑤𝑙𝑖𝑛 = ∑𝑁 𝑖=1 𝑤𝑙𝑖𝑛 (𝑖) ∙ (𝑥(𝑖) − 𝑦(𝑖)) , 𝑁
(34)
where 𝑖 = sample index, 𝑁 = chirp/sample length, 𝑤𝑙𝑖𝑛 = linear weight, 𝑥 = original signal, 𝑦 = chirp signal. The average Root Mean Squared Error 1
Avg. 𝑅𝑀𝑆𝐸 = ∑ √𝑀𝑆𝐸 𝑁 Fig. 4. Instantaneous frequency of the linear and/or exponential chirp phase model. The dashed line indicates a real fitted example.
𝜏 and a duration of 𝑇 is defined as:
(35)
of the whole data and the average linearly weighted root mean squared error Avg. 𝑅𝑀𝑆𝐸𝑤𝑙𝑖𝑛 =
1 𝑁
∑ √𝑀𝑆𝐸𝑤𝑙𝑖𝑛
(36)
(𝑡−𝜏)
𝜙𝑒𝑥𝑝𝑙𝑖𝑛 (𝑡) = 𝜆𝑓 (2𝜋𝑓0 𝑇 (
𝑘 𝑇 −1 ln 𝑘
)) + (1 − 𝜆𝑓 )
∙ (2𝜋𝑓0 (𝑡 − 𝜏) + 2𝜋
(𝑓1 −𝑓0 ) 2𝑇
(𝑡 − 𝜏)2 ),
(31)
where 𝑓0 = initial frequency, 𝑓1 = end frequency, 𝜆𝑓 = frequency nonlinearity coefficient,
and
𝑓
𝑘 = ( 1 ). 𝑓0
An
illustrative fitted instantaneous frequency example of the linear and exponential phase function (31) with associated parameters indicated is presented in Fig. 4. The final chirp waveform 𝑦(𝑡) using (30) and (31) is then defined as: 𝑦(𝑡) = 𝐴(𝑡) ∙ cos (𝜙0 + 𝜙𝑒𝑥𝑝𝑙𝑖𝑛 (𝑡)) .
(32)
are then used to evaluate which chirp model is the best model overall using all available 252 averaged SEP samples (i.e. 21 rats with baseline and 28 days post operation samples for both hind limbs in three anesthetic levels). An illustrative example of a fitted chirp using the linear and/or exponential phase chirp model (32) that combines the previous examples, dampening function (Fig. 3) and instantaneous frequency function (Fig. 4), into the final chirp model, associated parameters, and derived SEP features is presented in Fig. 5. C. SNR Testing SNR performance testing of the chirp model fitting was performed using all of the 252 averaged chirp recordings as clean signal templates. A bootstrap approach was used to
1534-4320 (c) 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TNSRE.2016.2525829, IEEE Transactions on Neural Systems and Rehabilitation Engineering
7
Fig. 5. An illustrative example of a fitted chirp model with the SEP features showed using the linear and/or exponential phase chirp model (31). An averaged SEP measurement x (solid gray line) is shown with the fitted chirp model y (solid black line) that also corresponds to the fitted example amplitude function A (Fig. 3) and the fitted example instantaneous frequency function in Fig. 4.
evaluate the stability of the averaged templates without added noise to evaluate reasonable lower limits for model fitting errors corresponding to the 150 averaged SEP signal fittings. 1,000 bootstrap means of 150 random overlapping selections of individual single sweep signals were taken from each 252 samples (i.e. a total of 252,000 bootstraps). Mean variances along all of the samples were calculated from the corresponding 1000 bootstrap means. The mean variances correspond to an unbiased estimators MSE and can then be used to produce estimates for the MSE and linearly weighed 𝑀𝑆𝐸𝑤𝑙𝑖𝑛 expected with the measurement process corresponding to the generation of the templates (i.e. an estimate of intrinsic error of the templates due to noises or other uncorrelated processes regardless of a SEP model). For the SNR testing, an array of Additive Gaussian White Noise (AWGN) measured at -12, -6, 0, 6, 12, 18, and 24 dB SNR levels (i.e. in relation to each template signal power to stabilize the SNR at each level) was then applied to produce test signals. The selected SNR dB levels correspond to SNR increases attained when using multiple trial averaging of 4 n signal sweeps between n dB levels (e.g. 4 1 = 4 averages from 0dB to 6dB and 43 = 64 averages from -6dB to 12dB), in the presence of AWGN. A Monte Carlo method using 10 independent AWGNs for each sample was implemented to produce stable SNR performance readings at each dB level. In the testing a total of 2520 independent noise samples were fitted with a chirp model for each selected SNR level using both the 𝑀𝑆𝐸 and the 𝑀𝑆𝐸𝑤𝑙𝑖𝑛 as measures for the goodness of fit when comparing to the corresponding clean template. The 𝑅𝑀𝑆𝐸 and the 𝑅𝑀𝑆𝐸𝑤𝑙𝑖𝑛 of each of the 252 samples were calculated for each SNR level tested. To evaluate how good the chirp fitting is for individual signals three error rate criterions using normalized errors were used (i.e. in order to measure erroneous waveform fitting rather than different amplitudes). First, For each 252 sample, the bootstrap noise estimate was normalized with the clean template signal estimate (using both 𝑅𝑀𝑆𝐸 and the
𝑅𝑀𝑆𝐸𝑤𝑙𝑖𝑛 ) to produce intrinsic normalized bootstrap error estimates (both normalized 𝑅𝑀𝑆𝐸 and normalized 𝑅𝑀𝑆𝐸𝑤𝑙𝑖𝑛 ) for the samples. Next, for each SNR level, chirp fitting errors were similarly normalized using the corresponding clean template signal estimates to produce normalized chirp model fitting errors. Then, the differences of normalized chirp model fitting errors at each SNR level and the normalized bootstrap error estimates were calculated to produce fitting Δ error distributions (both Δ normalized 𝑅𝑀𝑆𝐸 and Δ normalized 𝑅𝑀𝑆𝐸𝑤𝑙𝑖𝑛 ) at each SNR level. The three thresholds for fitting errors were chosen as Δ errors larger than the median (50%), 3rd quartile (75%), or 9th decile (90%) values of the normalized bootstrap error distribution (both 𝑅𝑀𝑆𝐸 and 𝑅𝑀𝑆𝐸𝑤𝑙𝑖𝑛 ). The choice of thresholds allows for a strict, average, and lax evaluation of bounds, correspondingly, for a chirp fitting that takes into account the estimated intrinsic noise in each template (with no weighting by using the 𝑅𝑀𝑆𝐸 or emphasizing the beginning by using the 𝑅𝑀𝑆𝐸𝑤𝑙𝑖𝑛 ). D. Feature Calculation The resulting chirp model parameters can be used as features describing the SEP recordings. Specifically, the onset delay (𝜏) and amplitude parameters (𝑎, 𝑏 and 𝜆𝑤 ) but also the frequency parameters (𝑓0 , 𝑓1 and 𝜆𝑓 ) are of interest. The goodness-of-fit measure gained from the optimization, i.e. the MSE measure, can also be used to calculate a normalized sum of square residual (𝑁𝑆𝑆𝑅𝑒𝑠 ) feature relating to the ratio of the model explained signal and the residual unexplained signal 𝑁𝑆𝑆𝑅𝑒𝑠 =
2
∑𝑁 𝑖=1(𝑥(𝑖)−𝑦(𝑖)) 2 ∑𝑁 𝑖=1 𝑦(𝑖)
.
(37)
The fitted model can also be used to find the initial peaks and peak latencies of the chirp-like signal (see Fig. 5), either from the model itself, which is a trivial task of finding the zeros, minimums, and maximums due to the model properties (e.g. the chirps have no riding waves) or, using the model peak locations as an initial estimate, from the original SEP signal (i.e. by local search and fitting). Parameters of the fitted chirp signals and the residual goodness-of-fit measure are presented in Table I. SEP features derived from the fitted chirp models are presented in Table II. E. Anesthesia level and SCI analysis Two testing groups were used to evaluate the parameters of the fitted chirp models and also features derived from the fitted chirp models. In the first group of tests, i.e. anesthesia level, all the 126 averaged baseline recordings (i.e. 21 rats with two hind limb recordings in three anesthetic levels) were used to evaluate the effect of changing anesthetic concentration on the SEP parameters and features. Statistical testing was performed to determine if the different anesthetic levels have a significant impact on the measured SEP recordings between the anesthetic concentration levels used: 1.5% vs. 2.0%, 2.0% vs. 2.5%, 1.5% vs. 2.5%, and 1.5% vs. 2.0% vs. 2.5% (i.e. all levels).
1534-4320 (c) 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TNSRE.2016.2525829, IEEE Transactions on Neural Systems and Rehabilitation Engineering
8 TABLE I CHIRP MODEL PARAMETERS Parameter
Quantity
𝑇
Chirp length
𝜏
Chirp onset delay
𝑎
Start amplitude
𝑏
End amplitude
𝜆𝑤
Amplitude nonlinearity coefficient
𝑓0
Start frequency
𝑓1
End frequency
𝜆𝑓
Frequency nonlinearity coefficient
𝑁𝑆𝑆𝑅𝑒𝑠
Normalized sum of squares residual
evaluate the generated p-values in order to counteract the problem of multiple comparisons. V. RESULTS
TABLE II SEP FEATURES DERIVED FROM THE CHIRP MODEL Feature
Quantity
st
1 peak delay
Delay of the first peak
1st peak
First peak amplitude (absolute)
2nd peak delay
Delay of the second peak
2nd peak
Second peak amplitude (absolute)
p2p delay
1st peak to 2nd peak delay
p2p
1st peak to 2nd peak amplitude
TABLE III AVERAGE MSES OF THE TESTED CHIRP MODELS 𝐴𝑣𝑔. 𝑅𝑀𝑆𝐸𝑤𝑙𝑖𝑛 (µV) Chirp model 𝐴𝑣𝑔. 𝑅𝑀𝑆𝐸 (µV) Linear
4.88
3.76
3rd degree polynomial
4.45
3.37
4 degree polynomial
3.79
2.78
Exponential
3.86
2.72
Linear & exponential
3.66
2.63
Bootstrap noise estimate
3.00
2.12
th
First, the fitting of all available SEP data (252 averaged signals) was performed using linear, 3rd and 4th degree polynomial, exponential, and both linear and exponential phase chirp functions with linear and/or exponential amplitude dampening functions. The best performing model, i.e. the model resulting in the overall lowest 𝐴𝑣𝑔. 𝑅𝑀𝑆𝐸 and linearly weighted 𝐴𝑣𝑔. 𝑅𝑀𝑆𝐸𝑤𝑙𝑖𝑛 , presented in Table III, was selected as the proposed chirp model and used in the further tests. The bootstrap analysis of the averaged SEP signal variances resulted in an 𝐴𝑣𝑔. 𝑅𝑀𝑆𝐸 = 3.00 µV corresponding to an average SNR of 15.69 dB for the 150 sample average templates with a bootstrap mean estimate for signal 𝐴𝑣𝑔. 𝑅𝑀𝑆𝐸 = 18.25 µV. The template average SNR also corresponds to an average SNR of -6.19 dB for a single trial. When considering the template SNR, the added AWGN SNR array of -12, -6, 0, 6, 12, 18, and 24 dB, used in the Monte Carlo runs, therefore corresponds to an array of -12.01, -6.06, -0.22, 5.27, 9.94, 13.17, and 14.82 dB estimated total average SNRs for the 150 averaged samples, respectively. A corresponding single trial SNR can be calculated by subtracting 21.76 dB from each estimate SNR (i.e. the effect of 150 average signals assuming AWGS). An 𝑅𝑀𝑆𝐸 correlation analysis between the bootstrap variances and the template fitting errors with the linear and/or exponential chirp model for each 252 samples revealed a correlation of 0.73 (p
1
Automatic parametrization of somatosensory evoked potentials with chirp modeling Eero Väyrynen, Kai Noponen, Ashwati Vipin, Member, IEEE, Thow Xin Yuan, Hasan Al-Nashash, Senior Member, IEEE, Jukka Kortelainen, Member, IEEE, and Angelo All Abstract—In this paper, an approach using polynomial phase chirp signals to model somatosensory evoked potentials (SEPs) is proposed. SEP waveforms are assumed as impulses undergoing group velocity dispersion while propagating along a multipath neural connection. Mathematical analysis of pulse dispersion resulting in chirp signals is performed. An automatic parameterization of SEPs is proposed using chirp models. A Particle Swarm Optimization algorithm is used to optimize the model parameters. Features describing the latencies and amplitudes of SEPs are automatically derived. Rat model is then used to evaluate the automatic parameterization of SEPs in two experimental cases, i.e. anesthesia level and spinal cord injury (SCI). Experimental results show that chirp based model parameters and the derived SEP features are significant in describing both anesthesia level and SCI changes. The proposed automatic optimization based approach for extracting chirp parameters offers potential for detailed SEP analysis in future studies. The method implementation in Matlab technical computing language is provided online. Index Terms—Anesthesia, biological system modeling, parameter estimation, particle swarm optimization, spinal cord injury
I. INTRODUCTION
S
OMATOSENSORY evoked potentials (SEPs) represent a widely used neurophysiological method when examining the functional integrity of ascending sensory pathways. Today, they are clinically used for various diagnostic purposes such as the assessment of brain injury after cardiac arrest and intraoperative monitoring during high risk surgeries [1-3]. For studies involving laboratory animals, SEPs provide a convenient minimally-invasive method to quantitatively assess the functionality of various portions of the somatosensory Manuscript received July 21, 2015. This work was supported in part by grants 2013-MSCRFII-0109-00 from the Maryland Stem Cell Research Fund as well as R-175-000-121-133, R-175-000-121-733, R-175-000-122-112, and R-711-201-026-133 from the National University of Singapore. E. Väyrynen, K. Noponen, and J. Kortelainen are with the Physiological Signal Analysis Team, Center for Machine Vision and Signal Analysis, University of Oulu, Finland. A. Vipin and X. Thow are with the Singapore Institute for Neurotechnology (SINAPSE), National University of Singapore, Singapore. H. Al-Nashash is with the Department of Electrical Engineering, American University of Sharjah, United Arab Emirates. A. All is with the Department of Orthopedic Surgery and Biomedical Engineering, National University of Singapore, Singapore and the Department of Neurology and Biomedical Engineering, Johns Hopkins University, USA. Address correspondence to Jukka Kortelainen (e-mail: [email protected]).
pathways. However, the analysis of SEPs still relies mainly on the determination of simple parameters such as amplitudes and latencies from the averaged signal waveform [4] even though several more sophisticated measures for the signal quantification have been proposed [5,6]. Amplitude and/or frequency modulated non-stationary signals are common in many engineering problems (e.g. radar, communications, and biomedical imaging) as well as in natural signals (e.g. echolocation with natural sonar). The parameter estimation of chirps and in general polynomialphase signals has therefore been widely studied [7,8]. The estimation is generally based on models aiming for continuous chirp signal analysis using methods such as Maximum Likelihood (ML) estimators [6], polynomial phase transform, [9], variations of the higher order ambiguity functions [10], polynomial Wigner-Ville distributions [11], and cubic phase functions [12]. Furthermore, hybrid methods [13] to increase performance and phase unwrapping approaches [14] have recently been studied. These methods yield theoretically very effective parameter estimates for continuous chirp signals with Signal-to-Noise Ratio (SNR) over 0 dB, i.e. approach the Cramer-Rao lower bound for Mean Squared Error (MSE). For SEP measurements, however, the SNR is usually very low, e.g. -15 – 0 dB. In order to counteract the poor SNR, multiple trial average measurements are typically needed to analyze the signals. A growing interest has been for single trial measures of SEPs [15] to enable faster evaluation, e.g. important in spinal surgery, and also access the variability of parameters lost in the averaging processes with the expense of SNR. To work in the low SNR method using spatial measures such as slope analysis and morphological features [4,16], spectral domain features such as spectral coherence [5], wavelet and/or independent component analysis [15], and matching pursuit time-frequency analysis [17] have been proposed. In this paper, an analysis approach that incorporates a chirp based modeling in the normal SEP analysis structure is considered. In order to achieve robust operation, the aim is not to perform perfect reconstructions of the SEP signals with chirps but rather to find and extract best fitting polynomialphase waveform components, and associated parameters, of the SEPs that are consistent with a pulse propagation and dispersion model (i.e. imposes specific time-frequency transition characteristics). First, we discuss the mathematical implications of pulse propagation in dispersive channels
1534-4320 (c) 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TNSRE.2016.2525829, IEEE Transactions on Neural Systems and Rehabilitation Engineering
2 (Section II). We assume that the aggregate cortical signal in the sensor, caused by impulse-like stimuli to the peripheral nerve, can be explained by pulses propagating in a multipath dispersive channel with amplitude dampening resulting in chirp signals. Second, a chirp model for SEP is constructed and a Particle Swarm Optimization (PSO) model fitting technique is used in order to provide automatic measures for the chirp model parameters and also to estimate a set of more traditional amplitude/latency parameters (Section IV; A, B, and C). Third, the resulting parametrization is then evaluated with experimental rat model data in which SEPs are affected by two different factors, anesthetic level and spinal cord injury (SCI) (Section IV; D). Finally, statistical testing is used to identify significant parameters and features in both of the experimental cases (Section IV; E). II. PULSE DISPERSION AND CHIRP SIGNALS Empirically, the general functional form of the SEPs resembles decaying sine waves with varying frequency. Generally, they all seem to be rather closely modelled by dampened chirp signals which are a certain class of frequency modulated sinusoids with amplitude envelope decay 𝑥(𝑡) = 𝐴(𝑡)𝑒 𝑗φ(𝑡) ,
(1)
where 𝐴(𝑡) is the envelope function and 𝜑(𝑡) is the phase function of the model. While the individual neural transmissions are not chirp signals (rather, the neural transmissions more closely resemble to short pulses), the aggregated electric fields in the EEG sensors can be explained by a multipath propagation of pulses undergoing dispersion and amplitude dampening. A mathematical evaluation is provided next to show exactly the connection of dispersed pulses and chirp signals. A. Pulse Dispersion Pulse propagation in dispersive media is a process in which different spectral components of the pulse travel at different velocities. This may happen, for example, due to multipath propagation through channels having different passbands. In the frequency domain, this dispersion is seen as phase shifts of the different spectral components of the pulse. Let us denote the dispersive phase shift at the angular frequency 𝜔 by 𝜙(𝜔). In the frequency domain, a phase shift can be modeled as a linear time-invariant (LTI) system that has the frequency response of 𝐺(𝜔) = 𝑒 −𝑗𝜙(𝜔) .
(2)
Since |𝐺(𝜔)| = 1, the amplitude spectrum will stay unmodified as the signal passes such a phase modifying system. 1) Approximating the dispersion system If the Fourier transform of any signal 𝑥 is 𝑋, and it has the inverse Fourier transform
𝑥(𝑡) =
1 ∞ ∫ 𝑋(𝜔)𝑒 𝑗𝜔𝑡 𝑑𝜔, 2𝜋 −∞
(3)
then the output of the system after phase modification is 𝑦(𝑡) =
1 ∞ ∫ 𝑋(𝜔)𝑒 −𝑗𝜙(𝜔) 𝑒 𝑗𝜔𝑡 𝑑𝜔 . 2𝜋 −∞
(4)
Assuming 𝜙 is well-behaved; it can be approximated by a polynomial function. For example, any continuous 𝜙 can be uniformly approximated as closely as desired by polynomials over a closed interval [𝜔1 , 𝜔2 ], which is assured by the Weierstrass Approximation Theorem [18,19]. What is more, if 𝜙 is smooth enough, it can be expanded e.g. as a truncated Taylor series around a point 𝜔0 [18]. More generally, 𝜙 may be expanded e.g. in terms of Chebyshev polynomials that provide near optimal construction for a given degree [19]. This process is similar to Fourier analysis, but with a polynomial basis instead of the usual trigonometric one. In the following, we consider the approximation 𝜙(𝜔) ≈ 𝑎 + 𝑏ω + 𝑐𝜔2 ,
(5)
where 𝑎, 𝑏, 𝑐 ∈ ℝ are constants specific to 𝜙. In other words, we assume that the phase delay function may be adequately approximated by a second degree polynomial. 2) Impulse response of dispersion system Let us consider the impulse response of the dispersion system. Hence, let the input 𝑥(𝑡) = 𝛿(𝑡) for which 𝑋(𝜔) = 1. The response is now 𝑦(𝑡) = ≈
1 ∞ ∫ 1𝑒 −𝑗𝜙(𝜔) 𝑒 𝑗𝜔𝑡 𝑑𝜔 2𝜋 −∞ 2 1 ∞ 𝑒 −𝑗𝑎 ∫−∞ 𝑒 𝑗𝜔(𝑡−𝑏)−𝑗𝑐𝜔 𝑑𝜔 , 2𝜋
(6)
where the polynomial approximation was used. To evaluate the integral, we follow [20] and consider the Fourier integral ∞
1
∫ 𝑒 𝑗𝜔𝑡−(𝛼+𝑗𝛽)𝜔 2𝜋 −∞
2 /2
𝑑𝜔 =
1 √2𝜋(𝛼+𝑗𝛽)
𝑒
𝑡2 2(𝛼+𝑗𝛽)
−
,
(7)
where 𝛼, 𝛽 ∈ ℝ, and 𝛼 ≥ 0 showing that the transform of a real or complex Gaussian is also a Gaussian; real or complex correspondingly. Continuing with the approximation, we do a change of variable 𝑡 ↦ 𝑡 − 𝑏, set 𝛼 = 0, 𝛽 = 2𝑐, and choose from the polar representation 𝛼 + 𝑗𝛽 = 𝑟𝑒 𝑗𝜃 the square root 1
𝜃
𝑟 2𝑒 𝑗2 . Now since 𝛼 = 0, 𝛽 = 2𝑐, we have 𝜋 2 2 𝑟 = √0 + (2𝑐) = 2|𝑐|, and 𝜃 = ± depending on the sign of 𝑐. Thus √𝛼 + 𝑗𝛽 = √2|𝑐|𝑒 𝑦(𝑡) ≈
1 2√𝜋|𝑐|
𝜃
±𝑗 𝜃 4
𝑒 −𝑗(𝑎−𝑠𝑖𝑔𝑛(𝑐)4 ) 𝑒 𝑗
2
, and
(𝑡−𝑏)2 4𝑐
,
(8)
which is a linear chirp signal with a delay of 𝑏, phase shift of 𝜃 1 𝑎 − 𝑠𝑖𝑔𝑛(𝑐) , and amplitude scaling of . 4
2√𝜋|𝑐|
1534-4320 (c) 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TNSRE.2016.2525829, IEEE Transactions on Neural Systems and Rehabilitation Engineering
3 3) Gaussian pulse response of the dispersion system The Fourier integral (7) gives us also the Fourier transform of Gaussian pulses 𝑥(𝑡) =
1 √2𝜋𝜏
𝑒
𝑡2 2𝜏
−
(9)
as 𝜔2
𝑋(𝜔) = 𝑒 −𝜏 2 ,
1
∞
∫ 𝑒 2𝜋 −∞
≈ 𝑒 −𝑗𝑎
1
−𝜏
𝜔2 2
∞
𝑒 −𝑗𝜙(𝜔) 𝑒 𝑗𝜔𝑡 𝑑𝜔
∫ 𝑒 2𝜋 −∞
𝜏 2
𝑗𝜔(𝑡−𝑏)−𝑗(𝑐+ )𝜔2
𝑑𝜔 ,
(11)
which, based on (7-8), we can see to actually be a linear chirp with a modified instantaneous frequency slope. 4) Amplitude modulation Exponentially decreasing amplitude can be expressed as ℎ(𝑡) = 𝑒 −𝜆|𝑡| .
(12)
It has the Fourier transform 𝐻(𝜔) =
2𝜆
(13)
𝜆2 +𝜔2
that decays rather rapidly with increasing frequency (low-pass type characteristics). A chirp signal 𝑓(𝑡) that undergoes decaying amplitude may then be expressed as the product 𝑥(𝑡) = 𝑓(𝑡)ℎ(𝑡).
(14)
Since a multiplication in the time domain is a convolution in the frequency domain, we see that ∞
𝑋(𝜔) = 𝐹(𝜔) ∗ 𝑋(𝜔) = ∫−∞ 𝐹(𝜓)𝑋(𝜔 − 𝜓)𝑑𝜓.
(15)
This causes a blurring of spectral peaks. We may also combine the causality of the chirp (𝑓(𝑡) = 0 when 𝑡 ≤ 0) by using a unit step function 0 𝑢(𝑡) = { 1
,𝑡 < 0 . ,𝑡 ≥ 0
.
(18)
𝑦(𝑡) =
∞ 1 ∞ ∫ ∫ 𝐻(𝜓)𝑒 −𝑗𝜙(𝜔−𝜓) 𝑋(𝜔 2𝜋 −∞ −∞
− 𝜓)𝑒 𝑗𝜔𝑡 𝑑𝜓 𝑑𝜔 ,(19)
where 𝐻(𝜔) is the spectral smearing function, for which ∡𝐻(𝜔) is constant and 𝐻(𝜔) decreases with increasing |𝜔|, 𝜙(𝜔) is the phase dispersion function (|𝑒 −𝑗𝜙(𝜔) | = 1, ∡𝑒 −𝑗𝜙(𝜔) = −𝜙(𝜔)), and X(𝜔) is the Fourier transform of the input signal. For a single frequency signal we can now analyze the response. A sinusoid 𝑓(𝑡) = sin 𝜔0 𝑡 has the Fourier transform 𝜋 𝐹(𝜔) = [𝛿(𝜔 − 𝜔0 ) − 𝛿(𝜔 + 𝜔0 )]. So the convolution with 𝑗
the spectral smearing basically gives just the spectral smearing function which is then phase distorted (or if dispersion is first, the sinusoids phase is just altered before the smearing). An important finding is that, considering narrower band input signals, we can interpret the spectral smearing here as introduction of new spectral components. I.e. the original input frequency will be accompanied by others near it with dispersed phases. Furthermore, the approximations for dispersion and amplitude modulation reveals that a chirp model accounting for both linear and exponential amplitude as well as frequency components can be used to accommodate both effects in a system model. Chirp models with the desired properties are explored next. C. Chirp signals To accommodate the system model we next consider some suitable linear and nonlinear chirp signals and envelope decay functions. 1) Linear chirp A so-called linear chirp is described with the quadratic phase function that can be produced by integration: 𝜙𝑙𝑖𝑛 (𝑡) = 𝜙0 + 2𝜋 ∫[𝑓0 + ∫ 𝑘 𝑑𝑡]𝑑𝑡 𝑘
= 𝜙0 + 2𝜋(𝑓0 𝑡 + 𝑡 2 ), 2
(20)
where 𝜙0 is the initial phase, 𝑓0 is the initial frequency, (16)
So the output we are considering is actually of the form 𝑢(𝑡)𝑓(𝑡). Hence, we may combine the unit step with the amplitude decay part as the multiplicative factor ℎ(𝑡) = 𝑢(𝑡)𝑒 −𝜆𝑡 ,
1 𝜆+𝑗𝜔
B. System model (based on the aforementioned) To summarize the findings above, if the system consists of phase distortion and spectral smearing (convolution with a low-pass like in the frequency domain), it can exhibit decaying chirp signal responses e.g. to impulses and Gaussian pulses. The system output can be modeled as
(10)
which can easily be seen by setting 𝛼 = 𝜏 and 𝛽 = 0. So the output for a Gaussian pulse input would be 𝑦(𝑡) =
𝐻(𝜔) =
(17)
where 𝑢(𝑡) is the unit step function. The Fourier transform of ℎ is now also a low-pass type function
𝑓 −𝑓
𝑘 = 1 0 defines the chirp rate, and 𝑇 is the chirp duration. 𝑇 The complex representation of a linear chirp is thus 𝑘 2 )]
𝑓𝑙𝑖𝑛 (𝑡) = 𝑒 𝑗[𝜙0 +2𝜋(𝑓0𝑡+2𝑡
.
(21)
In addition to this fundamental frequency, the frequency modulation gives also rise to additional frequency components. The utilization of complex valued signal is useful for computations, but it should be noted that the actual observed chirp is the real part of 𝑓, i.e.
1534-4320 (c) 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TNSRE.2016.2525829, IEEE Transactions on Neural Systems and Rehabilitation Engineering
4
𝑘
ℛℯ 𝑓𝑙𝑖𝑛 (𝑡) = cos [𝜙0 + 2𝜋 (𝑓0 𝑡 + 𝑡 2 )].
(22)
2
2) Polynomial nonlinear chirp A polynomial chirp is given by expansion of the linear chirp 𝑐𝑛 𝑛 phase function integral of the form ∑𝑁 𝑡 where 𝑐𝑛 are the 𝑛=0 𝑛!
polynomial coefficients. The first coefficient 𝑐0 corresponds to the initial phase 𝜙0 , 𝑐1 is the initial frequency 𝑓0 , and 𝑐2 = 𝑘 is the linear chirp rate. The linear chirp is thus the special case when 𝑁 = 3. Now with the chirp length 𝑇 and the linear part separated we get the complex representation 𝑓𝑝𝑜𝑙𝑦 (𝑡) = 𝑒
𝑗[𝜙0 +2𝜋(𝑓0 𝑡+
𝑐𝑛 𝑘 2 𝑛 𝑡 +∑𝑁 𝑛=3𝑛!∙𝑇𝑛−1 𝑡 )] 2𝑇
,
ℛℯ 𝑓𝑝𝑜𝑙𝑦 (𝑡) = cos [𝜙0 + 2𝜋 (𝑓0 𝑡 +
∑𝑁 𝑛=3
+ 𝑐𝑛
𝑛!∙𝑇 𝑛−1
𝑡 𝑛 )].
(24)
3) Exponential nonlinear chirp An exponential chirp of the form 𝑘 𝑡 can be expressed in an exponential form by noting that 𝑘 = 𝑒 ln 𝑘 so that 𝑘 𝑡 = 𝑡 (𝑒 ln 𝑘 ) = 𝑒 𝑡 ln 𝑘 . The exponential form can be written as a Maclaurian series 𝑒 𝑥 = ∑∞ 𝑛=0
𝑥𝑛 𝑛!
corresponding to a special (ln 𝑘)𝑛 𝑡 𝑛
case of infinite polynomial series 𝑒 𝑡 ln 𝑘 = ∑∞ where 𝑛=0 𝑛! the polynomial coefficients have a uniform structure 𝑐𝑛 = (ln 𝑘)𝑛 . Using an exponential instantaneous frequency function of 𝑓 the form 𝑓0 𝑘𝑡/𝑇 , where 𝑘 = 1 , we get, by integration and 𝑓0
substitution 𝑢 = 𝜏/𝑇 𝑡
𝜏
𝑡
∫0 𝑓0 𝑘 𝑇 𝑑𝜏 = 𝑓0 𝑇 ∫0 𝑘 𝑢 𝑑𝑢 = 𝑓0 𝑇
𝑘 𝑡/𝑇 −1 ln 𝑘
.
(25)
Now, the complex representation of the exponential chirp becomes 𝑓𝑒𝑥𝑝 (𝑡) = 𝑒
𝑘𝑡/𝑇 −1 )] ln 𝑘
𝑗[𝜙0 +2𝜋𝑓0 𝑇(
(26)
and the corresponding time-domain signal is then ℛℯ 𝑓𝑒𝑥𝑝 (𝑡) = cos [𝜙0 + 2𝜋𝑓0 𝑇 (
𝑘 𝑡/𝑇 −1 ln 𝑘
)].
(𝑎−𝑏) 𝑇
𝑡.
(28)
An exponential decay defined as 𝐴𝑒𝑥𝑝 (𝑡) = 𝐶𝑒 −𝜆𝑡 is a monotonically decreasing function, when 𝐶 is a positive constant and 𝜆 > 0 is the decay rate. Because 𝑒 0 = 1 and we can subtract 𝑒 −𝜆𝑇 to set the function to 0 at time 𝑇, the function can be expressed on the amplitude interval [𝑎, 𝑏] and (𝑎−𝑏) time interval [0, 𝑇] with a scale term 𝐶 = −𝜆∙0 −𝜆∙𝑇 = 𝑒
(𝑎−𝑏) 1−𝑒
−𝑒
−𝜆𝑇 . Now, with the scale term 𝐶, the end amplitude 𝑏, and
subtracting 𝑒 −𝜆𝑇 , the exponential dampening can be expressed as 𝐴𝑒𝑥𝑝 (𝑡) = 𝐶(𝑒 −𝜆𝑡 − 𝑒 −𝜆𝑇 ) + 𝑏 to get an exponential decay
(23)
where 𝑐𝑛 [𝑛 ≥ 3] are the coefficients for the higher nonlinear polynomials. The time-domain signal is then
𝑘 2 𝑡 2
𝐴𝑙𝑖𝑛 (𝑡) = 𝑎 −
(27)
4) Linear and exponential decay The observed chirps show a decay of signal. A linear decay 𝐴𝑙𝑖𝑛 (𝑡) = 𝐶 − 𝑟𝑡, where 𝐶 is a constant and 𝑟 is the rate of decay, can be defined in an amplitude interval [𝑎, 𝑏] at time interval [0, 𝑇] as
𝐴𝑒𝑥𝑝 (𝑡) = (𝑎 − 𝑏) [
𝑒 −𝜆𝑡 −𝑒 −𝜆𝑇 1−𝑒 −𝜆𝑇
] + 𝑏.
(29)
III. DATA A. Data collection For the SEP model creation and testing, experimental data were collected using a rat model with two varying factors, anesthetic level and SCI, both known to affect the signal [21,22] (Fig. 1). The experimental protocol was approved by the Institutional Animal Care and Use Committee (IACUC) at the National University of Singapore. The data were collected from 21 adult female SpragueDawley rats (200-220 g). Four screw electrodes were implanted to the cranium of rats above the somatosensory cortex of different limbs (Fig. 2). A reference electrode was positioned in the right parietal area close to lambda. Carboxylate dental cement (3M) was used in fixing the electrodes to the cranium. The animals were anesthetized with ketamine (75mg/kg) and xylazine (10 mg/kg) cocktail during the electrode placement. Baseline SEP recoding was carried out at least seven days after the electrode implantation. During the recording, rats were anesthetized with a mixture of isoflurane and 100% oxygen delivered at a rate of 1.3 L/min. A rodent-size anesthesia mask connected to a diaphragm with a C-Pram circuit designed to deliver and evacuate the gas was used for maintaining the anesthesia. To assess the SEPs at different anesthetic levels, the recording was carried out as the isoflurane dosage was increased in a step-like manner from 1.5% to 2.5% using 0.5% increments. The dosage was kept fixed for 5 min at each step before recording the SEPs to guarantee the equilibrium. A neuro-electrophysiology monitoring setup (Tucker-Davis Technology) was used for multilimb acquisition of SEP [23]. During the recording, four pairs of stainless-steel stimulating needle electrodes (Safelead F-E3-48, Grass Technologies) were placed in proximity to the tibial nerve in both hindlimbs and to the median nerve in both forelimbs. The needle electrodes were then connected to a stimulus generator (Digitimer Ltd.) and the cranial screw electrodes to an amplifier. The stimulus generator was used to generate 3.5 mA stimulating pulses with a pulse width of 200 µs. Pulses were generated at a rate of 1 Hz sequentially to each of the four limbs approximately 150 times per limb. During the
1534-4320 (c) 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TNSRE.2016.2525829, IEEE Transactions on Neural Systems and Rehabilitation Engineering
5 Furthermore, when assessing the effects of SCI, the signal analysis was restricted to the data of severe injury group.
Fig. 1. The effect of anesthesia and spinal cord injury on SEP. The figure shows the SEP of a single rat before (left column) and after (right column) severe spinal cord injury at three different anesthetic levels (different rows anesthetic effect increasing downwards).
B. Preprocessing In the preprocessing procedure applied, individual single sweep signals are averaged to provide a mean SEP signal. A Band Pass (BP) filter is applied to remove low frequency drifting of the mean signals and noise above and below the frequency band of interest. The corner frequencies of 𝑓𝑙 = 10 Hz and 𝑓ℎ = 280Hz were chosen to remove much of the noises, EEG channels from uncorrelated brain activities, and also all signal components not fitting in the analysis window (100 ms). A zero-phase IIR filtering method [26] with reversed and flipped signal adjustment conserving the DC value is used to preserve the absolute amplitude and phase information of the SEP pulses. Maximum absolute value is searched within a preset minimum and maximum delay envelope to identify most likely first SEP peak direction in order to limit the optimization search space and to increase the robustness of the model fitting. IV. METHODS
Fig. 2. Electrode locations for somatosensory evoked potential (SEP) measurement. The forelimb recording sites (circles) are located 0.2 mm posterior and 3.8 mm lateral to the bregma. The hindlimb recording sites (stars) are located 2.5 mm posterior and 2.8 mm lateral to the bregma. Reference electrode (square) is located 3 mm lateral to lambda.
stimulation, SEPs were recorded simultaneously from all four cortical screw electrodes via the amplifier and data acquisition setup. In the analysis, only the recordings from the corresponding electrode of the stimulated limb were used. Furthermore, for the current study, the signal analysis was restricted to the hindlimb data. After the baseline SEP recording, the animals were divided into three SCI groups (n = 5, 6, and 5 for mild, moderate, and severe injury, respectively) and one control group (n = 5). For the rats in the SCI groups, a controlled spinal cord contusion was performed using the NYU-MASCIS (New York University - Multicenter Animal Spinal Cord Injury Study) weight-drop impactor [24]. In this procedure, a 10 g rod with a flat circular impact surface was dropped from a certain height (6.25, 12.5, and 25 mm for mild, moderate, and severe injury groups, respectively) on the midline of the spinal cord exposed by removing the lamina at T8. A more detailed description of the controlled spinal cord contusion procedure can be found in e.g. [25]. For the rats in the control group, only the laminectomy without contusion was performed. After the operation, SEPs were recorded following the above described protocol weekly for 4 weeks. In this paper, only the data from baseline recordings and 28 days after the operation are used.
An automatic parameterization of SEPs is proposed using a chirp model where the resulting SEP waveforms are assumed as impulses undergoing group velocity dispersion while propagating along a multipath neural connection. The dispersion is assumed to cause changes to the signal phase spectra that manifests in the observations as chirp-like result waveforms. The chirp waveforms are assumed as non-periodic narrowband signals with strictly monotonic decay in time and frequency. The method for the automatic parametrization of SEPs with chirp modeling was implemented in Matlab™ technical computing language and is available online1. A. Proposed Chirp Models Chirp waveforms are modeled as monotonically strictly decreasing sinusoids using linear (22), 3rd or 4th degree polynomial (24) or exponential (27) instantaneous phase functions. A phase model with both linear and exponential component was also considered. The amplitude dampening weight function 𝐴(𝑡) based on the linear and/or exponential dampening (29), a unit step function (16), and added with an onset delay of 𝜏 is 𝐴(𝑡) = {
0 (𝑎 − 𝑏) [
,𝑡 < 𝜏 𝑒 −𝜆𝑤 (𝑡−𝜏) −𝑒 −𝜆𝑤 𝑇 1−𝑒 −𝜆𝑤 𝑇
]+𝑏
,𝑡 ≥ 𝜏
,
(30)
where 𝑎 = initial weight, 𝑏 = end weight, 𝑡 = time, 𝜏 = onset delay, 𝜆𝑤 = amplitude nonlinearity coefficient (exponential weight), 𝑇 = chirp duration. An illustrative fitted example dampening function (30) with parameters indicated is presented in Fig. 3. A chirp phase model using a linear (22) and an exponential chirp (27) component weighted by 𝜆𝑓 , and with an onset delay 1
Available: http://www.ee.oulu.fi/~jukortel/SEP.html
1534-4320 (c) 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TNSRE.2016.2525829, IEEE Transactions on Neural Systems and Rehabilitation Engineering
6 B. Model fitting and optimization Parameters (𝜏, 𝑎, 𝑏, 𝜆𝑤 , 𝑓0 , 𝑓1 , and 𝜆𝑓 ) are fitted for the preprocessed mean SEP signals using a PSO method [27]. The PSO Toolbox [28] was used to implement a Trelea type 1 PSO. Random initial seed positions for 150 particles with a 10000 epoch maximum iteration time was used for robust finding of the model parameters. The initial peak direction was fixed to positive side to narrow the search space into half. Search ranges for the parameters are: 𝜏 = 0.005 – 0.02 s, 𝑎 = ½·max|𝑥| – 1½·max|𝑥|, 𝑏 = 0 – 0.3·max|𝑥|, 𝜆𝑤 = 0.0001 – 5, 𝑓0 = 30 – 120 Hz, 𝑓1 = 0.01 – 15 Hz, 𝜆𝑓 = 0 – 1. Chirp duration 𝑇 = 80 ms was used. The goodness of fit (objective function to be minimized) used is the Mean Squared Error (MSE) that is linearly weighted (weight decreases towards the ends of the recordings) to emphasize fitting of the beginning parts of the chirps. The linear weight 𝑤𝑙𝑖𝑛 is defined as Fig. 3. Chirp amplitude model. The dashed line indicates a real fitted amplitude function A example.
𝑤𝑎
, 𝑡(𝑖) < 𝜏
𝑤𝑙𝑖𝑛 (𝑖) = {𝑤𝑎 − (𝑤𝑎 − 𝑤𝑏 ) 𝑤𝑏
𝑡(𝑖)−𝜏 𝑇
, 𝜏 ≤ 𝑡(𝑖) < 𝜏 + 𝑇,(33) , 𝑡(𝑖) > 𝜏 + 𝑇
where 𝑖 = sample index, 𝑤𝑎 = linear weight start, 𝑤𝑏 = linear weight stop. Linear weights 𝑤𝑎 = 1 and 𝑤𝑏 = 0.1 was used. The linearly weighted 𝑀𝑆𝐸𝑤 𝑙𝑖𝑛 is then defined as: 2
1
𝑀𝑆𝐸𝑤𝑙𝑖𝑛 = ∑𝑁 𝑖=1 𝑤𝑙𝑖𝑛 (𝑖) ∙ (𝑥(𝑖) − 𝑦(𝑖)) , 𝑁
(34)
where 𝑖 = sample index, 𝑁 = chirp/sample length, 𝑤𝑙𝑖𝑛 = linear weight, 𝑥 = original signal, 𝑦 = chirp signal. The average Root Mean Squared Error 1
Avg. 𝑅𝑀𝑆𝐸 = ∑ √𝑀𝑆𝐸 𝑁 Fig. 4. Instantaneous frequency of the linear and/or exponential chirp phase model. The dashed line indicates a real fitted example.
𝜏 and a duration of 𝑇 is defined as:
(35)
of the whole data and the average linearly weighted root mean squared error Avg. 𝑅𝑀𝑆𝐸𝑤𝑙𝑖𝑛 =
1 𝑁
∑ √𝑀𝑆𝐸𝑤𝑙𝑖𝑛
(36)
(𝑡−𝜏)
𝜙𝑒𝑥𝑝𝑙𝑖𝑛 (𝑡) = 𝜆𝑓 (2𝜋𝑓0 𝑇 (
𝑘 𝑇 −1 ln 𝑘
)) + (1 − 𝜆𝑓 )
∙ (2𝜋𝑓0 (𝑡 − 𝜏) + 2𝜋
(𝑓1 −𝑓0 ) 2𝑇
(𝑡 − 𝜏)2 ),
(31)
where 𝑓0 = initial frequency, 𝑓1 = end frequency, 𝜆𝑓 = frequency nonlinearity coefficient,
and
𝑓
𝑘 = ( 1 ). 𝑓0
An
illustrative fitted instantaneous frequency example of the linear and exponential phase function (31) with associated parameters indicated is presented in Fig. 4. The final chirp waveform 𝑦(𝑡) using (30) and (31) is then defined as: 𝑦(𝑡) = 𝐴(𝑡) ∙ cos (𝜙0 + 𝜙𝑒𝑥𝑝𝑙𝑖𝑛 (𝑡)) .
(32)
are then used to evaluate which chirp model is the best model overall using all available 252 averaged SEP samples (i.e. 21 rats with baseline and 28 days post operation samples for both hind limbs in three anesthetic levels). An illustrative example of a fitted chirp using the linear and/or exponential phase chirp model (32) that combines the previous examples, dampening function (Fig. 3) and instantaneous frequency function (Fig. 4), into the final chirp model, associated parameters, and derived SEP features is presented in Fig. 5. C. SNR Testing SNR performance testing of the chirp model fitting was performed using all of the 252 averaged chirp recordings as clean signal templates. A bootstrap approach was used to
1534-4320 (c) 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TNSRE.2016.2525829, IEEE Transactions on Neural Systems and Rehabilitation Engineering
7
Fig. 5. An illustrative example of a fitted chirp model with the SEP features showed using the linear and/or exponential phase chirp model (31). An averaged SEP measurement x (solid gray line) is shown with the fitted chirp model y (solid black line) that also corresponds to the fitted example amplitude function A (Fig. 3) and the fitted example instantaneous frequency function in Fig. 4.
evaluate the stability of the averaged templates without added noise to evaluate reasonable lower limits for model fitting errors corresponding to the 150 averaged SEP signal fittings. 1,000 bootstrap means of 150 random overlapping selections of individual single sweep signals were taken from each 252 samples (i.e. a total of 252,000 bootstraps). Mean variances along all of the samples were calculated from the corresponding 1000 bootstrap means. The mean variances correspond to an unbiased estimators MSE and can then be used to produce estimates for the MSE and linearly weighed 𝑀𝑆𝐸𝑤𝑙𝑖𝑛 expected with the measurement process corresponding to the generation of the templates (i.e. an estimate of intrinsic error of the templates due to noises or other uncorrelated processes regardless of a SEP model). For the SNR testing, an array of Additive Gaussian White Noise (AWGN) measured at -12, -6, 0, 6, 12, 18, and 24 dB SNR levels (i.e. in relation to each template signal power to stabilize the SNR at each level) was then applied to produce test signals. The selected SNR dB levels correspond to SNR increases attained when using multiple trial averaging of 4 n signal sweeps between n dB levels (e.g. 4 1 = 4 averages from 0dB to 6dB and 43 = 64 averages from -6dB to 12dB), in the presence of AWGN. A Monte Carlo method using 10 independent AWGNs for each sample was implemented to produce stable SNR performance readings at each dB level. In the testing a total of 2520 independent noise samples were fitted with a chirp model for each selected SNR level using both the 𝑀𝑆𝐸 and the 𝑀𝑆𝐸𝑤𝑙𝑖𝑛 as measures for the goodness of fit when comparing to the corresponding clean template. The 𝑅𝑀𝑆𝐸 and the 𝑅𝑀𝑆𝐸𝑤𝑙𝑖𝑛 of each of the 252 samples were calculated for each SNR level tested. To evaluate how good the chirp fitting is for individual signals three error rate criterions using normalized errors were used (i.e. in order to measure erroneous waveform fitting rather than different amplitudes). First, For each 252 sample, the bootstrap noise estimate was normalized with the clean template signal estimate (using both 𝑅𝑀𝑆𝐸 and the
𝑅𝑀𝑆𝐸𝑤𝑙𝑖𝑛 ) to produce intrinsic normalized bootstrap error estimates (both normalized 𝑅𝑀𝑆𝐸 and normalized 𝑅𝑀𝑆𝐸𝑤𝑙𝑖𝑛 ) for the samples. Next, for each SNR level, chirp fitting errors were similarly normalized using the corresponding clean template signal estimates to produce normalized chirp model fitting errors. Then, the differences of normalized chirp model fitting errors at each SNR level and the normalized bootstrap error estimates were calculated to produce fitting Δ error distributions (both Δ normalized 𝑅𝑀𝑆𝐸 and Δ normalized 𝑅𝑀𝑆𝐸𝑤𝑙𝑖𝑛 ) at each SNR level. The three thresholds for fitting errors were chosen as Δ errors larger than the median (50%), 3rd quartile (75%), or 9th decile (90%) values of the normalized bootstrap error distribution (both 𝑅𝑀𝑆𝐸 and 𝑅𝑀𝑆𝐸𝑤𝑙𝑖𝑛 ). The choice of thresholds allows for a strict, average, and lax evaluation of bounds, correspondingly, for a chirp fitting that takes into account the estimated intrinsic noise in each template (with no weighting by using the 𝑅𝑀𝑆𝐸 or emphasizing the beginning by using the 𝑅𝑀𝑆𝐸𝑤𝑙𝑖𝑛 ). D. Feature Calculation The resulting chirp model parameters can be used as features describing the SEP recordings. Specifically, the onset delay (𝜏) and amplitude parameters (𝑎, 𝑏 and 𝜆𝑤 ) but also the frequency parameters (𝑓0 , 𝑓1 and 𝜆𝑓 ) are of interest. The goodness-of-fit measure gained from the optimization, i.e. the MSE measure, can also be used to calculate a normalized sum of square residual (𝑁𝑆𝑆𝑅𝑒𝑠 ) feature relating to the ratio of the model explained signal and the residual unexplained signal 𝑁𝑆𝑆𝑅𝑒𝑠 =
2
∑𝑁 𝑖=1(𝑥(𝑖)−𝑦(𝑖)) 2 ∑𝑁 𝑖=1 𝑦(𝑖)
.
(37)
The fitted model can also be used to find the initial peaks and peak latencies of the chirp-like signal (see Fig. 5), either from the model itself, which is a trivial task of finding the zeros, minimums, and maximums due to the model properties (e.g. the chirps have no riding waves) or, using the model peak locations as an initial estimate, from the original SEP signal (i.e. by local search and fitting). Parameters of the fitted chirp signals and the residual goodness-of-fit measure are presented in Table I. SEP features derived from the fitted chirp models are presented in Table II. E. Anesthesia level and SCI analysis Two testing groups were used to evaluate the parameters of the fitted chirp models and also features derived from the fitted chirp models. In the first group of tests, i.e. anesthesia level, all the 126 averaged baseline recordings (i.e. 21 rats with two hind limb recordings in three anesthetic levels) were used to evaluate the effect of changing anesthetic concentration on the SEP parameters and features. Statistical testing was performed to determine if the different anesthetic levels have a significant impact on the measured SEP recordings between the anesthetic concentration levels used: 1.5% vs. 2.0%, 2.0% vs. 2.5%, 1.5% vs. 2.5%, and 1.5% vs. 2.0% vs. 2.5% (i.e. all levels).
1534-4320 (c) 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TNSRE.2016.2525829, IEEE Transactions on Neural Systems and Rehabilitation Engineering
8 TABLE I CHIRP MODEL PARAMETERS Parameter
Quantity
𝑇
Chirp length
𝜏
Chirp onset delay
𝑎
Start amplitude
𝑏
End amplitude
𝜆𝑤
Amplitude nonlinearity coefficient
𝑓0
Start frequency
𝑓1
End frequency
𝜆𝑓
Frequency nonlinearity coefficient
𝑁𝑆𝑆𝑅𝑒𝑠
Normalized sum of squares residual
evaluate the generated p-values in order to counteract the problem of multiple comparisons. V. RESULTS
TABLE II SEP FEATURES DERIVED FROM THE CHIRP MODEL Feature
Quantity
st
1 peak delay
Delay of the first peak
1st peak
First peak amplitude (absolute)
2nd peak delay
Delay of the second peak
2nd peak
Second peak amplitude (absolute)
p2p delay
1st peak to 2nd peak delay
p2p
1st peak to 2nd peak amplitude
TABLE III AVERAGE MSES OF THE TESTED CHIRP MODELS 𝐴𝑣𝑔. 𝑅𝑀𝑆𝐸𝑤𝑙𝑖𝑛 (µV) Chirp model 𝐴𝑣𝑔. 𝑅𝑀𝑆𝐸 (µV) Linear
4.88
3.76
3rd degree polynomial
4.45
3.37
4 degree polynomial
3.79
2.78
Exponential
3.86
2.72
Linear & exponential
3.66
2.63
Bootstrap noise estimate
3.00
2.12
th
First, the fitting of all available SEP data (252 averaged signals) was performed using linear, 3rd and 4th degree polynomial, exponential, and both linear and exponential phase chirp functions with linear and/or exponential amplitude dampening functions. The best performing model, i.e. the model resulting in the overall lowest 𝐴𝑣𝑔. 𝑅𝑀𝑆𝐸 and linearly weighted 𝐴𝑣𝑔. 𝑅𝑀𝑆𝐸𝑤𝑙𝑖𝑛 , presented in Table III, was selected as the proposed chirp model and used in the further tests. The bootstrap analysis of the averaged SEP signal variances resulted in an 𝐴𝑣𝑔. 𝑅𝑀𝑆𝐸 = 3.00 µV corresponding to an average SNR of 15.69 dB for the 150 sample average templates with a bootstrap mean estimate for signal 𝐴𝑣𝑔. 𝑅𝑀𝑆𝐸 = 18.25 µV. The template average SNR also corresponds to an average SNR of -6.19 dB for a single trial. When considering the template SNR, the added AWGN SNR array of -12, -6, 0, 6, 12, 18, and 24 dB, used in the Monte Carlo runs, therefore corresponds to an array of -12.01, -6.06, -0.22, 5.27, 9.94, 13.17, and 14.82 dB estimated total average SNRs for the 150 averaged samples, respectively. A corresponding single trial SNR can be calculated by subtracting 21.76 dB from each estimate SNR (i.e. the effect of 150 average signals assuming AWGS). An 𝑅𝑀𝑆𝐸 correlation analysis between the bootstrap variances and the template fitting errors with the linear and/or exponential chirp model for each 252 samples revealed a correlation of 0.73 (p
