The Empirical Mode Decomposition: a Must-have Tool in Speckle Interferometry?

In a wider and wider range of research and engineering activities, there is a growing interest for full-field techniques, featuring nanometric sensitivities, and able to be addressed to dynamic behaviors characterization. Speckle interferometry (SI) techniques are acknowledged as good candidates to tackle this challenge. To get rid of the intrinsic correlation length limitation and simplify the unwrapping step, a straightforward approach lies in the pixel history analysis. The need of increasing performances in terms of accuracy and computation speed is permanently demanding new efficient processing techniques. We propose in this paper a fast implementation of the Empirical Mode Decomposition (EMD) to put the SI pixel signal in an appropriate shape for accurate phase computation. As one of the best way to perform it, the analytic method based on the Hilbert transform (HT) of the so-transformed signal will then be reviewed. For short evaluation, a zero-crossing technique which exploits directly the extrema finding step of the EMD will be presented. We propose moreover a technique to discard the under-modulated pixels which yield wrong phase evaluation. This work is actually an attempt to elaborate a phase extraction procedure which exploits all the reliable information in 3D, – two space and one time coordinates –, to endeavor to make the most of SI raw data. Temporal phase-unwrapping algorithm for automated interferogram analysis, " Appl. Opt. 32, 3047-3052 (1993). 4. X. Colonna de Lega and P. Jacquot, " Deformation measurement with object-induced dynamic phase-shifting, " Appl. Phase demodulation method from a single fringe pattern based on correlation technique with a polynomial form, " Appl.Transform method of fringe-pattern analysis for computer-based topography and interferometry, " J. Fringe-pattern analysis using a 2D Fourier transform, " Appl. Natural demodulation of two-dimensional fringe patterns-Part I, " J. Regularization methods for processsing fringe-pattern images, " Appl. A simple and effective method for filtering speckle interferometric phase fringe patterns, " Opt. Dynamic speckle pattern interferometry (DESPI) phase analyses with temporal Hilbert transform, " Opt. " The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis, " Proc. R. A confidence limit for the empirical mode decomposition and Hilbert spectral analysis, " Proc. R. Evaluation of dynamic speckle activity using the empirical mode decomposition, " Opt. Phase measurement improvement in temporal speckle pattern interferometry using empirical mode decomposition, " Opt. Commun. 275, 38-41 (2007). 21. S. Equis and P. Jacquot, " Phase extraction in dynamic speckle …


