## Abstract

The phase demodulation method of adaptive windowed Fourier transform (AWFT) is proposed based on Hilbert-Huang transform (HHT). HHT is analyzed and performed on fringe pattern to obtain instantaneous frequencies firstly. These instantaneous frequencies are further analyzed based on the condition of AWFT to locate local stationary areas where the fundamental spectrum will not be interfered by high-order spectrum. Within each local stationary area, the fundamental spectrum can be extracted accurately and adaptively by using AWFT with the background, which has been determined previously with the presented criterion during HHT, being eliminated to remove the zero-spectrum. This method is adaptive and unconstrained by any precondition for the measured phase. Experiments demonstrate its robustness and effectiveness for measuring the object with discontinuities or complex surface.

©2012 Optical Society of America

## 1. Introduction

Phase demodulation is a critical step in 3D shape measurement. It is usually completed by using phase-shift method [1] or transform-based method. The latter is chosen frequently because it needs only one fringe pattern, which is helpful for the dynamic measurement. Fourier transform (FT) is a global operation which can just do the frequency analysis for signals, so it is more suitable for analyzing the global stationary signals [2]. Here a stationary signal means the signal in a fringe pattern varies slowly, in other words, the fundamental spectrum of a stationary signal will not be superposed by other order spectra and so the fundamental spectrum can be extracted completely [3]. However, signals of deformed fringe pattern sometimes are no longer stationary globally due to the height modulation of object, so good ability of local analysis for fringe pattern is necessary. Various methods of time-frequency analysis have been proposed to balance the time-frequency resolution for phase retrieval. Windowed Fourier transform is performed with fixed window width which cannot meet Heisenberg uncertainty theory [4]. The ridge of Wavelet transform (RWT) method [5, 6] or the ridge of S-transform method [7, 8] are both carried out by extracting the phase on the ridge of spectra. Adaptive windowed Fourier transform (AWFT) is implemented through adjusting window width with scale factors got by RWT directly [9]. The methods of multiscale windowed Fourier transform (MWFT) in [3, 10] control the window width with instantaneous frequencies got by RWT instead of just using scale factors directly. In these correlative researches, either scale factors or instantaneous frequencies are always obtained by using the ridge method which is based on the assumption of$\phi (x)=\phi (b)+\phi \text{'}(b)(x-b)$, that is, the high-order terms in Taylor series of the measured phase $\phi (x)$are neglected and $\phi \text{'}(x)$ is assumed as a slow variation. However, if the surface of object is complex or carries steep edge, the measured phase cannot meet the assumption. Then scale factors or instantaneous frequencies got by the ridge method will be smoothed subjectively leading to the meaningless change of window for local fringe analysis. 2D continuous wavelet transform is also limited by the derivation of rigorous equations [11], and meanwhile it also has another constraint including the adaptive selection for mother wavelet or parameters as the same problem as 1D RWT [12]. There are also many other methods pursuing the much robust and accurate phase retrieval method such as the ensemble empirical mode decomposition method [13] and the windowed Fourier ridge algorithm assisted with statistical analysis [14], etc.

The presented method aims to enhance the ability of time-frequency analysis for phase demodulation of single fringe pattern. The basic conditions for using AWFT are analyzed firstly. The instantaneous frequencies are obtained with Hilbert-Huang transform (HHT) and then selected according to the analyzed conditions to locate local stationary area with the presented steps, the local stationary area which is defined the signals in a Gaussian window whose fundamental spectrum can be extracted completely. Within each local stationary area in one line of fringe pattern, the fundamental spectrum will not be interfered by high-order spectrum because the Gaussian window can change adaptively during AWFT, and it also can avoid the spectrum aliasing caused by zero-order spectrum because the background is eliminated, the background which is determined with the presented criterion during the analysis of HHT in advance. The proposed method has good adaptiveness and is unrestricted by any precondition for the measured phase. Experiments show that it is of great robustness under noisy environment and it can extract the phase accurately even for the fringe pattern which deforms dramatically.

## 2. The key issue of using the AWFT method in phase demodulation of fringe pattern

1D signal of the deformed fringe pattern is usually described as

where*A*(

*x*) and

*B*(

*x*) represent the background and the amplitude modulation, respectively,

*f*is the fundamental frequency,$\varphi (x)$is the evaluated phase modulation and

_{0}*η*(

*x*) is the noise.

