 Research
 Open access
 Published:
Detecting fixation on a target using timefrequency distributions of a retinal birefringence scanning signal
BioMedical Engineering OnLine volume 12, Article number: 41 (2013)
Abstract
Background
The fovea, which is the most sensitive part of the retina, is known to have birefringent properties, i.e. it changes the polarization state of light upon reflection. Existing devices use this property to obtain information on the orientation of the fovea and the direction of gaze. Such devices employ specific frequency components that appear during moments of fixation on a target. To detect them, previous methods have used solely the power spectrum of the Fast Fourier Transform (FFT), which, unfortunately, is an integral method, and does not give information as to where exactly the events of interest occur. With very young patients who are not cooperative enough, this presents a problem, because central fixation may be present only during very shortlasting episodes, and can easily be missed by the FFT.
Method
This paper presents a method for detecting shortlasting moments of central fixation in existing devices for retinal birefringence scanning, with the goal of a reliable detection of eye alignment. Signal analysis is based on the Continuous Wavelet Transform (CWT), which reliably localizes such events in the timefrequency plane. Even though the characteristic frequencies are not always strongly expressed due to possible artifacts, simple topological analysis of the timefrequency distribution can detect fixation reliably.
Results
In all six subjects tested, the CWT allowed precise identification of both frequency components. Moreover, in four of these subjects, episodes of intermittent but definitely present central fixation were detectable, similar to those in Figure 4. A simple FFT is likely to treat them as borderline cases, or entirely miss them, depending on the thresholds used.
Conclusion
Joint timefrequency analysis is a powerful tool in the detection of eye alignment, even in a noisy environment. The method is applicable to similar situations, where shortlasting diagnostic events need to be detected in time series acquired by means of scanning some substrate along a specific path.
Introduction
The fovea, which is the most sensitive part of the retina, is known to have birefringent properties, i.e. it changes the polarization state of light upon reflection. The main cause for this birefringence are the Henle fibers surrounding the fovea, arranged in a radially symmetric pattern. When illuminated or scanned with polarized light, the Henle fiber layer produces the macular polarization cross or “bow tie”, which is a windmill shaped pattern centered about the fovea [1–3]. The strength and contrast of the polarization cross and the orientation of its bright parts depend on the orientation and degree of polarization of the light striking the retina, which is a function of both the instrumentation and the individual corneal properties [4, 5]. In recent years, the birefringent properties of the human eye have been used to identify the position of the fovea and the direction of gaze. This allows for one to check for eye alignment and strabismus (crosssightedness), a risk factor for amblyopia (called also “lazy eye”), which can potentially lead to a loss of sight in the affected eye. Screening techniques for amblyopia have been reported that are based on the birefringence signal derived from foveal circular scanning [6–8]. In this approach, a signal s(t) consisting of several frequency components (f_{ 1 }= k_{ 1 }*f_{ s }, f_{ 2 }= k_{ 2 }*f_{ s }, f_{ 3 }= k_{ 3 }*f_{ s }, etc.) is produced, where each frequency is a multiple of the scanning frequency f_{ s }. Some frequencies prevail during central fixation, while others appear at paracentral fixation. The constants k_{ 1 } , k_{ 2 }, k_{ 3, } etc. depend on the optomechanical design. Figure 1 shows the timefrequency distribution of such signals, produced by a system which generates two signals. Let us assume that the eye does not fixate on the target in the time interval [t_{ a } … t_{ b }], and fixates during [t_{ c } … t_{ d }]. Consequently, in the timefrequency plane, we observe concentration of power at frequency k_{ 1 }*f_{ s } in the timeframe [t_{ a } … t_{ b }], and then shift of the dominating frequency to k_{ 2 }*f_{ s } during [t_{ c } … t_{ d }]. At time t_{ e } the eye loses fixation again, and consequently, the emphasis is shifted back to f_{ 1 }= k_{ 1 }*f_{ s }. In the simplest case, f_{ 2 }= 2f_{ s } is produced during central fixation, while f_{ 1 }= f_{ s } prevails during offcentral fixation, i.e. k_{ 1 }= 1 and k_{ 2 }= 2 [6]. Existing instruments (i.e. [7]), acquire consecutive epochs of s(t), with gaps between them, during which an FFT is performed. A problem with this approach is that the FFT power spectrum is a global measure, i.e., it provides information on how much of f_{ 1 } and f_{ 2 } are represented in the whole epoch analyzed, but it does not provide information on exactly where these frequencies appear and for how long. With lesscooperative pediatric patients, important shortlasting moments of central fixation (f_{ 2 }) may easily be hidden behind large lowfrequency (f_{ 1 }) components, and missed. Analyzing short time intervals is desirable, in order to increase temporal resolution, but this is where the FFT becomes prone to noise and loses spectral resolution. There thus remains a need for improved methods for detecting central fixation.
Method
The method proposed here assumes that an optomechanical apparatus already exists, which enables the acquisition of a birefringence signal from a circular retinal scan, where the scanning circle is positioned around the expected position of the fovea during fixation on a stationary target, appearing optically in the center of the scanning circle. Such systems have been developed and described elsewhere [6–10]. The general design of such a device is shown on Figure 2. In short, the system uses a lowpower laser beam of vertically polarized nearinfrared light, to scan the retina around the presumed center of the fovea and the polarization cross. The returning light is diverted to a polarization sensor, consisting of a polarizing beamsplitter and two photodetectors – one for the vertically polarized component s , and one for the horizontally polarized component p. The difference of the two gives a measure of the polarization state of light, called also S_{ 1 } and being the second component of the full fourcomponent Stokes vector [11]. The signal is amplified and filtered, to remove frequency components outside the region of interest, and finally digitized. The background, comprising the optical noise, is removed prior to any signal processing. This paper focuses on a design producing f_{ 2 }= 2f_{ s } during central fixation and f_{ 1 }= f_{ s } during offcentral fixation (k_{ 1 }= 1, k_{ 2 }= 2) as reported in [7]. Scanning signals were recorded from six subjects after obtaining written informed consent, following a protocol approved by the Institutional Review Board (IRB) and in compliance with the Helsinki Declaration. Personal identification data was not collected, and the signals collected were not used for identification purposes. From each subject, 12 records were acquired, each 400 ms long, for both the right eye (RE) and the left eye (LE). The sampling frequency was 5 kHz, while the scanning frequency was 96 rounds per second.
The problem of detecting shortlasting episodes of central fixation, characterized by the appearance of f_{ 2 }= 2f_{ s }, can be successfully solved by using timefrequency distributions (TFD) obtained by means of the Continuous Wavelet Transform (CWT). It allows excellent localization in both time and frequency domains and is computationally efficient. If the signal to be analyzed is s(t) and g(t) is the analyzing wavelet, the CWT is defined as a convolution of the type:
where * denotes complex conjugate , a  the scale (dilation), and τ  a time shift. The wavelet g(t) and the distribution W(τ,a) are complexvalued [12–17]. The constant 1/ is used for energy normalization. The analyzing wavelet satisfies the following conditions:
where c_{ g } is a constant that depends only on g(t) and a is positive. For an analytic wavelet this constant should be positive and convergent
which in turn imposes an admissibility condition on g(t). For a realvalued wavelet, the integrals from both ∞ to 0 and 0 to +∞ should exist and be greater than zero. Admissible wavelets have no zero frequency contribution, or, what amounts to the same, they are of zero mean, or equivalently G(ω) = 0 for ω=0:

(a)
Belong to L ^{2} (R), i.e. be square integrable (be of finite energy),

(b)
Be analytic (G(ω) = 0 for ω < 0 ) and thus be complexvalued,

(c)
Be admissible. This condition was shown to enable invertibility of the transform:
An appropriate choice for an analyzing wavelet is the admissible complexvalued wavelet of Morlet [12, 18] comprising a modulated Gaussian function of optimal timefrequency concentration:
As an example, Figure 3 shows a surrogate signal s(t) containing two bursts of sinosoidal signals. The first burst is a mixture of two sine waves (25 Hz and 100 Hz) occurring simultaneously, whereas the second burst, shifted in time, contains a sine wave of 50 Hz. Each of the two bursts was modulated by a threeterm BlackmanHarris window:
where f_{ 1 }= 25 Hz, f_{ 2 } =100 Hz, f_{ 3 }= 50 Hz and each modulation window of duration T and starting times respectively t_{1} =0 ms and t_{2}= 600 ms, was calculated according to:
The time domain signal is shown in the upper panel of Figure 3, whereas the FFT of the same signal is shown on the middle panel, and the timefrequency distribution calculated with the CWT (equations 1 and 5) is shown on the bottom panel as a contour plot. Notably, while the FFT does detect all three frequencies correctly in the frequency domain, it does not localize the time moments of the three bursts. The CWT, conversely, does localize all three events in both time and frequency.
Results
For every epoch (T = 0.4 s) of the incoming signal, a timefrequency distribution was computed by means of the CWT (equations 1 and 5). Figure 4, top panel, shows the time domain scan signal s(t) (right eye of a 7 year old boy, f_{ s }= 96 Hz; k_{ 1 }= 1, k_{ 2 }= 2). The subject was trying to fixate on a target. The middle panel shows the FFT power of the same signal. Two main frequency components are present – one at 96 Hz (f_{ 1 }= k_{ 1* }f_{ s }= 96 Hz), characteristic of paracentral fixation, and one at 192 Hz (f_{ 2 }= k_{ 2 }*f_{ s }= 192 Hz), due to central fixation. Based on the FFT plot, one can easily decide that the prevailing frequency is f_{ 1 } and therefore paracentral fixation was prevailing during this record. With the CWT (bottom panel), episodes of f_{ 1 } and f_{ 2 } can easily be localized in time and in frequency dimensions simultaneously. A close look at the appearances of these components gives enough evidence to conclude that there were substantial shortlasting episodes of central fixation, manifested by intermittent appearance of components at f_{ 2 }= 192 Hz. They last between 30 ms and 100 ms each, and show up throughout the epoch. In the first 270 ms they exist sideby side with the paracentral fixation components (f_{ 1 }= 96 Hz), which at times are of higher amplitude. One can easily see why the paracentral components mask the central fixation components in the FFT. It should be mentioned here that sometimes remnants of the optical noise, not entirely eliminated by the background subtraction, appear at the scanning frequency f_{ s }= f_{ 1 }= 96 Hz, and can further mask the central fixation components in the FFT. Likewise, despite the polarizationkilling properties of the human skin in general, some polarized portion of the light gets reflected by the face and returns into the optical system, reaching the sensors and giving rise to additional interference at f_{ s }= f_{ 1 }= 96 Hz. This emphasizes the importance of analyzing the signals in the timefrequency plane, which gives a better chance of separating signal from noise.
In all six test subjects, the CWT allowed precise identification of both frequency components. Moreover, in four of these subjects, episodes of intermittent but definitely present central fixation were detectable, similar to those in Figure 4. A simple FFT is likely to treat them as borderline cases, or entirely miss them, depending on the discrimination thresholds used. In one subject with stable central fixation, the frequency component f_{ 2 }= k_{ 2 }*f_{ s }= 192 Hz was strong and persistent when the subject was fixating on the target, and therefore both the FFT and CWT revealed frequency doubling. But the vast majority of recordings were a constantly changing mixture of the two spectral components. Depending on the time interval analyzed and the duration of central fixation episodes, the FFT was able to detect many episodes, usually the longer lasting ones (longer than 150 ms) and in the absence of noise, whereas the CWT always detected central fixation, even short instants lasting less than 60 ms. Uncompensated portions of the background noise did sometimes affect the f_{ s }= f_{ 1 }= 96 Hz signal, but never the central fixation frequency f_{ 2 }= k_{ 2 }*f_{ s }= 192 Hz.
Discussion
This work is related mainly to methods aiming at detection of strabismus (crosssightedness), eye alignment and amblyopia (“lazy eye”) in young test subjects, where patient cooperation is problematic and the detection of shortlasting episodes is important and even crucial in some cases. The robust application of the CWT makes it applicable to any type of optical systems acquiring and analyzing signals distinguishable in the joint timefrequency domain, especially in systems with abundant optical noise of origin different from that of the signals measured. A condition would be that the noise does not overlap in time and frequency simultaneously with the frequency and temporal location of the event. In addition to amblyopia, the technique described can be used for the detection of frequent and shortlasting losses of fixation, possibly indicative of nystagmus, attentiondeficithyperactivity disorder (ADHD), autism, or other neuropsychologic disorders.
The CWT is not the only technique available for timefrequency expansion. The shorttime Fourier transform, STFT, (also known as the windowed Fourier transform) localizes the signal in time and frequency domain by modulating the time signal with a window function before performing the Fourier transform, to obtain the frequency content of the signal in the region of the window. As a rule, it is a compromise between time and frequency resolution; the wider the window, the higher the frequency resolution, at the cost of poorer time resolution, and vice versa [16]. Any attempt to increase the frequency resolution causes a larger window size and therefore a reduction in time resolution, and viceversa. Also, in order to be able to analyze transients, overlapping windows need to be used, which can slow down analysis considerably, and acts like a lowpass filter in the time domain.
Another alternative, most widely used two or three decades ago, is the WignerVille distribution (WVD). Its definition for timefrequency analysis is:
where i=\sqrt{1} is the imaginary unit, and * denotes complex conjugation [19–22]. In essence, the WVD is the Fourier transform of the input signal’s autocorrelation function, i.e. the Fourier spectrum of the product between the signal and its delayed, time reversed copy, as a function of the delay. Unlike the shorttime Fourier transform, the Wigner distribution function is not a linear transform. A cross term (“time beats”) occurs when there is more than one component in the input signal, analogous in time to frequency beats.
In order to reduce the cross term problem, many other transforms have been proposed, the most prominent one perhaps being the Cohen’s class distribution. The best known member of Cohen’s class distribution function is the Choi–Williams distribution function [22]. This distribution function adopts an exponential kernel to suppress the crossterm:
where
and the kernel function is
However, the kernel gain does not decrease along the η and τ axes in the ambiguity domain, and, consequently, the kernel function of the Choi–Williams distribution function can only filter out the crossterms resulting from components away from the η and τ axes and away from the origin.
In summary, the CWT appears to be the tool of choice when timefrequency distributions are needed in order to detect different frequency components appearing simultaneously, or in different moments in time. This is particularly true of high frequency events of short duration and low amplitude (small scales), or longerlasting lowfrequency oscillations of higher amplitude (large scales), as is the case here.
There are some limitations of the optical and optomechanical hardware involved in this technology. Mechanical vibrations, optical backreflections (reflections from reflective surfaces back to the sensors before actually the light reaches the eyes), multiple internal reflections not captured by the light traps, some birefringence caused by the beam splitters, and others, can all cause instrumental noise which translates into parasitic frequency components. Some of them may overlap with the central fixation frequency of interest and thus mask it, making a differentiation in the joint timefrequency domain difficult. Such artifacts can be avoided or moved away from the region of interest, by careful choice of the mechanical parameters (like scanning speed), optical design (i.e. spinning wave plates), or proper design of the analog electronics (i.e. filters).
Other limitations arise from the measurement principle. Because only about 1/1000 of the polarized light used for measuring is returned by the fovea, the distance between the measurement system and the eye cannot be increased much beyond 40 cm. Further, the room illumination should be dimmed, to prevent pupil constriction and reduction of the nearinfrared light going through it. This not only decreases the signsltonoise ratio, but can cause patient discomfort, and is one more reason to aim for fast and reliable tests.
Conclusion
The Continuous Wavelet Transform is superior to the FFT and most other TFD techniques in localizing fixation frequencies in both the time and frequency domains. It is an excellent tool for precisely identifying central fixation in an uninterrupted manner, thus improving device reliability and shortening test time. When done in a binocular manner, i.e. for both eyes simultaneously, it allows the detection of shortlasting moments of eye alignment or misalignment, or even lack of attention to the target. Using modern digital signal processing hardware like digital signal processors, field programmable gate array (FPGA) logic, for example, the CWT can be performed in real time, and is expected to improve significantly detection sensitivity when testing uncooperative young subjects.
Abbreviations
 FFT:

