1. Introduction

Signal models of physiological processes are commonplace [1]; however there has been very little research regarding models of dynamic ocular aberrations other than defocus. Roberts [2] asserted that measuring and understanding the dynamic performance of the eye is critical to the future development of advanced customised refractive surgeries, which is in itself important for generating an overall improvement of refractive surgery results. Therefore, any improvements in the measurement and modelling of dynamic aberrations could be beneficial in areas such as the development of customised ablation algorithms. Such models could also be of interest in the study of the relationship between dynamics and visual performance [3], the development of customised contact lenses [4], and in the design of ocular adaptive optics systems with real-time correction [5].

Many dynamic properties of ocular aberrations have been established by previous authors [6, 7]. In particular, the time evolutions of some Zernike aberration modes have been shown to be correlated with each other [8], as well as with the cardiopulmonary system [9, 10]. Previous work has shown that the power spectral density (PSD) of aberration signals tends to follow a power-law relationship for frequencies below approximately 10 Hz, when the subject is not actively accommodating [3,7]. The dynamics of Zernike modes have also been shown to exhibit features such as non-stationarity, harmonic distortion, and harmonic amplitude modulation [6, 11].

Iskander et al. [11] proposed a methodology for modelling the dynamics of higher-order ocular aberrations using a parametric, amplitude-modulated, frequency-modulated (AM-FM) approach. This study focused in particular on the Zernike spherical aberration and coma coefficients, however measurements of just 5 s in duration (without taking into account the effects of subjects’ blinking) were considered. Other work has demonstrated a potential application of a model of dynamic aberrations, by employing a deformable mirror to generate aberrations in real time [5].

In this work, we present a methodology for constructing power-law models of the temporal dynamics of ocular aberrations. The modelling approach is based on a numerical method for coloured-noise generation, developed by Billah and Shinozuka [12]. We describe the connection between power-law relationships in the frequency domain and temporal stationarity, and how this may be exploited in the modelling of dynamic aberrations. We attempt to validate the performance of our model by comparing real and simulated data.

2. Methods

2.1. Data collection

All data was collected with a custom-built aberrometer based on a Shack-Hartmann wavefront sensor. The subject was a 26 year-old male. For each measurement, data was collected from the subject’s dominant eye, whereas the non-dominant eye was covered with an eye patch. Eye dominance was assessed through an eye examination by a qualified optometrist prior to our experiments. Data was collected over 40 s, with an active pupil diameter of 5.4 mm. The sampling frequency was 100 Hz, giving a total of 4,000 samples per trial. The defocus term had the greatest magnitude, but there were significant contributions from astigmatism, coma, spherical aberration, and other higher-order terms. The system noise was white, with an RMS value of approximately 3 nm, as evaluated with an artificial eye (consisting of an 18 mm singlet coupled with a screen in the focal plane).

Analysis of dynamic ocular wavefront sensor measurements often takes the form of analysis of the Zernike coefficients as time series [3, 6–9, 13–17]. The temporal reconstruction of ocular wavefronts using Zernike polynomials can be represented as given in Eq. (1) [6]:

where W(ρ, θ;t) is the dynamic wave aberration function, m and n are the radial and azimuthal orders respectively, $Z_n^m\left(\rho ,\theta \right)$ is the Zernike polynomial expansion term, $c_n^m\left(t\right)$ gives the corresponding Zernike coefficients, and ε(ρ, θ;t) represents the modelling error. For our studies of aberration signals, we used an expansion of Zernike terms up to and including the 8^th radial order.

Measurements of dynamic ocular aberrations of more than a few seconds in duration are inevitably degraded by the subject’s blinking. This leads to discontinuities and/or missing values in the measured signal [6, 28]. Though the blink interval is typically on the order of 250 ms [6], the transient components of the measured aberration signal associated with the blink can last for up to a second or more [28]. For all measurements, we systematically removed any part of the signal that was related to the blink interval: tear film break-up and build-up, eye movement, and increased noise in the reconstruction of the Zernike modes of the wavefront due to a vignetted pupil. To incorporate these phenomena would significantly increase the complexity of any modelling approach. Therefore, as an approximation, we did not include them as part of the modelling objectives in this study.

2.2. Modelling approach

The modelling method is based on the assumption of a linear or piecewise linear power-law structure of the power spectral density (PSD) of the aberration signals. Synthesis of signals is based on a modification of an algorithm for the generation of power-law processes (or “coloured noise”) [12]. This algorithm has also been utilised as the basis for modelling approaches used in the field of EEG research [18], which deals with the study of signals that have many similarities to ocular signals in terms of their non-stationarity, power-law form of the spectral density, and dynamic range [19].

We initially assume that the power spectral density of the process to be modelled can be approximated by a power-law of the form given in Eq. (2):

