Introduction

Eye tracking is used in many fields of research and the examination of oculomotor behavior yields an abundance of interesting parameters (Yarbus, 1967). The analysis of fixation locations, for example, can reveal how we perceive art (Quiroga & Pedreira, 2011; Yarbus, 1967), faces (Dahl, Wallraven, Bülthoff, & Logothetis, 2009), or sexual characteristics (Dixson, Grimshaw, Linklater, & Dixson, 2011). Nevertheless, fixations are only a very small part of the set of oculomotor behaviors. Due to their close connection to attention (Shepherd, Findlay, & Hockey, 1986), the fast eye movements occurring between two phases of fixation, known as saccades, are another interesting object of research. These eye movements have the unique property of a fixed relationship between amplitude and peak velocity, known as the “main sequence” (Bahill, Clark, & Stark, 1975; Baloh, Sills, Kumley, & Honrubia, 1975). Thus, unlike other goal-directed movements (e.g., grasping), saccades of similar amplitudes always have similar peak velocities. This property arises from the saccade generators located in the brain stem (Robinson, 1975; Sparks, 2002). Although quite a robust trait, deviations from the normal main sequence exist. Voluntarily executed saccades reach lower peak velocities than reflexively-driven saccades (Smit, van Gisbergen, & Cools, 1987). The same holds true for saccades toward dimly lit targets (Becker & Fuchs, 1969). In addition, saccadic peak velocities differ among participants. Video game players, for instance, reach higher peak velocities than non-players (Mack & Ilg, 2014). In a clinical context, Yee et al. (1985, p. 938) noted that due to the

“maximal recruitment of motoneurons and extraocular muscle fibers […] during large saccades […] it would be expected that most disorders that impair the recruitment of motoneurons and muscle fibers or that otherwise interfere with the contraction of muscle fibers will decrease the peak velocity.”

Indeed, several pathologies are associated with alterations in saccadic peak velocity (Anderson & MacAskill, 2013), including progressive supranuclear palsy (Newman, Gay, Stroud, & Brooks 1970), Huntington’s disease (Starr, 1967), and spinocerebellar degeneration (Zee, Optican, Cook, Robinson, & Engel, 1976). Therefore, it is crucial that amplitude and peak velocity of saccades can be accurately measured and restored from the recorded data. Unfortunately, peak velocity is especially sensitive to the recording method (Bahill, Clark, et al., 1975; Boghen et al. 1974; Hess, Müri, & Meienberg, 1986; Juhola, Jäntti, & Pyykko, 1985; Kimmel, Mammo, & Newsome, 2012; Yee et al., 1985) and successive data processing (Inchingolo & Spanio, 1985; Jäntti et al., 1984; Juhola, 1986; but see also Schmitt, Muser, Lanz, Walz, & Schwarz, 2007).

Today, most systems record and store only the eye position, which creates the need to compute the velocity profiles by digital differentiation. It has been shown that a one-sample, two-point central difference (2pCD) method is ideal for this purpose (Bahill, Kallman, & Lieberman, 1982). Importantly, this method produces no lag in the velocity channel with respect to the position channel. In addition, it implicitly applies a lowpass filter that reduces the high frequency noise introduced by the differentiation process (Bahill et al., 1982; Inchingolo & Spanio, 1985). However, the computed peak velocity is also affected by the sampling rate of the eye tracker itself. Only sampling rates above 200–300 Hz allow the undistorted recovery of saccadic peak velocities (Inchingolo & Spanio, 1985; Juhola et al., 1985; Schmitt et al., 2007), which would imply using the highest possible sampling rate for saccade measurements. Besides differentiation, though, additional noise sources are present during recording and higher sampling rates pick up more high frequency noise. In general, two noise sources are the main problem: Instrumentation noise adds unstructured variations to the data and is assumed to be independent and of Gaussian nature (“additive Gaussian white noise”; AWGN). Using artificial eyes, these assumptions have been shown to be correct for several eye trackers (Coey, Wallot, Richardson, & Van Orden, 2012; Wang, Mulvey, Pelz, & Holmqvist, 2016). Most systems also pick up at least a small portion of mains noise introduced through an alternating current power source. This noise has a high harmonic content with a fundamental frequency depending on the local power line frequency (e.g., 50 or 60 Hz). Both noise sources cannot be completely suppressed and thus require digital lowpass filtering. In contrast to other biological signals like ECG, we do not have any prior knowledge about the exact spatiotemporal layout of eye movements. This has been framed well by Pettersson et al. (2013, p. 2):

“Denoising […] is challenging since eye movements are usually non-repetitive making the […] signal unpredictable. Consequently, methods that need structural or temporal knowledge about the expected signal cannot be used.”

Previous studies have analyzed the effects of different filters and their respective settings, but these studies have some limitations. Although Inchingolo and Spanio (1985) performed a thorough examination of the processing steps required during the analysis of eye movements, they used a rather small number of sampling rates and did not analyze the effect of mains frequency. Other studies (Bahill, Brockenbrough, & Troost, 1981; Jäntti et al., 1984) applied only a limited set of filters, or chose a focus on computational efficiency (Juhola, 1986), which is not a major concern in most (offline processing) setups nowadays. A thorough analysis of filter effects also suffers from the unalterable noise composition in the measured data used in many of these studies. In addition, such data create the need for saccade detection, which introduces further variation in the computed amplitudes and durations depending on the detection method. For example, an erroneous choice of the velocity threshold in a simple velocity-based procedure (a so-called “IVT” method; Salvucci & Goldberg, 2000) significantly influences the durations of saccades (Inchingolo & Spanio, 1985). Similarly, a flawed position criterion may include post-saccadic oscillations or overshoots in the saccade, thereby adding variability to the measured amplitude and duration (Bahill et al., 1981). An alternative to automatic and, in this context, error-prone saccade detection is the creation of a ground truth by manually labeling the raw data. This method requires skilled observers and is very time consuming. Hessels, Niehorster, Kemner, and Hooge (2016) for example, needed 3 h per expert for labeling 40 min of eye movement data. Hitherto, there has been only one attempt to do this on a larger scale (The COGAIN Association, 2016). However, measured data are still affected by a fixed noise composition and, since experts do not always produce identical labeling, it is not clear if “expert coders serve as a good gold standard for event detection” (Hessels et al., 2016, p. 16) at all.

To overcome these limitations, one can create a model of an “ideal” saccadic waveform. With such “an ideal saccade […] there is no need to consider noise, secondary saccades, or other normally present factors” (Juhola et al., 1985, p. 68). This ideal saccade can act as a ground truth for the effects of subsequent signal processing steps and omits the need for saccade detection since all important parameters are determined by the model. In addition, ideal saccades allow the overlay of different mains noise frequencies, noise levels, sampling rates, filter types, and filter settings. This can be used to reveal the best filter for a given condition – a procedure which is not possible with empirically measured saccades. Early work on the effects of filtering and sampling rates on peak velocities used synthesized saccades from artificial acceleration profiles (Jäntti et al., 1984; Juhola, 1986; Juhola et al., 1985). This method is very flexible and produces saccadic waveforms that closely resemble real saccades. Nevertheless, it is not analytical and does not directly incorporate measured main sequences. In addition, these studies used only a small number of amplitudes for their modelled saccades.

