Applying the Hilbert–Huang Decomposition to Horizontal Light Propagation C2 data n 1 2 Mark P. J. L. Chang, Erick A. Roura, Carlos O. Font , Charmaine Gilbreath and Eun Oh 1 Physics Department, University of Puerto Rico, Mayagu¨ez, Puerto Rico 00680 2 U.S. Naval Research Laboratory, Washington D.C. 20375 ABSTRACT The Hilbert Huang Transform is a new technique for the analysis of non–stationary signals. It comprises two distinct parts: Empirical Mode Decomposition (EMD) and the Hilbert Transform of each of the modes found from the first step to produce a Hilbert Spectrum. The EMD is an adaptive decomposition of the data, which results in the extraction of Intrinsic Mode Functions (IMFs). We discuss the application of the EMD to the calibrationoftwoopticalscintillometersthathavebeenusedtomeasureC2 overhorizontalpathsonabuilding n rooftop, and discuss the advantage of using the Marginal Hilbert Spectrum over the traditional Fourier Power Spectrum. Keywords: Empirical Mode Decomposition, Hilbert Transform, Strength of Turbulence, Scintillation 1. INTRODUCTION ThecommonpracticewhenstudyingtimeseriesdataistoinvokethetoolsofFourierspectralanalysis. Although extremely versatile and simple, the technique suffers from some stiff contraints that limit its usefulness when attempting to examine the effects of optical turbulence in the frequency domain. Namely, the system must be linear and the data must be strictly periodic or stationary. Strict stationarity is a constraintthat is impossible to satisfy simply on practical grounds, since no detector can cover all possible points in phase space. The linearity requirement is also not generally fulfilled, since turbulent processes are by definition non–linear. Fortunately a new technique that has come to be known as the Hilbert Huang Transform (HHT) has been developed,1 patented by NASA. This allows for the frequency space analysis of non–stationary, non–linear signals. The HHT is composed of two main algorithms for filtering and analyzing such data series. Firstly it employsanadaptive technique to decompose the signalinto a number ofIntrinsic Mode Functions (IMFs) that have well prescribedinstantaneous frequencies, defined as the first derivative of the phase of ananalytic signal. The second step is to convert these IMFs into an energy–time–frequency relationship, by means of the Hilbert Transform. Asides from overcomingthe problems associatedwith more traditional Fourier methods, the HHT makes it possibletovisualizetheenergyspreadbetweenavailablefrequencieslocallyintime,ratherlikewavelettransform methods. The advantage the HHT has over wavelet transforms is that it is of much higher resolution, since it does not a priori assume a basis; rather it ”lets the data do the talking”. 2. INSTANTANEOUS FREQUENCY Key to the HHT is the idea of instantaneous frequency, which we will sometimes refer to as simply ”the frequency”. The ideal instantaneous frequency is quite simply the frequency of the signal at a single time point. No knowledge is required of the signal at other times. Naturally such a statement leads to difficulties in definition; Huang et al2 take it to be the derivative of the phase of the analytic signal, found from the real and imaginary parts of the signal’s Hilbert Transform, which we follow. Furtherauthor information: (Sendcorrespondence toM.P.J.L.C.) M.P.J.L.C.: E-mail: [email protected],Telephone: 1 787 265 3844 Advances in Stellar Interferometry, edited by John D. Monnier, Markus Schöller, William C. Danchi, Proc. of SPIE Vol. 6268, 62683E, (2006) · 0277-786X/06/$15 · doi: 10.1117/12.672280 Proc. of SPIE Vol. 6268 62683E-1 Report Documentation Page Form Approved OMB No. 0704-0188 Public reporting burden for the collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection of information, including suggestions for reducing this burden, to Washington Headquarters Services, Directorate for Information Operations and Reports, 1215 Jefferson Davis Highway, Suite 1204, Arlington VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to a penalty for failing to comply with a collection of information if it does not display a currently valid OMB control number. 1. REPORT DATE 3. DATES COVERED 2006 2. REPORT TYPE 00-00-2006 to 00-00-2006 4. TITLE AND SUBTITLE 5a. CONTRACT NUMBER Applying the Hilbert-Huang Decomposition to Horizontal Light 5b. GRANT NUMBER Propagation C2n data 5c. PROGRAM ELEMENT NUMBER 6. AUTHOR(S) 5d. PROJECT NUMBER 5e. TASK NUMBER 5f. WORK UNIT NUMBER 7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 8. PERFORMING ORGANIZATION Naval Research Laboratory,Code 8123, Advanced Systems Technology REPORT NUMBER Branch,4555 Overlook Ave SW,Washington,DC,20375 9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSOR/MONITOR’S ACRONYM(S) 11. SPONSOR/MONITOR’S REPORT NUMBER(S) 12. DISTRIBUTION/AVAILABILITY STATEMENT Approved for public release; distribution unlimited 13. SUPPLEMENTARY NOTES 14. ABSTRACT The Hilbert Huang Transform is a new technique for the analysis of non?stationary signals. It comprises two distinct parts: Empirical Mode Decomposition (EMD) and the Hilbert Transform of each of the modes found from the first step to produce a Hilbert Spectrum. The EMD is an adaptive decomposition of the data, which results in the extraction of Intrinsic Mode Functions (IMFs). We discuss the application of the EMD to the calibration of two optical scintillometers that have been used to measure C2n over horizontal paths on a building rooftop, and discuss the advantage of using the Marginal Hilbert Spectrum over the traditional Fourier Power Spectrum. 15. SUBJECT TERMS 16. SECURITY CLASSIFICATION OF: 17. LIMITATION OF 18. NUMBER 19a. NAME OF ABSTRACT OF PAGES RESPONSIBLE PERSON a. REPORT b. ABSTRACT c. THIS PAGE Same as 8 unclassified unclassified unclassified Report (SAR) Standard Form 298 (Rev. 8-98) Prescribed by ANSI Std Z39-18 Theimmediateproblemindealingwithaphasesodefinedisthat,forthemostpart,theHilbertTransforms ofthe directsignalsarenotwellbehavedresultinginnegativeinstantaneousfrequencieswhichdo notrepresent physicaleffects. The methodbywhichthis iscircumventedistoensurethatthe inputtotheHilbertTransform obeys the following conditions: (a) The number of local extrema of the input and the number of its zero crossings must be either equal or differ at most by one. (b) At any point in time t, the mean value of the upper envelope (determined by the local maxima) and the lower envelope (determined by the local minima) is zero. The functions that obey these are considered the IMFs. 3. EMPIRICAL MODE DECOMPOSITION We have implemented an IMF filtering algorithm, known as Empirical Mode Decomposition (EMD), following Huang et al.1,3 The IMFs and the residual trend line thus obtained are verified to be complete by simply summing them to recreate the signal. The maximum relative error we have found is of the order 10−9 %. Empirical Mode Decomposition al n g si 1 mf i 2 mf i 3 mf i 4 mf i 5 mf i 6 mf i 7 mf i 8 mf i 9 mf i0 1 mf i1 1 mf i2 1 mf is. e r Figure 1. TheIMFsfoundfromatypicalinputsignal(takenon9-March-2006), withthesignalitselfshownatthetop. For convenience, we refer to the lowest order IMF as one with the fastest oscillation. The bottom–most graph is the residual after removing all theIMFs and represents theoverall trend. The IMFs show that in a very real sense the EMD method is acting as a filter bank, separating the more rapid oscillations from the slower oscillations. It seems that a subset of the individual IMFs may be added to determine the effect of physical variables, as suggested in Figure 2. In the absence of the major effects of Proc. of SPIE Vol. 6268 62683E-2 Addition of IMFs from slowest to fastest oscillations ginal gnal orisi 1 2 3 4 5 6 7 8 9 0 1 1 1 2 1 al ot T C2 vs 24−hour time x 10−13 n 14 C2 n IMFs+trend 12 10 8 2Cn6 4 2 0 −2 0 5 10 15 20 25 24−hour time Figure 2. (a) Stepwise summation of the 9-March-2006 IMFs and trend line to recreate the original input signal. (b) The sum of thetrendline and theslowest 3 IMFs superimposed on theinput signal. The overshoot into negative values of Cn2 is unphysical, and serves to demonstrate that the information content of the subset is incomplete. Nevertheless, thefit does suggest that theIMFs represent an underlyingphysical process (primarily solar insolation). Proc. of SPIE Vol. 6268 62683E-3 solar insolation, the HHT technique reveals that the majority contribution to the C2 signal lies in the highest n order (slowest oscillation) IMFs, as can be seen from the extremely faithful fit to the data composed of the trend line and the 3 highestorderIMFs shownin Figure 3. As a guide to the significanceof the variousmodes, x 10−14 Cn2 vs 24−hour time − Instrument A x 10−13 Cn2 vs 24−hour time 9 1.4 C2 C2 n n IMFs + trend IMFs + trend 8 1.2 7 1 6 0.8 2Cn5 2Cn 0.6 4 0.4 3 2 0.2 1 0 16 17 18 19 20 21 22 23 24 16 17 18 19 20 21 22 23 24 24−hour time 24−hour time Figure 3. The20-Feb-2006Cn2 andfitcomposed ofthetrendlineand3highestorderIMFsforbothinstruments(Aon theleft and B on theright). Sunset was at 18:31 local time. we examined the energy in the IMFs and compared the energy in each mode to the energy distribution of red noise. As aninitial na¨ıveestimator,we take IMF1 to be representativeof the noise containedin the data. This is unlikely to be completely correct; probably a better noise estimator would be IMF1 - <IMF1>, where the angle brackets signify the mean value overthe same temporalepoch (e.g. month or season). We do not do this simply because we do not have sufficient data. We define the red noise (random) time series to be an AR1 process r(tn)=σE(tn)+ρr(tn−1) (1) The terms are: σ := standard deviation of IMF1 (2) E := uniform distribution of random numbers between 1 and -1 tn := the nth timestep ρ := the autocorrelationat 1 lag step of IMF1 r := the random time series The IMFs of an ensemble of AR1 time series are generated and then this Monte Carlo is used to simulate the power distribution of the noise. The power of each C2 derived IMF is then compared to the noise power n distribution to determine the mode’s significance. Figure 4 shows that of the IMFs, the first three lie at or belowthemedianrednoisepower. IMF4toIMF9lieabovethemediannoisepower,withthehigherorderIMFs being most significant. We interpret this simulationto mean that IMF6-IMF9 are highly physically significant. 4. HILBERT TRANSFORM OF IMFS FollowingthedecompositionintoIMFsoftheoriginalsignal,thederivedcomponentscanbeHilberttransformed to produce a time–frequency map or spectrum. Figure 5 shows that during the hours when there is no solar insolation the C2 energy is distributed in the lowest frequencies. The highest frequencies sampled are only n reached after the Sun is contributing energy into the lower atmosphere. The discontinuous, filamentary aspect of the plot indicates a large number of phase dropouts which shows that the data are non–stationary. Proc. of SPIE Vol. 6268 62683E-4 Average Power per IMF 10−15 10−16 s nit u y ar bitr10−17 ar er, w o P 10−18 10−19 1 2 3 4 5 6 7 8 9 IMF order Figure 4. MonteCarlo simulations showthesignificance of the Cn2 (solid line) derivedIMFs for asingle instrument on theafternoon/evening of 20-Feb-2006, compared to theequivalent red noise power (open circles). Hilbert−Huang spectrum 0.05 0.045 0.04 0.035 0.03 Hz) ncy (0.025 e u q e fr 0.02 0.015 0.01 0.005 2 4 6 8 10 12 14 16 18 20 22 time Figure 5. The Hilbert spectrumof the9-March-2006 IMFs shown in Figure 1 plotted as a series of contour lines. Proc. of SPIE Vol. 6268 62683E-5 WearealsoabletofindaMarginalSpectrumbyintegratingtheHilbertspectrumacrosstime. TheMarginal SpectrumshowninFigure6clearlysuffersfromlessleakageintothehighfrequenciesthanthePowerSpectrum. The interpretation of both spectra are quite different: the Fourier Power Spectrum indicates that certain frequencies exist throughout the entire signal with a given squaredamplitude. The MarginalSpectrum, on the other hand, describes the probability that a frequency exists at some local time point in the signal. It is clear Marginal Hilbert Spectrum vs Fourier Power Spectrum 10−8 Marginal Spectrum Power Spectrum 10−9 10−10 m) u10−11 ctr e p S al 10−12 n gi ar 2C (Mn10−13 10−14 10−15 10−16 10−4 10−3 10−2 10−1 Frequency Figure6.TheHilbertMarginalspectrum(solidline)ofthe9-March-2006IMFscomparedtotheFourierPowerSpectrum (dottedline)ofthesamedata. ThePowerSpectrumhasbeenshiftedsothatitsmaximumfrequencycoincideswiththat of the Marginal Spectrum. Note also that theunits of the ordinate axis are arbitrary for the Fourier Power Spectrum. thatthetwospectrahaveadifferenceinthestandarddeviations: thelogarithmoftheMarginalSpectrumhasa standarddeviationof1.88comparedtothelogarithmoftheFourierPowerSpectrum,whosestandarddeviation is 2.24. ApplyingaKolmogorov–SmirnovtesttothetwospectrareturnsaP–valueof0andacumulativedistribution function distance of 1, indicating that the data sets represent different distributions, as we might guess from their different gradients. We may therefore state that the spectra are unrelated and we argue on the basis of non–stationarity that the Fourier Power Spectrum has little, if any, physical meaning. 5. INSTRUMENTS AND DATA REDUCTION The C2 data usedinthis study wascollectedduring 2006atthe UniversityofPuerto Rico,Mayagu¨ezCampus, n on the rooftop of the Physics Department. The data were obtained with two commercially available scintillometers (model LOA-004) from Optical Scientific Inc, co–located such that the transmitter of system 1 was next to the receiver of system 2. Each LOA-004 instrument comprises of a single modulated infrared transmitter whose output is detected by two single pixel detectors. For these data, the separation between transmitter and receiver was just under 100-m. The sample ratewas setto 10 seconds,sothat eachC2 point wasfound froma 10secondtime average. n The path integrated C2 measurements are determined by the LOA instruments by computation from the log– n amplitude scintillation(Cχ(r)) ofthe tworeceivingsignals.4,5 The algorithmfor relatingCχ(r) to Cn2 is based on an equation for the log–amplitude covariance function in Kolmogorovturbulence by Clifford et al.6 Proc. of SPIE Vol. 6268 62683E-6 The data was collected by dedicated PCs, one per instrument. During analysis, the data were smoothed by a 120 point (10 minute) boxcar rolling average. This value was chosen for future ease of comparison with local weather station data, sampled at one reading per 10 minutes. Figure 7 compares the extracted IMFs from a single day, from midnight to midnight. There are no data dropouts in the time signal for instrument A, while instrument B is 99.81% valid. A visual examination reveals that the measured C2 functions are very similar n and both instruments have 11 IMFs. Differences are to be found in the IMFs themselves. In Figure 8 we show Empirical Mode Decomposition Empirical Mode Decomposition signalmf1 signalmf1 i i mf2 mf2 i i mf3 mf3 i i mf4 mf4 i i mf5 mf5 i i mf6 mf6 i i mf7 mf7 i i mf8 mf8 i i mf9 mf9 imf10 imf10 imf11 imf11 ires. ires. Figure7.TheIMFsextractedfromeachinstrumentforthesameday(3-March-2006),beginningandendingatmidnight (A on theleft and B on the right). the Hilbert MarginalSpectra derivedfromthe IMFs togetherwith a KolmogorovPowerSpectrum trend scaled tostartcoincidentwiththeMarginalSpectra. BothMarginalSpectrafolloweachotherfairlywell,withsimilar frequency probabilities. Applying a Kolmogorov–Smirnov test to the two Marginal Spectra data sets gives the same means and standard deviations. This is indicative of a P–value of 1 and a cumulative distribution function distance of 0, so we may conclude that the instrument outputs have come from exactly the same distribution and are statistically identical. A further confirmation can be found by calculating the Hilbert phase difference between the two Marginal Spectra. Such a test displays phase synchronization, or lack thereof. In this case we find a phase difference of zero, so that the spectra are perfectly in phase. 6. CONCLUSIONS We have presented the results of applying the Hilbert Huang Transform to C2 time series data. When used n to compare the outputs of two of the same model of commercial scintillometer, we have been able to demon- strate that they provide identical output in terms of their Hilbert Marginal Spectra. It is clear that the HHT technique is a very useful toolin the analysisofnon–stationaryturbulence data andpromisesmuch in terms of understanding the nature of optical turbulence. ACKNOWLEDGMENTS MPJLC would like to thank Norden Huang for introducing him to the Hilbert Huang Transform. Thanks also are due to Sergio Restaino and Christopher Wilcox for making available the scintillometers and providing data processing software. Proc. of SPIE Vol. 6268 62683E-7 REFERENCES 1. N. E. Huang, Z. Shen, S. R. Long,M. C. Wu, H. H. Shih, Q. Zheng, N.-C. Yen, C. C. Tung, and H. H. Liu, “The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis,” Proc. R. Soc. Lond. Ser. A 454, pp. 903–995,1998. 2. N. E. Huang, S. R. Long, and Z. Shen, “A new view of water waves - the Hilbert spectrum,” Ann. Rev. Fluid Mech. 31, pp. 417–457,1999. 3. N. E. Huang, M.-L. C. Wu, S. R. Long, S. S. P. Shen, W. Qu, P. Gloersen, and K. L. Fan, “A confidence limit for the empirical mode decompositionand Hilbert spectralanalysis,” Proc. R. Soc. Lond. Ser. A 459, pp. 2317–2345,2003. 4. G. R. Ochs and T.-I. Wang, “Finite aperture optical scintillometer for profiling wind and C2,” Applied n Optics 17, pp. 3774–3778,1979. 5. T.-I. Wang, “Optical flow sensor.” USA Patent No. 6,369,881B1, April 2002. 6. S. F. Clifford, G. R. Ochs, and R. S. Lawrence, “Saturation of optical scintillation by strong turbulence,” Journal of the Optical Society of America 64, pp. 148–154,1974. Marginal Hilbert Spectrum 10−8 B A 10−9 10−10 frequency−5/3 10−11 2Cn 10−12 10−13 10−14 10−15 10−4 10−3 10−2 10−1 Frequency Figure 8. The 3-March-2006 Marginal Spectra of both instruments with a comparison frequency−5/3 (Kolmogorov Spectrum) trend line. Proc. of SPIE Vol. 6268 62683E-8