where f is the cyclic frequency, γ is the power-law exponent, and c is a constant. In general, for a finite-length, continuous-time signal x(t), the power spectral density can be written as shown in Eq. (3) [1]:

where the angular brackets denote ensemble average. X(f) is the Fourier transform of x(t), and can be written in the complex exponential form of Eq. (4):

where ϕ(f) denotes the phase spectrum. To synthesise a time-domain signal from X(f), we can specify the magnitude and phase spectra, and then take the inverse Fourier transform. We assume that the magnitude spectrum follows a power-law of slope γ/2.

In practice, the method described above can be implemented on a computer using an FFT algorithm. A detailed discussion of the discrete implementation of such an algorithm is given by Billah and Shinozuka [12]. Rankine et al. [18] suggested that the “roughness” of real-world power spectra (as opposed to the “smooth” power-law form of the model) can be approximated by performing the synthesis for several different realisations of ϕ(f), and then summing together the resultant signals in the time-domain.

2.3. Parameter estimation

For estimating the slope value to fit to the PSD, the well-known Hill estimator was used [22]. This method is simple to implement and gives an asymptotically normal and consistent estimate [23]. For the phase spectrum, ϕ(f), one may use the phase spectrum of a measured signal as a “surrogate” phase. Alternatively, Billah and Shinozuka [12] as well as Kasdin [21] asserted that the phase spectrum of a power-law process can be assumed to be composed of random phase angles uniformly distributed over the interval [0,2π]. Artificial dynamic aberration signals can thus be generated by randomly choosing phase angles from a uniform distribution, to go with a PSD slope value arbitrarily chosen or empirically selected based on estimates from real data.

2.4. Validation of model

To justify the use of the model, we must assess its validity. Unfortunately, for non-stationary signal models, no well-established validation methodology exists [26]. Rankine et al. [18] suggested using the Pearson product-moment correlation coefficient to validate simulated non-stationary data, in part due to its widespread familiarity. The application of this measure to processes that exhibit non-stationarity and non-linearity is however questionable. Muma et al. [20] described a method of using time-frequency coherence estimates to obtain a normalised estimate of the coupling between two non-stationary signals. We are grateful to the authors for providing us with scripts for the implementation of this method. In the absence of a more rigorous validation approach, we use the time-frequency coherence between real and simulated signals as a means to assess model performance. We do not claim this to be a conclusive or quantitative method of validation (our use of the computed phase spectrum of measured data as a model parameter means that a significant level of coupling is to be expected), but merely an indicator that the simulated data does indeed resemble measured data to a reasonable degree.

3. Results

Figure 1 (left) shows the measured time-evolution of several Zernike coefficients; Fig. 1 (right) shows the corresponding PSD estimates (via the Lomb-Scargle periodogram [24,25]). As mentioned earlier, the temporal behaviour of ocular aberrations is in general non-stationary. This is illustrated by Fig. 2, which gives a time-frequency representation of a measured Zernike astigmatism $(Z_2^{-2})$ signal. The time-frequency plot was obtained using the ZAM distribution, a member of Cohen’s class of distribution functions [32]. The plot shows low frequency content modulated in the < 0.5 Hz range over the full duration of the 40 s trial. There is significant power at approximately 1 Hz, which also varies in frequency and is likely to be related to cardiopulmonary activity [10, 20].

 

Fig. 1 Dynamics of selected Zernike aberrations (left) and the associated PSD estimates (right).

Download Full Size | PDF

 

Fig. 2 Time-frequency representation of Zernike astigmatism for subject CML (via the ZAM distribution).

Download Full Size | PDF

To assess the validity of our model, we generated signals from the model using parameters obtained directly from real data. From a measured signal, we determined the measured phase spectrum and an estimate $\widehat{\gamma }$ of the power-law slope, and used these as the model parameters ϕ(f) and γ respectively. We then synthesised a time-domain representation using these model parameters. For all measurements, our estimates of the power-law slope on a doubly logarithmic plot were greater than unity, and were similar to the PSD slope estimates of previous authors [3, 7]. Figure 3 compares the dynamics of measured Zernike spherical aberration to a simulated signal generated using our modelling method, while Fig. 4 shows the corresponding PSD estimates. Figure 5 is a plot of the measured phase spectrum ϕ(f) that is used as the surrogate phase.

 

Fig. 3 Comparison of measured and simulated Zernike spherical aberration signals in the time-domain.

Download Full Size | PDF

 

Fig. 4 Comparison of measured and simulated Zernike spherical aberration signals via their estimated power spectral densities.

Download Full Size | PDF

 

Fig. 5 Discrete Fourier phase spectrum computed from measured data and used as the “surrogate” phase to produce the simulated aberration signal of Fig. 3. This Fourier phase spectrum representation was computed using an interpolation algorithm to replace missing data points.

Download Full Size | PDF