We wanted to fill this gap with a simple analytical model that is capable of creating saccades of arbitrary amplitudes and also incorporates a main sequence. Three important constraints were made. First, since most eye trackers nowadays are based on videooculography, we omitted the option and effect of analog filtering and assumed the raw data from the eye tracker is fairly unfiltered. Second, to simplify signal processing, we used only one filtering step on the raw position data and a subsequent application of the mentioned 2pCD-method for the computation of the velocity profile. Third, we only considered offline processing and thus applied non-causal filters.

Methods

Through the methods sections, ranges of numbers will be given by “start:step:end”, where step is the increment from start to end.

The saccade model

We used the following hyperbolic tangent model to simulate saccadic waveforms:

$$ p(t)\ :={\displaystyle {\int}_{-\infty}^tv\left(\tau \right)\ \mathrm{d}\tau = \frac{a_{\infty }(a)}{2}}\cdot \left( \tanh \left(g(t)\right)+1\right) $$
(1)
$$ v(t)\ :={v}_{max}^{*}(a)\cdot {\mathrm{sech}}^2\left(g(t)\right) $$
(2)
$$ g(t)\ : = 2f(a)\cdot \left(t-\phi \right) $$
(3)

with the position profile p(t), the velocity profile v(t), the time scaling function g(t), the saccade amplitude a ≠ 0, and the time of occurrence ϕ of the amplitude-dependent peak velocity v * max (a). Given the amplitude, the peak velocity can be computed from any amplitude-velocity main sequence. Solving these equations with the constraint of a fixed velocity threshold v 0 > 0 at saccade on- and offset, yields the following shorthands:

$$ \lambda (a)=\sqrt{1-\frac{v_0}{\left|{v}_{max}^{*}(a)\right|}} $$
(4)
$$ {a}_{\infty }(a)=\frac{a}{\lambda (a)} $$
(5)
$$ f(a)=\frac{v_{max}^{*}(a)}{a_{\infty }(a)} $$
(6)

where f(a) is the amplitude-dependent time scale (“frequency”) of the saccade, a ∞ is the asymptotic amplitude for t → ∞, and λ(a) is a scaling function incorporating the velocity threshold and the main sequence. The model should yield real-valued non-singular profiles, and therefore we require λ(a) ∈ ℝ+ ⇔ |v * max (a)| > v 0. However, the peak velocity approaches the velocity threshold for small amplitudes, e.g., \( \left|{v}_{max}^{*}(a)\right|\underset{a\to 0}{\to }{v}_0 \) and, thus \( {a}_{\infty }(a)\underset{a\to 0}{\to}\infty \). Therefore, the local minimum

$$ {a}_{min}^{v_0}:=\underset{a}{\mathrm{argmin}}\kern.2em {a}_{\infty }(a),\kern1.5em a>0 $$
(7)

of a ∞(a) gives the lower bound on the amplitudes.

A detailed derivation of the model is given in the Supplementary Methods, Section 1. An implementation in Matlab (R2015b, MathWorks, Natick, MA, USA) can be found at Mack, Belfanti, and Schwarz (2016) or requested from the first author.

Since the shape of the hyperbolic functions is fully determined by specifying a and v 0, the duration d(a) of the simulated saccade cannot be specified by a measured main sequence but has to be computed with the following equation:

$$ d(a)=l(a)\cdot \frac{a}{v_{max}^{\ast }(a)},\kern1.5em l(a):=\frac{\mathrm{arctanh}\left(\lambda (a)\right)}{\lambda (a)} $$
(8)

This is not a linear function as is expected for measured amplitude-duration main sequences, due to the lambda-term l(a). The non-linear behavior strongly depends on the selection of the amplitude-velocity main sequence. Therefore, care should be taken when choosing the according function. It is, for example, very helpful if the analytical limit for λ(a) exists (see Realization of the saccade model, below), since then a bound on the error can be computed.

Creation of saccade profiles

Realization of the saccade model

We modelled saccades with 155 different amplitudes (0.5:0.1:15° and 20:5:60°). The profiles were created at a sampling rate of 20 kHz to avoid any aliasing artifacts in the ideal data. Each profile was simulated over the course of 1 s (“trial duration”) to allow all valid filters to be applied at each sampling rate. We chose a velocity threshold for the model of \( {v}_0=1\ \raisebox{1ex}{${}^{\circ}$}\!\left/ \!\raisebox{-1ex}{$s$}\right. \) and the following main sequence functions from Inchingolo and Spanio (1985; data from Fig. 7):

$$ {v}_{max,\ I85}^{*}(a):=\frac{1000\cdot a}{b_o+{b}_1\cdot \left|a\right|},\kern1.25em {b}_0=18\ ms,\kern0.5em {b}_1=1.2\ \raisebox{1ex}{$ ms$}\!\left/ \!\raisebox{-1ex}{${}^{\circ}$}\right. $$
(9)
$$ {d}_{I85}^{\ast }(a)\kern.2em :=\kern.1em 38\kern.1em ms+2\kern.2em \raisebox{1ex}{$ ms$}\!\left/ \!\raisebox{-1ex}{${}^{\circ}$}\right.\kern.1em \operatorname{}\cdot \left|a\right| $$
(10)

Using equation (9) to compute the durations from the model as given in equation (8), yields

$$ {d}_{I85}(a)={l}_{I85}(a)\cdot \left({b}_o+{b}_1\cdot \left|a\right|\right) $$
(11)

A careful choice of the amplitude-velocity main sequence is important, as mentioned above. This can be seen by the thoughtful formulation of Inchingolo and Spanio, for which the analytical limit of λ(a) exists:

$$ {\lambda}_{\infty, I85}:=\underset{a\to \infty }{ \lim\ }{\lambda}_{I85}(a)=\sqrt{1-\frac{b_1}{1000\cdot {v}_0}} $$
(12)

The limit can be used to compute the lambda-term l I85(a) in equation (11), resulting in a “linearized” duration function:

$$ {\widehat{d}}_{I85}(a):=\frac{\mathrm{arctanh}\left({\lambda}_{\infty, I85}\right)}{\lambda_{\infty, I85}}\cdot \left({b}_o+{b}_1\cdot \left|a\right|\right) $$
(13)

The linearized duration now allows the computation of a bound on the relative error. For amplitudes ≥9°, the relative error between the model durations and the linearized version is below 5%. Hence, the model durations can be considered linear from this point.

Finally, from the constraint on a ∞(a) in equation (7), the smallest possible saccade amplitude for this model is:

$$ {a}_{min,I85}^{v_0}=\frac{3}{2}\cdot \frac{b_0\cdot {v}_0}{1000-{b}_1\cdot {v}_0}\approx 0.03{}^{\circ} $$
(14)