For the Fourier based method, the core work of extracting phase is to extract the fundamental spectrum accurately. From the description in Fig. 1
, it can be seen two conditions should be met to extract the fundamental spectrum *Q*_{1} completely:

*f*is analogous to an instantaneous frequency and (

_{m}*f*)

_{m}_{min}is the minimum instantaneous frequency for the

*m*th spectrum component [15], (

*f*

_{1})

_{max}and (

*f*

_{1})

_{min}are the maximum and minimum instantaneous frequency for the fundamental spectrum component, respectively, and

*f*is the cut-off frequency for the zero-spectrum component.

_{b}AWFT is superior to FT because it can provide spectrum analysis for local signals in time domain while FT is just good at processing global stationary signals [9]. 1D AWFT performed on *I* (*x*) in Eq. (1) is expressed as

*b*and the window width changes narrow or wide with changing scale factor

*a*. The Gaussian function has the statistic characteristic such asSo it follows

As the above analysis, the key of using AWFT is also to extract the fundamental spectrum completely as the Gaussian window changed each time. Moreover, the width of Gaussian window can be controlled by

where*L*(

*x*) is defined as the full width at 1/

*e*

^{2}for the Gaussian window [10] and it also represents the length of a local stationary area in a line of fringe pattern. Therefore, Eq. (2) and Eq. (3) can also be thought as the conditions for using AWFT and locating local stationary areas in one line of fringe pattern turns to be the key of applying AWFT in phase demodulation.

## 3. HHT applied in the AWFT method

#### 3.1 HHT

HHT is a data-driven approach which is unconstrained by Heisenberg uncertainty theory and it needs no basic function to do time-frequency analysis. There are two contents contained in HHT [16, 17]: empirical mode decomposition (EMD) and Hilbert transform. Non-stationary signal is multi-component, thus it can be decomposed into a collection of intrinsic mode functions (IMFs) by EMD firstly. Each IMF is the mono-component and is constrained by two conditions: at any point, the mean value of the local maxima envelope and the local minima envelope is zero; the number of extreme points and zero crossings should be equal or differ at most by one. The process of EMD is carried out as

where*v*(

*x*) is the arbitrary time series, $v{(x)}_{\text{max\_envelope}}$ and $v{(x)}_{\text{min\_envelope}}$ stands for the upper envelope and lower envelope of the maxima and minima got by the cubic spline interpolation, respectively.

*h*(

*x*) is judged whether an IMF, if not, the process of Eq. (8) is repeated, or

*h*(

*x*) is expressed as m

_{1}(

*x*) as an IMF and then subtracted from

*v*(

*x*). Then the residual becomes a new

*v*(

*x*) which continues the process of Eq. (8). The sifting is continued until there is no more IMF contained. The result of decomposition is expressed aswhere ${m}_{n}(x)$ is the

*n*th IMF,

*N*represents the total number of IMFs and

*r*(

*x*) is the residual which contains no more IMF. The IMFs are ordered from high frequency to low and there is no loss of energy after decomposition. More issues for using EMD such as border effect, stopping criteria for sifting, interpolation etc. can refer to [18, 19].

For any IMF *m _{n}* (

*x*), its Hilbert transform

*m*′

*(*

_{n}*x*) is

*P*is the Cauchy principal value. Then the analytic signal of

*m*′

