Phase retrieval in optical fiber modal interference for structural health monitoring

Optics Communications 283 (2010) 1278–1284

Contents lists available at ScienceDirect

Optics Communications journal homepage: www.elsevier.com/locate/optcom

Phase retrieval in optical fiber modal interference for structural health monitoring S.K. Ghorai a,*, Soumya Sidhishwari a, S. Konar b a b

Department of Electronics and Communication Engineering, Birla Institute of Technology, Mesra, Ranchi 835215, Jharkhand, India Department of Applied Physics, Birla Institute of Technology, Mesra, Ranchi 835215, Jharkhand, India

a r t i c l e

i n f o

Article history: Received 25 May 2009 Received in revised form 12 November 2009 Accepted 4 December 2009

Keywords: Phase retrieval Data dependent system Modal interference Strain sensor

a b s t r a c t A method based on data dependent system (DDS) for extraction of phase in fiber modal interference is presented. The interference patterns of LP01 & LP11, LP01 & LP02 and LP06 &LP07 within the fiber have been recorded under different launching conditions. The patterns were characterized by means of autoregressive model and the self coherence functions of the corresponding interferogram were determined. It would provide the phase distribution of the pattern and the modulation of group delay due to the measurand. An application has been made for measuring strain in a simply supported beam under different loading conditions. Results are presented for the applied strain in the range of 270–1500 l strain. Ó 2009 Elsevier B.V. All rights reserved.

1. Introduction In recent years, there is a growing demand for health monitoring of civil structures as well as aircraft structures. Nondestructive evaluation of the health of engineering structures is the critical need for the deterioration of infrastructure elements. Structural health monitoring is an essential management tool for safely working of advanced structures. Advanced composite and concrete structures are now being widely used in modern vehicles, ships, civil infrastructures and aerospace industry. In those structures, in situ structural monitors are highly desirable to detect a decrease in performance or imminent failure due to the variation of relevant physical parameters such as strain, temperature, corrosion, vibration etc. A typical health monitoring system is composed of a network of sensors that measure the parameters relevant to the state of the structure and its environment. Optical fiber sensors will be ideally suited for the long term continuous structural health monitoring systems. These sensors can be developed to selectively detect the variation of physical parameters. In advanced structures, the location and extent of strain are the important information in order to understand their behavior under loading condition. Several different types of fiber optic sensors have been reported for measuring strain in different structures and beams. Fiber Bragg grating sensors have been considered as promising tools for measuring strain in composite structures and beams [1–3]. Distributed fiber sensors based on Brillouin Scattering have been the focus of * Corresponding author. Tel.: +91 651 2275750; fax: +91 651 2275401. E-mail addresses: [email protected] (S.K. Ghorai), [email protected] (S. Sidhishwari). 0030-4018/$ - see front matter Ó 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.optcom.2009.12.011

great attention for measuring strain distribution in large structures [4–6]. However in all these sensors, complexity arises during implementation and also the measurement error increases for large strain difference. Interferometric sensors are very attractive because of their high resolution and sensitivity. Different types of interferometric sensors based on Mach–Zehnder, Michelson, Fabry–Perot, Sagnac etc. have been reported for measuring strain in different structures [7–9]. But in these two beam interferometers, the reference and sensing arms are placed at different locations which can lead to error in measurement. Also the accuracy of measurement depends on the extraction of phase from the interferogram. In the present work, we develop two-mode fiber interferometer using birefringent fiber for measuring strain in a simply supported beam. Modal interference in optical fibers has been reported by several workers for its application in sensors [10–12]. In interference phenomenon, the determination of phase has been an active area of research, as it carries the key information of the physical parameters. Over the years, different methods have been demonstrated for phase retrieval from interferogram. Phase shifting techniques [13,14], Fourier transform method [15], Synchronous method [16], Quasi one frame algorithm method [17], Regularization technique with low pass filtering [18] were reported to determine phase from the fringe pattern. In recent years, more advanced techniques such as windowed Fourier transform method [19], parameter estimation method [20] and wavelet transform method [21] have been proposed for phase extraction. For higher phase resolution at a high temporal bandwidth, a technique based on sinusoidally modulated phase shifting interferometry has been reported [22]. The above phase retrieval methods may be either

S.K. Ghorai et al. / Optics Communications 283 (2010) 1278–1284