Figure 1 shows the difference between the main sequence functions and the model results.

Fig. 1
figure 1

Main sequence of modelled and measured saccades. Values are shown for the ideal data created at 20 kHz. Markers depict the parameters returned by the model, thin black lines show the main sequence functions from Inchingolo and Spanio (1985). Lower panels show the according gains (model/measured). Triangles give the mean of the gain over all amplitudes

Full size image

Choice of sampling rates

For the creation of a set of broadly applicable filter specifications, we determined the most commonly used sampling rates in oculomotor research. We gathered information on 94 commercially available eye trackers from 27 manufactures, inspired by the list of The COGAIN Association (2015). Eye trackers that are able to run at several sampling rates were counted for each rate. The eight most common sampling rates were 30, 50, 60, 120, 240, 250, 500, and 1 kHz. The complete list of eye trackers is given in Supplementary Table 1 and more information on the distribution in Supplementary Methods, Section 2.

Although our model permits the computation of saccadic position and velocity profiles at arbitrary sampling rates, we computed the position profiles for the different sampling rates by downsampling the ideal position data. This resembles natural conditions under which an ideal signal is sampled by the eye tracker. Care was taken that the time of occurrence of the peak velocity was not in between two samples. This is important since the downsampling procedure would otherwise have lowered the peak velocity by clipping the maximum through erroneous sample spacing.

To show the data in laboratory settings, we applied a simple method for saccade identification by a velocity threshold (IVT) of 1°/s to the filtered data for exemplary 10°-saccades. The IVT determined on- and offset of the saccade. From these values, the amplitude, duration and peak velocity (as maximum between on- and offset) of the saccade were computed. The resulting parameters serve as comparison to the simulated reference values.

The noise model

Noise was modelled as a combination of AWGN, reflecting instrumentation noise, and a sine wave with a uniformly random phase shift for mains noise. The AWGN component was created at 9 different amplitudes (0°/no noise, 0.01°, 0.04°, 0.07°, 0.1°, 0.3°, 0.5°, 1°, and 2°), as given by the standard deviation. Mains noise was added with amplitudes proportionally to the AWGN amplitude at three different contributions (0%/no mains noise, 5%, and 30%) and two power line frequencies (50 Hz and 60 Hz). The final noise profile was added to the ideal position profile and, accordingly, submitted to the downsampling procedure to mimic the implicit lowpass filter effect of the lower sampling rates. The noisy downsampled position profiles were then differentiated with the 2pCD-method (see Introduction) to obtain the velocity profiles and to reflect the need for digital differentiation under laboratory settings. Since white noise was added to the position profile, the according derivation for the velocity profiles created a purple noise component, e.g., a noise component which increases its power with frequency. Figure 2 shows one realization of the noise profile for the ideal data (Fig. 2a, upper panel) and the according derivation (Fig. 2a, lower panel). The derivation clearly causes the white noise to become purple (Fig. 2b, lower panel). In addition, the lowpass filter characteristic of the 2pCD at half the bandwidth of the original signal leads to a decline in amplitude above 5 kHz, instead of the expected further increase. It can also be seen that the mains noise component in the velocity profile is far less significant than in the position profile (Fig. 2b, upper panel).

Fig. 2
figure 2

(a) Exemplary noise realization and (b) according Fourier spectra. Ideal data created at 20 kHz is shown. AWGN amplitude was 2° with 30% (0.6°) mains noise at 60 Hz mains frequency (red lines). White lines in (b) show linear regressions for the lower and upper half of the frequency spectrum. The gray dotted line shows the expected increase in power for true purple noise which is dampened by the lowpass characteristic of the 2pCD-method. Fourier spectra are computed over the whole trial duration zero-padded to a frequency resolution of 10 mHz

Full size image

In total, ten independent AWGN components (“repetitions”) were computed for each noise amplitude. Each mains noise profile was analyzed with each of the AWGN components to obtain a solid estimate of the filter effects on peak velocity.

Filter settings

We chose two finite-impulse response (FIR) filters and one infinite-impulse response (IIR) filter. The first FIR filter was a simple moving average filter (MA), which stretches the signal in time, i.e. the transition between different signal levels is limited to a finite rise time. This rise time is directly proportional to the window size as illustrated in Fig. 3a. The larger the window gets, the more pronounced the reduction of the peak velocity becomes. This directly limits the peak velocities of the signal, which makes it a suboptimal filter in eye tracking applications. Nevertheless, it is a popular choice as it is easy to understand and implement. In total, we used 24 different window sizes (3:2:37 samples and 41:4:61 samples).

Fig. 3
figure 3

Exemplary step responses of (a) moving average, (b) Savitzky-Golay, and (c) Butterworth filters. The data are shown for a sampling rate of 1 kHz. The black line gives the original step function. The gray bars in the top show the two window sizes of the FIR filters centered on the step function. BW Butterworth, MA moving average, SG Savitzky-Golay, cXX cut-off frequency in Hz, oX filter order, sXX window size in samples (one sample corresponds to 1 ms at a 1-kHz sampling rate)

Full size image

An extension of the MA filter is the Savitzky-Golay filter (SG; Savitzky & Golay, 1964). SG filters fit a polynomial function of a given order over a moving window. They are optimal for polynomial or polynomial-like signals (Bromba & Ziegler, 1981), but show some stretching and some ringing at transients within half the window size (Fig. 3b). SG filters were applied for 23 window sizes (5:2:37 samples and 41:4:61 samples) and eight filter orders (2:1:9). In SG filters, larger window sizes and lower orders produce smoother but also more stretched-out signals. Higher orders increase the degree of the ringing at transients (compare the third- and fifth-order filters in Fig. 3b). Since the order of an SG filter must not be larger than its size, 12 of the 184 possible order-size combinations were not valid SG filter settings and thus only a total of 172 SG filters were applied.

Finally, the Butterworth filter (BW; Butterworth, 1930) was chosen as a commonly used IIR filter. These filters have very flat passbands and steep roll-offs, but introduce some stretching and long ringing artifacts at transients (thus the name infinite impulse response; see Fig. 3c). 175 BW filters were designed for seven different orders (3:1:9) and 25 cut-off frequencies (15 Hz, 25:5:45 Hz, and 46:1:64 Hz) at each of the eight sampling rates. Since BW filters do not have a linear phase response as do the FIR filters, zero-phase filtering was applied. This is important to avoid phase distortions in the signal. However, it effectively doubles the filters’ orders and squares the magnitudes of their transfer functions. In addition, it makes the filters non-causal, since the backward pass uses samples from the “future.” Similar to SG filters, higher orders increase the degree of the ringing artifacts (compare the third- and ninth-order filters at the 60-Hz cut-off frequency in Fig. 3c) but increase the steepness of the roll-off. Lower cut-off frequencies result in some stretching (contrast the curves for the 30-Hz and 60-Hz cut-offs in Fig. 3c).