*(*

_{n}*x*) is described aswhere ${\lambda}_{n}(x)={\{{m}_{n}{}^{2}(x)+{[{m}_{n}\text{'}(x)]}^{2}\}}^{1/2}$ is the amplitude of the analytic signal and ${\phi}_{n}(x)=\mathrm{arctan}[{m}_{n}\text{'}(x)/{m}_{n}(x)]$ is the phase. And instantaneous frequency is defined asAfter performing Hilbert transform on each IMF component, the original signal can be expresses as

#### 3.2 Locating local stationary areas with HHT

### 3.2.1 Determining the needed instantaneous frequencies through analysis of HHT

According to the analysis in Section 2, the Eq. (2) and Eq. (3) are used to locate local stationary areas for one line of fringe pattern directly. So it needs to get the needed instantaneous frequencies accurately in the two equations first.

The multi-component signal *I*(*x*) in Eq. (1) is decomposed into a collection of mono-component IMFs by EMD, so these IMFs are sorted into three groups such as

*m*(

_{n}*x*) is the high-frequency IMFs corresponding to

*η*(

*x*) in Eq. (1), the second group is the fundamental IMFs corresponding to$B(x)\mathrm{cos}[2\pi {f}_{0}x+\varphi (x)]$ and the third matches the background

*A*(

*x*). The three groups are separated by

*m*

_{k}_{1}(

*x*) and

*m*

_{k}_{2}(

*x*) which are the first fundamental IMF and the first background IMF, respectively. There is always an instantaneous frequency at a moment if performing Hilbert transform on each IMF component. But only one of the instantaneous frequencies reflects the real instantaneous change in phase. Therefore, the (

*f*

_{1})

_{max}in Eq. (2) must exist in

*m*

_{k}_{1}(

*x*) and (

*f*)

_{m}_{min}must exist in

*m*

_{(}

_{k}_{1-1)}(

*x*).

*m _{k}*

_{1}(

*x*) is defined as the first fundamental IMF which contains the most fundamental components and contains rarely high-frequency components. It can be determined with

*h*

_{max}is the maximum marginal spectrum in the 2D distribution of marginal spectrum and frequency, $Ins{f}_{n}^{{h}_{\mathrm{max}}}$is the corresponding instantaneous frequency at

*h*

_{max}for the

*n*th IMF,

*f*is the fundamental frequency and

_{0}**min**{…} denotes the operation for extracting the ordinal

*n*when$\gamma (n)$is the smallest. Equation (16) describes that the instantaneous frequency at the maximum marginal spectrum for

*m*

_{K}_{1}(

*x*) must be the closest with the fundamental frequency. This criterion is effective even if there is the mode-mixing problem between the IMFs, the problem which is caused by mixed different scales of extreme points for signals [20]. In Eq. (16),

*f*

_{0}is the fundamental frequency of the designed fringe pattern, which should be known in advance. However, it may be unknown in rare cases. Then it can be estimated by detecting the frequency at the biggest maximum marginal spectrum for all the IMFs, because the maximum marginal spectrum of the fundamental IMF should be the biggest among all the IMFs.

An example is used to describe the criterion clearly. Figure 2(a) shows the signal got by

where $\varphi (x)$ is the simulated phase and $\eta (x)$ is the random noise with amplitude 0~0.02. After performing EMD on*S*(

*x*),

*m*

_{1}(

*x*)~

*m*

_{7}(

*x*) are the IMFs arranged from high frequency to low frequency and

*r*(

*x*) is the residual, which are shown in Fig. 2(b). Figure 2(c) shows the corresponding instantaneous frequencies

*Insf*

_{1}(

*x*)~

*Insf*

_{7}(

*x*) for

*m*

_{1}(

*x*)~

*m*

_{7}(

*x*), and the corresponding marginal spectrum

*h*

_{1}(

*f*)~

*h*

_{7}(

*f*) is shown in Fig. 2(d). The corresponding $\gamma (n)$ listed in Table 1 is computed according to Eq. (16). It’s clear that$\gamma (2)$is the smallest value and thus

*K*1 is obtained as 2. As

*m*

_{2}(

*x*) is identified as the first fundamental IMF, (

*f*

_{1})

_{max}and (

*f*)

_{m}_{min}in Eq. (2) exist in

*Insf*

_{2}(

*x*) and

*Insf*

_{1}(

*x*), respectively.

In this example, the mode-mixing problem exists as shown in Fig. 2(b). *m*_{2}(*x*) loses some fundamental components while these losses occur at the corresponding positions in *m*_{1}(*x*) and *m*_{3}(*x*) because of the completeness character of EMD [16]. In Fig. 2(c), there are peak values at the mode-mixing positions in *Insf*_{2}(*x*), which are produced by the tiny values (close to zero) in *m*_{2}(*x*) at the corresponding positions when computing instantaneous frequencies. They are false instantaneous frequencies. However, the values at the corresponding positions in *Insf*_{1}(*x*) represent the actual instantaneous frequencies because they are obtained by the actual fundamental components. These actual instantaneous frequencies are much smaller than the peak values, and so the signals at the border of mode-mixing positions will be detected as the non-stationary signals according to Eq. (2). (Here a non-stationary signal is the signal which varies rapidly somewhere but slowly in other place, namely the signal whose fundamental spectrum is superposed by other order spectra. And just non-stationary signals can divide one line of fringe pattern into several local stationary areas, because when a signal is non-stationary relive to one of its two adjacent local stationary areas, it will be stationary relative to the other one.) Actually, the mode-mixing problem is usually caused by the irregular different scale of signals, so the signals at the border of the mode-mixing positions are really the non-stationary signals. Therefore, the criterion of determining the needed instantaneous frequencies is effective regardless of whether mode-mixing exists.

### 3.2.2 Locating local stationary signals with the selected instantaneous frequencies

After a line of fringe pattern being decomposed into numbers of IMFs with EMD, the first fundamental IMF *m*_{k1}(*x*) is determined firstly as analyzed above, and then the maxima and the minima of instantaneous frequencies can be found in *m*_{k1}(*x*) and *m*_{(k1-1)}(*x*), respectively, namely (*f*_{1})_{max} and (*f _{m}*)

_{min}in Eq. (2). (If

*K*1 = 1, the maxima and minima of instantaneous frequencies are both found in

*m*

_{k1}(

*x*).) Using these maxima and minima to locate local stationary signals is such a complex problem because they are not one-to-one correspondence. For example, as for the 1D signal analyzed in Fig. 2, Fig. 3(a) shows the instantaneous frequencies

*Insf*

_{1}(

*x*) for the last one high-frequency IMF

*m*

_{1}(

*x*) and Fig. 3(b) shows the instantaneous frequencies

*Insf*

_{2}(

*x*) for the first fundamental IMF

*m*

_{2}(

*x*). The positions of all the minima in Fig. 3(a) and the positions of all the maxima in Fig. 3(b) are marked by ‘*’. There are 194 minima and 93 maxima, respectively. The positions of these minima and maxima are not one-to-one correspondence, and meantime it’s unrealistic to select local signals in advance and then identify whether they meet Eq. (2).

Any two local stationary areas of signals are separated by a non-stationary signal whose instantaneous frequency is the maximum (or the minimum) but the maximum (or the minimum) doesn’t meet the Eq. (2) with its nearest minimum of instantaneous frequency (or its nearest maximum of instantaneous frequency). So if all the non-stationary signals can be found, the local stationary signals can be located easily. The following steps are given to find all the non-stationary signals.

**Step 1** Find all the minima and maxima of instantaneous frequencies for the last one high-frequency IMF *m*_{(}_{k}_{1-1)}(*x*) and the first fundamental IMF *m _{k}*

_{1}(

*x*), respectively. Denote the minima as the array${f}_{\mathrm{min}}[1,2,,q]$and the maxima as the array${f}_{\mathrm{max}}[1,2,,p]$, respectively. The coordinate positions of these minima and maxima relative to this line of original signal are recorded as$Q[1,2,,q]=[{x}_{1},{x}_{2},,{x}_{q}]$and$P[1,2,,p]=[{x}_{1}^{\text{'}},{x}_{2}^{\text{'}},,{x}_{p}^{\text{'}}]$, respectively.

**Step 2** If *p* ≤ *q*, let each element in${f}_{\mathrm{max}}[1,2,,p]$subtract the whole array${f}_{\mathrm{min}}[1,2,,q]$. Then a matrix$F[p][q]$is got, which is expressed as

If *p* > *q*, let each element of ${f}_{\mathrm{min}}[1,2,,q]$subtract the whole array${f}_{\mathrm{max}}[1,2,,p]$to get the matrix$F[q][p]$and set each positive value in$F[q][p]$to be zero.

**Step 3** In the matrix *F*[*p*][*q*] (or *F*[*q*][*p*]), find the only first line whose elements are all zero but its next line contains non-zero element, and write down the row-coordinate of this line in a new 1D array *l* and denote it as *l* [1] = ${x}_{{r}_{0}}^{\text{'}}$. If the first line is just the non-zero line, it means *r*_{0} = 0. (Fig. 4
is depicted to understand the **Step 3**~**Step 5**)

**Step 4** Beginning from the element$F({r}_{\text{0}}\text{+1,}1)$, look for the largest sub-matrix with all elements being zero from upper left to bottom right, and write down the coordinate denoted as $({x}_{{r}_{1}}^{\text{'}},{x}_{{c}_{1}})$of the sub-matrix’s last element$F({r}_{\text{1}}\text{,}{\text{c}}_{\text{1}})$ in the array *l* successively, that is *l* [2] = ${x}_{{r}_{1}}^{\text{'}}$ and *l* [3] = ${x}_{{c}_{1}}$. Then continue to look for the next biggest zero sub-matrix from $F({r}_{1}+1,{c}_{1}+1)$ along the direction from upper left to bottom right, get the coordinate$({x}_{{r}_{2}}^{\text{'}},{x}_{{c}_{2}})$of the last element for the second sub-matrix and also record$({x}_{{r}_{2}}^{\text{'}},{x}_{{c}_{2}})$in the array *l* successively. (During this process, if there is no zero sub-matrix but just single non-zero element along the direction, also record both of the element’s row-coordinate and column-coordinate in the array *l* successively.) Continue the same work until the last element of the sub-matrix is in the last row for the matrix *F*, and write down the row-coordinate${x}_{p}^{\text{'}}$ (or${x}_{q}^{\text{'}}$) and the column-coordinate as${x}_{{c}_{n}}$in the array *l* successively.

**Step 5** If the last element of the last zero sub-matrix in **Step 4** is not in the last column of the matrix *F*, go back to the (*r*_{0} + 1)th row of the matrix *F*, and continue the same work as the **Step 4** starting from the element$F({r}_{0}+1,{c}_{n}+1)$. The work is continued until the last element of a zero sub-matrix is at the last column of the matrix *F*, namely *F*(*r _{n}*,

*q*) (or

*F*(

*r*,

_{n}*p*)). And record the row-coordinate ${x}_{{r}_{n}}^{\text{'}}$and the column-coordinate

*x*(or

_{q}*x*).

_{p}**Step 6** The final array *l* is expressed as $l=[{x}_{{r}_{0}}^{\text{'}},{x}_{{r}_{1}}^{\text{'}},{x}_{{c}_{1}},{x}_{{r}_{2}}^{\text{'}},{x}_{{c}_{2}},,{x}_{p}^{\text{'}},{x}_{{c}_{n}},,{x}_{{r}_{n}}^{\text{'}},{x}_{q}]$where $0\le r\le p$and$1\le c\le q$(or$l=[{x}_{{r}_{0}}^{\text{'}},{x}_{{r}_{1}}^{\text{'}},{x}_{{c}_{1}},{x}_{{r}_{2}}^{\text{'}},{x}_{{c}_{2}},,{x}_{q}^{\text{'}},{x}_{{c}_{n}},,{x}_{{r}_{n}}^{\text{'}},{x}_{p}]$where$0\le r\le q$and$1\le c\le p$). Order the elements in *l* from small to big and merge the same ones. The final elements are the coordinates of the non-stationary signals which can divide the original signal sequence into several local stationary areas. If there are no elements in the final *l*, the whole line is thought as global stationary signals.

From Figs. 3(a) and 3(b), it is got that *q* and *p* in **Step 2** are 194 and 93, respectively. The local stationary areas located using the presented steps above is shown in Fig. 5
. Then the Gaussian window can be determined easily by the length of each local stationary area with Eq. (7). It is shown in Fig. 5 the changing Gaussian window centered at the 229th, the 445th and the 829th pixel, respectively. It’s clear the window width can change adaptively according to the changing of signals.

#### 3.3 Determining background through analysis of HHT

Another condition discussed in using AWFT is expressed in Eq. (3). This condition is set for avoiding the interference of zero-spectrum from the fundamental spectrum. However, the most effective means ought to be the way eliminating the background directly.

It has been illustrated that EMD is an effective method of eliminating background from fringe pattern in [21], but how to determine *K*2 in Eq. (15) is still a problem. We present another criterion to determine *K*2 as below:

**mean**{…} is the operation of computing the mean value of instantaneous frequencies for the

*n*th IMF and

**min**{…} is the operation of extracting the ordinal

*n*for the smallest mean value. The mean values of instantaneous frequencies should have been arranged from large to small as the arrangement of IMFs. However, some of them may become larger contrarily for the low-frequency IMFs. This situation is due to a certain amount of false peak instantaneous frequencies which are produced by the zero-intensity components. The more zero-intensity components, the more false instantaneous frequencies there are. If the zero-intensity components contained in a low-frequency IMF reach a certain level, the mean value of instantaneous frequencies becomes larger contrarily. And it means this IMF almost contains no fundamental components and mixes with some background components.

Table 2
shows the corresponding mean values of *Insf*_{1}(*x*) ~*Insf*_{7}(*x*) in Fig. 2(c). According to Eq. (18), it can be got *K*2 = 8 and it means just the residual *r*(*x*) is determined as the background even if there are different level of false instantaneous frequencies in *Insf*_{1}(*x*), *Insf*_{2}(*x*), *Insf*_{3}(*x*), *Insf*_{4}(*x*) and *Insf*_{6}(*x*), respectively. Figure 6(a)
is the result got by subtracting the background with the presented method from *S*(*x*), Fig. 6(b) is the standard values of simulation without background 0.5 and Fig. 6(c) is the error got by detecting the difference between Fig. 6(a) and Fig. 6(b). It shows the biggest error is no bigger than 0.015.

#### 3.4 The flow chart of the proposed method

Figure 7
shows the flow chart of the proposed method (the left part) and the method of locating local stationary areas namely **Step 1~5** in subsection 3.2.2 (the right part within the dotted box) for one line of signals. Then the whole phase map can be obtained by the proposed method line by line.

## 4. Simulation and experiments

#### 4.1 Simulation

The simulated fringe pattern (1020 × 1020) is produced by

where*f*= 0.05,$\eta (x,y)$is the random noise and $\varphi (x,y)$is the phase produced by peaks function:

_{0}*S*(

*x*) in Eq. (17) actually, and the red box shows a random local area.

Instantaneous frequencies can be obtained with the ridge method by detecting the scale factors at the ridge of wavelet transform or S-transform firstly and then taking the reciprocal of the scale factors. Here the RWT method is taken to do a comparison with our method. Still take *S*(*x*) as an example, the instantaneous frequencies are obtained by RWT with complex Morlet function whose bandwidth and center frequency is 0.8 and 1, respectively, and the scale factors of Morlet are set from 1 to 90 with step 0.2. The function and these parameters are chosen based on the better results with a large number of experiments. In [10], it is difficult to determine the local area [*x*-Δ*x*, *x* + Δ*x*] firstly and then to identify whether this local area being stationary. So the proposed steps in Subsection 3.2.2 in this manuscript is referred to locate the local stationary signals but just detecting the *f*_{min}[1,2,……, *q*] and *f*_{max}[1,2,……, *q*] both from the only array of instantaneous frequencies obtained by the RWT method.

Figure 9(a) shows the located local stationary signals by the method in [10]. It can be seen the local stationary signals are located too widely to perform detailed local analysis for the severe deformation areas. That is because the instantaneous frequencies got by the ridge method are smoothed subjectively due to the assumption that the phase of fringe pattern changes slowly or linearly. In Fig. 9(b), the local stationary areas are located in detail with the deformation of the signals. Figures 10(a) and 10(b) show the whole map of scale factors got by RWT and our method, respectively. If a whole line of signals are identified as global stationary, the scale factors are set zero and then the color appears pure black. In Fig. 10(a), many lines are pure black which means the corresponding lines of signals are determined global stationary but actually they are not.

Figure 11 shows the wrapped phase of the local area shown in Fig. 8, which is got by the FT method, the RWT method, the RWT-MWFT method in [10] and our method, respectively. It can be seen the quality of the first three phase maps is such bad while the last one is better.

To test the accuracy of the proposed method quantitatively, the quality-guided flood-fill algorithm [22] is used to unwrapped the obtained wrapped phase, which is selected because it needs just one fringe pattern and it is robust relatively. The unwrapped phase of fringe pattern for reference plane is subtracted from the unwrapped phase of the deformed fringe pattern to get the 3D phase map. Figure 12 shows the comparison of the true phase and the phase obtained by the four methods for the 500th line in Fig. 8. Computing the difference between the true phase and the result in Fig. 12, the error of each method can be obtained. Figure 13 shows the error comparison between our method and the other three methods. The error of our method is no bigger than 0.4 while there is always large local errors in the other method. In Fig. 13(c), the large error is produced at the position of the first large deformation, and then the error is propagated during the phase unwrapping.

Figure 14 shows the error comparison of phase retrieval for different fringe patterns. The wrapped phase got by each method is also unwrapped by the quality-guided flood-fill algorithm. And the error is obtained by subtracting the real phase of simulation from the retrieved phase got by each method. The experiments are carried out considering about different level of deformation and noise for fringe pattern by adjusting parameter$\alpha $in Eq. (20) and $\eta $in Eq. (19), respectively. The compared fringe patterns are shown in Figs. 14(a)–14(d).

The statistical errors are shown in Table 3
, where${\mu}_{\text{err}}$ is the mean of errors,${\sigma}_{err}$ is the standard deviation of errors and${\gamma}_{\text{err}}$ is the maximum error. When the level of deformation is small namely$\alpha \text{=1}\text{.5}$, the instantaneous frequencies are smoothed got by RWT, so each line is determined global stationary and is used with the FT method directly when applying RWT-MWFT. That is why the result got by using FT is the same with the result got by RWT-MWFT. When the level of deformation becomes serious namely$\alpha \text{=}7$and $\text{|}\eta \text{|=}0.02$, the errors become very large at the location of large deformation by using FT and RWT-MWFT, respectively, which is due to the improper selecting of window width leading to the spectra aliasing. When using RWT, large errors also occur at the larger gradient phase points. As the noise is added severely with$\text{|}\eta \text{|=}0.08$, large error propagation exists in the method of FT, RWT and RWT-MWFT. But our method is more accurate comparatively. That may benefit from the proper locating of local stationary areas and the effective eliminating background. For example, when the deformation level is $\alpha \text{=}7$and the noise level is$\text{|}\eta \text{|=}0.02$, the mean error and the standard deviation of determining background with the proposed criterion is just 4.63 × 10^{−4} and 9.1 × 10^{−3}, respectively.

#### 4.2 Experiments

By projecting a sinusoidal fringe pattern on the foam board which is with complex surface and carries steep edge as shown in Fig. 15(a) , the gray deformed fringe pattern (700pixels × 1000pixels) is obtained by CCD as shown in Fig. 15(b). The same four comparison methods are taken to obtain the demodulated phase. Moreover, the result of the four-step phase-shift method [23] is taken as a standard to assess the four methods. Figure 16 shows the demodulated phase for the 380th line, namely the red line in Fig. 15(b). And Fig. 17 shows the whole maps of the demodulated phase. In Fig. 17(c), it can be seen that the phase at some areas are the same with the result in Fig. 17(a), which is because some lines of signals are thought as global stationary and are processed by FT directly. There are large errors at the edge of the object in Figs. 17(a), 17(b) and 17(c). Figure 17(d) shows the demodulated phase got by using our method. Although there are some of phase errors, the result of our method shows our method can provide more accurate result for the phase demodulation of deformed fringe pattern.

## 5. Conclusion

An accurate adaptive method of AWFT based on HHT is presented for phase demodulation of fringe pattern. By analyzing the maximum amplitude of marginal spectrum for each IMF got by HHT, the criterion is presented to determine the first fundamental IMF adaptively. Then the instantaneous frequencies of the determined IMF and its adjacent high-frequency IMF are selected to locate local stationary regions in order to obtain the corresponding scale factors further for using AWFT. The proposed criteria are effective no matter whether the mode-mixing problem exists during the processing of EMD. As shown in simulation, the instantaneous frequencies got by RWT are smoothed by the limitation of the ridge method and thus the derived window widths are useless for extracting the detailed phase accurately. But the instantaneous frequencies got by our method describe the detailed changes of deformed fringe pattern better.

The presented steps for locating local stationary areas are simple to be carried out without iterative calculation. So feasibility and efficiency of the AWFT method is increased greatly. With these steps, lengths of local stationary areas are determined adaptively and objectively. However, frequency resolution is poor leading to the severe spectrum aliasing which is caused by zero-spectrum for the small length of local stationary area. Therefore, the background components, which have been determined by analyzing the mean value of the instantaneous frequencies for each IMF during HHT, are eliminated to tackle this problem. This processing is very important for measuring the detailed phase in local region accurately especially when the phase change is complex or serious. It is worth mentioning that minority over-segmentation may exist if the mode-mixing problem appears during EMD and it will lead to the very small size of window (such as the right window in Fig. 5), but the eliminating of background as the above analysis will reduce the error effectively if such phenomenon appears. The time consumed by the presented method is longer than FT method but near to RWT-MWFT. But its robustness, adaptive ability and accuracy is much better, and it can give good result even if the fringe pattern contains discontinuities.

## Acknowledgments

This research is supported by the National Natural Science Foundation of P. R. China (51175081, 61107001), Natural Science Foundation of Jiangsu Province (BK2010058) and the Scientific Research Foundation of Graduate School of Southeast University.

## References and links

**1. **S. Gai and F. Da, “A novel phase-shifting method based on strip marker,” Opt. Lasers Eng. **48**(2), 205–211 (2010). [CrossRef]

**2. **M. Takeda, H. Ina, and S. Kobayashi, “Fourier-transform method of fringe-pattern analysis for computer-based topography and interferometry,” J. Opt. Soc. Am. **72**(1), 156–160 (1982). [CrossRef]

**3. **J. Zhong and H. Zeng, “Multiscale windowed Fourier transform for phase extraction of fringe patterns,” Appl. Opt. **46**(14), 2670–2675 (2007). [CrossRef] [PubMed]

**4. **Q. Kemao, “Windowed Fourier transform for fringe pattern analysis,” Appl. Opt. **43**(13), 2695–2702 (2004). [CrossRef] [PubMed]

**5. **J. Zhong and J. Weng, “Phase retrieval of optical fringe patterns from the ridge of a wavelet transform,” Opt. Lett. **30**(19), 2560–2562 (2005). [CrossRef] [PubMed]

**6. **S. Li, X. Su, and W. Chen, “Wavelet ridge techniques in optical fringe pattern analysis,” J. Opt. Soc. Am. A **27**(6), 1245–1254 (2010). [CrossRef] [PubMed]

**7. **S. Özder, Ö. Kocahan, E. Coşkun, and H. Göktaş, “Optical phase distribution evaluation by using an S-transform,” Opt. Lett. **32**(6), 591–593 (2007). [CrossRef] [PubMed]

**8. **M. Zhong, W. Chen, and M. Jiang, “Application of S-transform profilometry in eliminating nonlinearity in fringe pattern,” Appl. Opt. **51**(5), 577–587 (2012). [CrossRef] [PubMed]

**9. **S. Zheng, W. Chen, and X. Su, “Adaptive Windowed Fourier transform in 3-D shape measurement,” Opt. Eng. **45**(6), 063601 (2006). [CrossRef]

**10. **J. Zhong and J. Weng, “Generalized Fourier analysis for phase retrieval of fringe pattern,” Opt. Express **18**(26), 26806–26820 (2010). [CrossRef] [PubMed]

**11. **Z. Wang, J. Ma, and M. Vo, “Recent progress in two-dimensional continuous wavelet transform technique for fringe pattern analysis,” Opt. Lasers Eng. **50**(8), 1052–1058 (2012). [CrossRef]

**12. **S. Fernandez, M. A. Gdeisat, J. Salvi, and D. Burton, “Automatic window size selection in Windowed Fourier Transform for 3D reconstruction using adapted mother wavelets,” Opt. Commun. **284**(12), 2797–2807 (2011). [CrossRef]

**13. **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]

**14. **W. Gao and Q. Kemao, “Statistical analysis for windowed Fourier ridge algorithm in fringe pattern analysis,” Appl. Opt. **51**(3), 328–337 (2011). [CrossRef] [PubMed]

**15. **M. Takeda and K. Mutoh, “Fourier transform profilometry for the automatic measurement of 3-D object shapes,” Appl. Opt. **22**(24), 3977–3982 (1983). [CrossRef] [PubMed]

**16. **N. E. Huang, Z. Shen, S. R. Long, M. C. Wu, H. H. Shih, Q. N. 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. A **454**(1971), 903–995 (1998). [CrossRef]

**17. **S. Equis and P. Jacquot, “The empirical mode decomposition: a must-have tool in speckle interferometry?” Opt. Express **17**(2), 611–623 (2009). [CrossRef] [PubMed]

**18. **G. Rilling, P. Flandrin, and P. Goncalves, “On empirical mode decomposition and its algorithms,” in IEEE-EURASIP Workshop on Nonlinear signal and Image Processing, NSTP-03, GRADO (I) (2003).

**19. **G. Rilling and P. Flandrin, “One or two frequencies? The empirical mode decomposition answers,” IEEE Trans. Signal Process. **56**(1), 85–95 (2008). [CrossRef]

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

**21. **S. Li, X. Su, W. Chen, and L. Xiang, “Eliminating the zero spectrum in Fourier transform profilometry using empirical mode decomposition,” J. Opt. Soc. Am. A **26**(5), 1195–1201 (2009). [CrossRef]

**22. **S. Zhang, X. Li, and S. T. Yau, “Multilevel quality-guided phase unwrapping algorithm for real-time three-dimensional shape reconstruction,” Appl. Opt. **46**(1), 50–57 (2007). [CrossRef] [PubMed]

**23. **S. Gai and F. Da, “Fringe image analysis based on the amplitude modulation method,” Opt. Express **18**(10), 10704–10719 (2010). [CrossRef] [PubMed]