temporal or spatial. In the temporal techniques, three or more phase steps are required for reducing the errors. In the spatial techniques, a single interference pattern is used, but it requires a large number of fringes in the pattern. All of these techniques require global processing to find the phase at a particular point i.e. full length of data set is utilized for phase determination. For example, in FTM method, by Fourier transform, the point wise information is globally distributed over the whole pattern. In our work, we have used the data dependent system (DDS) method for phase extraction from modal interference pattern. The DDS methodology analyses discrete data to develop a statistically adequate mathematical model for a real system. There is no information loss between the data and the model, if the model is adequate. The method utilizes the correlation between the neighbouring data points of a profile. The model takes into account the dynamic dependence between the neighbouring data points that represent a real system. The dynamic dependence is expressed by a regression model [23]. The method has been widely used in manufacturing and process design [24].The method has also been used for phase extraction from noisy interferogram [25]. In our earlier work, we have used DDS method for extraction of phase in a Mach–Zehnder type fiber interferometer [26]. In the present context, we extract phase from the beat signals produced due to the interference of LP01 & LP11, LP01 & LP02 and LP06 &LP07 modes in single mode birefringent fiber and multimode fiber under loading conditions in a beam. Higher order modes have more coupling coefficients than lower ones [11]. A set of data has been selected from the interferogram and characterized by means of autoregressive (AR) model. The autoregressive parameters determine the self coherence function of the pattern and it would provide the phase information due to change in strain. The group delay obtained from the phase contour is utilized to determine the strain of the beam under loading conditions.

1279

change in propagation constants of LP modes and Dh = h1  h2. Whenever there is a change in external disturbance, cos (Db  Dh) will change due to the variation of DbL, and accordingly intensity pattern will change, otherwise it remains constant. A small perturbation due to strain at one point of the fiber causes a coupling of light to the other mode. The phase delay due to the different mode velocities leads to a beat frequency. The phase difference between the two modes may be written as

DU ¼ ðb1  b2 ÞL;

ð2Þ

where L is the distance of light propagation along the fiber. Under weakly guiding approximation, the propagation constant bi (i = 1,2) can be expressed in terms of normalized propagation constant bi as

  2 b2i ¼ k0 n22 þ bi n21  n22 ;

ð3Þ

where k0, n1 and n2 are the free space wave vector, core and cladding refractive index, respectively. Assuming the fiber as homogeneous material in elasticity, the refractive index variation under strain can be written as [12]:

oni n3i ½p  rðp11 þ p12 Þ; ¼ oL 2L 12

ð4Þ

where p11, p12 are Strain Optic Coefficients (for fused silica 0.12 and 0.27, respectively), r is the Poisson’s ratio (0.17). In our work, we consider a simply supported beam where a concentrated load is applied at the center of the beam. If the fiber is attached at the distance ‘Y’ from the neutral axis in the direction of the beam length, then the strain distribution within elastic limit is given by

WY Z; 4EI

2. Modal interference

SðZÞ ¼

Intermodal interference within a single mode fiber is based on the principle of modal Mach–Zehnder interferometer. In single mode fiber, although fundamental mode is guided, higher order modes can be propagated for wavelength above the cut off. In a single mode fiber, if the ‘V’ parameter is set in the range 2.408– 3.8317, then LP01 & LP11 modes are excited and in the range 3.8317–5.1356, LP01 & LP02 modes are excited with the launching of axially symmetric light beam at the input end. The fundamental mode LP01 is two fold degenerate and LP11 mode is four fold degenerate with x or y polarizations. The modes that can propagate withe;y o;y o;x in a two-mode fiber are LPx01 ; LPy01 ; LPe;x 11 ; LP11 ; LP11 and LP11 . In low birefringence fiber, both LP01 mode and LP11 mode have their same propagation constants in x and y polarizations. However, in high birefringence fiber, degeneracy is lifted and different polarization modes will have different propagation constants [27]. Similar phenomenon occurs in LP02 mode, but in conventional communication fibers, propagation constants are equal in x and y polarization. In the present analysis, intermodal beating has been considered, e.g. LPx01 and LPx11 or LPy01 and LPy11 . Higher order modes e.g. LP06 & LP07 are excited by splicing a single mode fiber with a multimode fiber with no axial offset between them. The intensity pattern due to the interference of two LP modes may be written as