Fast fourier transform
 CWT:

Continuous wavelet transform
 s(t):

A signal consisting of several frequency components, produced by retinal scanning
 fs:

Scanning frequency
 f1:

,f2, f3: Frequency components, produced by retinal scanning
 k1:

k2, k3: Constants, dependingf on the optomechanical design
 [ta … tb]:

[tc … td]: Time intervals
 RE:

Right eye
 LE:

Left eye
 TFD:

Timefrequency distribution(s)
 g(t):

An analyzing wavelet
 *:

Denotes complex conjugate
 a:

The scale (dilation), used with the continuous wavelet transform
 τ:

A time shift), used with the continuous wavelet transform
 W(τ:

a): A timefrequency distribution, calculated with the wavelet transform (complexvalued)
 L2 (R):

The set of square integrable numbers (being of finite energy)
 cg:

A constant that depends on g(t) and a
 G(ω):

The fourier image of the analyzing wavelet g(t)
 ω:

The circular frequency
 WBH:

A blackmanharris window
 T:

Duration of the blackmanharris window
 ADHD:

Attentiondeficithyperactivity disorder
 FPGA:

Field programmable gate array
 STFT:

The shorttime fourier transform
 WVD:

Wignerville distribution
 Ws(t:

f): A timefrequency distribution, calculated with the wignerville distribution
 Cs(t:

f): A timefrequency distribution, calculated with the Choi–Williams distribution
 Φ(η:

τ): Fernel function of the Choi–Williams distribution.
References
Klein Brink HB, van Blokland GJ: Birefringence of the human foveal area assessed in vivo with Muellermatrix ellipsometry. J Opt Soc Am A 1988, 5(1):49–57. 10.1364/JOSAA.5.000049
Dreher AW, Reiter K, Weinreb RN: Spatially resolved birefringence of the retinal nerve fiber layer assessed with a retinal laser ellipsometer. Appl Opt 1992, 31: 3730–3735. 10.1364/AO.31.003730
Weinreb RN, Shakiba S, Zangwill L: Scanning laser polarimetry to measure the nerve fiber layer of normal and glaucomatous eyes. Am J Ophthalmol 1995, 119(5):627–636.
Elsner AE, Weber A, Cheney MC, Vannasdale DA: Spatial distribution of macular birefringence associated with the Henle fibers. Vision Res 2008, 48(26):2578–2585. 10.1016/j.visres.2008.04.031
Irsch K, Gramatikov B, Wu YK, Guyton D: Modeling and minimizing interference from corneal birefringence in retinal birefringence scanning for foveal fixation detection. Biomed Opt Express 2011, 2(7):1955–1968. 10.1364/BOE.2.001955
Guyton DL, Hunter DG, Patel SN, Sandruck JC, Fry RL: Eye Fixation Monitor and Tracker. U.S. Patent No 6,027,216; 2000.
Hunter DG, Nassif DS, Piskun NV, Winsor R, Gramatikov BI, Guyton DL: Pediatric vision screener 1: instrument design and operation. J Biomed Opt 2004, 9(6):1363–1368. 10.1117/1.1805560
Hunter DG, Patel SN, Guyton DL: Automated detection of foveal fixation by use of retinal birefringence scanning. Appl Optics 1999, 38(7):1273–1279. 10.1364/AO.38.001273
Hunter DG, Sandruck JC, Sau S, Patel SN, Guyton DL: Mathematical modeling of retinal birefringence scanning. J Opt Soc Am A 1999, 16(9):2103–2111. 10.1364/JOSAA.16.002103
Hunter DG, Shah AS, Sau S, Nassif D, Guyton DL: Automated detection of ocular alignment with binocular retinal birefringence scanning. Appl Opt 2003, 42(16):3047–3053. 10.1364/AO.42.003047
Born M, Wiolf E: Principles of Optics. 7th edition. Cambridge: Cambridge University Press; 1999.
KronlandMartinet R, Morlet J, Grossmann A: Analysis of sound patterns through wavelet transforms. Intern J Pattern Recognit Artif Intell 1987, 1: 273–302. 10.1142/S0218001487000205
Rioul O, Vetterli M: Wavelets and Signal Processing. IEEE Signal Processing Magazine 1991, 14–38.
Gramatikov B, Thakor N: Wavelet analysis of coronary artery occlusion related changes in ECG. San Diego, CA: 15th Annual International Conference of the IEEE Engineering in Medicine and Biology Society; 1993:731.
Chui CK Series: Wavelet Analysis and Its Applications, Ser. Vol. 2. In Wavelets: A tutorial in theory and applications. Academic press, Inc; 1992:723.
Gramatikov B, Georgiev I: Wavelets as an alternative to STFT in signalaveraged electrocardiography. Med Biol Eng Comput 1995, 33(3):482–487. 10.1007/BF02510534
Gramatikov B, Brinker J, Yichun S, Thakor NV: Wavelet analysis and timefrequency distributions of the body surface ECG before and after angioplasty. Comput Methods Programs Biomed 2000, 62(2):87–98. 10.1016/S01692607(00)000602
Thakor NV, Gramatikov BI, Sherman D, Bronzino J: Wavelet (TimeScale) Analysis in Biomedical Signal Processing. In The Biomedical Engineering Handbook. Volume 1 2nd edition. Edited by: CRC Press LLC, IEEE Press. 1999, 56–51. 56–27
Claasen TACM, Mecklenbrauker WFG: The wigner distribution  a tool for timefrequency signal analysis. Phillips J Res 1980, 35: 217–250. 276–300; 1067–1072
Cohen L: Wigner distribution for finite duration or bandlimited signals and limited cases. IEEE Trans Acoustics, Speech, Signal Processing 1987, 35: 796–806. 10.1109/TASSP.1987.1165201
Velez EF, Garudadri H: Speech analysis based on smoothed WignerVille distribution. In TimeFrequency Signal Analysis  Methods and Applications. Edited by: Boashash B. New York, Toronto, Chichester: Halsted Press and John Wiley & Sons; 1992:351–374.
Qian S, Chen D: Joint TimeFrequency Analysis  Methods and Applications. Upper Saddle River, NJ: Prentice Hall; 1996.
Acknowledgements
This work was supported by an Individual Biomedical Research Award from The Hartwell Foundation, April 2010March 2013.
Author information
Authors and Affiliations
Corresponding author
Additional information
Competing interests
The author declares that he has no competing interests.
Authors’ contribution
The device with which the signals were measured was built several years ago by a team of researchers, of which the author was a part, being responsible for the optoelectronics (lasers and sensors), analog and digital electronics, software for data acquisition, analysis, user interface, archival etc. This device has been extensively evaluated, and the results have been published in different publications, as seen from the references. The timefrequency analysis, subject of this paper, was developed and tested by BG alone. All authors read and approved the final manuscript.
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.
Rights and permissions
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
About this article
Cite this article
Gramatikov, B. Detecting fixation on a target using timefrequency distributions of a retinal birefringence scanning signal. BioMed Eng OnLine 12, 41 (2013). https://doi.org/10.1186/1475925X1241
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/1475925X1241