Open Access
Methods  |   August 2017
A nonlinear generalization of the Savitzky-Golay filter and the quantitative analysis of saccades
Author Affiliations
Journal of Vision August 2017, Vol.17, 10. doi:https://doi.org/10.1167/17.9.10
  • Views
  • PDF
  • Share
  • Tools
    • Alerts
      ×
      This feature is available to authenticated users only.
      Sign In or Create an Account ×
    • Get Citation

      Weiwei Dai, Ivan Selesnick, John-Ross Rizzo, Janet Rucker, Todd Hudson; A nonlinear generalization of the Savitzky-Golay filter and the quantitative analysis of saccades. Journal of Vision 2017;17(9):10. https://doi.org/10.1167/17.9.10.

      Download citation file:


      © ARVO (1962-2015); The Authors (2016-present)

      ×
  • Supplements
Abstract

The Savitzky-Golay (SG) filter is widely used to smooth and differentiate time series, especially biomedical data. However, time series that exhibit abrupt departures from their typical trends, such as sharp waves or steps, which are of physiological interest, tend to be oversmoothed by the SG filter. Hence, the SG filter tends to systematically underestimate physiological parameters in certain situations. This article proposes a generalization of the SG filter to more accurately track abrupt deviations in time series, leading to more accurate parameter estimates (e.g., peak velocity of saccadic eye movements). The proposed filtering methodology models a time series as the sum of two component time series: a low-frequency time series for which the conventional SG filter is well suited, and a second time series that exhibits instantaneous deviations (e.g., sharp waves, steps, or more generally, discontinuities in a higher order derivative). The generalized SG filter is then applied to the quantitative analysis of saccadic eye movements. It is demonstrated that (a) the conventional SG filter underestimates the peak velocity of saccades, especially those of small amplitude, and (b) the generalized SG filter estimates peak saccadic velocity more accurately than the conventional filter.

Introduction
The smoothing and differentiation of time series is important in numerous fields. The Savitzky-Golay (SG) filter (Savitzky & Golay, 1964) is widely used for this purpose, especially for biomedical data. For example, the SG filter has been advocated for electroencephalography and electrocardiography (Azami, Mohammadi, & Bozorgtabar, 2012), elastography (Luo, Bai, He, & Ying, 2004), near infrared spectroscopy (Schneider & Kovar, 2003), functional magnetic resonance imaging (Geissler et al., 2007), speech enhancement (Shajeesh, Kumar, Pravena, & Soman, 2012), and eye movement analysis (Nyström & Holmqvist, 2010; Quaia, Joiner, FitzGibbon, Optican, & Smith, 2010; Shaikh, Xu-Wilson, Grill, & Zee, 2010). 
Despite its widespread use, the SG filter obeys an unavoidable trade-off between noise suppression and signal distortion, as does any linear filter. Correspondingly, the SG filter tends to oversmooth time series that exhibit abrupt deviations such as sharp waves or steps, yet these are often of physiological interest (e.g., a saccade is an abrupt change in eye position). Consequently, the SG filter can lead to physiological parameters being systematically underestimated (e.g., the peak angular velocity of the eye during a saccade). Accurate estimation of such parameters is of interest, since they are used in clinical research and practice. 
Saccades are rapid eye movements redirecting sight from one object of interest to another (Leigh & Zee, 2015). These movements are brief, rarely lasting more than 80 ms (Leigh & Zee, 2015). Advances in the study of anatomy and pathology of the brain and eye movements have reinforced the substantial utility of saccades in clinical practice. Abnormalities of saccades offer important clues in the diagnosis of a number of neurological disorders (Ramat, Leigh, Zee, & Optican, 2007). Quantitative saccade parameters include latency, amplitude, duration, peak velocity, and peak acceleration, etc. Among these measures, the peak velocity of saccades is of particular interest. Slow saccadic eye movements may suggest a brainstem lesion in excitatory saccade burst neurons (Horn & Büttner-Ennever, 1998; Horn, Büttner-Ennever, Suzuki, & Henn, 1995; Ramat, Leigh, Zee, Shaikh, & Optican, 2008). Examples of conditions that typically cause saccadic slowing include progressive supranuclear palsy (Chen et al., 2010; Leigh & Riley, 2000) or spinocerebellar ataxia type 2 (Garbutt, Harwood, Kumar, Han, & Leigh, 2003; Wadia & Swami, 1971). Additional factors, such as saccade adaption and visual salience, may also affect saccade peak velocity (Ethier, Zee, & Shadmehr, 2008; Schütz, Braun, & Gegenfurtner, 2011); even subtle reductions in the peak velocity of otherwise normal saccades have been observed in states of mental fatigue (Di Stasi et al., 2012). Thus, accurate determination of saccade peak velocity is critical to the understanding of normal and pathological saccade behavior. 
We propose a generalization of the SG filter that can track abrupt deviations in time series more accurately than the conventional SG filter. The generalized SG filter leads to more accurate parameter estimates when it is desirable to measure abrupt changes, such as the peak velocity of saccadic eye movements. The proposed filtering methodology models a time series as the sum of two component time series: a low-frequency time series for which the conventional SG filter is well suited, and a second time series that exhibits instantaneous deviations (e.g., sharp waves, steps, or more generally, discontinuities in a higher-order derivative). We formulate the generalized SG filter via sparse signal modeling. The output of the filter is expressed in terms of the solution to a sparse-regularized linear inverse problem, and is calculated via convex optimization. The generalized SG filter is thus nonlinear and can overcome the limitations of the conventional (linear) SG filter. 
We model an abrupt deviation in a time series as a discontinuity either in the underlying signal, in the derivative of the underlying signal, or in a higher order derivative of the underlying signal. In turn, the next-higher order derivative of the time series exhibits impulses (isolated values of high amplitude). Such a time series therefore has a sparse component (a time series is said to be sparse if it is mostly zero in value, except for isolated values of high amplitude). Accordingly, in the proposed generalization of the SG filter, we adopt tools from sparse signal processing (Elad, 2010). An early instance of filtering based on a sparse signal model is total variation denoising (Rudin, Osher, & Fatemi, 1992), which models a step signal as having a sparse first-order derivative. 
The proposed generalization of the SG filter is closely related to the sparsity-assisted signal smoothing (SASS) method (Selesnick, 2015, 2017; Selesnick, Graber, Pfeil, & Barbour, 2014). The SASS method can be considered a generalization of the digital Butterworth filter (a particular recursive filter), while the filter proposed in this work constitutes a generalization of the SG filter (a particular nonrecursive filter). Additionally, SASS was developed only for smoothing, while the filter proposed here is also developed for differentiation, which is important for the estimation of velocity time series and associated parameters (e.g., peak velocity). 
In this article, we demonstrate that the proposed generalization of the SG filter can estimate peak saccadic velocity from eye position time series more accurately than the conventional SG filter. The evaluation of the proposed filter uses (a) simulated data based on a recently developed parametric model for saccadic eye movements (Dai, Selesnick, Rizzo, Rucker, & Hudson, 2016), (b) physiologic data of saccadic eye movements that have been objectively recorded and manually identified by two experts, Nyström and Andersson (2017), and (c) data recorded in our lab (Rizzo et al., 2016). 
Methods
The SG smoothing filter
The SG smoothing filter is a particular low-pass filter defined by two parameters that we shall denote K and M (Orfanidis, 1995; Savitzky & Golay, 1964; Schafer, 2011b). The SG filter can be defined and implemented as a weighted moving average, i.e., a finite impulse response (FIR) filter. We denote by x(n) the time series to be estimated, and we suppose the observed time series is given by y(n) = x(n) + w(n) where w(n) represents additive white Gaussian noise. Applied to time series y(n), the SG filter produces an output time series Display Formula\(\def\upalpha{\unicode[Times]{x3B1}}\)\(\def\upbeta{\unicode[Times]{x3B2}}\)\(\def\upgamma{\unicode[Times]{x3B3}}\)\(\def\updelta{\unicode[Times]{x3B4}}\)\(\def\upvarepsilon{\unicode[Times]{x3B5}}\)\(\def\upzeta{\unicode[Times]{x3B6}}\)\(\def\upeta{\unicode[Times]{x3B7}}\)\(\def\uptheta{\unicode[Times]{x3B8}}\)\(\def\upiota{\unicode[Times]{x3B9}}\)\(\def\upkappa{\unicode[Times]{x3BA}}\)\(\def\uplambda{\unicode[Times]{x3BB}}\)\(\def\upmu{\unicode[Times]{x3BC}}\)\(\def\upnu{\unicode[Times]{x3BD}}\)\(\def\upxi{\unicode[Times]{x3BE}}\)\(\def\upomicron{\unicode[Times]{x3BF}}\)\(\def\uppi{\unicode[Times]{x3C0}}\)\(\def\uprho{\unicode[Times]{x3C1}}\)\(\def\upsigma{\unicode[Times]{x3C3}}\)\(\def\uptau{\unicode[Times]{x3C4}}\)\(\def\upupsilon{\unicode[Times]{x3C5}}\)\(\def\upphi{\unicode[Times]{x3C6}}\)\(\def\upchi{\unicode[Times]{x3C7}}\)\(\def\uppsy{\unicode[Times]{x3C8}}\)\(\def\upomega{\unicode[Times]{x3C9}}\)\(\def\bialpha{\boldsymbol{\alpha}}\)\(\def\bibeta{\boldsymbol{\beta}}\)\(\def\bigamma{\boldsymbol{\gamma}}\)\(\def\bidelta{\boldsymbol{\delta}}\)\(\def\bivarepsilon{\boldsymbol{\varepsilon}}\)\(\def\bizeta{\boldsymbol{\zeta}}\)\(\def\bieta{\boldsymbol{\eta}}\)\(\def\bitheta{\boldsymbol{\theta}}\)\(\def\biiota{\boldsymbol{\iota}}\)\(\def\bikappa{\boldsymbol{\kappa}}\)\(\def\bilambda{\boldsymbol{\lambda}}\)\(\def\bimu{\boldsymbol{\mu}}\)\(\def\binu{\boldsymbol{\nu}}\)\(\def\bixi{\boldsymbol{\xi}}\)\(\def\biomicron{\boldsymbol{\micron}}\)\(\def\bipi{\boldsymbol{\pi}}\)\(\def\birho{\boldsymbol{\rho}}\)\(\def\bisigma{\boldsymbol{\sigma}}\)\(\def\bitau{\boldsymbol{\tau}}\)\(\def\biupsilon{\boldsymbol{\upsilon}}\)\(\def\biphi{\boldsymbol{\phi}}\)\(\def\bichi{\boldsymbol{\chi}}\)\(\def\bipsy{\boldsymbol{\psy}}\)\(\def\biomega{\boldsymbol{\omega}}\)\(\hat x(n)\) given by  
\begin{equation}\tag{1}\hat x(n) = \sum\limits_{k = - M}^M h (k)\,y(n - k)\end{equation}
where M is a parameter of the filter and h(n) is the impulse response of the filter (defined for Display Formula\(|n|\, \le \,M)\).  
The SG filter is defined via a least-squares polynomial approximation problem. Specifically, the SG filter defines Display Formula\(\hat x(0)\) (i.e., the output at time n = 0) to be the value of the coefficient of a polynomial of order K that best fits the time series data y(n) over the interval Display Formula\(|n|\, \le \,M\). Given the time series y(n) and the parameters K and M, define the polynomial  
\begin{equation}\tag{2}p(n) = \sum\limits_{k = 0}^K {{c_k}} {n^k}\end{equation}
as the polynomial minimizing the square error E,  
\begin{equation}\tag{3}E = \sum\limits_{n = - M}^M ( y(n) - p(n){)^2}.\end{equation}
The output value Display Formula\(\hat x(0)\) of the SG filter at time n = 0 is defined to be p(0)—that is, the value of the best-fitting polynomial at n = 0. Hence, for the SG filter, we have Display Formula\(\hat x(0)\) = p(0) = c0 where c0 is a coefficient of the best-fitting polynomial. The output values Display Formula\(\hat x(n)\) at time instants other than n = 0 are similarly defined. This process wholly determines the impulse response h(n) of the filter (Orfanidis, 1995). We note that the impulse response of the SG filter is symmetric—that is, h(–n) = h(n). Given M, the parameter K must be chosen so that Display Formula\(0 \le K \le 2M.\)  
The transfer function of a filter with impulse response h(n) is defined as  
\begin{equation}\tag{4}H(z) = \sum\limits_n h (n)\,{z^{ - n}}\end{equation}
and the frequency response is given by Display Formula\(H({{\it{e}}^{{\it{j}}\omega }})\). Since the impulse response of the SG filter is symmetric, its frequency response is real-valued, and it is a zero-phase filter. Figure 1a shows the frequency response of the SG smoothing filter for several parameter values. The cut-off frequency of the SG filter depends on the parameters K and M, as described by Schafer (2011a). We highlight a property of the SG smoothing filter we use below.  
Figure 1
 
