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 OSA
Phase demodulation is a critical step in 3D shape measurement. It is usually completed by using phase-shift method  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 . 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 . 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 . 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 . 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, that is, the high-order terms in Taylor series of the measured phase are neglected and 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 , and meanwhile it also has another constraint including the adaptive selection for mother wavelet or parameters as the same problem as 1D RWT . There are also many other methods pursuing the much robust and accurate phase retrieval method such as the ensemble empirical mode decomposition method  and the windowed Fourier ridge algorithm assisted with statistical analysis , 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
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 Q1 completely:15], (f1)max and (f1)min are the maximum and minimum instantaneous frequency for the fundamental spectrum component, respectively, and fb is the cut-off frequency for the zero-spectrum component.
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 . 1D AWFT performed on I (x) in Eq. (1) is expressed as
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 by10] 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
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 asEq. (8) is repeated, or h(x) is expressed as m1(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 as18, 19].
For any IMF mn (x), its Hilbert transform m ′n (x) is16], thus it means the probability that certain frequency energy exists in an IMF. The largerthe larger possibility certain frequency exists.
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 asEq. (1), the second group is the fundamental IMFs corresponding to and the third matches the background A(x). The three groups are separated by mk1(x) and mk2(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 (f1)max in Eq. (2) must exist in mk1(x) and (fm)min must exist in m(k1-1)(x).
mk1(x) is defined as the first fundamental IMF which contains the most fundamental components and contains rarely high-frequency components. It can be determined withEquation (16) describes that the instantaneous frequency at the maximum marginal spectrum for mK1(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 . In Eq. (16), f0 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 byFig. 2(b). Figure 2(c) shows the corresponding instantaneous frequencies Insf1(x)~Insf7(x) for m1(x)~m7(x), and the corresponding marginal spectrum h1(f)~h7(f) is shown in Fig. 2(d). The corresponding listed in Table 1 is computed according to Eq. (16). It’s clear thatis the smallest value and thus K1 is obtained as 2. As m2(x) is identified as the first fundamental IMF, (f1)max and (fm)min in Eq. (2) exist in Insf2(x) and Insf1(x), respectively.
In this example, the mode-mixing problem exists as shown in Fig. 2(b). m2(x) loses some fundamental components while these losses occur at the corresponding positions in m1(x) and m3(x) because of the completeness character of EMD . In Fig. 2(c), there are peak values at the mode-mixing positions in Insf2(x), which are produced by the tiny values (close to zero) in m2(x) at the corresponding positions when computing instantaneous frequencies. They are false instantaneous frequencies. However, the values at the corresponding positions in Insf1(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 mk1(x) is determined firstly as analyzed above, and then the maxima and the minima of instantaneous frequencies can be found in mk1(x) and m(k1-1)(x), respectively, namely (f1)max and (fm)min in Eq. (2). (If K1 = 1, the maxima and minima of instantaneous frequencies are both found in mk1(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 Insf1(x) for the last one high-frequency IMF m1(x) and Fig. 3(b) shows the instantaneous frequencies Insf2(x) for the first fundamental IMF m2(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(k1-1)(x) and the first fundamental IMF mk1(x), respectively. Denote the minima as the arrayand the maxima as the array, respectively. The coordinate positions of these minima and maxima relative to this line of original signal are recorded asand, respectively.
Step 2 If p ≤ q, let each element insubtract the whole array. Then a matrixis got, which is expressed as
If p > q, let each element of subtract the whole arrayto get the matrixand set each positive value into 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  = . If the first line is just the non-zero line, it means r0 = 0. (Fig. 4 is depicted to understand the Step 3~Step 5)
Step 4 Beginning from the element, look for the largest sub-matrix with all elements being zero from upper left to bottom right, and write down the coordinate denoted as of the sub-matrix’s last element in the array l successively, that is l  = and l  = . Then continue to look for the next biggest zero sub-matrix from along the direction from upper left to bottom right, get the coordinateof the last element for the second sub-matrix and also recordin 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 (or) and the column-coordinate asin 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 (r0 + 1)th row of the matrix F, and continue the same work as the Step 4 starting from the element. The work is continued until the last element of a zero sub-matrix is at the last column of the matrix F, namely F(rn, q) (or F(rn, p)). And record the row-coordinate and the column-coordinate xq (or xp).
Step 6 The final array l is expressed as where and(orwhereand). 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 , but how to determine K2 in Eq. (15) is still a problem. We present another criterion to determine K2 as below:
Table 2 shows the corresponding mean values of Insf1(x) ~Insf7(x) in Fig. 2(c). According to Eq. (18), it can be got K2 = 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 Insf1(x), Insf2(x), Insf3(x), Insf4(x) and Insf6(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
The simulated fringe pattern (1020 × 1020) is produced byFig. 8 is modulated bywith and added the random noise with amplitude 0~0.02. The red line denotes the 500th row of signals which is just the signal sequence 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 , 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 fmin[1,2,……, q] and fmax[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 . 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  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  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 parameterin Eq. (20) and 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 is the mean of errors, is the standard deviation of errors and is the maximum error. When the level of deformation is small namely, 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 namelyand , 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, 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 and the noise level is, 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.
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  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.
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.
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]
9. S. Zheng, W. Chen, and X. Su, “Adaptive Windowed Fourier transform in 3-D shape measurement,” Opt. Eng. 45(6), 063601 (2006). [CrossRef]
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]
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]
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]