## Abstract

Digital speckle pattern interferometry (DSPI) fringes contain low spatial information degraded with speckle noise and background intensity. The denoising technique proposed recently based on bi-dimensional empirical mode decomposition (BEMD) could implement noise reduction adaptively. However, the major drawback of BEMD, called mode mixing, has affected its practical application. With noise-assisted data analysis (NADA) method, bi-dimensional ensemble empirical mode decomposition (BEEMD) was proposed, which has solved the problem of mode mixing. The denoising approach based on BEEMD will be presented, compared with other classic denoising methods and evaluated both qualitatively and quantitatively using computer-simulated and experimental DSPI fringes.

© 2011 OSA

## 1. Introduction

Digital speckle pattern interferometry (DSPI) is a full field, non-contact and real time testing method to measure the difference in out-of-plane displacement. The characteristic of real time makes non-specialists be able to visualize the measurement results as rapidly and readily as possible. DSPI is faster in operation and less sensitive to environmental perturbations than holographic interferometry. However, this method cannot attain high quality image as that of holographic interferometry due to the low resolution of video camera system. Consequently, the development of noise reduction techniques is such a critical subject in DSPI. The Fourier methods don’t preserve details of the object due to the reason that it expands the original function in terms of orthonormal basis functions of sine and cosine waves of infinite duration. The wavelet transforms depend on the decomposition of DSPI fringes in fixed basis functions, requiring the intervention of an external operator. To reduce noise in DSPI fringes adaptively, the denoising method based on BEMD was proposed by M. B. Bernini *et al* [1]. Although this BEMD-based method has proved to be effective to some extent, it still leaves some annoying difficulties for further research.

The major drawback of this denoising method is the frequent appearance of leaky wave, which is caused by 2D mode mixing in BEMD. Leaky wave can be explained that some elements originally belonging to one IMF may leak into other IMF components. This phenomenon makes BEMD-based denoising method lose its power in reducing noise and preserving the fundamental component of DSPI fringes. To overcome the problem of mode mixing generally existing in BEMD without introducing a subjective intermittence test, a new noise-assisted data analysis (NADA) method, called bi-dimensional ensemble empirical mode decomposition (BEEMD), was proposed. With the assisted white noise, 2D mode mixing can be eliminated almost completely and the denoising method based on BEEMD presents better performance since the phenomenon of leaky wave has been restrained obviously.

## 2. Bi-dimensional ensemble empirical mode decomposition

#### 2.1 The phenomenon of mode mixing in bi-dimensional space

There generally exist some unresolved drawbacks as long as BEMD algorithm is used. The major one is 2D mode mixing, which is defined [2] as a single IMF either consisting of signals of widely disparate scales, or a signal of a similar scale residing in different IMF components. Mode mixing is consequence of signal intermittency. Intermittency cannot only cause serious aliasing in the time-frequency distribution, but also make the physical meaning of each individual IMF unclear in the mode-mixing area.

A typical example will be presented to illustrate the problem of 2D mode mixing in the BEMD method. Lena has been a popular image in the image processing fields and has served as the bench-mark image for testing various methods.

Through careful examination of the second orders IMFs component in Fig. 1(c) , one can clearly see several zones with deep color, for example, at the edge of face, nose, arm and hat, considered as 2D mode mixing, which is so serious as to make the fundamental components unclear. The similar conclusions can be obtained from other IMFs. This is the so-called 2D mode mixing.

#### 2.2 Bi-dimensional ensemble empirical mode decomposition

To overcome the problem of 2D mode mixing, a new noise-added signal analysis method was proposed, called BEEMD [3]. The components are defined as the mean of an ensemble. The subject in each trial consists of original signal and white noise with finite amplitude. One can separate scales more precisely through BEEMD without any subjective criterion. This ensemble approach takes full advantage of statistical characteristic of white noise. Since the white noise is different in each trial, it is canceled out in the ensemble mean of sufficient trails and the ensemble mean is treated as the true answer. The BEEMD algorithm can be described as follows: (1) add a 2D white noise signal to the original signal; (2) decompose the compound signal into IMFs; (3) repeat steps 1 and 2 with disparate white noise signals each time; (4) obtain the ensemble means of corresponding IMFs of all the decompositions as the final result.

However, the algorithm for BEEMD in some references is not the real BEEMD, usually called pseudo-BEEMD. The main reason is that the sifting algorithm for BEEMD still utilizes the 1D method. In procedure, the decomposition subject is a 1D signal, a row or a column extracted from image matrix and the decomposition method is EMD. The decomposition goes on row (column) after row (column). Cubic spline interpolation function “spline” is often used, while in this paper one can detect the extreme points with 4-connected neighbors technique and treat a signal with radial basis interpolation function toolbox “FastRBF”, with which the 2D envelope can connect every extreme point and the second derivative is continuous so that the surface is smooth enough. This is the real BEEMD.