Figure 6 shows the time-frequency coherence between the real and modelled signals of Fig. 1. The value of the time-frequency coherence function is considerably close to unity for most values of t and ω, which indicates that the simulated signal has spectral content that is close to that of the original measured signal from which the model parameters were estimated.

 

Fig. 6 Time-frequency coherence between real and simulated signals for Zernike spherical aberration.

Download Full Size | PDF

4. Discussion

4.1. Performance of models

Rankine et al. [18] suggested that greater flexibility in the modelling of non-stationary data can be achieved by allowing the estimated power-law exponent $\widehat{\gamma }$ to be time-varying, i.e. $\widehat{\gamma }$ = $\widehat{\gamma }$(t). We have observed that $\widehat{\gamma }$ can vary significantly over a single epoch of data. However, by allowing $\widehat{\gamma }$ to vary over time in our model, we observed only marginal improvement in the model performance at best. This could be partly due to the difficulties in validating model performance. However, we feel that the potential benefits in allowing $\widehat{\gamma }$ to be time-varying must be weighed against some drawbacks. Firstly, the time-scales upon which γ varies significantly are difficult to determine. Secondly, by attempting to estimate γ over shorter time-scales, one must take into account the additional bias in the estimate of γ.

4.2. Two-slope model

The shape of the PSD varies with the level of accommodative effort [27, 28]. When the eye actively accommodates, the shape of the PSD in this range often exhibits two distinct slopes with a “breaking frequency” between them [28]. This effect is not observed when the eye is in its unaccommodated state, such as in Fig. 1. For the cases where the estimated PSD is better fitted by two slopes rather than one, the model can be adapted accordingly. A bilinear fit can be performed by adopting a regression approach and imposing a constraint on the regression function. In this case, $\widehat{\gamma }$ = $\widehat{\gamma }$(f), with $\widehat{\gamma }$(f) given by Eq. (5):

where γ ₁ and γ ₂ are constants representing estimates of the PSD slope for lower and higher frequencies respectively. The two slope regions are separated by the breaking frequency f_br.

4.3. Relationship Between Power-Law Slope and Stationarity

Kasdin [21] stipulated that for power-law processes, the stationarity of the process is related to the slope γ of the PSD. In particular, the slope of the PSD at low frequencies is indicative of whether the process is stationary or not. This can be summarised as follows:

We illustrate this notion in Fig. 7, with some examples of signals synthesised using the two-slope model. The higher frequency slope was fixed to γ ₂ = 2.0. We then generated signals using five different values of γ ₁. In each case, an identical random phase spectrum (composed of uniform random variates distributed over the interval [0,2π]) was used. To employ a random phase spectrum (rather than a measured one), is potentially more useful as a general approach to modelling the temporal behaviour of the eye’s aberrations. Both methods can be used to generate either stationary or non-stationary signals with frequency content shaped by $\widehat{\gamma }$(f), as required. As the value of γ ₁ exceeds unity, the effect on the stationarity of the signals becomes clearly visible (for example, the mean of the process wanders significantly). We can consider the value of γ ₁ to relate to the “statistical stability” of the signal, rather than to stationarity in a more rigorous sense. It is useful in this context to think of stationarity in a relative sense, as truly stationary signals are mathematical constructs that do no not typically occur in real-world data [29].

 

Fig. 7 Illustration of the relationship between statistical stability about a mean and the PSD slope at low frequencies. All signals were generated using the two-slope model, with an identical phase spectrum generated from a uniform distribution.

Download Full Size | PDF

4.4. Application of models

The model for the dynamics of ocular aberrations was originally conceived as a way to extend static models of aberrations over a large population, such as the model developed by Thibos et al. [30]. If the dynamic model parameters could be similarly estimated over a large population, the authors’ idea of creating a database of “virtual eyes” could be extended to include the dynamics of aberrations.

A notable feature of the aberration signals is their dependence on the mean accommodative effort [28]. The modelling technique described in this work could similarly be used to model dynamic accommodation, as measures of accommodation can be computed from Zernike terms [31]. A model of dynamic accommodation could be adapted to include the state-dependence by modelling the relationship between the power-law slope values and the accommodative effort.

It should be noted that in our modelling efforts we have not attempted to compensate for the movements of the eye during measurement. Given that the ultimate aim of this work is to simulate the temporal behaviour of the eye’s intrinsic aberrations (rather than those measured by the wavefront sensor), it follows that these movements should ideally be compensated for, e.g., by using some active tracking of the pupil.

Another property that we have not directly addressed in our modelling approach is the non-linearity of aberration and accommodation signals. We previously discussed the presence of harmonic distortion in our measured signals, which suggests the presence of non-linearity. This is not surprising for a physiological process. Non-linearity has been directly addressed in EEG signal models [26], and we feel that proper study of non-linearity in ocular signals is also warranted (e.g., using bispectrum or bicoherence methods) for the improvement of the modelling effort.