Figure 4 shows examples of the transfer functions of some of the applied filters. As can be seen, the cut-off frequencies of MA filters decrease with window size (Fig. 4a; 26 Hz at 17 samples vs. 7 Hz at 61 samples). Where this is also true of SG filters, the cut-off frequencies are also affected by the filter order (Fig. 4b; 63 Hz with order 3 vs. 101 Hz with order 5 at 17 samples and 17 Hz with order 3 vs. 28 Hz with order 5 at 61 samples). In case of the BW filters, the cut-off frequency is a design parameter and the order has only negligible effects on it (Fig. 4c); instead, it affects the roll-off of the transfer function.

Fig. 4
figure 4

Exemplary transfer functions of (a) moving average, (b) Savitzky-Golay, and (c) Butterworth filters. The data are shown for a sampling rate of 1 kHz. The thin gray line gives the -3 dB-point for the cut-off frequencies (dotted lines and markers). The transfer functions of the Butterworth filters are shown for zero-phase filtering which results in lower actual cut-off frequencies than specified in the design parameters. BW Butterworth, MA moving average, SG Savitzky-Golay, cXX cut-off frequency in Hz, cXX cut-off frequency in Hz, oX filter order, sXX window size in samples (one sample corresponds to 1 ms at a 1-kHz sampling rate)

Full size image

As not all filters are applicable for all sampling rates (e.g., too few samples, too low sampling rate), only the valid settings have been used at the according sampling rates.

Determining the best filter

After applying the filters to the noisy data, peak velocities were determined at the time of occurrence of the ideal peak velocity.

The velocity divergence measure

Although peak velocity was the most important parameter for our simulation, simply using this as a quality measure fell short of real-life conditions. For example, high-order BW filters with low cut-off frequencies might produce peak velocities close to the original value but will introduce strong ringing artifacts at saccade on- and offset. Therefore, we needed a more complete measure for a filter’s performance. We were inspired by the discrete Kullback-Leibler divergence (Kullback & Leibler, 1951), which measures the information loss when assuming some empirical distribution P instead of the true distribution P*, and is given by:

$$ {D}_{KL}\left(P,{P}^{*}\right)=-{\displaystyle {\sum}_{i,\ {P}_i\ne 0,{P}_i^{*}\ne 0}{P}_i^{*}}\cdot { \log}_2\left(\frac{P_i}{P_i^{*}}\right),\kern2em {D}_{KL}\left(P,{P}^{*}\right)\ge 0 $$
(15)

Starting from this formulation, we derived a velocity divergence measure to rate a filter’s performance. Note that we only used saccades with positive amplitudes for the simulation, resulting in non-negative, non-constant velocity profiles for the noise-free data. This is important for the following derivation to hold (although the adaptation for negative amplitudes is quite straight forward) and allows us to omit the sign since we are only interested in absolutely low values.

For a specific sampling rate f s , let \( {V}_{f_s}^{*} \) be the corresponding noise-free velocity profile. \( {V}_{f_s}^{*} \) will serve as a reference since it is the best result a filter will achieve at this particular sampling rate. Let