where E and I are the elastic modulus and the moment of inertia with respect to the neutral axis, W is the load applied. As the strain is symmetric with respect to the center of the beam, the Eq. (5) provides the strain distribution for one half of the beam, maximum at the center while zero at the fixed end. When the strain is applied to the beams, mode coupling occurs and the phase difference between the two modes leads to a change in interferogram pattern. The interferogram is captured by a CCD camera and the intensity distribution in one dimension can be expressed as

Iðr; /Þ ¼ A21 f12 ðrÞ þ A22 f22 ðrÞ cos2 / þ 2A1 A2 f1 ðrÞf2 ðrÞ  cos / cosðDbL  DhÞ;

ð1Þ

where A1 and A2 are the amplitude coefficients of the two LP modes, and f1 & f2 are their radial distribution functions. In the third term of Eq. (1), cos (DbL  Dh) contributes towards the change in phase due to external perturbations, where Db = b1  b2, represents the

ð5Þ

IðxÞ ¼ 2I0 ðxÞ þ CðxÞ þ C ðxÞ;

ð6Þ

where C(x) denotes the self coherence function and C*(x) is its conjugate, I0(x) is the background intensity and x is the pixel number. The argument of this coherence function will provide the required phase information due to measurand (i.e. strain, temperature, pressure etc.). When the phase is expanded around a centre optical frequency x0, it may be expressed as,

 1 /ðxÞ ¼ L bðx0 Þ þ b1 ðx0 Þðx  x0 Þ þ b2 ðx0 Þðx  x0 Þ2 2  1 3 þ b3 ðx0 Þðx  x0 Þ þ    6  m    where bm ¼ ddxmb 

x¼x0

ð7Þ

; m ¼ 0; 1; 2; 3; . . . represent the imbalances in

the propagation constant b and its derivatives. Using DDS method, the phase is recovered from the interferogram and the values of the coefficients of Eq. (7) are determined. They would provide the values of the measurand producing the modulation of group delay.

1280

S.K. Ghorai et al. / Optics Communications 283 (2010) 1278–1284

BIREFRINGENT FIBER POLARIZER LASER

POLARIZER

CM

W CP PC

Fig. 1. Experimental setup for recording interference pattern in a simply supported beam. CM – camera; CP – optical coupler; PC – computer.

Fig. 2. Intensity patterns of LP01  LP11 modal interference under different loading conditions.

Fig. 4. Intensity patterns of LP06  LP07 modal interference under different loading conditions.

tion x in terms of its predecessors (Ix1, Ix2, . . .) and the current stimulus ax by autoregressive model of order n, AR (n) [28]:

Ix ¼ /1 Ix1 þ /2 Ix2 þ /3 Ix3 þ    þ /x Ixn þ ax ;

ð8Þ

where /i are the autoregressive parameters, and ax is Gaussian homogeneous noise. The linear least square algorithm has been used to estimate /i of the AR model. The AR (n) model can be arranged in a transfer function form and written in terms of the Green’s function Gj, where Gj has the form:

Gj ¼ g 1 kj1 þ    þ g n kjn ; . . . ;

ð9Þ

where ki are the roots of the autoregressive polynomial given by Fig. 3. Intensity patterns of LP01  LP02 modal interference under different loading conditions.

kn  /1 kn1  /2 kn2      /n ¼ 0; and the weighting factor of each

3. Phase extraction In our work, we have used the DDS method to extract phase from the interferogram. The DDS methodology analyzes discrete data to derive a mathematical model of a real system. It exploits the correlation between the neighboring data points of a profile. The intensity pattern obtained from the interferometer may be modulated to get the dependence of a current pixel value Ix at loca-

gi ¼

kn1 i ; Pðki  kj Þ

kji

ð10Þ is given by

ð11Þ

1 6 j – i 6 n: The Green’s function characterizes the deterministic components of the process [28]. The convolution of Gj with the residuals ax produces the real data Ix, obtained during experiment. For each of the residuals ax, there will be responses like Gjaxj and the sum

S.K. Ghorai et al. / Optics Communications 283 (2010) 1278–1284



ImðCx Þ UðxÞ ¼ tan1 : ReðCx Þ

Storing an image of size m x n in a file.