The frequency response of Savitzky-Golay (a) smoothing and (b) differentiation filters for several parameter values.
Figure 1
 
The frequency response of Savitzky-Golay (a) smoothing and (b) differentiation filters for several parameter values.
Property 1
Let H(z) be the transfer function of a SG smoothing filter with polynomial order parameter K. Then H(z) can be expressed as  
\begin{equation}\tag{5}H(z) = 1 - {(1 - {z^{ - 1}})^{K + 1}}Q(z)\end{equation}
where Q(z) is the transfer function of another FIR filter. This property implies that the frequency response of the SG smoothing filter is flat (of order K) at ω = 0 (Orfanidis, 1995). The flatness of the frequency response can be seen in Figure 1a.  
The SG differentiation filter
SG filters are used not only for smoothing but also for differentiation of time series. Like the SG smoothing filter, the SG differentiation filter is a particular zero-phase low-pass differentiation filter defined by two parameters, K and M, and it can be implemented as an FIR filter. The SG differentiation filter is defined via the differentiation of the polynomial that best fits the data, determined already for the SG smoothing filter. Differentiation of the polynomial p in Equation 2 yields  
\begin{equation}\tag{6}p^{\prime} (n) = \sum\limits_{k = 1}^K {{{\bar c}_k}} \,{n^{k - 1}}\end{equation}
where Display Formula\({\bar c_k} = k{c_k}\). Like the SG smoothing filter, the impulse response Display Formula\(\bar h(n)\) of the SG differentiation filter is determined by parameters K and M. Figure 1b shows the frequency response of the SG differentiation filter for several parameter values. We highlight a property of the SG differentiation filter we use below.  
Property 2
Let Display Formula\(\bar H(z)\) be the transfer function of a SG differentiation filter with polynomial order K. Then Display Formula\(\bar H(z)\) can be expressed as  
\begin{equation}\tag{7}\bar H(z) = (1 - {z^{ - 1}}) + {(1 - {z^{ - 1}})^{K + 1}}U(z)\end{equation}
where U(z) is the transfer function of another FIR filter. This property implies that the frequency response of the SG differentiation filter is tangent (of order K) to a line at Display Formula\(\omega {\rm{\ }} = {\rm{\ }}0\) (Luo, Ying, He, & Bai, 2005) The high-order tangency of the frequency response can be seen in Figure 1b.  
A generalization of the smoothing filter
In this section, we derive a framework for extending the SG smoothing filter. We first express the conventional SG filter using linear algebra and define related filters. We then formulate the proposed filter using sparse-regularized optimization. 
Let H(z) be a SG smoothing filter. Since the SG filter is linear, we express the output time series Display Formula\(\hat x\left( n \right)\) in Equation 1 via matrix-vector multiplication,  
\begin{equation}\tag{8}\hat x = Hy\end{equation}
where y is a vector of time series values y(n) and H is a Toeplitz matrix. For example, if M = 1, then the matrix H has the form  
\begin{equation}\tag{9}H = \left[ {\matrix{ {h(0)}&{h(1)}&{}&{}&{} \cr {h( - 1)}&{h(0)}&{h(1)}&{}&{} \cr {}&\ddots&\ddots&\ddots&{} \cr {}&{}&{h( - 1)}&{h(0)}&{h(1)} \cr {}&{}&{}&{h( - 1)}&{h(0)} \cr } } \right].\end{equation}
Since H is a zero-phase low-pass filter, the filter G defined as  
\begin{equation}\tag{10}G = I - H,\end{equation}
is a zero-phase high-pass filter (I is an identity matrix, equivalently the identity filter). Its transfer function is G(z) = 1 – H(z).  
We define DS as the S-order difference matrix. The matrix D1 has the form  
\begin{equation}\tag{11}{D_1} = \left[ {\matrix{ { - 1}&1&{}&{}&{} \cr {}&{ - 1}&1&{}&{} \cr {}&{}&\ddots&\ddots&{} \cr {}&{}&{}&{ - 1}&1 \cr } } \right].\end{equation}
The matrix D1 is a discrete form of first-order differentiation, and D1x is a discrete first-order derivative of the time series x. The matrix D3 has the form  
\begin{equation}\tag{12}{D_3} = \left[ {\matrix{ { - 1}&3&{ - 3}&1&{}&{}&{} \cr {}&{ - 1}&3&{ - 3}&1&{}&{} \cr {}&{}&\ddots&\ddots&\ddots&\ddots&{} \cr {}&{}&{}&{ - 1}&3&{ - 3}&1 \cr } } \right].\end{equation}
The transfer function of the S-order difference filter is given by  
\begin{equation}\tag{13}{D_S}(z) = {(1 - {z^{ - 1}})^S}.\end{equation}
Using Equations 5 and 10, we write the transfer function of the high-pass filter G as  
\begin{equation}\tag{14}G(z) = 1 - H(z)\end{equation}
 