$$ \Delta \left({X}_{f_s}(t)\right):=\left\{\begin{array}{cc}\hfill 0\hfill & \hfill {V}_{f_s}^{*}(t)\le {V}_{f_s\ min}^{*}\hfill \\ {}\hfill {X}_{f_s}(t)-{V}_{f_s\ min}^{*}\hfill & \hfill otherwise\hfill \end{array}\right. $$
(16)

be any velocity \( {X}_{f_s}(t) \) at that sampling rate centered on the minimum noise-free velocity \( {V}_{f_s\ min}^{*} \) found at saccade onset t0 or offset t1. This can be used to define the normalization function with respect to \( {V}_{f_s}^{*}(t) \):

$$ \Theta \left({X}_{f_s}(t)\right)\ :=\frac{\Delta \left({X}_{f_s}(t)\ \right)}{{\displaystyle {\sum}_{t_i={t}_0}^{t_1}}\Delta \left({X}_{f_s}\left({t}_i\right)\ \right)} $$
(17)

The normalized, noise-free velocity \( \Theta \left({V}_{f_s}^{*}(t)\right) \) now resembles a probability density function. Nevertheless, the normalized noisy velocity \( \Theta \left({V}_{f_s}(t)\right) \) still can take any real value due to the contribution of the noise. Consequently, the argument of the log-term in equation (15) is not assured to be strictly positive. However, this term can be approximated by the absolute gain deviation:

$$ \gamma (t):=\left\{\begin{array}{cc}\hfill 1\hfill & \hfill {V}_{f_s}^{*}(t)\le {V}_{f_s\ min}^{*}\hfill \\ {}\hfill \left|\frac{\Theta \left({V}_{f_s}(t)\right)}{\Theta \left({V}_{f_s}^{*}(t)\right)}-1\right|+1\hfill & \hfill otherwise\hfill \end{array}\right.,\kern1em \gamma (t)\ge 1 $$
(18)

Since the number of samples in the velocity profiles differs depending on the duration of the saccade, the final velocity divergence measure has to be normalized additionally by the number of samplesFootnote 1 and is given by:

$$ \widehat{D}\left(V,{V}_{f_S}^{\ast}\right)\kern.2em :=\frac{1}{n}{\displaystyle \sum_{t={t}_0}^{t_1}}\Theta \kern.1em \left({V}_{f_s}^{\ast }(t)\right)\cdot { \log}_2\left(\gamma (t)\right),\kern1em \widehat{D}\left(V,{V}_{f_S}^{\ast}\right)\ge 0 $$
(19)

\( \widehat{D}\left(V,{V}_{f_S}^{\ast}\right) \) is measured in bits per sample. It can be conceived as the weighted mean of the log-gain deviation with most weight put on the peak velocity and the equality:

$$ \widehat{D}\left(V,{V}_{f_S}^{\ast}\right)=0\kern.8em \Longleftrightarrow \kern.8em V={V}_{f_S}^{\ast } $$
(21)

Evaluation of filter performances

The velocity divergence of all saccades was determined for each filter at each noise setting and all repetitions. The mean and its 95% confidence interval (CI) were computed over all repetitions for each saccade. The mean was chosen since the divergence is an average over samples and thus the values are expected to be normally distributed.

The final measure of a filter’s performance \( {\overline{D}}_{{}_{Filt}} \) was the median of the upper CI of the velocity divergence over all amplitudes. The median accounts for the non-normal distribution of the peak velocities, and the upper CI provides a worst-case estimate incorporating any residual noise after filtering. Filters were rated according to \( {\overline{D}}_{{}_{Filt}} \) with lower values corresponding to better performance. Since several filters may perform similarly well, we computed the filter performance score \( {P}_{Filt}:=100\cdot {\overline{D}}_{\min }/{\overline{D}}_{{}_{Filt}} \) at a given noise setting and sampling rate with respect to the divergence of the best filter \( {\overline{D}}_{\min } \) at that setting. Filters with a performance score of P Filt  ≥ 95% were considered similar to the best filter. It has to be noted that at very low noise levels applying even the best filter increases the divergence and thus the unfiltered data will have performance scores larger than 100%. Any filtered or unfiltered data was ranked according to their divergence (“best-best” filter = 1, “worst-best” filter = k). Among these similar filters, the filter which had the least effect on the signal was determined. For this purpose we defined the “filter softness” as not filtering the data at all (softest), using a moving average (minimum size; soft), or a Savitzky-Golay filter (minimum size and maximum order; harder), or ultimately applying a BW filter (maximum cut-off frequency, minimum order; hardest).

Finally, as a more graspable quality measure, we will report the median \( {\overline{\phi}}_{Filt} \) of the velocity score \( \phi (a):=100\cdot \left(1-\left|\frac{{\overset{\smile }{v}}_{\max }(a)}{v_{f_s\ max}^{*}(a)}-1\right|\right) \) with \( {\overset{\smile }{v}}_{\max }(a) \) as the lower CI on the mean of the peak velocity in the filtered data for a given amplitude a over all repetitions. ϕ(a) is equal to the velocity gain \( \Gamma (a):=100\cdot \frac{{\overset{\smile }{v}}_{\max }(a)}{v_{f_s\ max}^{*}(a)} \), when both peak velocities have the same sign and \( \left|{\overset{\smile }{v}}_{\max }(a)\right|\le \left|{v}_{f_s\ max}^{*}(a)\right| \). Similar to \( {\overline{D}}_{Filt} \), \( {\overline{\phi}}_{Filt} \) is a worst-case estimate. It measures how close the filtered peak velocity resembles the noise-free peak velocity \( {v}_{f_s\ max}^{*} \) at that sample rate. We explicitly did not compute the velocity score with respect to the ideal peak velocities, since this would have led to values below 100% for noise-free data in cases where the sampling rate already affects the peak velocity. Filters achieving velocity scores of ≥95% can be considered as restoring the peak velocity reasonably well within the limits of the sampling rate.

Due to the varying signal-to-noise ratio of saccades of different amplitudes at a fixed noise composition, we also determined the best and softest filters for three amplitude ranges (micro-/small saccades: 0.5–5°, normal saccades: 5–15°, and large saccades 15–60°). The velocity divergence was computed separately for the amplitudes in these ranges.

Frequency analysis

For the analysis of the frequency content of the saccades, position and velocity profiles were clipped at on- and offset with a 25-ms margin similar to Harris, Wallman, and Scudder (1990). All saccades in the following figures are shown over this time course. These clipped saccades were then zero-padded to achieve a frequency resolution of 10 mHz and submitted to a fast Fourier transform to get the power spectra.

Results

Modeling saccades

Several real saccades of a healthy subject and a Parkinson’s patient are shown in comparison to matched simulated saccades in Fig. 5; whereas the simulated velocity profiles are symmetric, this is in general not the case for the measured data (Fig. 5a, middle panel). Only the visually-guided prosaccades of the healthy subject are close to the simulated waveforms. The antisaccade of the healthy subject and, particularly, the data from the Parkinson’s patient, however, are more skewed. Nevertheless, the Fourier spectra of the position profiles (Fig. 5b) are quite similar. The small deviations of the simulated from the measured data above 10 Hz are a result of the asymmetrical shape in the measured data and reflect the fine structure of the velocity profiles.

Fig. 5
figure 5

Examples of (a) simulated and real saccades and (b) the according Fourier spectra of the position profiles. Saccades from a young, healthy subject (male, age 35 years) and a Parkinson’s patient (male, age 71 years) are shown. Measured data (solid lines) were recorded using a search coil (sampling rate: 1 kHz; zero-phase filtered at 48 Hz with a seventh-order Butterworth filter with respect to the highest noise level). The simulated saccades (dashed lines) were based on the Inchingolo model with individually fitted main sequence parameters for each subject and saccade type. The data were created using the same sampling rate, noise amplitude and filter conditions as the real data. The angles given in the legend depict target eccentricity. The spectra of measured and simulated saccades in (b) did not differ outside the displayed frequency range. Note the logarithmic scaling of the x-axes in (b). AS antisaccade (healthy subject), PS visually-guided prosaccade (healthy subject), Park visually-guided prosaccade (Parkinson’s patient)

Full size image

The effect of sampling rate

Figure 6 shows a noise-free 10°-saccade at all eight sampling rates. In each case, the final position (e.g., the amplitude) is reached at the end of the trace. However, the position profile is affected by strong aliasing at the lowest sampling rate (Fig. 6a, upper panel). This becomes even more visible in the velocity profile where the peak velocity drops below 50% of the ideal value at a sampling rate of 30 Hz (Fig. 6a, lower panel). Only data sampled at least with 240 Hz allows a reasonably good reconstruction of the peak velocity (inset in Fig. 6a, lower panel). This is paralleled in the Fourier spectra where the deformation of the velocity profile at the lower sampling rates is reflected by an earlier onset of the drop-off (Fig. 6b). In addition, the occurrence of the characteristic local minimum in the spectra of the position profiles is shifted to lower frequencies at the lower sampling rates. This shift is caused by the prolongation of the saccade at these sampling rates.

Fig. 6
figure 6

(a) Effect of sampling rate on a 10°-saccade and (b) according Fourier spectra. The inset in the lower panel of (a) gives a close-up of the peak velocities for sampling rates ≥240 Hz. Markers in (a) show sampling intervals. The spectra in (b) are constant below 2 Hz. Note the logarithmic scaling of the x-axes in (b)

Full size image

The effect of sampling rates on peak velocities for saccades of different amplitudes is given in Fig. 7a. As in the case of the 10°-example, the peak velocities are significantly reduced for sampling rates below 240 Hz. A comparison of the simulated data from Juhola et al. (Juhola et al. 1985; data from Fig. 3), who used a different model and higher peak velocities, shows a similar effect but with a stronger drop-off (dashed yellow line in Fig. 7a). We also subjected the visually-guided prosaccade of the healthy subject (see Fig. 5) to the downsampling procedure and computed the peak velocities at all sampling rates (dashed green line in Fig. 7a). The resulting peak velocities similarly show a reduction with sampling rate and closely resemble the results from the simulated saccade. In general, the drop-off is not the same for different amplitudes, indicating a differential effect of sampling rate on peak velocities depending on the size of the saccade. To get the full picture over the whole amplitude range, we computed the velocity gain for all sampling rates as shown in Fig. 7b. It can be seen that the velocity gain generally increases with amplitude. To see how well each sampling rate recovered the ideal peak velocities, we determined the minimal amplitude for which the velocity gain was ≥95% (circles in Fig. 7b). Only the highest sampling rates (500 Hz and 1 kHz) were able to restore the peak velocity over the whole amplitude range. At the other end, the lowest sampling rates (30 Hz and 50 Hz) reduced the peak velocity even at the largest amplitudes. Although the detrimental effect of sampling rate gets better with increasing frequencies, peak velocities sampled at 240 Hz still do not reach 95% of the ideal value for saccades with amplitudes below 3°.

Fig. 7
figure 7

(a) Peak velocities of several saccades at different sampling rates and (b) peak velocity gains for all modelled amplitudes and sampling rates. The dashed lines in (a) show the results from Juhola et al. (1985) for an artificial 20°-saccade (yellow) and the downsampled measured 10°-prosaccade of the healthy subject from Fig. 5 (green). Peak velocity gain in (b) is computed with respect to the ideal data. Markers in (b) show the minimal amplitude needed to reach a gain of 95% (circles; 30 Hz and 50 Hz did not reach the threshold; 500 Hz and 1 kHz were always above) and the medium gain over all amplitudes (triangles). Note the logarithmic scaling of the x-axes

Full size image

Some eye-tracker manufacturers upsample their data to a higher sampling rate for further analysis. To quickly see if this reduces the effect of sampling rate on the peak velocities, we used a simple spline interpolation on the noise-free position profiles of the sampling rates which did not reach 95% velocity gain over the whole amplitude range (30–250 Hz). The resulting profiles were interpolated to resemble a sampling rate of 1 kHz. Figure 8 shows the results for a 10°-saccade. It can be seen that the interpolation indeed increases the peak velocity but also introduces some ringing artefacts at the lowest sampling rates.

Fig. 8
figure 8

Effect of interpolation on a 10°-saccade. The profiles are interpolated to resemble a sampling rate of 1 kHz. Markers depict the peak velocities achieved with the original sampling rate (see Fig. 6a). XX us sampling rate of the original data

Full size image

To see the full picture over all saccade amplitudes, we again computed the velocity gain for the interpolated profiles and determined the according minimal amplitudes. Tbl. 1 shows the result in comparison to the original data.

The improvement in the minimal amplitude is quite pronounced, especially for 120 Hz, where the interpolation reduces the minimal amplitude such that also small saccades can be restored reasonably well.

The effect of filtering

Figure 9 shows the velocity divergence over the whole range of amplitudes for different sampling rates, noise levels and filters. In Fig. 9a it can be seen that the sampling rates do not yield the same values at the same amplitudes. However, they all follow the same trend. The decrease in the velocity divergence at higher amplitudes is due to the increasing signal-to-noise ratio, with higher amplitudes for fixed noise levels. At the lower sampling rates, slight rounding artifacts are visible resulting from the quantized durations at these sampling rates. Figure 8b shows the effect of applying a filter on the velocity divergence. A good filter simply removes the noise without changing the velocity profile. This is visible in the rather constant offset of the best filters with respect to the noisy data. A bad filter, on the other hand, alters the shape of the velocity profile either by residual noise or by introducing ringing or stretching artifacts. This can be clearly seen in the filters which yielded the highest velocity divergence values (“worst” filters).

Fig. 9
figure 9

Velocity divergence at (a) different noise levels and sampling rates and (b) filtered with the worst and best filters of each type. Data in (b) are simulated at 1 kHz and 0.1° noise amplitude. Note the logarithmic scaling of the axes. BW Butterworth, MA moving average, SG Savitzky-Golay, UF unfiltered, cXX cut-off frequency in Hz, oX filter order, sXX window size in samples (one sample corresponds to 1 ms at a 1-kHz sampling rate)

Full size image

Figure 10a shows the effect of applying the worst filter of each type (\( {\overline{D}}_{\mathrm{MA}} \) = 5.35 mbits/sample, \( {\overline{D}}_{\mathrm{SG}} \) = 2.2 mbits/sample, \( {\overline{D}}_{\mathrm{BW}} \) = 3.56 mbits/sample) to a 10°-saccade simulated with a sampling rate of 1 kHz. It can be seen that stretching occurs for all filter types, leading to profoundly lower peak velocities (MA 159°/s, SG 268°/s, BW 239°/s) compared to the unfiltered data (UF 333°/s). In addition, ringing is present for both the SG and the BW filter. Especially in the latter case, the ringing affects the shape of the velocity profile. Using the simple IVT method to detect a saccade under such conditions will alter the measured durations (UF 108 ms, MA 138 ms, SG 82 ms, BW 80 ms). Interestingly, the stretching introduced by the MA filter prolongs the saccade whereas this is not the case for the SG filter. The ringing introduced by the BW filter shortens the saccade and slightly increases its detected size (UF/MA/SG 10°, BW 10.6°).

Fig. 10
figure 10

Effects of worst-case filtering on (a) a simulated 10°-saccade with no noise and (b) a visually-guided prosaccade towards a 10°-target. Markers show time of peak velocity, and on- and offset of the saccades as detected with the IVT method. The insets give close-ups of saccade onsets. Data in (a) ares simulated at a 1-kHz sampling rate. Measured data in (b) correspond to the unfiltered 10°-example of the healthy subject in Fig. 5. Unfiltered raw data (0.018° AWGN amplitude, no mains noise contribution; gray line) are shown together with the best filter applied (black line). The unfiltered data is omitted for clarity in the close-ups. BW Butterworth, MA moving average, SG Savitzky-Golay, UF unfiltered, cXX cut-off frequency in Hz, oX filter order, sXX window size in samples (one sample corresponds to 1 ms at a 1-kHz sampling rate)

Full size image

Figure 10b shows a measured visually-guided prosaccade towards a 10°-target and the same worst filters applied. In addition, the best filter for the according noise level has been determined. The amplification of the noise through differentiation can be seen in the unfiltered velocity profiles (Fig. 10b, lower panel). The best filter (black line) simply removes the noise from the data and has no effect on the peak velocity or the shape of the profile. In contrast, the worst filters introduce similar artifacts as in the case of the simulated saccades.

The best filter settings

Unsurprisingly, refraining from using a filter for noise-free data is the best option at all sampling rates. Thus, we will not consider this case any further in the remainder of this section.

Figure 11 shows the decision tree for some of the highest ranking filters with a performance score ≥95% at two common sampling rates for the three amplitude ranges. It can be seen that the actual number of valid filters is rather small (60 Hz: 14 unique filters, 1 kHz: 31 unique filters) compared to the total number of 162 cases.

Fig. 11
figure 11

Decision tree of the filter choice at (a) a 60-Hz and (b) a 1-kHz sampling rates Node labels in the first layer give saccade amplitude ranges in degrees. The second layer shows the noise amplitudes in degrees. The third layer depicts the mains noise frequencies in Hz. Edge labels between the third layer and the leaf nodes show the mains noise proportion in percentages. Leaves are labeled with the according best or similar-to-best filter. Omitted nodes or edge labels indicate that the according factor has no influence on filter choice. For brevity, labels have been shortened: “x|y” denotes a set of values and “~x” indicates everything but the given value. Filters that are valid for several settings are depicted in bold face. BW Butterworth, SG Savitzky-Golay, UF unfiltered, cXX cut-off frequency in Hz, oX filter order, sXX window size in samples (one sample corresponds to 17 ms at a 60-Hz sampling rate and to 1 ms at a 1-kHz sampling rate)

Full size image

It can be seen that not filtering the data is still the best choice at low noise levels at 60 Hz (Fig. 11a) independent of mains noise composition in the normal saccade range. In contrast, at a 1-kHz sampling rate (Fig. 11b) it is better to apply a filter even under such conditions. Furthermore, at a 60-Hz sampling rate, BW filters only occur at the highest noise-levels, due to the limitation of their application at that sampling rate where only very strong filters (e.g., very low cut-off frequencies) can be used. At 1 kHz, however, this filter type is always the best choice.

Table 2 shows a comparison of the values suggested by Inchingolo and Spanio (1985) and the filter settings obtained by us. Whereas they used a third-order Butterworth filter with a cut-off frequency ≥50 Hz under all conditions, our results were located at higher orders and lower cut-off frequencies.

The full list of the best and softest filters for each condition is given in Supplementary Table 2 or can be requested from the first author.

Discussion

The present study developed a flexible and easy-to-use model for the generation of saccades of different amplitudes rooted in real main sequence data. In addition, we derived a similarity measure that can be used to compare velocity profiles of different saccades. We used the model and this measure to examine the effects of lowpass filtering on saccadic peak velocities at several noise compositions and sampling rates. Finally, the according best filters with the least influence on the velocity profiles were identified.

The saccade model

The proposed model should not be understood as an alternative to elaborated neuronal models for saccades (e.g., Meeter, Van der Stigchel, & Theeuwes, 2010). Rather, it is intended as a simple method for the generation of saccade-like profiles during the examination of processing procedures. In the development of saccade detection algorithms, this model can serve as a quickly accessible ground truth for performance tests under different conditions such as noise and sampling rate. This reduces the need for manually labeling real measurements, which can become especially unwieldy when different noise compositions should be considered. Our model not only paves the way for the comparison of a single algorithm at different settings, but also provides a means for the comparison of different algorithms or processing procedures as we did with the filters.

An obvious disadvantage of the model are the prolonged durations of the generated saccades. These are much longer than expected from measured main sequences due to the symmetry and long settling phases of the underlying hyperbolic functions. Real saccades are skewed with longer deceleration than acceleration phases. In general, this skewness is rather low at small amplitudes but becomes stronger as saccade size increases (Baloh et al., 1975), which can already be observed in the velocity profiles of the measured saccades shown in Fig. 5. Moreover, skewness is affected by other factors like fatigue, medication, or variations in peak velocity (van Opstal & van Gisbergen, 1987). Saccades measured under such conditions can be more readily fit by an altered Gamma function (Smit et al., 1987). This function simultaneously yields an elegant and reliable measure of the skewness of the velocity profiles (van Opstal & van Gisbergen, 1987). However, in the context of modeling saccades, the Gamma function has several drawbacks. First, it provides no simple solution for the inverse problem, e.g., finding saccade on- and offset, or any other point in time for a given velocity. Second, velocity profiles created with the Gamma function cannot be concatenated directly in a differentiable fashion. This makes it hard to compute biologically plausible position and acceleration profiles of sequences of saccades. Both drawbacks are not present in our model and despite the lack of skewness, the simulated saccades possess similar frequency spectra as real saccades. The position spectra of the simulated saccades also show a duration-dependent minimum and a 20 dB/decade roll-off (see Fig. 5a) as is expected from real saccades (Harris et al., 1990; Zuber, Semmlow, & Stark, 1968). Similarly, the velocity spectra show a larger bandwidth than the position spectra and the power is not 40 dB down before 50 Hz (see Fig. 6b, lower panel; Bahill et al., 1981), emphasizing the model’s validity for the analysis of filtering effects.

The effect of sampling rate

Sampling rate affects saccadic peak velocities, and previous studies recommended using at least 200 Hz (Inchingolo & Spanio, 1985), 250 Hz (Schmitt et al., 2007), or even 300 Hz (Juhola et al., 1985). Replicating these results, we found that peak velocities can be restored with an accuracy of 95% or more at sampling rates ≥250 Hz for saccades ranging from 1.6° to 60°. Interestingly, increasing the sampling rate from 500 Hz to 1 kHz does not substantially improve accuracy any further. On the other end of the spectrum, peak velocities were significantly reduced for sampling rates of 120 Hz and below – which is good in agreement with Juhola et al. (1985). This drop-off, however, was neither constant nor linear over all amplitudes and its shape became more curved at lower sampling rates. Such a non-linear effect makes it hard to compare results between studies using different sampling rates, since the peak velocities cannot be adjusted simply by adding an offset. This non-linearity is even worse within studies comparing saccades of different amplitudes, as the acquired main sequences will be significantly distorted. What might reasonably look like the velocity profile of a 10°-saccade at a sampling rate of 55 Hz (Gibaldi, Vanegas, Bex, & Maiello, 2016) might not be true for other amplitudes. Finally, the reduction in peak velocities reported here is a best-case estimate as we made sure that the actual time of the peak velocity coincides with a sample. Since this time is a one-point measure, it can easily be missed by an eye tracker during a real measurement (Andersson, Nystrom, & Holmqvist, 2010), resulting in further lowered peak velocities (Bahill et al., 1982). Therefore the reduction in peak velocities for lower sampling rates reported here might even be worse for real data.

For these reasons, studies intending to analyze saccadic peak velocities should either refrain from using low-sampling eye trackers, or apply more elaborated signal processing methods to restore the peak velocities to their original values (e.g.,Wierts, Janssen, & Kingma, 2008). For this purpose, upsampling (colloquially including interpolation) is a promising candidate. We found that simple spline interpolation improves the restoration accuracy for all the lower sampling rates (see Table 1). This improvement is most remarkable for the 120-Hz data. Whereas a minimal amplitude of 19° is outside the range of most naturally occurring saccades and experiments (Bahill, Adler, & Stark, 1975), 2.5° covers this range quite well and thus makes a 120-Hz eye tracker a more valid choice. However, these results should not be understood as a complimentary ticket to use upsampling as a universal cure for sampling rate related peak velocity reduction. For example, when using interpolation, the underlying function has to be at least twice differentiable and should be monotonous between samples. The latter is not the case for spline interpolation as can be seen from the ringing artifacts at 30 Hz and 50 Hz in Fig. 8. Therefore, other functions should be used (e.g., Hessels et al., 2016). On the other hand, when using upsampling in its strict sense, an additional filter step is required which might again adversely affect the peak velocity. Finally, we used the noise-free data for the interpolation which is not available in real life. Therefore, further research is needed to explore the effects of noise and upsample method on peak velocity restoration before drawing any ultimate conclusions.

Table 1 Effect of interpolation on the minimal amplitude needed to reach a velocity gain of 95%. The original data are sampled at the depicted sampling rates. The interpolated data resemble a sampling rate of 1 kHz
Full size table

The effect of filtering

Our results show that choosing a filter is not a straight forward task but requires a thorough analysis of the raw data. At very low noise levels it might be better to refrain from using a filter at all. At very high noise levels, the application of a strong filter induces additional artifacts that change the shape of the signal (see Fig. 10). Besides the noise amplitude itself, the sampling rate of the eye tracker, the mains noise composition, and the amplitudes of the measured saccades are important factors to consider. Previous studies have provided guiding values based on the frequency content of saccades and considerations of computational efficiency (Bahill et al., 1981; Inchingolo & Spanio, 1985; Jäntti et al., 1984; Juhola, 1986), but many filter choices still seem to be “guided by heuristics, or ‘rules of thumb’ motivated by visual inspection of the data” (Holmqvist et al., 2011, p. 48). The measured data in Fig. 10b shows that such visual inspection can be quite erroneous if not performed carefully. The worst SG and BW filters both produce smooth and rather inconspicuous position profiles. However, the resulting peak velocities are quite off in comparison to the best filter. Therefore, the extension of the previous findings to a wide range of filter settings, sampling rates, noise compositions, and amplitude ranges and the resulting filter list should help to improve the filter choice. In comparison to the guiding values provided by Inchingolo and Spanio, the best filters identified here are slightly harder (see Table 2). However, they differ quite a lot over the noise range specified by Inchingolo and Spanio. These authors also used only third-order BW filters which have a flatter roll-off than higher order filters and therefore require higher cut-off frequencies to not affect the frequency component of the signal. That fact further emphasizes the complexity of the correct filter choice.

Table 2 Best filter settings compared to Inchingolo and Spanio (1985). All filters are Butterworth. Mains noise frequency was 60 Hz. Settings are shown for the whole amplitude range (0.5–60°) as the original data used saccades ranging from 2–70°
Full size table

A general observation from the best filter list is the increasing prevalence of BW filters at higher sampling rates, ending in a total absence of other filter types at 1 kHz (see Supplementary Methods, Section 3). This can be explained by considering the smoothness of the signal. At higher sampling rates more noise is present in the data. Such high frequency noise can be more efficiently suppressed by the steeper roll-off of the BW filters compared to the two FIR filters, resulting in a smoother signal. On the other end of the spectrum, the number of valid BW filters becomes smaller as sampling rates decrease, due to the cut-off frequency as a design parameter. The latter can obviously not be higher than the sampling rate, which limits the possible selection to BW filters with increasingly lower cut-off frequencies. This, however, is not needed since the signal is already smoothed by the lowpass effect of the lower sampling rates. Applying such a BW filter to a low-sampled signal leads to strong ringing artifacts, whereas the left-over noise from the FIR filters might be acceptable.

Signal smoothness is also an important factor during event detection. Threshold-based methods heavily rely on smooth data to avoid false positives from noise oscillations (Salvucci & Goldberg, 2000). Using such a method after applying a hard filter will result in unintuitive outcomes. Whereas an MA filter will lead to prolonged saccades due to the introduced stretching in time, a same-sized SG filter will actually shorten their duration (see Fig. 10a). In addition, the FIR filters most often fail to completely suppress the noise, while BW filters produce very smooth signals but introduce long ringing artifacts. Both points are less problematic for dispersion-based algorithms that are more robust to noise, but have a poorer temporal resolution (Nystrom & Holmqvist, 2010). It would be interesting to see which filters are most suitable for the different methods and what amount of left-over noise or introduced ringing is tolerable. This could be a future application of our saccade model, further extending the finding by Andersen, Larsson, Holmqvist, Stridh, and Nystöm (2016).

Methodological considerations

The prolongation of the modelled saccades’ durations is most pronounced at very large amplitudes, since small saccades are approximately symmetrical (Baloh et al., 1975). However, the frequency content of our saccades is similar to fairly skewed saccades compared to the Gamma function-based model (Smit et al., 1987; a 10°-saccade is equal to γ = 9; see Supplementary Methods, Section 4). Therefore, we remain confident that our results are applicable to real data. Nevertheless, the optimal filter settings presented here are a best-case estimate and should be used as lower bounds. Under conditions where extremely skewed saccades are expected (e.g., fatigue, drugs, senescence), we recommend using the softest filter instead and carefully examining the results.

The longer durations and the symmetric velocity profiles should also be kept in mind when testing an algorithm for saccade detection, if duration limits or a differential treatment of on- and offset are applied. Care should be taken not to adapt the algorithm too much to these idealized saccades to retain good detection results with real saccades.

Despite the main sequence relationship, there is still a considerable amount of subject- and experiment-related variability in peak velocities as mentioned in the Introduction (Becker & Fuchs, 1969; Mack & Ilg, 2014; Smit et al., 1987). Therefore, fitting a main sequence to experimental data will always yield slightly different model parameters and thus produce somewhat dissimilar peak velocities. Since we wanted to provide widely applicable filter settings we chose a main sequence model which provides a solid basis. For this reason, the model by Inchingolo and Spanio and its parameters are outstanding as the model is well formulated, computed with a carefully selected set of visually-guided saccades over a large range of amplitudes and it apparently “provides more accurate fitting of the experimental data than the exponential laws” (Inchingolo & Spanio, 1985, p. 686). However, when comparing modelled saccades to a particular set of saccades (as we did in Fig. 5) one should compute specific main sequence parameters for the dataset to obtain similar peak velocities. The Matlab function of our saccade model (Mack et al., 2016) additionally provides a power law for the main sequence with two sets of parameters (Bahill et al., 1981; Baloh et al., 1975). It also allows the use of completely arbitrary parameters for all supported models.

Besides saccadic peak velocity, other oculomotor parameters also depend on the filter settings. For example, ocular drift speed is affected by the cut-off frequency of the filter (Ko, Snodderly, & Poletti, 2016), and fractal analysis relies on an unaltered correlational structure of the data (Coey et al., 2012) – which is not preserved by many filters. As our results are based on the saccadic velocity profile, they should not be blindly used in these cases.

Finally, we are aware that the unfiltered raw data are not always accessible since manufacturers might implicitly apply possibly unspecified filters. Additionally, in more complex eye-tracking setups, it might be convenient to convert the digital raw data back to an analog signal for tight experimental control or synchronized recording with physiological data. This may involve another analog filtering step. Under such circumstances, our assumption of one-step filtering does not hold anymore. Furthermore, a cascade of filters might perform better than just one filter. However, we tried to avoid these more complex scenarios to keep our results as easily and widely applicable as possible. Future work should explore other processing configurations and more complex filter types (e.g.,Erer, 2007; Tibshirani, 2014).

Conclusion

The simplicity and the resulting ease of application of the proposed saccade model can pave the way for many interesting comparisons of saccade detection algorithms and their settings. The list of the best filters should bolster future filters choices by providing a more data-driven foundation.