Characterization of Laser Propagation Through Turbulent Media by Quantifiers Based on the Wavelet Transform^{*}^{*}footnotemark: *
Abstract
The propagation of a laser beam through turbulent media is modeled as a fractional Brownian motion (fBm). Time series corresponding to the center position of the laser spot (coordinates and ) after traveling across air in turbulent motion, with different strength, are analyzed by the wavelet theory. Two quantifiers are calculated, the Hurst exponent, , and the mean Normalized Total Wavelet Entropy, . It is verified that both quantifiers give complementary information about the turbulence state.
Departamento de Física,
Facultad de Ciencias Exactas,
Universidad Nacional de La Plata (UNLP) y Centro de Investigaciones Ópticas, CC. 124 Correo Central, 1900 La Plata, Argentina.
Instituto de Cálculo,
Facultad de Ciencias Exactas y Naturales,
Universidad de Buenos Aires (UBA).
Pabellón II, Ciudad Universitaria.
1428 Ciudad de Buenos Aires, Argentina.
1 Introduction
Waveletsbased tools have been shown to be wellsuited to fractal processes and their analysis. This is mainly due to the fact that the wavelet transform incorporates in its definition two basic features, time and scale, which are of primary importance for fractal processes. Another remarks concerns the structure of the wavelet transform which, by construction, builds a signal by successive refinements, starting from a coarse approximation and adding finer and finer details at each step. Such a procedure is of course reminiscent of basic fractal constructions, thus suggesting we should make use of wavelets to synthesize a fractal process [1, 2, 3].
Since the earlier leading work of Mandelbrot[4, 5] there is a lot of evidence that several facets of fully developed turbulent flows are fractals. Most of the applications of fractals to turbulence have been devoted to the study of subsets of the region occupied by turbulent flow where a given property is satisfied. For example, the turbulent/nonturbulent interface, some constant property surfaces (such as the isovelocity and isoconcentration surfaces) and the structure of spatial distribution of dissipation rates have been characterized as fractals or multifractals and the dimension has been measured [6, 7, 8]. Then, it is natural to apply wavelet tools for analyzing atmospheric turbulence [9, 10, 12, 11].
The fractional Brownian motion (fBm) was introduced by Mandelbrot and Van Ness [13] as an example of a process which contains an infinite domain of dependence with the intention of explaining the results reported by Hurst in 1951 [14]. fBm as a model for turbulence is not new[12, 11, 15]. In particular, the fBm can be introduced as an alternative model for the turbulent index of refraction, and it can be shown that these processes allow to reconstruct most of the index properties[15].
This work reports on the basic characteristics of the centroid position of a laser spot after the light has traveled through turbulent media. It is analyzed as a stochastic process, within a mathematically defined theoretical model: the fractional Brownian motion.
That is, given a component of the centroid position as a time function, the geometrical nature of the graph is studied. Afterwards, the wavelet theory is used to characterize this centroid. Two quantifiers are obtained: the Hurst exponent, , and the mean Normalized Total Wavelet Entropy, . Their behaviors are compared; the analysis shows they describe different properties of the turbulence.
2 Fractional Brownian motion
The fractional Brownian motion of exponent (Hurst exponent), with , is a zeromean Gaussian process with such that:

.

