The Continuous Wavelet Transform in the detection of short lasting moments of central fixation

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 (cross-sightedness), 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₁=k₁*f_s, f₂=k₂*f_s, f₃=k₃*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 para-central fixation.  The constants k₁ , k₂,  k_3,  etc. depend on the opto-mechanical design. Figure 1 shows the time-frequency 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 time-frequency plane, we observe concentration of power at frequency k₁*f_s  in the timeframe [t_a … t_b], and then shift of the dominating frequency to k₂*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₁=k₁*f_s. In the simplest case, f₂=2f_s is produced during central fixation, while f₁=f_s  prevails during off-central fixation, i.e. k₁=1 and  k₂=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₁ and  f₂ are represented in the whole epoch analyzed, but it does not provide information on exactly where these frequencies appear and for how long.  With less-cooperative pediatric patients, important short-lasting moments of central fixation (f₂) may easily be hidden behind large low-frequency (f₁) 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.

Figure 1.   Time-frequency distribution of signals, produced by a retinal birefringence scanning system. Central fixation is characterized by the appearance of frequency component k₂*f_{s .}

Method

The method proposed here assumes that an opto-mechanical 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 low-power laser beam of vertically polarized near-infrared 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₁ and being the second component of the full four-component 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₂=2f_s  during central fixation and  f₁=f_s  during off-central fixation (k₁=1, k₂=2) as reported in [7].  Data from six subjects, properly consented, were acquired (12 records, each 400 ms long, both for the right eye (RE) and the left eye (LE). The sampling frequency was 5 kHz, while the scanning frequency was 96 rounds per second.

Figure 2.   A simplified block diagram of a device for retinal birefringence scanning. The system measures changes in the polarization state of light caused by the retina.

The problem of detecting short-lasting episodes of central fixation, characterized by the appearance of f₂=2f_s, can be successfully solved by using time-frequency 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:

                                                  (1)

where  *  denotes complex conjugate ,  a - the scale (dilation),  and  t  - a time shift.   The wavelet g(t) and the distribution W(t,a) are complex-valued. [12] [13] [14] [15] [16] [17]  The constant 1/ is used for energy normalization. The analyzing wavelet satisfies the following conditions: 

(a) Belong to L² (R), i.e. be square integrable (be of finite energy),

(b) Be analytic (G(w) = 0   for w < 0 ) and thus be complex-valued,

(c) Be admissible. This condition was shown to enable invertibility of the transform:

                                                                                                                    (2)

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

      (3)

which in turn imposes an admissibility condition on g(t).  For a real-valued 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(w )=0 for w=0:

(4)

An appropriate choice for an analyzing wavelet is the admissible complex-valued wavelet of Morlet [12] [18] comprising a modulated Gaussian function of optimal time-frequency concentration:

                                                         (5)

Figure 3.   Analysis of 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.  FFT and CWT are performed on the same time signal.

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 three-term Blackman-Harris window:

               (6)

where f₁ = 25 Hz,  f₂ =100 Hz,  f₃ = 50 Hz and each modulation window of duration T and starting times respectively t₁ =0 ms and t₂=600 ms, was calculated according to:

for  0 < t < T              (7)

 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 time-frequency 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

Figure 4.  Time-frequency distribution of a real patient data (right eye of a 7 year old boy). The top panel shows the time domain scan signal s(t). The middle panel shows the FFT power of the same signal. The CWT (bottom panel) reveals short-lasting episodes of central fixation, with f_{2 =} 192 Hz.

For every epoch (T=0.4s) of the incoming signal, a time-frequency 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, k₂=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₁ =  k_1* f_s = 96 Hz), characteristic of para-central fixation, and one at 192 Hz (f₂=k₂*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 para-central fixation was prevailing during this record.  With the CWT (bottom panel), episodes of   f₁  and f₂  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 short-lasting 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 side-by side with the para-central fixation components (f_{1 =} 96 Hz), which at times are of higher amplitude. One can easily see why the para-central 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₁ = 96 Hz, and can further mask the central fixation components in the FFT. Likewise, despite the polarization-killing 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₁ = 96 Hz. This emphasizes the importance of analyzing the signals in the time-frequency 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₂=k₂*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₁ = 96 Hz signal, but never the central fixation frequency f₂=k₂*f_s = 192 Hz.

This work is related mainly to methods aiming at detection of strabismus (cross-sightedness), eye alignment and amblyopia (“lazy eye”) in young test subjects, where patient cooperation is problematic and the detection of short-lasting 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 time-frequency 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 short-lasting losses of fixation, possibly indicative of nystagmus, attention-deficit-hyperactivity disorder (ADHD), autism, or other neuropsychologic disorders. 

Figure 5. Time-Frequency distributions vs. Fast Fourier Transform. The former allows localization in time, while the latter gives a generalized information, pertaining to the whole epoch.

The CWT is not the only technique available for time-frequency expansion. The short-time 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 vice-versa.  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 low-pass filter in the time domain.

Another alternative, most widely used two or three decades ago, is the Wigner-Ville distribution (WVD). Its definition for time-frequency analysis is:

                                    (8)

where i= 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 short-time 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 cross-term:

                                         (9)

where

                                               (10)

and the kernel function is

                                                          (11)

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 cross-terms 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 time-frequency 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 longer-lasting low-frequency oscillations of higher amplitude (large scales), as is the case here.

 There are some limitations of the optical and opto-mechanical hardware involved in this technology. Mechanical vibrations, optical back-reflections (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 time-frequency 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 near-infrared light going through it. This not only decreases the signal-to-noise ratio, but can cause patient discomfort, and is one more reason to aim for fast and reliable tests.

In 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 short-lasting 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.

References

1.            klein Brink HB, van Blokland GJ: Birefringence of the human foveal area assessed in vivo with Mueller-matrix ellipsometry. J Opt Soc Am A 1988, 5(1):49-57.

2.            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.

3.            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.

4.            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.

5.            Irsch K, Gramatikov B, Wu YK, Guyton D: Modeling and minimizing interference from corneal birefringence in retinal birefringence scanning for foveal fixation detection. Biomedical optics express 2011, 2(7):1955-1968.

6.            Guyton DL, Hunter DG, Patel SN, Sandruck JC, Fry RL: Eye Fixation Monitor and Tracker. In US Patent No 6,027,216. 2000.

7.            Hunter DG, Nassif DS, Piskun NV, Winsor R, Gramatikov BI, Guyton DL: Pediatric Vision Screener 1: instrument design and operation. Journal of biomedical optics 2004, 9(6):1363-1368.

8.            Hunter DG, Patel SN, Guyton DL: Automated detection of foveal fixation by use of retinal birefringence scanning. Appl Optics 1999, 38(7):1273-1279.

9.            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.          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.

11.          Born M, Wiolf E: Principles of Optics. 7th (expanded) edn. Cambridge: Cambridge University Press; 1999.

12.          Kronland-Martinet R, Morlet J, Grossmann A: Analysis of sound patterns through wavelet  transforms. Intern J of Pattern Recognition and Artificial Intelligence 1987, 1:273-302.

13.          Rioul O, Vetterli M: Wavelets and Signal Processing. IEEE Signal Processing         Magazine 1991(October):14 -38.

14.          Gramatikov B, Thakor N: Wavelet analysis of coronary artery occlusion related changes in ECG. In 15th Annual International Conference of the IEEE Engineering in Medicine and Biology Society; San Diego, CA. 1993:731.

15.          Chui CK: Wavelets: A tutorial in theory and applications. Academic press.; 1992.

16.          Gramatikov B, Georgiev I: Wavelets as an Alternative to STFT in signal-averaged electrocardiography. Medical and Biological Engineering and Computing 1995, 33(3):482-487.

17.          Gramatikov B, Brinker J, Yi-chun S, Thakor NV: Wavelet analysis and time-frequency distributions of the body surface ECG before and after angioplasty. Comput Methods Programs Biomed 2000, 62(2):87-98.

18.          Thakor NV, Gramatikov BI, Sherman D: Wavelet (Time-Scale) Analysis in Biomedical Signal Processing. In The Biomedical Engineering Handbook. Volume 1. 2 edition. Edited by Bronzino J: CRC Press LLC, IEEE Press; 1999: 56-51 - 56-27.

19.          Claasen TACM, Mecklenbrauker WFG: The Wigner Distribution - A tool for time-frequency signal analysis. Phillips J Res 1980, 35:217-250; 276-300; 1067-1072.

20.          Cohen L: Wigner distribution for finite duration or band-limited signals and limited cases. IEEE Trans Acoustics, Speech, Signal Processing 1987, 35:796-806.

21.          Velez EF, Garudadri H: Speech analysis based on smoothed Wigner-Ville distribution. In Time-Frequency Signal Analysis - Methods and Applications. Edited by Boashash B. New York, Toronto, Chichester: Halsted Press and John Wiley & Sons; 1992: 351-374.

22.          Qian S, Chen D: Joint Time-Frequency Analysis - Methods and Applications. Upper Saddle River, NJ: Prentice Hall; 1996.

See also:

· Gramatikov, B. Detecting fixation on a target using time-frequency distributions of a retinal birefringence scanning signal. BioMedical Engineering OnLine 2013, 12:4. Published on May 13, 2013. doi:10.1186/1475-925X-12-41

http://www.biomedical-engineering-online.com/content/12/1/41.

· Nitish V. Thakor, Boris Gramatikov and David Sherman. Wavelet (Time-Scale) Analysis in Biomedical Signal Processing. In The Biomedical Engineering Handbook, Ed. Joseph Bronzino, 2^nd Edition., Volume I, Section VI (Biomedical Signal Analysis), Chapter 56, pp. 56-1 – 56-27. CRC Press LLC, IEEE Press, December 1999

https://books.google.co.in/books/about/Biomedical_Engineering_Handbook.html?id=ijhP6BZxCMAC

http://www.crcnetbase.com/doi/abs/10.1201/9781420049510.ch56 (1999)

http://www.crcnetbase.com/doi/abs/10.1201/9781420003864.ch5 (2006)

· Irsch,K., Gramatikov,B.I., Wu, Y. K., Guyton,D. L. A new pediatric vision screener employing polarization-modulated, retinal-birefringence-scanning-based strabismus detection and bull’s eye focus detection with an improved target system: Opto-mechanical design and operation. Journal of Biomedical Optics, 2014, Jun 1;19(6):67004. doi: 10.1117/1.JBO.19.6.067004.

http://biomedicaloptics.spiedigitallibrary.org/article.aspx?articleid=1881172