OSA's Digital Library

Optics Express

Optics Express

  • Editor: C. Martijn de Sterke
  • Vol. 20, Iss. 16 — Jul. 30, 2012
  • pp: 18459–18477
« Show journal navigation

Phase demodulation using adaptive windowed Fourier transform based on Hilbert-Huang transform

Chenxing Wang and Feipeng Da  »View Author Affiliations


Optics Express, Vol. 20, Issue 16, pp. 18459-18477 (2012)
http://dx.doi.org/10.1364/OE.20.018459


View Full Text Article

Acrobat PDF (2078 KB)





Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

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 OSA

1. Introduction

Phase demodulation is a critical step in 3D shape measurement. It is usually completed by using phase-shift method [1

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

] 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

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]

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

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]

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

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

]. The ridge of Wavelet transform (RWT) method [5

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

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]

] or the ridge of S-transform method [7

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

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]

] 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

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

]. The methods of multiscale windowed Fourier transform (MWFT) in [3

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]

, 10

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

] 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φ(x)=φ(b)+φ'(b)(xb), that is, the high-order terms in Taylor series of the measured phase φ(x)are neglected and φ'(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

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]

], and meanwhile it also has another constraint including the adaptive selection for mother wavelet or parameters as the same problem as 1D RWT [12

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]

]. There are also many other methods pursuing the much robust and accurate phase retrieval method such as the ensemble empirical mode decomposition method [13

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]

] and the windowed Fourier ridge algorithm assisted with statistical analysis [14

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]

], 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
I(x)=A(x)+B(x)cos[2πf0x+f(x)]+η(x),
(1)
where A(x) and B(x) represent the background and the amplitude modulation, respectively, f0 is the fundamental frequency,ϕ(x)is the evaluated phase modulation and η(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
Fig. 1 Condition for extracting the fundamental spectrum completely.
, it can be seen two conditions should be met to extract the fundamental spectrum Q1 completely:
f0<(f1)max<(fm)min,(m=2,3,...)
(2)
andfb<(f1)min<f0,
(3)
where the local spatial frequency fm is analogous to an instantaneous frequency and (fm)min is the minimum instantaneous frequency for the mth spectrum component [15

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]

], (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 [9

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

]. 1D AWFT performed on I (x) in Eq. (1) is expressed as
AWFTab(f)=[I(x)gab(x)exp(i2πfx)]dx,
(4)
where gab(x)=1ag(xba)=12πaexp[12(xba)2] is the Gaussian window selected for its best statistical properties in either time domain or frequency domain. To balance the time-frequency resolution, the center of Gaussian window changes with moving shift factor b and the window width changes narrow or wide with changing scale factor a. The Gaussian function has the statistic characteristic such as
gab(x)=1ag(xba)dx=1.
(5)
So it follows
AWFTab(f)db=[I(x)exp(i2πfx)ga(xb)db]dx=FT(f).
(6)
It means the sum of the fundamental spectrum in each Gaussian window equals to the fundamental spectrum just got by Fourier transform. So the phase of fringe pattern can be obtained by performing inverse Fourier transform on the sum of the fundamental spectrum.

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
a(x)=L(x)4,
(7)
where L(x) is defined as the full width at 1/e2 for the Gaussian window [10

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

] 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

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

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]

]: 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
h(x)=v(x)v(x)max_envelope+v(x)min_envelope2,
(8)
where v(x) is the arbitrary time series, v(x)max_envelope and v(x)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 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 as
v(x)=n=1Nmn(x)+r(x),(1nN)
(9)
where mn(x) is the nth 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

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

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]

].