The difference has finite dimensional normal distribution.
This means that their increments are stationary Gaussian processes whose variance is
(1) 
where is a parameter— is, of course, the average. The nonstationary character of the fBm is evidenced by its covariance function given by
(2) 
The covariance of future increments with past ones is:
(3)  
Note that is independet of and the parameter determine the correlation of the increments.
For the correlation of past and future increments vanishes for all , as is required for an independent random process. However, for one has . For this quantity is positive and the process is called persistent[17]. In this case, an increasing trend in the past implies on the average an increasing trend in the future and, conversely, a decreasing trend in the past implies on the average a continued decrease in the future. For the process is called antipersistent. Now, an increasing trend in the past implies on the average a decreasing trend in the future, while a decreasing trend in the past makes on the average an increasing trend in the future—see Fig. 1.
K. Helland and Van Atta[18] were the first to study the Hurst exponent in grid generated turbulence by applying a rescaledrange analysis—an usual measure of longterm persistence in geophysical time series—to turbulence velocity records measured within a wind tunel. They showed that there are some deviations from —referred as the Hurst phenomenon.
As a nonstationary process, the fBm does not have a spectrum defined in the usual sense; however, it is possible to define an empirical power spectrum of the form:
(4) 
This equation is not a valid power spectrum in the theory of stationary processes since it is a nonintegrable function, but it could be considered as a generalized spectrum. Through this interpretation of Eq. (4) a selfsimilarity relation can be shown for the fBm. That is for with and parameters, one has that have the same finite dimensional distributions for all and . The fractal dimension of sample functions of these processes is .
3 TimeFrequency Analysis
3.1 Wavelet transform
First introduced by Dennis Gabor [19], wavelet analysis is a method which relies on the introduction of an appropriate basis and a characterization of the signal by the distribution of amplitude in this basis. If the basis is required to be a proper orthogonal basis, any arbitrary function can be uniquely decomposed and the decomposition can be inverted [20, 21, 22, 23]. Wavelet analysis is a suitable tool for detecting and characterizing specific phenomena in time and frequency planes.
The wavelet is a simple and quickly vanishing (compactly supported) oscillating function. Unlike sine and cosine of Fourier analysis, which are precisely localized in frequency but extend infinitely in time, wavelets are relatively localized in both time and frequency. Furthemore, wavelets are bandlimited; they are composed of not one but a relatively limited range of several frequencies.
A wavelet family is the set of elementary functions generated by dilations and translations of a unique admissible mother wavelet :
(5) 
where , are the scale and translation parameters respectively, and is the time. As increases, the wavelet becomes narrower. Thus, one have a unique analytic pattern and its replications at different scales and with variable time localization.
The continuous wavelet transform (CWT) of a signal (the space of real square summable functions) is defined as the correlation between the function with the family wavelet for each and :
(6) 
For special election of the mother wavelet function and for the discrete set of parameters, and , with (the set of integers) the family
(7) 
constitutes an orthonormal basis of the Hilbert space consisting of finiteenergy signals.
The correlated decimated discrete wavelet transform (DWT) provides a nonredundant representation of the signal, and the values constitute the coefficients in a wavelet series. These wavelet coefficients provide relevant information in a simple way and a direct estimation of local energies at the different scales. Moreover, the information can be organized in a hierarchical scheme of nested subspaces called multiresolution analysis in . In the present work, we employ orthogonal cubic spline functions as mother wavelets. Among several alternatives, cubic spline functions are symmetric and combine in a suitable proportion smoothness with numerical advantages. They have become a recommendable tool for representing natural signals [24, 25]—figure 2 shows the cubic spline wavelet function used in this work.
In what follows, the signal is assumed to be given by the sampled values , corresponding to an uniform time grid with sampling time . For simplicity, the sampling rate is taken as . If the decomposition is carried out over all resolutions levels the wavelet expansion is:
(8) 
where and the wavelet coefficients can be interpreted as the local residual errors between successive signal approximations at scales and , and is the residual signal at scale . It contains the information of the signal corresponding to the frequencies with the sampling frequency.
3.2 Relative Wavelet Energy
Since the family is an orthonormal basis for , the concept of energy is linked with the usual notions derived from Fourier’s theory. The wavelet coefficients are given by and the corresponding asociated energy will be its square. The energy at each resolution , will be
(9) 
where represents the number of wavelet coefficients at resolution . The total energy can be obtained in the fashion
(10) 
Finally, we define the normalized values, which represent the Relative Wavelet Energy (RWE) by
(11) 
for the resolution levels . The yield, at different scales, the probability distribution for the energy. Clearly, and the distribution can be considered as a timescale density that constitutes a suitable tool for detecting and characterizing specific phenomena in both the time and the frequency planes.
3.3 Wavelet entropy
The Shannon entropy[26] gives a useful criterion for analyzing and comparing probability distribution. It provides a measure of the information contained in any distribution. We define the Normalized Total Wavelet Entropy (NTWS) [27, 28, 29] as
(12) 
where
(13) 
The NTWS appears to be a measure of the degree of orderdisorder of the signal. It provides useful information about the underlying dynamical process associated with the signal[29]. Indeed, a very ordered process can be represented by a periodic monofrequency signal (signal with a narrow band spectrum). A wavelet representation of such a signal will be resolved at one unique wavelet resolution level, i. e., all RWE will be (almost) zero except at the wavelet resolution level which includes the representative signal frequency. For this special level the RWE will be (in our chosen energy units) almost equal to one. As a consequence, the NTWS will acquire a very small, vanishing value. A signal generated by a totally random process or chaotic one can be taken as representative of a very disordered behavior. This kind of signal will have a wavelet representation with significant contributions coming from all frequency bands. Moreover, one could expect that all contributions will be of the same order. Consequently, the RWE will be almost equal at all resolutions levels, and the NTWS will acquire its maximum possible value.
Figure 3 presents two different relative wavelet energy (probability) distribution corresponding to five wavelet resolution levels (). It is clear from the figure that distribution presents broad band spectrum. In contrast, distribution shows a clear dominance of the resolution level . According to the description given above, for the NTWS the following relations can be expected: NTWS NTWS. This is observed in the numerical values given at the figure.
3.4 Wavelet quantifiers time evolution
In order to follow the temporal evolution of the quantifiers defined above, RWE and NTWS, the analyzed signal is divided into nonoverlapping temporal windows with length (, with ). Afterwards, appropriate signalvalues for these quantifiers are assigned to the middle point of each time window. In the case of a diadic wavelet decomposition, the number of wavelet coefficients at resolution level is two times smaller than at the previous, , one. The minimum length of the temporal window will therefore include at least one wavelet coefficient at each level.
The wavelet energy at resolution level for the time window is given by
(14) 
where represents the number of wavelet coefficients at resolution level corresponding to the time window ; while the total energy in this time window will be
(15) 
The time evolution of RWE and NTWS will be given by:
(16) 
(17) 
In order to obtain a quantifier for the whole time period under analysis [28] the temporal average is evaluated. The temporal average of NTWS is given by
(18) 
and for the wavelet energy at resolution level
(19) 
then the total wavelet energy temporal average is defined as
(20) 
In consequence, a mean probability distribution representative for the whole time interval (the complete signal) can be defined as
(21) 
with and the corresponding mean NTWS as
(22) 
4 Fractional Brownian motion and Wavelet Transform
A relevant property of the wavelet based multiresolution analysis is the stationary character of the wavelet coeficient series corresponding to each level of resolution[30]. Another important property is that the reconstruction of the original time series from the stationary series of wavelet coefficients reproduces the original signal with small error.
In relation with fractional Brownian motion it can be shown that:

fBm exhibits positive longrange correlation in the range but wavelet coefficients have a correlation which can be arbitrarily reduced;
All the above properties evidence that wavelet analysis is naturally wellsuited to fBm. Each of them provides in fact a key ingredient for a problem of major importance when analyzing fBm: the estimation of the Hurst exponent or of the related spectral exponent [3, 30]. Starting from the above mentioned observations that, in the wavelet transform, variance progression follows a power law across scales
(24) 
one can simply make use of the empirical variance estimators at scale (based on coefficients for a sample of total length )
(25) 
This is made possible because the fBm, when decomposed via the wavelet transform, becomes stationary at each scale.
In this way an estimator of the parameter can be obtained by: a) estimating the variance of the wavelet coefficients with Eq. (25); b) plotting versus and fitting a minimum square line. The slope of the line give the estimator of .
Waveletbased estimators dedicated to fBm can be viewed as versatile generalizations of previous techniques. An important feature of fBm is that its increments are stationary and such that the structure function is proportional to —see Eq. (1)—thus suggesting that this variance should be estimated in order to find . Although feasible, this approach is faced with a difficulty due to longrange dependence. Classical (empirical) variance estimators are obviously poor estimators in such a context and specific estimators have to be designed[33, 34].
5 Experimental Setup and Data Adquisition
Time series corresponding to the fluctuations in the position of a laser beam’s spot (wandering) over a screen, after propagation through laboratory generated convective turbulence, were recorded with a position sensitive detector located as screen at the end of the path. Twenty records of 5000 spot beam coordinates measurements every five minutes were obtained. Temperature along the laser beam path, for each record, were also measured and stored.
The experimental measures were performed in the laboratory by producing convective turbulence over a path of lenght m with two electric heaters in a row and covering the path. We use a mW continuous wave HeNe laser (Melles Griot Model 05LHP991), wavelenght laser beam nm (red), with pointing stability, after 15 minutes of warmup, less than m and a beam divergence of mrad.
Each electric heater was cm long and W of power. The height of the laser beam’s propagation path was m above the electrical heaters and a box of expanded polyurethane was used as a thermal protection for the equipment. Three thermometers were allocated on the top of the box, cm above the laser beam. They were used to determine the temperature in three positions along the propagation path, see Fig. 4. The geometrical experimental arrangement was similar to that used by Consortini et al. [35] and, Consortini and O’Donnell[36].
The position sensor has a relative accuracy of the order of m so that very small position fluctuation can be measured. It was interfaced to a computer which allowed to measure at a rate of about samples per second ( laser spot coordinates in approximately s). Thus, with these coordinates stored on a hard disk the Hurst exponent and the mean NTWS were computed off line.
Three different intensity levels of convective turbulence were generated by changing the amount of heat dissipated for each electrical heater:

Normal turbulence, the electrical heaters were off;

Soft turbulence, each electrical heaters dissipated half of its available power;

Hard turbulence, each electrical heaters dissipated at the maximum available power.
Figure 5 shows the temperature for the three turbulence levels over the twenty records. Note that there are temperature gradients over the twenty records for the three turbulence levels and that the temperature (, and , see Fig. 4) is not homogeneously distributed. So, an uniform flux of warm air and uniformity of turbulence along the path was not present.
6 Results and discussion
We transformed the twenty time records corresponding to the laser beam and coodinates (), under the three turbulence conditions, to the timefrequency domain by means of orthogonal discrete wavelet transform (ODWT)[20, 21, 22, 23] obtaining in this way the corresponding wavelet coefficients series. In the present work, we consider eight wavelet resolution levels and cubic spline mother wavelet[27, 28, 29]. Note that the wavelet coefficients were nonoverlapping for each scale.
The wavelet coefficients were squared to obtain the associted wavelet energy and the estimator of wavelet coefficient variance . To evaluate the Hurst exponent the procedure described at the end of Sec. 4 is used. Figure 6 presents versus the resolution level , for the coordinates and for the three turbulence conditions. The vertical lines represent the scaling region used for the square fitting. The values of the coefficient obtained were averaged over the twenty records at each turbulence condition for the coordinates respectively.
On the other hand, each signal under analysis was divided in time windows of length and for each one the probability wavelet energy distribution was evaluated. With these values the mean wavelet energy distribution and the mean NTWS, , were obtained. This quantifier is taken as representative of the orderdisorder of the whole signal. As before, the obtained values were averaged over the twenty records at each turbulence condition for the coordinates respectively.
A comparison between average Hurst exponent and mean NTWS is given in Fig. 7. It is known the entropy is a measure of the order of a given system. In this case, the mean NTWS shows the same behavior for both coordinates. As the turbulence increases the mean NTWS does the same. But, the Hurst exponent discriminates between coordinates: for the axis no noticeable change is observed, and for the axis the Hurst exponent decreases (increasing the roughness, see Fig. 1) with the increasing turbulence.
Thus, it is observed that the Hurst exponent is sensitive to the mean flow of the warm air. It distinguishes the anisotropy characteristic of the convective turbulence. Because the entropy measures the order obviously it does not detect the mean flow.
In any case, the value of the obtained quantities indicates that memory as well as selfsimilarity and scale invariance are significant property of these time series.
In relation with the Hurst exponent a useful generalization consists in allowing the singularity exponent to become timedependent , thus generating a new process such that
(26) 
In such a situation the increment process is no longer stationary. Then, it is impossible to globally apply the technique mentioned above for estimating . Provided that variations of are smooth enough, the time series can be divided in time windows where this requirement is satisfied and for each one the similar techniques based on timescale energy distributions can be applied locally[37]. In any case, we must emphasize that for each one an scaling region must be defined. The time evolution of NTWS could be easily implemented[28, 29]. Moreover, the NTWS is capable of detecting changes in a nonstationary signal due to the localization characteristic of the wavelet transform, the computational time is significantly shorter since the algorithm involves the use of wavelet transform in a multiresolution framework and, the NTWS is parameterfree and scaling region is not necessary for its evaluation.
In the next future the intention of our project is to make experiments varying the intensity of turbulence by adjusting the voltage applied to the heating element using a voltage controller, in order to study the quantifiers temporal evolution and characterize in a quantitative way the dynamics of the process. In another hand, the same study will be doing in outdoor experiments at different moments of the day to characterize the ground atmospheric turbulence.
References
 [1] H. Aubert and D. L. Jaggard, “Fractal superlattices and their wavelet analyses,” Opt. Commun. 149, 207–212 (1998).
 [2] M. E. Degaudenzi and C. M. Arizmendi, “Waveletbased fractal analysis of electrical power demand,” Fractals 8, 3, 239–245 (2000).
 [3] G. Frantziskonis and A. Hansen, “Waveletbased multiscaling in selfaffine random media,” Fractals 8, 4, 403–411 (2000).
 [4] B. B. Mandelbrot, “Intermittent turbulence in selfsimilar cascades: Divergence of high moments and dimension of the carrier,” J. Fluid Mech. 62, 331–358 (1974).
 [5] B. B. Mandelbrot, “On the geometry of homogeneous turbulence, with stress on the fractal dimension of the isosurfaces of scalars,” J. Fluid Mech. 72, 401–416 (1975).
 [6] K. Sreenivasan and C. Meneveau, “The fractal facets of turbulence,” J. Fluid Mech. 173, 357–386 (1986).
 [7] C. Meneveau and K. Sreenivasan, “The multifractal nature of the turbulent energy dissipation,” J. Fluid Mech. 224, 429–484 (1991).
 [8] A. Scotti, C. Meneveau and S. Saddoughi, “Fractal dimension of velocity signals in highReynoldsnumber hydrodynamic turbulence,” Phys. Rev. E 51, 6, 5594–5608 (1995).
 [9] L. Hudgins, C. Friehe and M. Mayer, “Wavelet transforms and atmospheric turbulence,” Phys. Rev. Lett. 71, 20, 3279–3282 (1993).
 [10] Qian Minwei and Liu Shida, “Scale exponents of atmospheric turbulence under different stratifications,” Fractals 4, 1, 45–48 (1996).
 [11] G. Katul, B. Vidakovic and J. Albertson, “Estimating global and local scaling exponents in turbulent flows using wavelet transformations,” Phys. Fluids 13, 241–250 (2001).
 [12] G. Papanicolaou and K. Sølna, “Wavelet based estimation of local Kolmogorov turbulence,” in Longrange dependence. Theory and Applications, edited by P. Doukhan, G. Oppenheim and M. S. Taqqu. Birkhauser (2002), pp. 473–505.
 [13] B. B. Mandelbrot and J. Van Ness, “Fractional Brownian motion, fractional noises and applications,” SIAM Rev. 10, 4, 422–437 (1968).
 [14] H. Hurst, “Longterm storage capacity of reservoirs,” Trans. Am. Soc. Civil Eng. 116, 770–808 (1951).
 [15] D. G. Pérez, Propagación de Luz en Medios Turbulentos, (in Spanish) Tesis de la Universidad Nacional de La Plata, Argentina, http://arxiv.org/abs/physics/0307144 (2003).
 [16] J. LévyVehel, “FracLab: A fractal analysis toolbox for signal and image processing,” http://fractales.inria.fr/, (2001).
 [17] J. Feder, Fractals (Plenum Press, New York, 1988), pp. 170–171.
 [18] K. Helland and C. W. Van Atta, “The ‘Hurst phenomenon’ in grid turbulence,” J. Fluid Mech. 85, 573–589 (1978).
 [19] D. Gabor, “Theory of communication,” J. IEE (London) 93, 429–457 (1946).
 [20] I. Daubechies, Ten Lectures on Wavelets, (SIAM, Philadelphia, 1992).
 [21] A. Aldroubi and M. Unser (Eds), Wavelets in Medicine and Biology, (CRC Press, Boca Raton, 1996).
 [22] S. Mallat, A Wavelet Tour of Signal Processing, 2nd. Edition, (Academic Press, San Diego, 1999).
 [23] V. Samar, A. Bopardikar, R. Rao and K. Swartz, “Wavelet analysis of neuroelectric waveforms: a conceptual tutorial,” Brain and Language 66, 7–60 (1999).
 [24] P. Thévenaz, T. Blue and M. Unser, “Interpolation revisited,” IEEE Trans. on Medical Imaging 19, 739–758 (2000).
 [25] M. Unser, “Spline: a perfect fit for signal and image processing, ” IEEE Signal Processing Magazine 16, 22–38 (1999).
 [26] C. E. Shannon, “A mathematical theory of communication,” Bell System Technol. J. 27, 379–423; 623–656 (1948).
 [27] S. Blanco, A. Figliola, R. Quian Quiroga, O. A. Rosso and E. Serrano. “Timefrequency analysis of electroencephalogram series (III): Wavelet Packets and Information Cost Function,” Phys. Rev. E 57, 932–940 (1998).
 [28] O. A. Rosso, S. Blanco, J. Yordanova, V. Kolev, A. Figliola, M. Schürmann, and E. Başar, “Wavelet entropy: a new tool for analysis of short duration brain electrical signals, ” J. Neuroscience Method 105, 65–75 (2001).
 [29] O. A. Rosso and M. L. Mairal. “Characterization of time dynamical evolution of electroencephalographic epileptic records,” Physica A 312, 469–504 (2002).
 [30] A. Peréz, C. E. D’Attellis, M. Rapacioli, G. A. Hirchoren and V. Flores, “Analyzing blood cell concentration as a stochastic process,” IEEE Ingineering in Med. and biology, November/December, 170–175 (2001).
 [31] P. Flandrin, “On the spectrum of fractional Brownian motions,” IEEE Trans. Inf. Theory IT35, 1, 197–199 (1989).
 [32] P. Flandrin, “Wavelet analysis and synthesis of fractional Brownian motion,” IEEE Trans. Inf. Theory, IT38, 2, 910–917 (1992).
 [33] J. Beran, Statistics for longmemory processes (Chapman and Hall, New York, 1994), pp. 81–120.
 [34] M. Taqqu, V. Teverovsky and W. Willinger, “Estimators for longrange dependence: an empirical study,” Fractals 3, 4, 785–798 (1995).
 [35] A. Consortini, Sun Yi Yi, Li Zhi Ping and G. Conforti, “A mixed method for measuring the inner scale of atmospheric turbulence,” J. Modern Optics 37, 10, 15551560 (1990).
 [36] A. Consortini and K. O’Donnell, “Beam wandering of thin parallel beams through atmospheric turbulence,” Waves in Random Media 3, S11S28 (1991).
 [37] R. F. Peltier and J. LévyVehel, Multifractional Brownian motion: definition and preliminary results, (Research Report No. RR2645, INRIA, 1995).