Introduction
In full-field techniques, each pixel of the photo-detector probes the surface as an independent sensor.For interferometric experiment, the phase is most likely to vary at a different rate along time from one pixel to another.Hence, the processing task consists in handling in parallel a large quantity of 1D non-stationary signals at a temporal sampling rate given by the frame frequency of the camera.To really take advantage of those 3D (2D space and 1D time) data, efficient signal processing procedures are absolutely mandatory.Although we focus here on SI signals -their large intensity and phase fluctuations making them the worst case to be treated -the method exposed here is not limited to SI signals analysis only.It is applicable to any non-stationary fluctuating signal as soon as there are several oscillations in the whole data set.This method is thus well adapted to deal with any full-field technique signals, and also to any punctual transducer signals.
Fringe analysis has been -and is still -the matter of many researches [1].Phase-shifting techniques have been massively used and substantially improved [2].But, though generally acknowledged as the choicest techniques in terms of measurement and spatial resolutions, they are not in essence good candidates for dynamic problems, despite attempts to adapt it to real time situations [3][4].Resorting to single-frame based techniques is then the natural trend.The mainstream techniques are here the local interpolation of the intensity pattern using splines [5][6], the Fourier transform (FT) method combined with spatial carrier fringes [7][8], the spiral phase quadrature transform method [9], the latter being basically a convenient 2D isotropic Hilbert Transform (HT), and regularization based methods [10].Most of the singleframe techniques need a 2D unwrapping step [11], and the task is not trivial with noisy wrapped phase maps, even if iterative filtering techniques have proven to be efficient and easy to implement [12].Besides the sign ambiguity that cannot be removed without a priori knowledge, single-image techniques depend on the quality of the correlation fringes.With standard apertures of usual systems, the correlation lengths are as small as a few tens of microns both in out-of-plane and in-plane measurements, and the fringes are not exploitable beyond this range.Leaving aside computational considerations for the moment, a fruitful approach lies in the pixel history analysis.It offers not only means to get rid of the intrinsic limit of correlation, but it reduces also greatly the complexity of the unwrapping step, as one deals from now on with 1D signals.The sign ambiguity is moreover easily cleared up by adding a temporal carrier, readily subtracted afterwards [4].The methods devoted to solve the 1D temporal problem include the Morlet wavelet transform [13] and the HT [14].
In order to develop faster and more reliable methods, we proposed in [15] to implement the EMD algorithm for the purpose of reshaping the SI pixel signal to conduct accurate phase computation in conjunction with the HT.The EMD has been successfully used in several domains, dealing with strongly non-stationary signals [16][17][18][19], and has been introduced lately in SI [15,[20][21].The main contribution of this work lies firstly in vindicating the use of the EMD in dynamic SI.A fast implementation of EMD will follow from this standpoint and will be shown to provide well-conditioned signals for an accurate subsequent phase computation, whichever computation method is actually chosen.Being a priori a good candidate for further processing the EMD-transformed signals, the analytic method will be reviewed.The EMD lends itself well also to the application of a zero-crossing search, providing an overall fast evaluation of the phase history.In a next to last section, we will present one way to manage the under-modulated pixels which lead to wrong phase computation.The purpose of this section is not to thoroughly study the behavior of singular pixels and their neighbors.It is rather to propose a solution which relies only on well-modulated pixels and discard the ones with a too low modulation depth, which surely be the scene of complex phenomena accompanied by large phase errors.To the best of our knowledge, this issue has been very seldom tackled [13] and it is definitely worth trying to use the data both in space and in time, to make the most of the huge quantity of readily available data.Some experimental results will conclude the paper.

SI signals: an inappropriate form for straightforward phase computation
SI signals obey the well-known 2-beam interference equation, whether it stands for spatially resolved or integrated regimes: where α is the mean intensity, β the fringe modulation, ψ s the speckle phase, ψ opd the optical path difference phase, i.e. the quantity of interest, and x, y, t are respectively the spatial and temporal coordinates.α, β and ψ s are random variables, and the pixel signals are thus strongly non-stationary, making the phase computation tricky.As an example, a genuine temporal SI signal is depicted in Fig. 1.It is well-known that the three sine qua non conditions to achieve a meaningful phase extraction from a real-valued signal are [21][22][23][24][25]: i) the amplitude and phase modulations spectra have to be well separated, ii) the mean has to be locally zero and iii) the signal spectrum has to be narrow-band.In SI, the spectra separation and the narrow-band conditions are actually intrinsically fulfilled.As can be shown from Fig. 1, the modulation depth variations (the term β in Eq.( 1)) are much slower than the local period of the cosine term.Modulation depth variations are ruled by the statistical spatial properties of the speckle fields and decorrelation effects, while the phase variations depend on the sensitivity of the interferometer.With adequately selected system parameters (the aperture of the recording system, the pixel size, the interferometer sensitivity, the frame frequency, and the rate of phase change), pixel signals take ipso facto the form of a random, slowly varying modulation depth, enveloping the rapid oscillations of the cosine term.So, for the kind of signals we are focused on, accurate phase extraction is ensured as soon as the zero local mean condition is fulfilled.We will show in the next section in what ways the EMD is a neat candidate to clear up this issue.

Standard algorithm
The EMD has been introduced by Huang et al in the late 90's [16] to process non-stationary data.This technique actually decomposes any non-stationary signal s(t) into its intrinsic oscillation modes, acting basically like a filtering process from higher to lower frequencies, but with self-adaptive band-filters.The starting idea is to consider the signal constituted by a detail part (local high frequency) d(t), and a residue part (local low frequency) m(t).The detail part is sifted out from the raw signal by removing the mean envelope, whose computation is based on a cubic spline fitting between the signal extrema.The residue is then considered itself as a signal to process and thus split into a detail and a residue part as well.The modes, namely the intrinsic mode functions (IMF), that could be non-stationary too, have to satisfy two conditions: i) in the whole data set, the number of extrema and the number of zeros differ from each other at most by one, and ii) the mean envelope is zero.The final algorithm [27] is represented in Fig. 2. By construction, the IMFs have a well-behaved HT [16], and more generally they allow good phase extraction.We get in fine the following final decomposition at the rank K: (2) where the d k are the IMFs and m K is the final residue.