We will still use the Lena to present the potential in resolving the problem of 2D mode mixing. The decomposition results are displayed in Fig. 2 when BEEMD is used to decompose Lena image with the same stop criterion as BEMD previously mentioned in this paper, 100 trails and the standard deviation of white noise signal set to one fifth of that of the original signal.

Through careful examination of Fig. 2, we will obtain the conclusion: we can scarcely see the so-called zones with deep color and the fundamental components are uncontaminated by the intermittence and much more visible than the corresponding one in Fig. 1. We can conclude that traditional BEMD cannot realize the ideal effect of BEMD. From the experimental results, BEEMD can resolve the problem of 2D mode mixing successfully and can be applied in 2D space.

## 3. Adaptive noise reduction method based on BEEMD for DSPI fringes

DSPI fringes contain low spatial information degraded with speckle noise and background intensity. To avoid the application of a phase retrieval algorithm [4], several methods based on the Fourier and the wavelet transforms [5,6] have been proposed to reduce noise in DSPI fringes. However, the intervention of an external operator, subjective basis functions, is required in these methods. As an adaptive method of data analysis, BEEMD is introduced to reduce noise for DSPI fringes in this paper.

#### 3.1 The relationship between DSPI fringes and IMF components

The intensity distribution of the interference pattern of the vibrating object can be expressed as:

The Eq. (1) cannot well reflect the relationship between DSPI fringes and IMF components. Therefore, it has to be converted into another expression [7] as fellows by some mathematical transforms.

#### 3.2 The denoising scheme for DSPI fringes

According to Eq. (3), it is concluded that *K* is a critical coefficient to determine which IMFs ${\overline{c}}_{j}\left(x,y\right)$ stand for speckle noise and which IMFs ${\overline{c}}_{j}\left(x,y\right)$ correspond to the fundamental component. Therefore, the BEEMD-based denoising approach can be summarized for retaining the terms $\sum _{j=K+1}^{N}{\overline{c}}_{j}}\left(x,y\right)$ and removing the terms $\sum _{j=1}^{K}{\overline{c}}_{j}}\left(x,y\right)$ and $\overline{r}\left(x,y\right)$.

Hence finding a proper coefficient *K* is an important problem to be tackled. Generally speaking, the methods for determinating coefficient *K* may be divided into three types: (1) just let $K=1$. In other word, simply consider the first orders IMFs as the speckle noise. (2) use the trial method to find proper *K* manually. (3) design an algorithm to choose the suitable value of *K* automatically. Obviously, the third method is more scientific and convenient. Denoting the first and the second terms of Eq. (3) as ${n}_{K}\left(x,y\right)$ and ${c}_{K}\left(x,y\right)$ respectively, we calculate the quantity [7]

*N*with an increment of 1. Since the autocorrelation degree of the fundamental component must be an abrupt change in the curve of $X\left(K\right)$, which represents the critical power ratio of the autocorrelation functions between the fundamental component and the noise. To find this critical point, the waveform $Y\left(K\right)=X\left(K\right)/X\left(K+1\right)$ is calculated. The exact value of

*K*at the abrupt point is obtained when $Y\left(K\right)$ reaches its maximum. For example, in Fig. 3(a) we can calculate the value of $X\left(K\right)$ and plot its curve. Then the value of $Y\left(K\right)$ is computed and its curve is plotted in Fig. 3(b). It is clearly observed that

*K*should be equal to 2 and it means that the first two orders IMFs stand for speckle noise. In practice, one can easily implement this algorithm and then find the appropriate value of

*K*after obtaining all the IMFs.

Besides, there are two processing techniques for selective use:

- 1. It is well known that DSPI fringes contain high density of extreme points. In order to detect extreme points more easily and implement the BEEMD algorithm faster, we need apply an average filter to DSPI fringes before denoising. One cannot only remove the highest frequencies of speckle noise but also preserve the desirable details effectively if using an average filter with an appropriate kernel size, usually$5\times 5$. This technique is used to reduce the computational cost of the BEEMD algorithm.
- 2. Generally speaking, an oversmoothing effect may be produced as long as the BEEMD-based denoising method is used. Therefore, the oversmoothed fringe has to be normalized after denoising. The normalization algorithm [8] uses two orthogonal bandpass filters. One removes the negative frequencies along the horizontal direction and the other removes the negative frequencies along the vertical direction. Then, the respective phases of two filters are combined to yield an estimation of the normalized modulation intensity. So this technique is a post-processing approach used to avoid oversmoothing.

#### 3.3 Denoising for computer-simulated DSPI fringes