Representing each character by the gray level value of a particular pixel of the image.

4. Experimental setup

Storing all the pixel gray level values of the image in a 2-D array.

Displaying some specific pixel gray level values, which were stored in the array. Fig. 5. Flowchart of image processing algorithm.

of all these responses yields the real data Ix. With the help of Green’s function, the mathematical equivalence between the DDS and the interferogram expression can be established as follows: x X

Gj axj ;

ð12Þ

j¼0

where n(x) is the noise present in the interferogram. When the variance component of the carrier frequency dominates, the self coherence function can be written as

Cx ¼

x X j¼0

Gj axj ¼

x X

g i kji axj

ð14Þ

This phase is wrapped within the range p to p and it is then unwrapped with suitable algorithm. In our program, we have used Matlab’s ‘unwrap’ function which is based on the assumption that the phase jumps by more than p radians must have been wrapped i.e. multiples of 2p are added or subtracted to each phase input to restore original phase values. For a matrix it operates column wise. A program may be developed to calculate the Taylor coefficients of the unwrapped phase and these would provide the required group delay and dispersion due to the measurand.

Reading the file information character by character.

Ix ¼ 2Io ðxÞ þ CðxÞ þ C ðxÞ þ nðxÞ ¼

1281

ð13Þ

j¼0

The phase of the interferogram can be recovered from the argument of the coherence function Cx:

Experiments have been performed for the verification of the proposed method. The experimental setup used for recording interference pattern is shown in Fig. 1 for an acrylic beam. Single mode birefringent fibers (HB-800 and HB-1300, Oxford Electronics) were attached to a simply supported beam (Acrylic rod). The dimension of the beam was 51.0  3.0  0.6 cm for LP01 & LP11, and also for LP01 & LP02 modal interferences. For LP06 & LP07 modal interference the beam dimension was taken as 38.8  3.8  0.5 cm and a portion of multimode fiber (F-MLD) was spliced between two single mode fibers. Optical couplers were used at both ends of the fiber and at the input end, light from a He–Ne laser (k = 630 nm) was coupled to the birefringent fiber through the polarizer. The loads were applied at the centre of the beam. For LP01 & LP11 modes, the fiber taken was characterized by cut off wavelength of 800 nm. To excite the LP01 mode, the laser beam was focused at the input end of the fiber. At the operating wavelength of 630 nm, the normalized frequency ‘V’ of the fiber becomes 3.32. Thus the fiber can support only two lower modes LP01 & LP11. Initially the fundamental mode LP01 was excited by a suitable launching condition. When the loads are applied to the beam, power exchange takes place between the LP modes. For LP01 & LP02 modes excitation, the fiber was characterized by cut off wavelength of 1300 nm. He–Ne laser was replaced by a laser diode (Ando-Electric, AQ-1318, 850 nm). At the operating wavelength of 850 nm, ’V’ becomes 4.4308, and fiber supports four modes. By introducing proper longitudinal offset, both LP01 and LP02 were excited. The output end of the fiber was imaged/observed on a nearby

Fig. 6. Phase distribution obtained after removing carrier frequency for intensity patterns of LP01  LP11, for different loads.

1282

S.K. Ghorai et al. / Optics Communications 283 (2010) 1278–1284

Fig. 7. Phase distribution obtained after removing carrier frequency for intensity patterns of LP01  LP02 for loads: (a) 300 l strain, (b) 1200 l strain (DDS method).

Fig. 8. Phase distribution obtained after removing carrier frequency for intensity patterns of LP06  LP07 for loads: (a) 270 l strain, (b) 1100 l strain (DDS method).

difference between the two modes changes the interference patterns. The patterns were recorded using CCD camera and the corresponding phase fields were extracted using DDS technique, which provides the information of the measurand. 5. Results and discussions In our experiment, three different types of fiber modal interferences (LP01 & LP11, LP01 & LP02 and LP06 &LP07) were recorded using CCD camera. Figs. 2–4 show the intensity patterns recorded for the simply supported beam under different loading conditions. An image processing algorithm has been developed to digitize the interferogram data and thus to obtain the gray level intensity value for each pixel. The flowchart of the algorithm is shown in Fig. 5. An autoregressive model has been selected of the order (15) i.e. AR (15), n = 15, using the coding method described in article 3. For one dimension analysis, the current pixel value Ix may be written from Eq. (8) as [28]: Fig. 9. Unwrapped phase of Fig. 3a with the carrier frequency removed.