Assets and weaknesses of the method
One of the main assets of the EMD is its sparseness: an arbitrary signal is decomposed in fewer components than with classical Fourier or even wavelet analysis.Another important asset is the completeness: by design, the algorithm guarantees a lossless decomposition.However, one of the main drawbacks is the non uniqueness of the final decomposition.It is indeed strongly dependent on the different algorithm parameters and choices, like the sifting ending criterion, the boundaries ending technique (signal continuation), the interpolation method, etc.Some details will be given on the sifting process and on the sampling influence in the coming subsections.The boundary ending is also a sensitive part of the mean estimation: an interpolation kernel is chosen to link the extrema with smooth curves, and the cubic spline fitting requires also extrapolation at the edges of the dataset [16,28].
Bad choices for the aforementioned issues could lead to errors in the entire decomposition, like over-decomposition but also mode-mixing (for instance, higher frequencies oscillations are not caught locally by a given mode but by a successive one that should contain lower frequencies oscillations).It simply follows from the non uniqueness of the decomposition.From the mathematical point of view, this is the most annoying issue, but in practice the previous errors can be straightforwardly identified by simple visual check of the decomposition or by more involved means, as e.g. the quantification of the orthogonality between modes [15].Some readjustments can then be brought to the algorithm, to avoid the errors or compensate them [15,[17][18].At the end, a trade-off has to be found between imposing signal-dependent safeguards to give more robustness to the algorithm and preserving the self-adaptiveness of the method.The purpose is, here, to obtain the most sparse decomposition where each mode leads to an extracted phase with physical meaning.
The outcome of a well-controlled EMD algorithm is actually multiple, ranging from signal detrending to specific features extraction.In addition of the centered signal on which phase extraction, estimation or tracking can be carried out, the location of the signal extrema can provide a rough estimation of the "instantaneous" phase -rather the mean of the instantaneous phase over every half-period.Moreover, the envelopes computation gives also the modulation depth which is helpful to appraise the reliability of the phase computation, and can also be of use to phase tracking techniques (for instance Digital Phase Lock Loop).

Sampling considerations
EMD outcomes depend strongly on sampling conditions [29][30][31].The authors in [29] have shown that a slight misalignment of the extrema with the sampling points results in a wrong estimation of the amplitude and then in over-decomposition and mode leakage.The error is actually maximum when the extremum of the continuous time function is equally spaced from two samples.It emerges from this study that the amplitude difference between the continuous time signal and its sampled version becomes below 5% as soon as the signal is sampled 10 times finer than required by the Shannon sampling theorem.Their criterion provides indeed an indication on the wrong estimation of the amplitude, but the connection with a possible error on the extracted phase has not been addressed.We might lose indeed, the orthogonality of the extracted modes, leading to an energy spreading over successive modes, but the key point for our application is definitely the phase conservation.In [30][31], the authors consider a single tone and define the error as the deviation of the EMD outcome from the original signal.Here, with an original signal defined in Eq. (3), we quantify the phase error due to sampling with the averaging over ϕ (a uniform distribution of random phase pertaining to the speckle effect) of the standard deviation of the difference between the original phase and the phase φ imf extracted from the 1 st IMF: f s is the signal frequency and (2πf s k+ϕ) is the theoretical phase.The phase extraction is carried out through the analytic method as explained in section 4. To reduce boundaries errors due to the discrete FT, a Gaussian window is used to pre-filter the dataset.Simulations have been conducted over 100 realizations of the signal x fs,ϕ in Eq. ( 3) and the dependence of the phase error with the signal frequency 2π/2 r (r being a real number) is shown in Fig. 3.It appears clearly that the EMD behaves better when 2 r is an even integer greater than 3, and especially a power of 2. The largest error is made when 2 r is an odd integer.In the latter case the odd number of samples per period leads indeed to a detrimental asymmetry for the upper and lower envelopes computation.The maximum standard deviation is in average about one hundredth of period and an upper bound varying with the square of the signal frequency gives the global trend of the error (see Fig. 3).Finally, the sampling error is much lower than what we could have expected from the study in [29], and our results are actually closer to the ones presented in [30][31].We end this section with the comforting conclusion that sampling has indeed an influence on the phase extraction but that we can neglect the error as long as the sampling frequency is at least 3 times higher than the Nyquist frequency.

The sifting process influence
Another key step of EMD algorithm is the sifting process and the criterion to stop it.The authors in [32] consider the following sum of two tones: where α is a normalised amplitude ratio ranging from say 10 -2 up to 10 2 , ν is a normalised frequency ratio ranging from say 10 -2 to 1 and ϕ is a random phase uniformly distributed in [-π,π].Their performance criterion measures the deviation of the outcome of the EMD algorithm for the 1 st IMF from the expected mode, i.e. the high frequency term in Eq.( 5): where d i 1 stands for the 1 st IMF obtained after i iterations of sifting process, and || || l2 stands for the l 2 norm (square root of the sum of the square values).When the two tones frequencies are close to each other, i.e. ν≈1, we guess the separation will fail and the EMD outcome will be the signal itself leading to a performance criterion close to 1.However, if the two tones frequencies are well separated, i.e. ν<<1, the EMD will perform a perfect separation, and thus  Those results actually provide a wealth of insights.We observe two very distinct areas basically for α < 1 and for α > 1.It has been shown [32] that the area where α <1 can be modelled by a filtering process with the kernel used for the interpolation − the cubic spline one being acknowledged to perform the best in most cases.However, the filter model is no longer appropriate for the area where α >1.It appears indeed that the sharp boundary between no separation at all and complete separation does not depend on the number of sifting process iterations.Before drawing any conclusion we must study the non-stationary case, but it is a good omen for a fast implementation for SI temporal signals.Indeed, as our signals have necessarily an amplitude ratio greater than 1 and well-separated spectra, we are most likely to be in an area where perfect separation will be achieved in one single iteration.As previously said, the signals we are interested in are strongly non-stationary, and it would be valuable to define a criterion similar to the one in Eq. ( 6) for SI signals defined in Eq. (1).To this aim, we simulate temporal SI signals with the model described in [33].Our model actually lies on the convolution equation between a given impulse response and a complex field of random argument and unitary amplitude, according to both the diffraction theory in the Fresnel approximation and the mechanism of speckle formation.The statistics of the simulated interferometric signals in question are in close agreement with theory [26].Our purpose here is to show that for SI signals, the separation is actually perfect between the AM-FM term and the varying bias in Eq. ( 1) after one single iteration.For each temporal signal we compute the following criterion: ( ) A simulated signal is shown in Fig. 5, while the histograms of the criterion δ 1 i computed for 4096 temporal signals are depicted for i = 1,5 and 10 in Fig. 6.We impose a non-linear digital instantaneous frequency (IF) evolution distributed in the [2π/16, 2π/8] range.It finally turns out from the simulations that whatever is the number of sifting process iterations the criterion defined in Eq. ( 7) is below 0.015: the similarity between the EMDprocessed temporal signal and the true sought-after AM-FM signal is better than 98.5%.We observe that the sifting process actually degrades rather than improves the components separation.It is likely to be due to the boundary choice and to the fact that the sifting process makes the boundaries errors propagating within the dataset.

A fast and accurate implementation of EMD applied to SI
For a pixel signal with enough modulation and as soon as the illumination is well adapted to the sensor dynamic, the relevant information is carried by the high frequency part of the signal (the cosine term in Eq. ( 1)).This information is thus contained in the 1 st IMF, and it becomes useless to proceed further in the modes extraction.We have moreover seen in details that the EMD is effective after one single iteration of the sifting process with non-stationary, narrowband signals containing two well-separated IF.We can extract the 1 st IMF (or pseudo-IMF for the sake of rigorousness) through a single iteration of the EMD algorithm.At this point, it is worth mentioning the method developed by Vikhagen [34] and improved by Carlsson and Wei [35] for deformation measurement in dynamic SI experiments.The phase evaluation method [34] consists in scanning the pixel history signal within a local oscillation to detect a maximum and a minimum value: (α+β) and (α-β) with the notations of Eq. ( 1).There is only one unknown left, the phase ψ that is finally easily computed modulo 2π at each instant using again Eq. ( 1).The improvement of the method [35] consists in a better evaluation of the initial speckle phase, i.e. before deformation, and in a least-squares estimation of the phase during the deformation allowing at the same time the resolution of the sign ambiguity.Even if the methods come from very different starting points, we arrive to a quite similar technique, likely to be nonetheless much faster in our case.Indeed, the EMD does not need to analyze the signal at the oscillation level.With discrete time signals, especially in the non-stationary case, the extrema are very unlikely to be located at sampling points, adding thus some quantization noise to the computation.The higher the FM is, the better the temporal resolution will be, as usual to the detriment of the frequency resolution and, in corollary, of the quantization noise.Not withstanding the strong limitations of this method, we will use is as it is, either to provide a fast evaluation of the phase in long experiments for instance, or to provide an initial guess for phase tracking methods.

How to discard pixels with low modulation: 3D piecewise processing (3DPP)
Due to the intrinsic randomness of speckle, the temporal pixel signal will experience strong fluctuations of modulation depth.When this modulation depth drops to zero, the apparent IF does not reflect the underlying physical phenomenon with a possible frequency doubling.In addition, in the vicinity of those areas of null modulation depth, noise is preeminent and the EMD and, as a consequence, the phase extraction become very inaccurate, as any other method would be by the way.Due to the 1D unwrapping operation, the error will propagate and corrupt the whole phase dataset.For a long experiment (several correlation lengths), time intervals of low modulation will surely occur in the temporal history of each pixel.This new unavoidable issue can be overcome by using the spatial information contained in neighboring pixels.Between the areas where the pixel is classified as irrelevant, the phase will be obtained up to an additive constant, but the discrete IF will be valid.The basic idea is thus to identify the areas where the pixel signal has a modulation depth lower than a certain threshold empirically set (chosen equal to 30 gray levels on Fig. 7-a  Each pixel being independent, when we gather again the temporal sets to form frames, we obtain arrays of discrete IF sampled on a fluctuating non-uniform grid.An interpolation step has to take place to obtain arrays of discrete IF uniformly sampled, which are summed up afterwards to recover the phase map at each instant.The interpolation method will not be detailed here but it is based on the Delaunay triangulation implemented using the Quickhull algorithm [36].Due to the interpolation process, an unavoidable loss of spatial resolution is introduced.But at one point of the probed surface, the spatial resolution, or in other words the density of relevant pixels, changes at each frame and except for areas with almost no activity, a pixel always features enough modulation somewhere in its history to allow meaningful IF extraction.As for any interpolation process, some care has to be put near the dataset boundaries.So, we simply make sure that the area of interest (AOI) is always surrounded by pixels referenced as valid with smooth values of IF (either imposed by a priori knowledge or computed by averaging the AOI border).We compare in the next section our technique with the median filter, which is acknowledged to remove efficiently "pepper and salt" noise and to be edge-preserving.We will see that our technique clearly outperforms the median filter.

Experimental results
A SI setup with an in-plane sensitivity (Leendertz configuration) has been built to follow the rotation of a diffusing metallic plate.The object is illuminated by two divergent laser beams of equal intensity, each of them making an angle θ with the normal to the object surface.The smooth in-plane rotation of 0.25° is achieved through the use of a DC motor.During the motion, 4640 images (1000×1000 encoded on 10 bits) are taken at a rate of 48 fps, and processed by blocks of 1024 frames with an overlapping of 100 frames between them to reduce the boundaries errors.In Fig. 8, S 1 and S 2 designate the two unit vectors of the illumination directions, while the observation direction coincides with the object surface normal.The overall phase change between the two arms when a displacement L occurs is given by: ( ) where S and L x are respectively the sensitivity of the interferometer and the projection of the motion on the x-axis.Due to its symmetry, the sensitivity is actually simply related to the physical parameters of the setup by the following relation: In our experiment, the sensitivity S is equal to 8.5 rad/µm.For off-axis points, due to the divergence of the illumination beams, the interferometer has slight sensitivity to L y , and L z .Those sensitivities are more than one order of magnitude lower than S, and will be thus neglected.The divergent illumination leads also to variations of the sensitivity S within the field that appear to be negligible (less than 0,2 %).We limit the AOI to the half the rotating plate (see Fig. 8).In this experiment, we do not use an optical temporal carrier -mandatory to find out the direction of the deformation -to have the largest measurement bandwidth, which otherwise would have to be shared with the carrier frequency.It is thus enough to process the temporal pixels of half the plate.In Fig. 9(a), the movie of the phase computed through the use of the fast EMD followed by the analytic method is shown.The result is very noisy and very quickly unexploitable for pixels experiencing less and less fringes of displacement (toward the central part of the plate).This area is denoted as the "blindness" area as the algorithm does not allow measuring accurately what happens inside.There are actually two main cases to distinguish: the case where there are few fringes of displacement and the case where there is almost none.The latter one corresponds to the very central part of the metal plate and cannot be solved without a temporal carrier.When there is very low activity during the pixel history, the presence of noise and the sparsity of fringe extrema make the extrema finding step very inaccurate.The outcome of the EMD algorithm contains then long areas of only noise.Applying the analytic method to such signals leads to a result where the useful phase information is overwhelmed by large and meaningless phase values.Filtering the raw phase maps with a box 5×5 median filter leads to the result shown in Fig. 9-b, while the phase coming from the 3DPP is shown in Fig. 9(c).The benefit of this latter technique is obvious and greatly extends the measurement bandwidth.The cross sections of the total final phase map in each case are compared in Fig. 10.Not surprisingly, errors exist even when the pixel experienced a lot of activity.The pixel size in the object plane is 70×70 µm 2 , i.e. the fourth of the largest displacement on the metal plate.So decorrelation will surely occur for all pixels outside a circle four times smaller than the plate.The good point is that with the 3DPP technique, only a strip of 5 mm wide in the center of the plate cannot be characterized.This "blindness" strip is 4 times wider with the fast-EMD & HT technique and still 3 times wider after the filtering step with a 5×5 median filter.Adding optically a temporal carrier would suppress this blind zone, at the expense of a reduction of the measurement bandwidth as said before.

Conclusion and outlooks
This paper proposed a new efficient and flexible processing tool to handle SI and more generally any full-field technique signals in dynamic regimes.We have shown that the combination of the EMD, the HT and the 3DPP is perfectly able to tackle the challenge of dynamic behaviours characterization even for displacements and/or deformations beyond the classical limit of the correlation length.

Fig. 1 .
Fig. 1.A typical temporal pixel signal experimentally obtained in SI (gray levels in ordinate and time measured in frame number in abscissa).

Fig. 3 .
Fig. 3. Standard deviation of the difference between the extracted phase from the 1 st IMF and the theoretical phase with an upper bound.
$15.00 USD Received 6 Nov 2008; revised 22 Dec 2008; accepted 23 Dec 2008; published 7 Jan 2009 (C) 2009 OSAlead to a criterion close to 0. For intermediate cases, the guess is not trivial.We report a result from[31]  in Fig.4-a, where the average criterion over 100 realisations of ϕ is mapped for 10 iterations of sifting process.

Fig. 4 .
Fig. 4. (a) averaged criterion δ for 10 iterations of sifting process, (b) EMD filter model (cubic spline kernel) with simulations in solid line and predictions in dashed line.
signals leading to a certain separation criterion #103798 -$15.00USD Received 6 Nov 2008; revised 22 Dec 2008; accepted 23 Dec 2008; published 7 Jan 2009 (C) 2009 OSA ), discard the IF values within those areas (highlighted by gray hatching on Fig.7-a and b; pixel validity indicator equal to 1), while keeping them unchanged outside (pixel validity indicator equal to 0).

Fig. 7 .
Fig. 7. (a) Pseudo-IMF obtained after fast EMD with areas identified as invalid (1) or valid (0) (gray levels in ordinate and sample number in abscissa), and (b) its discrete extracted IF.