In most cases encountered in DSPI fringes, noise $n\left(x,y\right)$ is intermittent and will cause 2D mode mixing in BEMD and leaky wave in BEMD-based denoising method. The phenomenon of 2D mode mixing has been presented in this paper. To tackle the problem of leaky wave, the BEEMD-based noise reduction method is introduced into DSPI fringe preprocessing.

In order to make a comparison of the denoised fringes with the noise-free ones, the fringes are simulated with a technique [9] that reproduces the generation of DSPI fringes of vibration mode, with a resolution of $256\times 256$ pixels and 256 gray levels.

In order to evaluate the performance of the BEEMD method numerically, we used image quality index *Q* [10] defined as

*O*is the noise-free image,

*E*is the filtered fringe pattern, $\overline{E}$, ${\sigma}_{E}$ and $\overline{O}$, ${\sigma}_{O}$ are the means and the standard deviations of the filtered fringe pattern and the noise-free image, respectively, and ${\sigma}_{EO}$ is the covariance between

*E*and

*O*. The

*Q*index allows to model distortions of any kind as a combination of three factors: luminance distortion, loss of correlation and contrast distortion. The dynamic range of the index is [-1,1], being $Q=1$ the best value. The index is calculated using a $5\times 5$ window that is displaced pixel by pixel across the horizontal and vertical direction of the image and the final result was average of all previously computed values.

Processed results using two kinds of denoising approaches are displayed in Fig. 4 . Figure 4(b) shows a fringe pattern denoised with the BEEMD method while Fig. 4(c) displays the result obtained on the same fringe pattern but using the BEMD approach.

Table 1
presents a numerical comparison of quality index *Q* between the results obtained with the proposed BEEMD denoising method and those given by the BEMD approach.

Figure 5 shows the noise, the fundamental frequency component and the background intensity of the simulated fringe pattern by applying the BEEMD-based denoising method.

#### 3.4 Denoising for experimental DSPI fringes

Equation of motion for free vibration of plate of uniform thickness made of homogeneous isotropic material is given by [11]

where*ω*is the out-of-plane displacement, $D=E{h}^{3}/12\left(1-{\upsilon}^{2}\right)$ the plate stiffness,

*E*the modulus of elasticity,

*h*the thickness of the plate surface,

*ρ*the mass density per unit area of plate surface,

*t*the time, ${\nabla}^{4}$the biharmonic differential operator, and ${\nabla}^{2}=\left({\partial}^{2}/\partial {x}^{2}\right)+\left({\partial}^{2}/\partial {y}^{2}\right)$ in rectangular coordinate.

In order to validate this denoising method in practice, we implement a simple, general but classic experiment. The experimental subject is a cantilever beam (i.e. aspect ratio = 1, Young’s modulus of elasticity = 70GPa, mass density = 270kg/m^{3}, and Poisson’s ratio = 0.3). Function generator regulates the frequency and magnitude of the force of the exciter and was set to generate the sinusoidal signal at the center of the cantilever. One of the edges was fixed and the boundary conditions for fixed edge are $\omega =0$ and $\partial w/\partial x=0$.

Denoising methods are tested for effectiveness to improve signal to noise ratio (SNR) of DSPI fringes. A parameter to be calculated for testing denoising method is speckle index. The speckle index [12] can be expressed as fellows:

*C*is the speckle index,

*σ*is the standard deviation, and

*m*is the mean.

The SNR is a reciprocal of the speckle index [13]

We use six filtering schemes to denoise the experimental DSPI fringes in Fig. 6 in order to see the potential of each one for the reduction in speckle index and improvement in SNR. The processed results are displayed in Fig. 7 and the corresponding values of speckle index and SNR are presented in Table 2 . The six denoising schemes are described as follows:

## 4. Conclusion

Although denoising method based on either BEMD or BEEMD posses the same scheme, it still can be clearly observed that Fig. 4(c) contains much more noise than Fig. 4(b). The main reason is that the elements originally belonging to a certain IMFs component leak into the others due to 2D mode mixing. This may make noise reduction incomplete or fundamental components lack some important elements. From the simulated results, this phenomenon of leaky wave can be restrained by the BEEMD algorithm. Therefore, the denoising method has been greatly improved when using BEEMD instead of BEMD. The similar conclusion can be precisely obtained by numerical evaluation in Table.1, in which the *Q* value is higher when using BEEMD-based denoising method.

In order to verify this denoising method in practical experiment, we also use other five classic methods to improve the SNR of experimental DSPI fringes. The processed results are presented in Fig. 7. It is obvious that the best SNR is obtained in Table.2 when using BEEMD-based denoising method. The BEEMD-based denoising method presents better performance than other congeneric methods (EEMD, BEMD) and is also superior to different kinds of methods (Fourier transforms, Daubechies wavelet and Symlet wavelet). Moreover, the method based on BEEMD proposed in this paper is an adaptive denoising technique.