For any IMF mn (x), its Hilbert transform mn (x) is
m'n(x)=1πPmn(τ)xτdτ,
(10)
where P is the Cauchy principal value. Then the analytic signal of mn (x) is described as
zn(x)=mn(x)+jm'n(x)=λn(x)ejφn(x),
(11)
where λn(x)={mn2(x)+[mn'(x)]2}1/2 is the amplitude of the analytic signal and φn(x)=arctan[mn'(x)/mn(x)] is the phase. And instantaneous frequency is defined as
Insfn(x)=12πdφn(x)dx.
(12)
After performing Hilbert transform on each IMF component, the original signal can be expresses as
v(x)=n=1Nλn(x)exp{j2π[Insfn(x)]dx}.
(13)
Here the residual is not considered temporarily. The time-frequency distribution of the amplitude is designated as Hilbert spectrum, and then the marginal spectrum for an IMF component is obtained by
hn(f)=H(f,x)dx=|λn(x)exp{j2π[Insfn(x)]dx}|dx,
(14)
whereH(f,x) represents the Hilbert spectrum of an IMF. The marginal spectrum offers a measure of total amplitude contribution from each instantaneous frequency [16

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]

], thus it means the probability that certain frequency energy exists in an IMF. The largerhn(f)the 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 as
I(x)=n=1k11mn(x)+n=k1k21mn(x)+[n=k2Nmn(x)+r(x)],(1nN)
(15)
where the first group of mn(x) is the high-frequency IMFs corresponding to η(x) in Eq. (1), the second group is the fundamental IMFs corresponding toB(x)cos[2πf0x+ϕ(x)] 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 with
k1=min{γ(n)}=min{|Insfnhmaxf0|},(1nN)
(16)
where hmax is the maximum marginal spectrum in the 2D distribution of marginal spectrum and frequency, Insfnhmaxis the corresponding instantaneous frequency at hmax for the nth IMF, f0 is the fundamental frequency and min{…} denotes the operation for extracting the ordinal n whenγ(n)is the smallest. Equation (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 [20

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]

]. 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)
Fig. 2 (a) Simulated signal S(x). (b) The IMFs m1(x) ~m7(x) and the residual r(x) after EMD. (c) Instantaneous frequencies Insf1(x) ~Insf7(x) for each IMF. (d) Marginal spectra h1(f) ~h7(f) for each IMF.
shows the signal got by
S(x)=0.5+0.25cos[2π×0.05x+ϕ(x)]+η(x),(0x1020)
(17)
where ϕ(x) is the simulated phase and η(x) is the random noise with amplitude 0~0.02. After performing EMD on S(x), m1(x)~m7(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 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 γ(n) listed in Table 1

Table 1. Corresponding γ for each IMF

table-icon
View This Table
| View All Tables
is computed according to Eq. (16). It’s clear thatγ(2)is 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 [16

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]

]. 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)
Fig. 3 (a) The minima in Insf1(x). (b) The maxima in Insf2(x).
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 arrayfmin[1,2,,q]and the maxima as the arrayfmax[1,2,,p], respectively. The coordinate positions of these minima and maxima relative to this line of original signal are recorded asQ[1,2,,q]=[x1,x2,,xq]andP[1,2,,p]=[x1',x2',,xp'], respectively.

Step 2 If pq, let each element infmax[1,2,,p]subtract the whole arrayfmin[1,2,,q]. Then a matrixF[p][q]is got, which is expressed as
F[p][q]=[fmax[1]fmin[1],fmax[1]fmin[2],,fmax[1]fmin[q]fmax[2]fmin[1],fmax[2]fmin[2],,fmax[2]fmin[q]fmax[p]fmin[1],fmax[p]fmin[2],,fmax[p]fmin[q]].
Set each negative value in the matrix to be zero.

If p > q, let each element of fmin[1,2,,q]subtract the whole arrayfmax[1,2,,p]to get the matrixF[q][p]and set each positive value inF[q][p]to be zero.

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 elementF(r0+1,cn+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(rn, q) (or F(rn, p)). And record the row-coordinate xrn'and the column-coordinate xq (or xp).

Step 6 The final array l is expressed as l=[xr0',xr1',xc1,xr2',xc2,,xp',xcn,,xrn',xq]where 0rpand1cq(orl=[xr0',xr1',xc1,xr2',xc2,,xq',xcn,,xrn',xp]where0rqand1cp). 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
Fig. 5 The located local stationary areas for S(x) in Fig. 2(a) and any Gaussian window centered at different pixel.
. 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

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]

], but how to determine K2 in Eq. (15) is still a problem. We present another criterion to determine K2 as below:
K2=min{mean{Insfn(x)}}+1,
(18)
where mean{…} is the operation of computing the mean value of instantaneous frequencies for the nth 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