screen/on a computer and it enabled the identification of lobe patterns of LP modes. When load is applied to the beam, the phase

Ix ¼ /1 Ix1 þ /2 Ix2 þ /3 Ix3 þ    /15 Ix15 þ ax ;

ð15Þ

for x = 16, 17, . . ., m Using fifteen pixel intensities, the roots of the autoregressive polynomial were determined, and those were used to obtain the

S.K. Ghorai et al. / Optics Communications 283 (2010) 1278–1284

self coherence function of the interferogram which would provide the phase field. Figs. 6–8 show the phase distribution for the corresponding loads. The unwrapped phase obtained from the intensity pattern of Fig. 3a is shown in Fig. 9 after carrier frequency removal. Figs. 10, 11 show the phase distribution obtained from Fourier transform method. The variations of group delay for different values of relative strain are shown in Fig. 12. Results are compared with those obtained by the Fourier transform method. It may be seen that more accurate phase extraction has been obtained using DDS method as compared to Fourier transform method. This is due to the fact that the noise present in the interferogram is also being transformed due to the use of FFT which is a part of the Fourier transform method and it would remain in the frequency domain as before. On the contrary, in DDS method, embedded noises arising from different sources or measurements can be separated out with the help of Green’s function’s decomposition and it is possible to isolate the noise free coherence term of the interference pattern. The DDS method requires minimum data processing skill and real time processing is possible with proper interfacing. However, inaccurate model selection can lead to errors in phase extraction. The modal interference pattern changes due to

1283

phase variation under loading condition. It may be seen that better contrast of the fringe pattern can be obtained using higher modal interference because of higher coupling than the lower order modes and hence it will be more useful to determine the strain of the beam. The sensitivity of our sensor system depends on the term (DbL). The birefringent fiber used in our experiment has the value of Db(=2pn1D/k*ob) = 0.0339 rad lm1 at 630 nm, for LP01 & LP11 modal interference (values of the parameters were as follows: r.i. of the core n1 = 1.464, fractional index difference D = 0.0055, normalized propagation constant for LP01 mode b0 = 0.6970 and normalized propagation constant for LP11 mode b1 = 0.2727), for ‘V’ no. = 3.32. Fig. 13 shows the phase variation for different values of strain. It may be seen a good agreement between experimental  results and theoretical one. The strain sensitivity 1L oðoDeUÞ for LP01  LP11 modal interference was 0.052 rad le1m1, and for LP01  LP02 it was 0.059 rad le1 m1. With this fiber, we have obtained strain resolution 20 l strain which is less than FBG or LPG and greater than Brillouin Scattering based strain sensor [29,30]. The proposed sensor system can be applied to monitor strain in large structures with no need of high resolution. In FBG, strain

Fig. 10. Phase distribution obtained by Fourier transform method for intensity patterns of LP01  LP11 for different loads.

Fig. 11. Phase distribution obtained by Fourier transform method for intensity patterns of LP01  LP02 for loads: (a) 300 l strain, (b) 1200 l strain.

1284

S.K. Ghorai et al. / Optics Communications 283 (2010) 1278–1284

the phase variation using DDS method. The phase distribution would provide the group delay and with proper calibration, the value of the measurand can be determined for structural health monitoring. 6. Conclusion

Fig. 12. Variation of group delay with strain for LP01  LP11modal interference.

We have demonstrated the data dependent system method for extraction of phase in fiber modal interferometers. The method has been used to determine strain in a simply supported beam. Measurement sensitivity would depend upon the mode coupling within the birefringent fiber as well as on appropriate model selection. Interference of higher order modes provides better contrast of fringe pattern due to higher mode coupling than lower order modes. The present technique requires minimum data processing skill and would provide real time processing with proper interfacing. The method provides better result in presence of high noise and insensitive to errors associated with other techniques e.g. FFT which is used in Fourier transform method. Another advantage is that it does not require a full length of data set for evaluating the phase at a particular point, only certain number of data is required depending upon the model. The method would be very much useful to determine the physical parameters in structural health monitoring. Acknowledgements The authors deeply acknowledge the institute sponsored project for carrying out the research work. The research work was partially supported by ISRO (VSSC) under RESPOND scheme. References