However, it has to be pointed out that such an ensemble BEMD approach requires larger computational resource, proportional to the number of the trials as the regular BEMD. A typical implementation of BEEMD is usually 100 independent trials in the ensemble. If processing a $256\times 256$ DSPI fringe using BEEMD-based denoising method, it will take about 70 seconds in each trial with the CPU configuration Intel(R) Core(TM)2 Quad CPU Q8400 @ 2.66GHz 2.67GHz and the memory 2GB. Pseudo-BEEMD method may require much more time for calculation. The large demand of computational resource by far has imposed a barrier of developing BEEMD.

## Acknowledgment

The authors are grateful for the research support received from the National Natural Science Foundation of China (NSFC No. 10972137).

## References and links

**1. **M. B. Bernini, A. Federico, and G. H. Kaufmann, “Noise reduction in digital speckle pattern interferometry using bidimensional empirical mode decomposition,” Appl. Opt. **47**(14), 2592–2598 (2008). [CrossRef] [PubMed]

**2. **Z. Wu and N. E. Huang, “Ensemble empirical mode decomposition: a noise assisted data analysis method,” Adv. Adapt. Data Anal. **1**(01), 1–41 (2009). [CrossRef]

**3. **Z. Wu, N. E. Huang, and X. Chen, “The multi-dimensional ensemble empirical mode decomposition,” Adv. Adapt. Data Anal. **1**(03), 339–372 (2009). [CrossRef]

**4. **A. Federico and G. H. Kaufmann, “Phase retrieval in digital speckle pattern interferometry by application of two-dimensional active contours called snakes,” Appl. Opt. **45**(9), 1909–1916 (2006). [CrossRef] [PubMed]

**5. **R. Kumar, D. P. Jena, and C. Shakher, “Application of wavelet transform and image morphology in processing vibration speckle interferogram for automatic analysis,” Proc. SPIE **8082**, 80821Y, 80821Y-5 (2011). [CrossRef]

**6. **A. Federico and G. H. Kaufmann, “Comparative study of wavelet thresholding methods for denoising electronic speckle pattern interferometry fringes,” Opt. Eng. **40**(11), 2598–2604 (2001). [CrossRef]

**7. **X. Zhou, H. Zhao, and T. Jiang, “Adaptive analysis of optical fringe patterns using ensemble empirical mode decomposition algorithm,” Opt. Lett. **34**(13), 2033–2035 (2009). [CrossRef] [PubMed]

**8. **A. Federico and G. H. Kaufmann, “Local denoising of digital speckle pattern interferometry fringes by multiplicative correlation and weighted smoothing splines,” Appl. Opt. **44**(14), 2728–2735 (2005). [CrossRef] [PubMed]

**9. **P. D. Ruiz and G. H. Kaufmann, “Evaluation of a scale-space filter for speckle noise reduction in electronic speckle pattern interferometry,” Opt. Eng. **37**(8), 2395 (1998). [CrossRef]

**10. **Z. Wang and A. C. Bovik, “A universal image quality index,” IEEE Signal Process. Lett. **9**(3), 81–84 (2002). [CrossRef]

**11. **A. W. Leissa, “The free vibration of rectangular plates,” J. Sound Vibrat. **31**(3), 257–293 (1973). [CrossRef]

**12. **C. Loizou, C. Christodoulou, C. S. Pattichis, R. Istepanian, M. Pantziaris, and A. Nicolaides, “Speckle reduction in ultrasound images of atherosclerotic carotid plaque,” *Proc. IEEE 14th Intl. Conf. Digital Signal Process*, 525–528 (2002).

**13. **T. R. Crimmins, “Geometric filter for speckle reduction,” Appl. Opt. **24**(10), 1438–1443 (1985). [CrossRef] [PubMed]

**14. **C. Shakher, R. Kumar, S. K. Singh, and S. A. Kazmi, “Application of wavelet filtering for vibration analysis using digital speckle pattern interferometry,” Opt. Eng. **41**(1), 176 (2002). [CrossRef]

**15. **R. Kumar, I. P. Singh, and C. Shakher, “Measurement of out-of-plane static and dynamic deformations by processing digital speckle pattern interferometry fringes using wavelet transform,” Opt. Lasers Eng. **41**(1), 81–93 (2004). [CrossRef]

**16. **R. Kumar, “Wavelet filtering applied to time-average digital speckle pattern interferometry fringes,” Opt. Laser Technol. **33**(8), 567–571 (2001). [CrossRef]

**17. **K. Creath, “Temporal phase method,” in *Interferogram Analysis*, D. Robinson, and G. Reid, eds. (Institute of Physics, 1993), pp. 94–140.

**18. **Y. Lei, Z. He, and Y. Zi, “Application of the EEMD method to rotor fault diagnosis of rotating machinery,” Mech. Syst. Signal Process. **23**(4), 1327–1338 (2009). [CrossRef]