\begin{equation}\tag{15} = {(1 - {z^{ - 1}})^{K + 1}}Q(z).\end{equation}
If Display Formula\(S \le K + 1\), then using Equations 13 and 15, we may factor DS(z) out of G(z),  
\begin{equation}\tag{16}G(z) = {(1 - {z^{ - 1}})^{K + 1 - S}}Q(z){D_S}(z)\end{equation}
 
\begin{equation}\tag{17} = R(z){D_S}(z),\end{equation}
where R(z) is the transfer function of another FIR filter. We also write G = RDS where R is a matrix correspondingly defined.  
The two-component model
We assume the time series data x to be estimated can be modeled as  
\begin{equation}\tag{18}x = {x_1} + {x_2}\end{equation}
where x1 and x2 are two component time series. In this model, we assume the time series x1 has the property that DSx1 is sparse, and that x2 is a low-frequency time series for which conventional SG filters are well suited. We assume the time series data is given by  
\begin{equation}\tag{19}y = x + w\end{equation}
 
\begin{equation}\tag{20} = {x_1} + {x_2} + w\end{equation}
where w is additive white Gaussian noise.  
We assumed x2 is a time series for which the SG filter is well suited. By this, we mean the SG filter H preserves the time series x2 and suppresses white noise w, (i.e., we have x2H(x2 + w). Using Equation 20, we have  
\begin{equation}\tag{21}{x_2}{\rm{\ }} \approx H({x_2} + w)\end{equation}
 
\begin{equation}\tag{22} = H({y - {x_1}}).\end{equation}
Hence, we propose to estimate x2 in terms of an estimate Display Formula\({\hat x_1},\)  
\begin{equation}\tag{23}{\hat x_2} = H(y - {\hat x_1}).\end{equation}
In turn, we express an estimate of the unknown time series x as  
\begin{equation}\tag{24}\hat x = {\hat x_1} + {\hat x_2}\end{equation}
 
\begin{equation}\tag{25} = {\hat x_1} + H\left( {y - {{\hat x}_1}} \right)\end{equation}
 
\begin{equation}\tag{26} = (I - H){\hat x_1} + Hy\end{equation}
 
\begin{equation}\tag{27} = G{\hat x_1} + Hy\end{equation}
 
\begin{equation}\tag{28} = R{D_S}{\hat x_1} + Hy\end{equation}
where G is the high-pass filter (Equation 10), which admits the factorization (Equation 17). By assumption, Display Formula\({D_S}{\hat x_1}\) is sparse. Hence, we define û = Display Formula\({D_S}{\hat x_1}\) and write  
\begin{equation}\tag{29}\hat x = R\hat u + Hy.\end{equation}
Since we model û as sparse, we estimate û via standard sparse-regularized least squares (Elad, 2010). That is, û is calculated as the solution to the optimization problem  
\begin{equation}\tag{30}\hat u = \arg \ \mathop {\min }\limits_u \left\{ {{1 \over 2}||y - (Ru + Hy)||_2^2 + \lambda ||u|{|_1}} \right\}\end{equation}
where  
\begin{equation}\tag{31}||x||_2^2: = \sum\limits_n x {(n)^2}\end{equation}
 
\begin{equation}\tag{32}||x|{|_1}: = \sum\limits_n {|x(n)|{\rm{\ }}} \end{equation}
and λ is a positive parameter. The parameter λ should be chosen according to the noise level. Several algorithms are available to solve Equation 30, e.g., the forward-backward splitting algorithm (Combettes & Pesquet, 2011). Once Display Formula\(\hat u\) is obtained, the estimated time series Display Formula\(\hat x\) is then given by Equation 29.  
A generalization of the differentiation filter
In this section, we use the framework developed in the previous section to extend the SG differentiation filter. 
We write the derivative of the time series x as x′ and we express it as  
\begin{equation}\tag{33}x^{\prime} = \bar Fx\end{equation}
where Display Formula\(\bar F\) is a zero-phase full-band differentiation filter (i.e., its impulse response is symmetric and it accurately differentiates high-frequency sinusoids). The ideal full-band differentiation filter cannot be realized by an FIR filter, but it can be well approximated by one (Al-Alaoui, 2007; Kumar, Roy, & Shah, 1992; Pei & Wang, 2001; Rabiner & Schafer, 1974; Selesnick, 2002). In fact, the SG differentiation filter with parameters K = 2M is an approximate full-band differentiation filter. Therefore, we let Display Formula\(\bar F\) be a zero-phase full-band differentiation FIR filter. A problem in practice with full-band differentiators is that they severely amplify noise. Hence, by a differentiation filter, we mean by default a low-pass differentiation filter. The SG filter with parameters K < 2M is a low-pass differentiation filter.  
Let Display Formula\(\bar h(n)\) be the impulse response of the SG differentiation filter with parameters K and M. If the input of the filter is the time series y, then the output of the filter is given by the weighted moving average,  
\begin{equation}\tag{34}\hat x^{\prime} (n) = \sum\limits_{k = - M}^M {\bar h(k)y(n - k).} \end{equation}
As in Equation 8, we express the output via matrix-vector multiplication,  
\begin{equation}\tag{35}\hat x^{\prime} = \bar Hy\end{equation}
where Display Formula\(\bar H\) is a matrix determined by the SG differentiation filter.  
In the previous section, we subtracted the zero-phase low-pass filter H from the identity, I, to define a zero-phase high-pass filter G (i.e., G = IH). Here, we subtract the zero-phase low-pass differentiation filter Display Formula\(\bar H\) from a zero-phase full-band differentiation filter Display Formula\(\bar F\), to define  
\begin{equation}\tag{36}\bar G = \bar F - \bar H.\end{equation}
The filter Display Formula\(\bar G\) is a zero-phase high-pass differentiation filter. Here, we let Display Formula\(\bar F\) be a SG differentiation filter with parameters K = 2M. Then the transfer function of Display Formula\(\bar F\) can be written as  
\begin{equation}\tag{37}\bar F(z) = (1 - {z^{ - 1}}) + {(1 - {z^{ - 1}})^{2M + 1}}V(z)\end{equation}
where V(z) is the transfer function of another FIR filter. Using Equations 7 and 37, we write the transfer function of filter Display Formula\(\bar G\) as  
\begin{equation}\tag{38}\bar G(z) = \bar F(z) - \bar H(z)\end{equation}
 
\begin{equation}\tag{39} = (1 - {z^{ - 1}}) + {(1 - {z^{ - 1}})^{2M + 1}}V(z) - (1 - {z^{ - 1}}) - {(1 - {z^{ - 1}})^{K + 1}}U(z)\end{equation}
 
\begin{equation}\tag{40} = {(1 - {z^{ - 1}})^{2M + 1}}V(z) - {(1 - {z^{ - 1}})^{K + 1}}U(z).\end{equation}
If K ≤ 2M, then  
\begin{equation}\tag{41}\bar G(z) = {\left( {1 - {z^{ - 1}}} \right)^{K + 1}}W(z)\end{equation}
where W(z) = (1 – z–1)2MKV(z) – U(z).  
If SK + 1, then using Equations 13 and 41, we may factor DS(z) out of Display Formula\(\bar G(z)\),  
\begin{equation}\tag{42}\bar G(z) = {(1 - {z^{ - 1}})^{K + 1}}W(z)\end{equation}
 
\begin{equation}\tag{43} = {(1 - {z^{ - 1}})^{K + 1 - S}}W(z){D_S}(z)\end{equation}
 
\begin{equation}\tag{44} = \bar R(z){D_S}(z)\end{equation}
where Display Formula\(\bar R(z) = {(1 - {z^{ - 1}})^{K + 1 - S}}W(z).\)  
The two-component model
We assume the time series x to be differentiated can be modeled as x = x1 + x2 as in the previous section. Given the noisy time series y = x + w = x1+ x2 + w, we seek an estimate Display Formula\(\hat x^{\prime} \) of the unknown differentiated time series x′. Since we assumed x2 is a signal for which conventional SG filters are well suited, its derivative can be well estimated by the SG differentiation filter Display Formula\(\bar H\)—that is, Display Formula\({x^{\prime} _2} \approx \bar H({x_2} + w).\) Consequently,  
\begin{equation}\tag{45}{x^{\prime} _2} \approx \bar H({x_2} + w)\end{equation}
 
\begin{equation}\tag{46} = \bar H(y - {x_1}).\end{equation}
Hence, we propose to estimate Display Formula\({x^{\prime} _2}\) in terms of an estimate Display Formula\({\hat x_1}\),  
\begin{equation}\tag{47}{\hat x^{\prime} _2} = \bar H(y - {\hat x_1}).\end{equation}
In turn, an estimate of the unknown time series Display Formula\(x^{\prime} \) can be expressed as  
\begin{equation}\tag{48}\hat x^{\prime} = {\hat x^{\prime} _1} + {\hat x^{\prime} _2}\end{equation}
 
\begin{equation}\tag{49} = \bar F{\hat x_1} + \bar H(y - {\hat x_1})\end{equation}
 
\begin{equation}\tag{50} = (\bar F - \bar H){\hat x_1} + \bar Hy\end{equation}
 
\begin{equation}\tag{51} = \bar G{\hat x_1} + \bar Hy\end{equation}
 
\begin{equation}\tag{52} = \bar R{D_S}{\hat x_1} + \bar Hy\end{equation}
 
\begin{equation}\tag{53} = \bar R\hat u + \bar Hy\end{equation}
where Display Formula\(\bar G\) is the high-pass filter (Equation 36, which admits the factorization in Equation 44) and Display Formula\(\hat u = {D_S}{\hat x_1}.\) Once Display Formula\(\hat u\) is obtained as the solution to optimization problem (Equation 30), the estimated time series Display Formula\(\hat x^{\prime} \) is given by Equation 53. The framework developed here generalizes SG filters for higher order differentiation as well.  
The generalized SG filter: Summary
We summarize the procedure for the proposed generalized SG filter. Let y denote the time series data to be filtered. 
  1.  
    Set the parameters, K and M of the conventional SG filter. They must be positive integers with K < 2M.
  2.  
    Set the additional parameters S and λ. Parameter S must be a positive integer with SK + 1, which determines the type of abrupt deviation in the time series. Parameter λ must be a positive real number, which is nominally set proportional to the noise level.
  3.  
    Define linear filters:
     
    •  
      Define H to be the conventional SG smoothing filter.
    •  
      Define R(z) = [1 – H(z)] / (1 – z–1)S.
    •  
      Define Display Formula\(\bar H\) to be the conventional SG differentiation filter.
    •  
      Define Display Formula\(\bar R(z) = [\bar F(z) - \bar H(z)]{\rm{/}}{(1 - {z^{ - 1}})^S}.\)
  4.  
    Solve problem (Equation 30) to obtain Display Formula\(\hat u\).
  5.  
    Set Display Formula\(\hat x = R\hat u + Hy\) (smoothing).
  6.  
    Set Display Formula\(\hat x^{\prime} = \bar R\hat u + \bar Hy\) (differentiation).
The filters H, Display Formula\(\bar H\) and Display Formula\(\bar F\) can be determined using any realization of the conventional SG filter. For the entire calculation, Matlab software is available from the authors. 
Results
In this section, we apply the proposed filtering method for the problem of peak saccadic velocity estimation, and we evaluate its accuracy. First, we simulate saccadic eye movements (horizontal angular displacement) using a recently developed parametric model (Dai et al., 2016). This saccade model has three parameter: η, c, and A. The formula for the saccade model is  
\begin{equation}\tag{54}s(t) = cf(\eta t{\rm{/}}c) - cf(\eta t{\rm{/}}c - A{\rm{/}}c)\end{equation}
where f is defined as  
\begin{equation}\tag{55}f(t) = \left\{ \matrix{ t + 0.25{{\rm{e}}^{ - 2t}},\,\,\,\,\,t \ge 0 \hfill \cr 0.25{{\rm{e}}^{2t}},\,\,\,\,\,\,\,\,\,\,\,\,\,t \le 0.\, \hfill \cr} \right.\end{equation}
It has been shown that this saccade model can well approximate physiologic saccades (Dai et al., 2016). In this model, the saccade s(t) has an amplitude of A and a peak velocity of  
\begin{equation}\tag{56}{V_p} = \eta (1 - {{\rm{e}}^{ - A{\rm{/}}c}}).\end{equation}
We note that the relation between peak velocity Vp and amplitude A in Equation 56 is the “main sequence” formula proposed by Baloh, Sills, Kumley, and Honrubia (1975).  
The use of simulated saccades is useful in this work because true saccade velocities are unknown for human subjects. However, by simulating saccades with known velocity, the error of a velocity estimation method can be quantitatively evaluated. The simulation consists of horizontal saccades made to the right and then back to the central angular position with amplitude 5°, 10°, and 15°, as shown in Figure 2. We simulate data at a sampling rate of 500 samples/s to match the recording condition in our lab. (We used the EyeLink 1000 Plus eye tracker.) While the eight most common sampling rates for eye trackers are 30, 50, 60, 120, 240, 250, 500, and 1000 Hz (Mack, Belfanti, & Schwarz, 2017), eye trackers with low sampling rates are not suitable for the study of saccade dynamics, such as saccade peak velocity. A sampling rate of 333 Hz is considered sufficient for studying position, velocity, and acceleration of human saccadic eye movements (Bahill, Kallman, & Lieberman, 1982; Inchingolo & Spanio, 1985). To this simulated time series, we add a zero-mean white Gaussian random process 𝒩(0, 0.2) to model noise. We use both the conventional and proposed generalized SG filters for smoothing and differentiation. For the conventional SG filter, we use parameters K = 2, M = 5, as recommended by Nyström and Holmqvist (2010) for eye movement analysis. For the generalized SG filter, we use parameters K = 3, M = 10, S = 4. 
Figure 2
 
Smoothing (a) and differentiation (b) of simulated saccades of 5°, 10°, and 15° amplitude. Unlike the conventional SG filter, the generalized SG filter does not underestimate the saccade peak velocity.
Figure 2
 
Smoothing (a) and differentiation (b) of simulated saccades of 5°, 10°, and 15° amplitude. Unlike the conventional SG filter, the generalized SG filter does not underestimate the saccade peak velocity.
The generalized SG differentiation filter estimates the peak saccadic velocity more accurately than the conventional SG differentiation filter (Figure 2). The conventional SG filter underestimates the peak velocity, especially for small-amplitude saccades. We also evaluate the filters in terms of the root mean square error (RMSE) between the noiseless and filtered time series (Figure 2). The RMSE values further indicate that the generalized SG filter performs smoothing and differentiation more accurately than the conventional SG filter. 
To evaluate the generalized SG filter across a range of noise levels, we simulate eye movement data as in Figure 2, but vary the noise standard deviation, σ. We set λ in Equation 30 proportional to σ, because greater λ results in stronger smoothing. Across a range of noise levels, the generalized SG filter estimates position, velocity, acceleration, and saccade peak velocity more accurately (i.e., with lower RMSE) than the conventional SG filter (Figure 3). 
Figure 3
 
Evaluation of conventional and generalized SG filter across noise level (σ). The generalized SG filter has better RMSE than the conventional SG filter.
Figure 3
 
Evaluation of conventional and generalized SG filter across noise level (σ). The generalized SG filter has better RMSE than the conventional SG filter.
To further investigate the accuracy of peak saccadic velocity estimation using the conventional and generalized SG differentiation filters, we simulate 500 saccades of 5° and 15° at a sampling rate of 500 samples/s and add zero-mean white Gaussian random process 𝒩(0, 0.2) to model noise. The distribution of estimated saccade peak velocity is shown as histograms in Figure 4 for each case. The histograms reveal that the conventional SG filter is biased; it systematically underestimates the saccade peak velocity. The histograms show that the peak saccadic velocity, as estimated by the generalized SG filter, is generally closer to the true value (indicated by the dashed line). 
Figure 4
 
Histogram of peak saccadic velocity of 5° and 15° saccades. (a–b) Estimated using conventional SG filter. (c–d) Estimated using generalized SG filter. Dashed lines indicate the true peak velocities.
Figure 4
 
Histogram of peak saccadic velocity of 5° and 15° saccades. (a–b) Estimated using conventional SG filter. (c–d) Estimated using generalized SG filter. Dashed lines indicate the true peak velocities.
Systematic underestimation of peak saccadic velocity leads to distortion of the main sequence relationship, an important diagnostic tool for clinicians (Bahill, Clark, & Stark, 1975). To investigate the relative impact of the filter, we simulate eye movement data comprising 50 saccades of various amplitudes and with random intersaccadic intervals and add a zero-mean white Gaussian random process 𝒩(0, 0.1) to model noise. Figure 5a shows the true main sequence curve and its 95% prediction interval (the true peak velocities and amplitudes are known, since the saccades are simulated). The main sequence curve follows an exponential equation (Equation 56) where Vp is the saccade peak velocity, A is the saccade amplitude, and η and c are parameters to be estimated (Baloh et al., 1975). Figure 5b and c show the main sequence curves calculated using the conventional SG filter and generalized SG filter, respectively. It can be observed that the conventional SG filter results in a distorted main sequence, and the best-fitting main sequence curve differs substantially from the true main sequence curve. This is because the peak velocities of small-amplitude saccades are underestimated. The generalized SG filter results in a main sequence curve that more accurately estimates the true one. 
Figure 5
 
Main sequence relationship between peak velocity and amplitude of simulated saccades. (a) True values. (b) Calculated using conventional SG filter. (c) Calculated using generalized SG filter. The generalized SG filter yields a more accurate main sequence than the conventional SG filter.
Figure 5
 
Main sequence relationship between peak velocity and amplitude of simulated saccades. (a) True values. (b) Calculated using conventional SG filter. (c) Calculated using generalized SG filter. The generalized SG filter yields a more accurate main sequence than the conventional SG filter.
To examine the generalized SG filter for different individuals, we simulate eye movements as in Figure 5a, but vary the main sequence parameter η or c. We vary the main sequence parameters because different individuals have different main sequence relations. The shape of the simulated saccade (Equation 54) depends on the parameters η and c. Therefore, we vary either parameter η or c to represent saccades made by 20 individuals. For each individual, we simulate eye movement data comprising 50 saccades of various amplitudes and with random intersaccadic intervals and add zero-mean white Gaussian noise 𝒩 (0, 0.1). The main sequence parameters calculated using the conventional SG filter and generalized SG filter are listed in Table 1. It can be observed that the generalized SG filter estimates the main sequence parameters more accurately than the conventional SG filter. 
Table 1
 
Main sequence parameters, η and c, of simulated saccades estimated using the conventional and generalized SG filters.
Table 1
 
Main sequence parameters, η and c, of simulated saccades estimated using the conventional and generalized SG filters.
The oversmoothing of eye movement data also leads to prolonged estimates of saccadic duration. Although this phenomenon is not as prominent as the estimation of peak saccadic velocity, it affects the empirical relationship between saccade duration and amplitude. Here, we define the onset and end of a saccade using a velocity threshold of 30°/s. As shown in Figure 6b, the main sequence relationship calculated using the conventional SG filter is different from the reference one, shown in Figure 6a. In comparison, the generalized SG filter yields a more accurate estimate of the true main sequence relationship as shown in Figure 6c
Figure 6
 
Main sequence relationship between duration and amplitude of simulated saccades. (a) True values. (b) Calculated using conventional SG filter. (c) Calculated using generalized SG filter. The generalized SG filter estimates the true relationship more accurately than the conventional SG filter.
Figure 6
 
Main sequence relationship between duration and amplitude of simulated saccades. (a) True values. (b) Calculated using conventional SG filter. (c) Calculated using generalized SG filter. The generalized SG filter estimates the true relationship more accurately than the conventional SG filter.
We also compare the conventional and generalized SG filters as applied to real eye movement data. Although no “ground truth” for the real saccade peak velocity or amplitude is available, it is nevertheless informative to observe that the two filters yield different results, demonstrating the impact of the utilized filter. Figure 7 shows part of an eye movement time series recorded at a sampling rate of 500 samples/s while the participant views an image. This time series is part of a publicly available dataset for which eye movements have been manually annotated (Larsson, Nyström, & Stridh, 2013). Many automatic saccade detection algorithms have been developed and evaluated (Andersson, Larson, Holmqvist, Stridh, & Nyström, 2017; Behrens, MacKeben, & Schröder-Preikschat, 2010; Dai et al., 2016; Engbert & Kliegl, 2003; Komogortsev, Gobert, Jayarathna, Koh, & Gowda, 2010; König & Buffalo, 2014; Nyström & Holmqvist, 2010; Otero-Millan, Castro, Macknik, & Martinez-Conde, 2014). Since our focus is on methods to accurately estimate saccade parameters rather than methods to reliably detect saccades, we use the manual annotation instead of an automatic saccade detection algorithm to identify saccades. The horizontal and vertical components of eye movement time series are processed separately. The velocity and amplitude of each eye movement is then calculated by taking the square root of the sum of squares of the horizontal and vertical components (Nyström & Holmqvist, 2010). Using the conventional and generalized SG filters to calculate saccade peak velocity and amplitude, we obtain the values illustrated in Figure 8. The best-fitting exponential curves are calculated based on the respective estimated saccade parameters. It can be observed that the two filters yield distinctly different results. In particular, the main sequence curve obtained using the generalized SG filter is substantially steeper than that obtained using the conventional SG filter. This is consistent with the results using simulated saccades: The conventional SG filter tends to underestimate the peak velocity of small-amplitude saccades. 
Figure 7
 
Eye movement of one participant viewing an image. Saccades identified by experts are indicated by shaded gray areas. The annotated data is publicly available online.
Figure 7
 
Eye movement of one participant viewing an image. Saccades identified by experts are indicated by shaded gray areas. The annotated data is publicly available online.
Figure 8
 
Main sequence relationship between peak velocity and amplitude of real saccades. (a) Calculated using conventional SG filter. (b) Calculated using generalized SG filter.
Figure 8
 
Main sequence relationship between peak velocity and amplitude of real saccades. (a) Calculated using conventional SG filter. (b) Calculated using generalized SG filter.
We further investigate the conventional and generalized SG filter as applied to real eye movement data of 20 individuals. The data was recorded using the EyeLink 1000 Plus eye tracker and sampled at a rate of 500 samples/s while the participants read numbers on a test card from left to right and top to bottom (Rizzo et al., 2016). Saccades are manually annotated by author WD for this study. The main sequence parameters calculated using the conventional SG filter and generalized SG filter are listed in Table 2. The parameter c estimated by the conventional SG filter is larger than that estimated by the generalized SG filter (8.10 ± 2.07 vs. 5.76 ± 2.17, p = 0.001, two-sample t test). The parameter η estimated using either the conventional or generalized SG filter is similar. The same observation is found when the filter is applied to different data: The conventional SG filter tends to underestimate the peak velocity of small-amplitude saccades. Therefore, we conclude that the main sequence estimated using the generalized SG filter is more accurate than the one estimated using the conventional SG filter. 
Table 2
 
Main sequence parameters, η and c, of real saccades estimated using the conventional and generalized SG filters.
Table 2
 
Main sequence parameters, η and c, of real saccades estimated using the conventional and generalized SG filters.
Discussion
The SG filter is a particular (linear) low-pass filter widely used for smoothing and differentiation of time series. Other commonly used filters for eye movement data is the moving average filter and the Butterworth filter (Mack et al., 2017). They are used for data smoothing and usually followed by a two-point central difference filter for data differentiation. All three filters are linear filters and hence obey an unavoidable trade-off between noise suppression and signal distortion. It should also be recognized that linear low-pass filters tend to oversmooth abrupt deviations from a smooth trend. Therefore, quantitative measures of signal extremes (e.g., the peak value of a signal derivative) might be underestimated when calculated using the SG filter, or any other linear filter. The proposed nonlinear generalization of the SG filter aims to mitigate this issue, by modeling a signal component having discontinuities in its order-S derivative of some prescribed order S
Historically, the relationship between saccade peak velocity and amplitude has been a useful metric in clinical practice for the detection and assessment of various disease states; however, saccades of different amplitudes have different power spectra (Harris, Wallman, & Scudder, 1990). Hence, the use of a conventional low-pass filter with a fixed cut-off frequency may be suboptimal for the analysis of eye movement data comprising saccades of various sizes. One comprehensive study used a saccade model to evaluate the three commonly used linear filters (moving average, Savitzky-Golay, and Butterworth) as applied to eye movement data (Mack et al., 2017). Their results suggest various filters should be used for saccades of different sizes. However, eye movement data usually contains saccades of various sizes. For example, we make small intraline saccades and large interline saccades during reading. Recent studies showed that chronic concussion patients and multiple sclerosis patients made more saccades to complete a rapid number naming test, and these extra saccades appear to be small (Hainline et al., 2017; Rizzo et al., 2016). As shown in the Results section, the conventional SG filter tends to underestimate the peak velocity of saccades, especially small-amplitude saccades. The proposed nonlinear generalization of the SG filter may more accurately estimate saccade peak velocity and other saccade parameters. 
In the Results section, we used white Gaussian random processes to model noise. Studies of artificial eyes for the measurement of noise in eye trackers showed that the noise is white (it has the same energy at all frequencies) when artificial eyes are used (Coey, Wallot, Richardson, & Van Orden, 2012; Wang, Mulvey, Pelz, & Holmqvist, 2017). Both studies point out that human data produce pink noise (the energy decreases as the frequency increases). This may be because the human eye is never still and always makes micro movements, which can be similar in amplitude to the noise. We conducted numerical experiments adding pink noise instead of white noise in the simulations. The generalized SG filter was found to more accurately estimate peak velocity in pink noise as well as white noise. 
Conclusions
In the quantitative analysis of saccades, the limitations of linear low-pass filters should be taken into account, given their potential to negatively impact objective metrics. This article demonstrates through simulation using a recently developed parametric saccade model that the SG filter, a widely used linear low-pass filter, may systematically underestimate the peak velocity of saccades. This article proposes a nonlinear generalization of the SG filter that can more accurately track deviations from the smooth background of an underlying signal. The generalized SG filter was shown to estimate saccade peak velocity more accurately than the conventional SG filter, which is commonly used for this purpose. In turn, the generalized SG filter leads to a more accurate estimation of the main sequence (i.e., the relationship between the saccade peak velocity and amplitude). Given the clinical significance of quantitative saccade analysis in disease diagnosis, the generalized SG filters may be useful in the study of eye movement data, particularly saccades. 
Acknowledgments
This work was supported in part by the National Science Foundation under Grant No. 1525398. 
Commercial relationships: none. 
Corresponding author: Weiwei Dai. 
Address: Department of Electrical and Computer Engineering, Tandon School of Engineering, New York University, Brooklyn, NY, USA. 
References
Al-Alaoui, M. A. (2007). Linear phase low-pass IIR digital differentiators. IEEE Transactions on Signal Processing, 55 (2), 697–706.
Andersson, R., Larsson, L., Holmqvist, K., Stridh, M., & Nyström, M. (2017). One algorithm to rule them all? An evaluation and discussion of ten eye movement event-detection algorithms. Behavior Research Methods, 49 (2), 616–637.
Azami, H., Mohammadi, K., & Bozorgtabar, B. (2012). An improved signal segmentation using moving average and Savitzky-Golay filter. Journal of Signal and Information Processing, 3 (1), 39–44.
Bahill, A. T., Clark, M. R., & Stark, L. (1975). The main sequence, a tool for studying human eye movements. Mathematical Biosciences, 24 (3), 191–204.
Bahill, A. T., Kallman, J. S., & Lieberman, J. E. (1982). Frequency limitations of the two-point central difference differentiation algorithm. Biological Cybernetics, 45 (1), 1–4.
Baloh, R. W., Sills, A. W., Kumley, W. E., & Honrubia, V. (1975). Quantitative measurement of saccade amplitude, duration, and velocity. Neurology, 25 (11), 1065–1070.
Behrens, F., MacKeben, M., & Schröder-Preikschat, W. (2010). An improved algorithm for automatic detection of saccades in eye movement data and for calculating saccade parameters. Behavior Research Methods, 42 (3), 701–708.
Chen, A. L., Riley, D. E., King, S. A., Joshi, A. C., Serra, A., Liao, K.,… Leigh, R. J. (2010). The disturbance of gaze in progressive supranuclear palsy: Implications for pathogenesis. Frontiers in Neurology, 1, 147.
Coey, C. A., Wallot, S., Richardson, M. J., & Van Orden, G. (2012). On the structure of measurement noise in eye-tracking. Journal of Eye Movement Research, 5 (4), 1–10.
Combettes, P. L., & Pesquet, J.-C. (2011). Proximal splitting methods in signal processing. In Bauschke, H. H. Burachik, R. S. Combettes, P. L. Elser, V. Luke, D. R. & Holkowicz H. (Eds.), Fixed-point algorithms for inverse problems in science and engineering (pp. 185–212). New York: Springer-Verlag.
Dai, W., Selesnick, I., Rizzo, J.-R., Rucker, J. C., & Hudson, T. E. (2016, December). A parameteric model for saccadic eye movement. In IEEE Signal Processing in Medicine and Biology Symposium (pp. 1–6), doi:10.1109/SPMB.2016.7846860. New York: IEEE.
Di Stasi, L. L., Renner, R., Catena, A., Cañas, J. J., Velichkovsky, B. M., & Pannasch, S. (2012). Towards a driver fatigue test based on the saccadic main sequence: A partial validation by subjective report data. Transportation Research Part C: Emerging Technologies, 21 (1), 122–133.
Elad, M. (2010). Sparse and redundant representations: From theory to applications in signal and image processing. New York: Springer.
Engbert, R., & Kliegl, R. (2003). Microsaccades uncover the orientation of covert attention. Vision Research, 43 (9), 1035–1045.
Ethier, V., Zee, D. S., & Shadmehr, R. (2008). Changes in control of saccades during gain adaptation. Journal of Neuroscience, 28 (51), 13929–13937.
Garbutt, S., Harwood, M. R., Kumar, A. N., Han, Y. H., & Leigh, R. (2003). Evaluating small eye movements in patients with saccadic palsies. Annales of the New York Academy of Sciences, 1004 (1), 337–346.
Geissler, A., Gartus, A., Foki, T., Tahamtan, A. R., Beisteiner, R., & Barth, M. (2007). Contrast-to-noise ratio (CNR) as a quality parameter in fMRI. Journal of Magnetic Resonance Imaging, 25 (6), 1263–1270.
Hainline, C., Rizzo, J.-R., Hudson, T. E., Dai, W., Birkemeier, J., Raynowska, J.,… Rucker, J. C. (2017). Capturing saccades in multiple sclerosis with a digitized test of rapid number naming. Journal of Neurology, 264 (5), 989–998.
Harris, C. M., Wallman, J., & Scudder, C. A. (1990). Fourier analysis of saccades in monkeys and humans. Journal of Neurophysiology, 63 (4), 877–886.
Horn, A. K. E., & Büttner-Ennever, J. A. (1998). Premotor neurons for vertical eye movements in the rostral mesencephalon of monkey and human: Histologic identification by parvalbumin immunostaining. Journal of Comparative Neurology, 392 (4), 413–427.
Horn, A. K. E., Büttner-Ennever, J. A., Suzuki, Y., & Henn, V. (1995). Histological identification of premotor neurons for horizontal saccades in monkey and man by parvalbumin immunostaining. Journal of Comparative Neurology, 359 (2), 350–363.
Inchingolo, P., & Spanio, M. (1985). On the identification and analysis of saccadic eye movements—A quantitative study of the processing procedures. IEEE Transactions on Biomedical Engineering, 9, 683–695.
Komogortsev, O. V., Gobert, D. V., Jayarathna, S., Koh, D. H., & Gowda, S. M. (2010). Standardization of automated analyses of oculomotor fixation and saccadic behaviors. IEEE Transactions on Biomedical Engineering, 57 (11), 2635–2645.
König, S. D., & Buffalo, E. A. (2014). A nonparametric method for detecting fixations and saccades using cluster analysis: Removing the need for arbitrary thresholds. Journal of Neuroscience Methods, 227, 121–131.
Kumar, B., Roy, S. C. D., & Shah, H. (1992). On the design of FIR digital differentiators which are maximally linear at the frequency π/p,p ∈ {positive integers}. IEEE Transactions on Signal Processing, 40 (9), 2334–2338.
Larsson, L., Nyström, M., & Stridh, M. (2013). Detection of saccades and postsaccadic oscillations in the presence of smooth pursuit. IEEE Transactions on Biomedical Engineering, 60 (9), 2484–2493.
Leigh, R. J., & Riley, D. E. (2000, March). Eye movements in parkinsonism: It's saccadic speed that counts. Neurology, 54 (5), 1020–1021.
Leigh, R. J., & Zee, D. S. (2015). The neurology of eye movements. Cambridge, UK: Oxford University Press.
Luo, J.-W., Bai, J., He, P., & Ying, K. (2004). Axial strain calculation using a low-pass digital differentiator in ultrasound elastography. IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, 51 (9), 1119–1127.
Luo, J.-W., Ying, K., He, P., & Bai, J. (2005). Properties of Savitzky-Golay digital differentiators. Digital Signal Processing, 15 (2), 122–136.
Mack, D. J., Belfanti, S., & Schwarz, U. (2017). The effect of sampling rate and lowpass filters on saccades—A modeling approach. Behavior Research Methods, E-pub ahead of print, doi:10.3758/s13428-016-0848-4.
Nyström, M., & Holmqvist, K. (2010). An adaptive algorithm for fixation, saccade, and glissade detection in eyetracking data. Behavior Research Methods, 42 (1), 188–204.
Orfanidis, S. J. (1995). Introduction to signal processing. Upper Saddle River, NJ: Prentice-Hall.
Otero-Millan, J., Castro, J. L. A., Macknik, S. L., & Martinez-Conde, S. (2014). Unsupervised clustering method to detect microsaccades. Journal of Vision, 14 (2): 18, 1–17, doi:10.1167/14.2.18. [PubMed] [Article]
Pei, S.-C., & Wang, P.-H. (2001). Closed-form design of maximally flat FIR Hilbert transformers, differentiators, and fractional delayers by power series expansion. IEEE Transactions on Circuits and Systems I, 48 (4), 389–398.
Quaia, C., Joiner, W. M., FitzGibbon, E. J., Optican, L. M., & Smith, M. A. (2010). Eye movement sequence generation in humans: Motor or goal updating? Journal of Vision, 10 (14): 28, 1–31, doi:10.1167/10.14.28. [PubMed] [Article]
Rabiner, L. R., & Schafer, R. W. (1974). On the behavior of minimax relative error FIR digital differentiators. Bell Labs Technical Journal, 53 (2), 333–361.
Ramat, S., Leigh, R. J., Zee, D. S., & Optican, L. M. (2007). What clinical disorders tell us about the neural control of saccadic eye movements. Brain, 130 (1), 10–35.
Ramat, S., Leigh, R. J., Zee, D. S., Shaikh, A. G., & Optican, L. M. (2008). Applying saccade models to account for oscillations. Progress in Brain Research, 171, 123–130.
Rizzo, J.-R., Hudson, T. E., Dai, W., Birkemeier, J., Pasculli, R. M., Selesnick, I.,… Rucker, J. C. (2016). Rapid number naming in chronic concussion: Eye movements in the King-Devick test. Annals of Clinical and Translational Neurology, 3 (10), 801–811.
Rudin, L. I., Osher, S., & Fatemi, E. (1992). Nonlinear total variation based noise removal algorithms. Physica D: Nonlinear Phenomena, 60 (1–4), 259–268.
Savitzky, A., & Golay, M. J. E. (1964). Smoothing and differentiation of data by simplified least squares procedures. Analytical Chemistry, 36 (8), 1627–1639.
Schafer, R. W. (2011a). On the frequency-domain properties of Savitzky-Golay filters. In Digital Signal Processing Workshop and IEEE Signal Processing Education Workshop (pp. 54–59). New York: IEEE.
Schafer, R. W. (2011b). What is a Savitzky-Golay filter? IEEE Signal Processing Magazine, 28 (4), 111–117.
Schneider, R. C., & Kovar, K.-A. (2003). Analysis of ecstasy tablets: Comparison of reflectance and transmittance near infrared spectroscopy. Forensic Science International, 134 (2), 187–195.
Schütz, A. C., Braun, D. I., & Gegenfurtner, K. R. (2011). Eye movements and perception: A selective review. Journal of Vision, 11 (5): 9, 1–30, doi:10.1167/11.5.9. [PubMed] [Article]
Selesnick, I. (2002). Maximally flat low-pass digital differentiators. IEEE Transactions on Circuits and Systems II, 49 (3), 219–223.
Selesnick, I. (2015). Sparsity-assisted signal smoothing. In Excursions in harmonic analysis (Vol. 4, pp. 149–176). Basel, Switzerland: Birkhäuser (Springer).
Selesnick, I. (2017, March). Sparsity-assisted signal smoothing (revisited). In IEEE International Conference on Acoustics, Speech, Signal Processing (pp. 4546-4550), doi:10.1109/ICASSP.2017.7953017. New York: IEEE.
Selesnick, I., Graber, H. L., Pfeil, D. S., & Barbour, R. L. (2014). Simultaneous low-pass filtering and total variation denoising. IEEE Transactions on Signal Processing, 62 (5), 1109–1124.
Shaikh, A. G., Xu-Wilson, M., Grill, S., & Zee, D. S. (2010). ‘Staircase' square-wave jerks in early Parkinson's disease. British Journal of Ophthalmology, 95 (5), 705–709, doi:10.1136/bjo.2010.179630.
Shajeesh, K. U., Kumar, S., Pravena, D., & Soman, K. P. (2012). Speech enhancement based on Savitzky-Golay smoothing filter. International Journal of Computer Applications, 57 (21), 39–44.
Wadia, N. H., & Swami, R. K. (1971). A new form of heredo-familial spinocerebellar degeneration with slow eye movements (nine families). Brain, 94 (2), 359–374.
Wang, D., Mulvey, F. B., Pelz, J. B., & Holmqvist, K. (2017). A study of artificial eyes for the measurement of precision in eye-trackers. Behavior Research Methods, 49 (3), 947–959.
Figure 1
 
The frequency response of Savitzky-Golay (a) smoothing and (b) differentiation filters for several parameter values.
Figure 1
 
The frequency response of Savitzky-Golay (a) smoothing and (b) differentiation filters for several parameter values.
Figure 2
 
Smoothing (a) and differentiation (b) of simulated saccades of 5°, 10°, and 15° amplitude. Unlike the conventional SG filter, the generalized SG filter does not underestimate the saccade peak velocity.
Figure 2
 
Smoothing (a) and differentiation (b) of simulated saccades of 5°, 10°, and 15° amplitude. Unlike the conventional SG filter, the generalized SG filter does not underestimate the saccade peak velocity.
Figure 3
 
Evaluation of conventional and generalized SG filter across noise level (σ). The generalized SG filter has better RMSE than the conventional SG filter.
Figure 3
 
Evaluation of conventional and generalized SG filter across noise level (σ). The generalized SG filter has better RMSE than the conventional SG filter.
Figure 4
 
Histogram of peak saccadic velocity of 5° and 15° saccades. (a–b) Estimated using conventional SG filter. (c–d) Estimated using generalized SG filter. Dashed lines indicate the true peak velocities.
Figure 4
 
Histogram of peak saccadic velocity of 5° and 15° saccades. (a–b) Estimated using conventional SG filter. (c–d) Estimated using generalized SG filter. Dashed lines indicate the true peak velocities.
Figure 5
 
Main sequence relationship between peak velocity and amplitude of simulated saccades. (a) True values. (b) Calculated using conventional SG filter. (c) Calculated using generalized SG filter. The generalized SG filter yields a more accurate main sequence than the conventional SG filter.
Figure 5
 
Main sequence relationship between peak velocity and amplitude of simulated saccades. (a) True values. (b) Calculated using conventional SG filter. (c) Calculated using generalized SG filter. The generalized SG filter yields a more accurate main sequence than the conventional SG filter.
Figure 6
 
Main sequence relationship between duration and amplitude of simulated saccades. (a) True values. (b) Calculated using conventional SG filter. (c) Calculated using generalized SG filter. The generalized SG filter estimates the true relationship more accurately than the conventional SG filter.
Figure 6
 
Main sequence relationship between duration and amplitude of simulated saccades. (a) True values. (b) Calculated using conventional SG filter. (c) Calculated using generalized SG filter. The generalized SG filter estimates the true relationship more accurately than the conventional SG filter.
Figure 7
 
Eye movement of one participant viewing an image. Saccades identified by experts are indicated by shaded gray areas. The annotated data is publicly available online.
Figure 7
 
Eye movement of one participant viewing an image. Saccades identified by experts are indicated by shaded gray areas. The annotated data is publicly available online.
Figure 8
 
Main sequence relationship between peak velocity and amplitude of real saccades. (a) Calculated using conventional SG filter. (b) Calculated using generalized SG filter.
Figure 8
 
Main sequence relationship between peak velocity and amplitude of real saccades. (a) Calculated using conventional SG filter. (b) Calculated using generalized SG filter.
Table 1
 
Main sequence parameters, η and c, of simulated saccades estimated using the conventional and generalized SG filters.
Table 1
 
Main sequence parameters, η and c, of simulated saccades estimated using the conventional and generalized SG filters.
Table 2
 
Main sequence parameters, η and c, of real saccades estimated using the conventional and generalized SG filters.
Table 2
 
Main sequence parameters, η and c, of real saccades estimated using the conventional and generalized SG filters.
×
×

This PDF is available to Subscribers Only

Sign in or purchase a subscription to access this content. ×

You must be signed into an individual account to use this feature.

×