Fig. 13. Phase variations vs. strain for LP01  LP11 modal interference.

measurement depends on the Bragg wavelength shift which is generally measured in costlier equipments like optical spectrum analyzer (OSA) etc. and complex arrangement is required for interrogation of FBGs. Our sensor system can be used for both distributed and quasi-distributed measurement system, whereas FBGs can be used only in quasi-distributed sensor system i.e. in predetermined points. The phase variation  also results due to temperature. The temDUÞ comes out to be 0.21 rad C1 m1 for perature sensitivity 1L oðoT LP01  LP11 modal interference. However, the sensitivity has been observed nearly constant for the range of strains in the experiment.It can be reduced by adjusting the thermal expansion coefficient and thermo-optic coefficient, or by selecting the wavelength of light when wavelength dispersion is minimum. In other way, the present sensor system can be used as a temperature sensor by selecting a fiber of large residual thermal stress. Efforts are also being made to record the modal interference pattern due to temperature variation and subsequently to extract

[1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30]

H.Y. Ling, K.T. Lau, L. Cheng, W. Jin, Measurement 39 (2006) 328. V.M. Murukeshan, P.Y. Chan, L.S. Ong, A. Asundi, Appl. Opt. 40 (2001) 145. H.Y. Ling, K.T. Lau, L. Cheng, Compos. Struct. 66 (2004) 317. K. Hotate, L.S.S. Ong, IEEE Photonic Technol. Lett. 15 (2003) 272. H. Naruse, M. Tateda, H. Ohno, A. Shimada, Appl. Opt. 41 (2002) 7212. M. DeMerchant, A. Brown, X. Bao, T. Bremner, Appl. Opt. 38 (1999) 2755. L.B. Yuan, L.M. Zhou, W. Jin, Opt. Lett. 25 (2000) 1074. L. Yuan, J. Yang, Sensor Actuat. A 105 (2003) 40. S.C. Huang, W.W. Lin, M.T. Tsai, M.H. Chen, Sensor Actuat. A 135 (2007) 570. T.J. Chen, Opt. Commun. 261 (2006) 43. E. Li, X. Wang, C. Zhang, Appl. Phys. Lett. 89 (2006) 091119-1. D. Kumar, S. Sengupta, S.K. Ghorai, Meas. Sci. Technol. 19 (2008) 065201-8. K. Creath, in: E. Wolf XXVI (Ed.), Progress in Optics, Elsevier, Amsterdam, 1988, p. 350. J.A. Ferrari, E.M. Frins, W. Dultz, J. Lightwave Technol. 15 (1997) 968. X. Su, W. Chen, Opt. Laser Eng. 35 (2001) 263. K.H. Womack, Opt. Eng. 23 (1984) 391. M.J. Huang, Opt. Commun. 210 (2002) 187. M. Servin, J.L. Marroquin, F.J. Cuevas, J. Opt. Soc. Am. A 18 (2001) 689. Q. Kemao, Opt. Laser Eng. 45 (2007) 304. A. Patil, R. Langoju, P. Rastogi, Opt. Laser Eng. 45 (2007) 293. L.R. Watkins, Opt. Laser Eng. 45 (2007) 298. K. Falaggis, D.P. Towers, C.E. Towers, J. Opt. A: Pure Appl. Opt. 11 (2009) 054008-11. S.M. Pandit, N. Jordache, G.A. Joshi, J. Opt. Soc. Am. A 11 (1994) 2584. G. Subhash, J.E. Loukus, S.M. Pandit, Mech. Mater. 34 (2002) 25. S.M. Pandit, N. Jordache, Appl. Opt. 34 (1995) 5945. S.K. Ghorai, S. Narayan, Opt. Laser Technol. 36 (2004) 571. T.A. Eftimov, W.J. Bock, J. Lightwave Technol. 11 (1993) 1925. S.M. Pandit, D.P. Chan, Appl. Opt. 38 (1999) 4095. K.T.V. Grattan, Dr. T. Sun, Sensor Actuat. A 82 (2000) 40. O. Frazao, L.A. Ferreira, F.M. Araujo, J.L. Santos, Fiber Integrated Opt. 24 (2005) 227.