Table 2. The mean values of Insf1(x) ~Insf7(x) in Fig. 2(c)

table-icon
View This Table
| View All Tables
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)
Fig. 6 (a)The result got by eliminating the background by the presented method. (b) The standard value of simulation without background. (c) The error of the presented method compared with the standard value.
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
Fig. 7 The flow chart of the proposed method and the method of locating local stationary areas.
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
I(x,y)=0.5+0.25cos[2πf0x+ϕ(x,y)]+η(x,y),
(19)
where f0 = 0.05,η(x,y)is the random noise and ϕ(x,y)is the phase produced by peaks function:
ϕ(x,y)=α{3(1-x)2exp[-x2-(y+1)2]-10(x5-x3-y5)exp(-x2-y2)-13exp[-(x+1)2-y2]},
(20)
where α is the parameter controlling the level of deformation for fringe pattern. The deformed fringe pattern shown in Fig. 8
Fig. 8 Simulated fringe pattern.
is modulated byϕ(x,y)with α=7 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 [10

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

], it is difficult to determine the local area [xx, 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)
Fig. 9 The located local stationary areas of S(x) got by: (a) the RWT method; (b) our method.
shows the located local stationary signals by the method in [10

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

]. 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)
Fig. 10 The respective maps of scale factors got by: (a) the RWT method; (b) our method.
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
Fig. 11 The wrapped phase of the local area got by: (a) the FT method, (b) the RWT method, (c) the RWT-MWFT method, (d) our method.
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

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

] 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

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]

] 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
Fig. 12 The comparison between the true phase and the retrieved phase for S(x) got by the FT method, the RWT method, the RWT-MWFT method and our method.
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
Fig. 13 Errors for phase demodulation of S(x) got by: (a) the FT method and our method; (b) the RWT method and our method; (c) the RWT-MWFT method and our method.
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
Fig. 14 Comparison of errors got by FT, RWT, RWT-MWFT and our method: (a)~(d): simulated fringe pattern with different level of deformation and noise; (a1)~(a4): errors for (a) with four methods respectively; (b1)~(b4): errors for (b) respectively; (c1)~(c4): errors for (c) respectively; (d1)~(d4): errors for (d) respectively.
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αin 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

Table 3. Error statistics for different methods

table-icon
View This Table
| View All Tables
, whereμerr is the mean of errors,σerr is the standard deviation of errors andγerr is the maximum error. When the level of deformation is small namelyα=1.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α=7and |η|=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|η|=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 α=7and the noise level is|η|=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)
Fig. 15 (a) The measured object. (b) The obtained fringe pattern.
, 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

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

] is taken as a standard to assess the four methods. Figure 16
Fig. 16 The retrieved phase for the 380th line of fringe pattern got by using five methods, respectively.
shows the demodulated phase for the 380th line, namely the red line in Fig. 15(b). And Fig. 17
Fig. 17 Restored 3D phase map got by: (a) FT; (b) RWT; (c) RWT-MWFT; (d) our method; (e) phase-shift method.
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]

OCIS Codes
(100.5070) Image processing : Phase retrieval
(110.2960) Imaging systems : Image analysis
(120.2650) Instrumentation, measurement, and metrology : Fringe analysis

ToC Category:
Image Processing

History
Original Manuscript: May 30, 2012
Revised Manuscript: July 18, 2012
Manuscript Accepted: July 19, 2012
Published: July 27, 2012

Citation
Chenxing Wang and Feipeng Da, "Phase demodulation using adaptive windowed Fourier transform based on Hilbert-Huang transform," Opt. Express 20, 18459-18477 (2012)
http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-20-16-18459


Sort:  Author  |  Year  |  Journal  |  Reset  

References

  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. A27(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. Express18(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. A454(1971), 903–995 (1998). [CrossRef]
  17. S. Equis and P. Jacquot, “The empirical mode decomposition: a must-have tool in speckle interferometry?” Opt. Express17(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. A26(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. Express18(10), 10704–10719 (2010). [CrossRef] [PubMed]

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.


« Previous Article  |  Next Article »

OSA is a member of CrossRef.

CrossCheck Deposited