Carrier peak isolation from single interferogram using spectrum shift technique Satoshi Tomioka,1,* Shusuke Nishiyama,1 and Samia Heshmat1,2 1

Faculty of Engineering, Hokkaido University, Kita 13 Nishi 8, Kita-ku, Sapporo, Hokkaido 060-8628, Japan 2

Faculty of Engineering, Aswan University, Aswan 81542, Egypt *Corresponding author: [email protected]

Received 7 April 2014; revised 23 July 2014; accepted 24 July 2014; posted 28 July 2014 (Doc. ID 209581); published 22 August 2014

This paper presents a new method to obtain a wrapped phase distribution from a single interferogram with a spatial carrier modulation. The Fourier transform of the interferogram has three peaks: one is a dc peak around the origin in the Fourier domain, and the other two are carrier peaks that have information of phase modulation by an object placed in the interferometer. Since the wrapped phase can be evaluated by one of the two carrier peaks, the dc peak and the adjoint peak that is the other peak of two carrier peaks should be removed by filters. The proposed filtering process consists of two stages: dc peak filtering and adjoint peak filtering. A spectrum shift filter based on symmetrical characteristics of the spectrum is applied in both stages as a basic filter that can remove most of the undesired spectrum. An additional two filters are applied to remove the remaining spectrum. The new method can automatically isolate the carrier peak, even when the boundary of peaks is not very clear. Numerical evaluations of simulation data and experimental data demonstrate that the proposed method can successfully isolate the carrier peak. © 2014 Optical Society of America OCIS codes: (100.2650) Fringe analysis; (100.5070) Phase retrieval; (070.2615) Frequency filtering; (070.4790) Spectrum analysis; (100.5088) Phase unwrapping; (120.2650) Fringe analysis. http://dx.doi.org/10.1364/AO.53.005620

1. Introduction

Phase measurement of optical waves has been widely applied for the purpose of quantitative measurement of the two-dimensional optical distance distribution by interferometers. The phase information cannot be obtained directly, because the recorded interferogram involves the phase as a cosine function in which the argument is the phase (cos
ϕ). In order to retrieve ϕ from cos
ϕ, there are two difficult issues. The first is the problem of a multivalued function. Since the inverse function of a cosine function is a multivalued function, the result of inversion has the ambiguity of a multiple of 2π, of which principal value is called the wrapped phase. This problem can be solved by phase unwrapping techniques that 1559-128X/14/255620-12$15.00/0 © 2014 Optical Society of America 5620

APPLIED OPTICS / Vol. 53, No. 25 / 1 September 2014

compensate phase jumps between neighboring points. The second problem is that the sign of the phase cannot be determined by the local phase distribution. Many experimental configurations were proposed, e.g., [1–12]. From the use of a simple experimental setup and from the applicability of measurements for fast phenomena, it can be said that the Fourier transform method (FTM) [1–5] is appropriate. The FTM uses only a single fringe pattern, but the phase distribution of the object is superposed onto the carrier phase that generates a stripe-shaped background fringe pattern. The carrier phase can be added with a small angle between the reference wave and the object wave in the interferometers. The fringe pattern, i
r, in the FTM is given as i
r  a
r  b
r cos ϕs
r;

(1)

ϕs
r  ϕ
r  s · r;

(2)

where a
r and b
r depend only on the profile of the reference wave and the object wave, and s is the spatial carrier frequency generated by the tilted angle between the two waves. By introducing a complexvalued function, cˆ
r, that has the information of the phase distribution, ϕ
r, as 1 cˆ
r  b
rejϕ
r ; 2

(3)

Eq. (1) can be rewritten as i
r  a
r  cˆ
rejs·r  cˆ 
re−js·r ;

(4)

where the caret denotes a complex-valued function in an object domain to distinguish the real-valued function. When we define a Fourier transform as F
k ≜ F ffˆ
rg 

Z

fˆ
re−jk·r dr

(5)

the Fourier transform of Eq. (4) is written as I
k  A
k  C
k − s  C
−k − s;

(6)

where the asterisk denotes a complex conjugate. It should be noted that C
k represents the complex conjugate of the Fourier transform of cˆ
r; i.e., C
k  F fˆc
rg ≠ F fˆc
rg. The first term on the right-hand side is distributed around the origin in k space, and the other two are distributed around k  s. In later discussion, we name those three distributions “dc peak”, “carrier peak”, and “adjoint peak”, respectively. The phase information is included in the carrier peak and the adjoint peak. Once we can obtain the carrier peak, C
k − s, by filtering in k space, unwrapped continuous phase can be evaluated by applying operators of Fourier inverse transform, complex argument evaluation, carrier phase shifting, and phase unwrapping as follows: 1 F −1 fC
k − sg  cˆ
rejs·r  b
rej
ϕs·r ; 2 ϕsw
r ≜ arg
F −1 fC
k − sg  Wfϕ
r  s · rg;

(7) (8)

ϕw
r ≜ Wfϕsw
r − s · rg;

(9)

ϕ
r  Ufϕw
rg;

(10)

where Wfg and Ufg denote phase wrapping and unwrapping operators, respectively, and the carrier frequency, s, in the evaluation of Eq. (9) is estimated from an interferogram without the object. The key issue of this sequence of processes is the accurate extraction of the carrier peak by filtering

out both the dc peak and the adjoint peak. Although the width of the dc peak is broadened in the case in which the region of the obtained interferogram is limited [13–17], a
r is generally smooth and the spectrum, A
k, does not have a high frequency component. The other two terms, C
k − s and C
−k − s, are the shifted spectrum by the carrier frequency, s, with and without axis inversion of C
k. When b
r has a similar profile to a
r, the spread of C
k is determined by F fejϕ
r g. Therefore, the condition of peak isolation is given by the relation between the gradient of the phase and the spatial carrier frequency as follows: jsj > maxfj∇ϕ
rjg:

(11)

Since the maximum gradient corresponds to the highest local frequency, the spectrum spread around s is represented by the maximum gradient. The carrier frequency, s, can be controlled by adjusting the tilted angle between two waves in interferometers to satisfy Eq. (11). In some cases, the adjustment of the angle is difficult due to changes in the experimental configuration during the measurement, e.g., phase tomography techniques to reconstruct the three-dimensional refractive index distribution [18] and low coherence interferometry to measure the depth distribution of reflection [19]. When Eq. (11) is satisfied, the spectrum of the dc peak can be removed by a simple Gaussian filter, since the three peaks are not overlapping each other. And a half-plane filter is applicable to remove the adjoint peak in which the boundary of the half-plane is determined by Bone et al.’s method [3]. In the case jsj ≳ j∇ϕj, the tails of the peaks overlap. For such a case, the filtering window with various shapes [20–24] was applied instead of the half-plane filter. However, selecting the right window size can be a drawback to these windows. Knowledge of the phase profile is essential to choosing an appropriate window size or the shape of the window, since it affects the retrieved phase. In some worse cases, closed fringes may appear. The proposed method is applicable without critical errors for the partially closed fringe pattern that cannot be analyzed by the past spectral filters; e.g., the dc peak does not overlap the others but the carrier peak and the adjoint peak overlap each other as shown in Fig. 4, and the dc peak is included in the overlapped region of the other peaks as shown in Fig. 5. In the case of s  0, which is not categorized into the FTM, the dc peak and twin peaks carrying the phase information are completely overlapping, and most of the fringes are closed. Hence, the phase cannot be retrieved by using the standard Fourier transform. This problem can be solved by using the methods proposed in the past [25–34]. From an analytical point of view, these methods are superior to the standard FTM because they can be solved even when the fringes are closed. However, most of them require several tests for different parameters 1 September 2014 / Vol. 53, No. 25 / APPLIED OPTICS

5621

(e.g., window size, the number of polynomials, and the number of terms for continuity of derivatives), or many evaluations of cost functions. Furthermore, the experimental adjustment to achieve s  0 is difficult, i.e., s ≠ 0 with general optical configurations. Therefore, the FTM-based approaches to isolate the carrier peak for jsj ≳ j∇ϕj are still useful. The filter proposed in this paper is applied in the k space similarly to the half-plane filter and the other spectral domain filters. However, knowledge of neither the size nor the shape of the window in the spectrum domain is necessary, since the window is computed automatically. The filtering consists of two stages; the first stage is to remove the dc peak, and the second stage is to remove the adjoint peak. In both stages, a frequency shifting of the spatial spectrum is applied as a basic filter. Since this filter cannot remove the undesired component completely, other filters are also applied. These filters are described in Section 2. After that, numerical examples to demonstrate the validity of these filters are shown in Section 3, and the conclusion is shown in Section 4. 2. Spectral Peak Isolation by Spectrum Shift A.

Basic Idea

As shown in Eq. (6), the Fourier transform of the observed fringe pattern, I
k, is given by the sum of three distributions of A
k, C
k − s, and C
−k − s, of which the centers are located at k  0, s, and −s, respectively. The known function is only I
k, and the others are unknown. When we ignore products between different terms such as jA
kC
k − sj, the power spectrum, P
k, is approximated as P
k ≜ jI
kj2 ≃ PA
k  PC
k − s  PC
−k − s; (12) where PA and PC are the power spectra of unknown functions of A and C, respectively. By replacing k with k ∓ 2s, the shifted power spectra are written as P
k − 2s ≃ PA
k − 2s  PC
k − 3s  PC
−k  s; (13)

PC
−ks  PC
k  s:

(17)

Evaluation of each term in Eq. (15) with substitutions of Eqs. (12)–(14) and the above symmetrical nature gives D1
k ≃ PA
k − PA
k − 2s − PA
k  2s − PC
k − 3s − PC
−k − 3s:

(18)

Each term except PA
k in the right-hand side is negative, since the power spectra are positive functions. Hence, taking positive value of Eq. (18) as PosfD1
kg ≃ PA
k;

(19)

Posff g ≜ maxff ; 0g;

(20)

we can obtain PA
k. Figure 1 shows the schematic of the above process. Subtracting PA
k from P
k, we can evaluate the twin peaks that have the information of the phase modulation: Ptwin
k  PC
k − s  PC
−k − s. By employing a similar procedure, a single carrier peak, PC
k − s, can be isolated. The treatment of the spectral domain of the shifted spectrum should be noted. In the discrete Fourier transform, the spectral domain is limited to the finite size determined by the sampling theorem. By applying the spectrum shift, some area is shifted to outside the limited domain, and some area is shifted into the limited domain. Since the signal outside the domain can be considered as zero, the spectrum at the area shifted into the domain can be set as zero. Meanwhile, the spectrum at the area shifted outside the domain is not required, because the evaluation point, k, in Eq. (15) is still located in the original limited domain. In addition, we must consider the discrete grid in the spectral domain because the carrier frequency, s, is not located on the grid. In such a case, the shifted spectra shown in Eqs. (13) and (14) can be evaluated from i
r by using Fourier transform as P
k  2s  jF fi
re2js·r gj2 :

(21)

P
k  2s ≃ PA
k  2s  PC
k  s  PC
−k − 3s: (14) These functions include a term of which peak is located at k  s. Here, in order to isolate PA
k we define the difference, D1
k, between P
k and the sum of the shifted power spectra: D1
k ≜ P
k −
P
k − 2s  P
k  2s:

(15)

Consider a simple case in which PA
k and PC
k are symmetrical functions with respect to the origin; namely, PA
−k  PA
k; 5622

(16)

APPLIED OPTICS / Vol. 53, No. 25 / 1 September 2014

Fig. 1. Scheme of dc peak isolation by spectrum shift technique.

The description in this subsection is based on the assumption that the spectra have a symmetrical nature as in Eqs. (16) and (17). Otherwise, this evaluation contains an error. B.

Using these expressions, the power spectra of PC
k and PC
−k have the following relations:

Symmetrical Property of Spectrum

In the case in which a function, f
r, is a real-valued function, the Fourier transform at the symmetrical point −k is evaluated by replacing k with −k in the definition of Fourier transform shown in Eq. (5). By using the relation of f 
r  f
r, the following symmetrical properties are obtained: F
−k  F 
k;

(23)

Since PA
k at the first term in the right-hand side of Eq. (12) is a power spectrum of the real-valued function, a
r, it has these symmetrical properties. The sum of the second term and the third term, Ptwin
k  PC
k − s  PC
−k − s, is also symmetric with respect to the origin. However, PC
k − s themselves are not symmetric with respect to s as shown below, because C
k is the Fourier transform of a complex-valued function, cˆ
r. A complex-valued function, cˆ
r, can be represented by a sum of two real-valued functions that correspond to a real and an imaginary part, respectively, as cˆ
r  cr
r  jci
r


ˆc ∈ C; cr ; ci ∈ R:

(24)

By denoting the Fourier transforms of each term as Cr
k ≜ F fcr
rg;

Ci
k ≜ F fci
rg;

(25)

we obtain the following relations: C
k  Cr
k  jCi
k;

(26)

C
k  Cr
k − jCi
k;

(27)

C
−k  Cr
k  jCi
k:

(28)

The Fourier transform of the symmetrical point, C
−k, is not equal to either C
k or its complex conjugate, C
k. The power of PC
k, therefore, is an asymmetric function: PC
−k ≠ PC
k:

(29)

The Fourier transforms of Cr
k and Ci
k in a polar coordinate system are written as jΘr
k

Cr
k  jCr
kje

;

Ci
k  jCi
kjejΘi
k :

(30)

(31)

PC
−k  PeC
k − PoC
k;

(32)

where PeC and PoC are defined as PeC
k ≜ jCr
kj2  jCi
kj2 ≥ 0; PoC
k ≜ 2jCr
kjjCi
kj sin
Θr
k − Θi
k:

(22)

PF
−k  PF
k:

PC
k  PeC
k  PoC
k;

(33) (34)

Note that PoC may have negative values. In Eqs. (33) and (34), both jCr j and jCi j are even functions, because both Cr and Ci are Fourier transforms of real-valued functions of cr and ci , respectively. Furthermore, we can find that the functions Θr
k and Θi
k are odd functions by considering Eqs. (22) and (30). As a result, the functions PeC and PoC are the even component and the odd component of PC, respectively: PeC
−k  PeC
k; PoC
−k  −PoC
k: C.

(35)

Isolation of dc Peak

1. Rough Estimation of dc Peak by Spectrum Shifting Technique From the symmetrical properties shown in the previous subsection, the distribution of the power spectrum, P
k, in Eq. (12) can be illustrated as shown in Fig. 2(a). It consists of three blocks: the dc power spectrum around the origin, the carrier power spectrum depicted by the three boxes at kx > 0, and the adjoint spectrum in kx < 0. The carrier frequency, s, is at the center of the middle box of the carrier spectrum. The carrier spectrum is not a symmetrical function with respect to s. Using the even and the odd spectra defined in Eqs. (31) and (32), the power around s can be represented as follows: PC
ks  PeC
ks  PoC
ks;

(36)

PC
−k  s  PeC
ks − PoC
ks:

(37)

Furthermore, the symmetry of PeC and the antisymmetry of PoC shown in Eq. (35) with respect to k  s can be represented as PeC
−k  s  PeC
ks;

(38)

PoC
−k  s  −PoC
ks:

(39)

By applying these relations to Eqs. (12)–(14), and computing D1
k defined in Eq. (15), we obtain 1 September 2014 / Vol. 53, No. 25 / APPLIED OPTICS

5623

(a)

(b)

ky

ky

narrower. Therefore, it is required to prepare an additional filter shown in the next subsection. The rough estimation of the dc component, which will be modified by the additional filter in the later stage, can be evaluated as

kx

kx

A†
k  (c)

ky

(d)

ky

kx

kx

Fig. 2. Scheme of two-dimensional dc peak isolation: (a) P
k, (b) D1
k by spectrum shifting, (c) PosfD1
kg, and (d) jA0
kj2 by a single connected support region filter. The number of small circles represents the spectrum density. The filled small circles and open ones express positive and negative spectrum, respectively. The enclosed circles represent dc peaks. Triplets of squares represent asymmetrical carrier peaks in which the carrier frequency is positioned at the center of each triplet.

− PC
−k − 3s2PoC
k − s  2PoC
−k − s: (40) It is found that the terms of PoC
k − s are added to the symmetrical case shown in Eq. (18). Unfortunately, it cannot be confirmed that these added terms are the negative functions. Even if PoC
k − s < 0, the values at the symmetrical point around s satisfy PoC
−k  s > 0 due to Eq. (39). This means that there are no cases in which the odd function around s takes a negative value at every point. Therefore, the positive component of D1
k includes the positive component of PoC around s: PosfD1
kg ≃ PA
k  2PosfPoC
k − sg 2PosfPoC
−k − sg:

(41)

Figures 2(b) and 2(c) illustrate these operations. Two boxes including single dots in Fig. 2(c) correspond to the second and the third terms in Eq. (41). In order to remove the undesired terms in Eq. (41), we should apply further filtering. One idea on the further filtering is the consideration of a weight of the subtracted terms. By denoting the weight of reduced terms as γ, the definition of the difference is modified as follows: Dγ
k  P
k − γP
k − 2s  P
k  2s:

(42)

When γ > 1, the spectrum of the second and the third terms in Eq. (41) becomes smaller. However, if γ is too large, the width of the estimated spectrum becomes 5624

APPLIED OPTICS / Vol. 53, No. 25 / 1 September 2014

(43)

where † shows the rough estimation. The factor multiplying I
k takes a value in the range [0, 1], which can be derived from Eq. (42). This filtering is equivalent to the so-called Wiener filter [35]. In the Wiener filter, an observed spectrum, O
k, is represented as O
k  S
k  N
k, where S
k and N
k represent a true unknown signal and the noise of its power spectrum, jN
kj2 , is known. The filtered sig~ nal, S
k, is evaluated by the following equation: ~ S
k  Ψ
kO
k;

Ψ
k 

(44)

jS
kj2 jS
kj2  jN
kj2 ≃

D1
k ≃ PA
k − PA
k − 2s − PA
k  2s − PC
k − 3s

PosfDγ
kg I
k; P
k

PosfjO
kj2 − jN
kj2 g ≜ ΨfjO
kj2 ; jN
kj2 g: jO
kj2 (45)

In the estimation of A†
k in Eq. (43), PosfDγ
kg corresponds to the power spectrum after subtracting the power spectrum of the noise, and P
k corresponds to the power spectrum of the observed signal. 2. Single Connected Support Region Filter As a further filtering, we propose a method to determine the dc peak region by controlling the threshold level. The roughly estimated spectrum, A†
k, consists of a few major peaks: the dc peak and the undesired residuals of the carrier peaks. In addition, the white noise spectrum spreads randomly over the whole spectral space as an uneven distribution with low intensity. When the threshold level is sufficiently high, the uneven background goes under the threshold level, and the major peaks are isolated like islands; however, the tail of the dc peak also goes under the threshold level, and the width of the dc peak becomes narrow. In contrast, if the threshold level is too low, the major peaks are connected to a single region. In order to obtain the isolated dc peak without narrowing of the tail, the major peaks except the dc peak should be reduced. In a comparison between Figs. 2(a) and 2(c), it can be understood that the spectrum at just s is completely removed, and PosfDγ
kg in the adjoining box to the box including s is smaller than the original one. The reason is that the even component, PeC , is completely removed. Therefore, the roughly estimated spectrum, A†
k,

enables us to set a lower threshold level than a threshold for the original spectrum, I
k. The threshold level, which is the average of the background of A†
k, can be obtained by an iterative method as shown below. In this process, the dc peak itself is considered as an outlier to estimate the background. The set of samples to estimate the background is denoted as L
j b , where j is an iteration number. An initial set of samples, L
0 b , is taken as the set of k with jA†
kj ≠ 0. The average, I b , and the standard deviations, σ b , are computed from the set of the samples: † I
j b  hjA
kjik∈L
j ; b

σ
j b 

r















































2 hjA†
kj2 ik∈L
j −
I
j b  ; b

(46)

(47)

(48)

After several iterations, the background average, I b , can be obtained. Among the points k with jA†
kj > I b , only the points connected to the origin constitute the set of points belonging to the dc peak, L0 . The points k ∈ L0 are obtained by using a flood fill algorithm used in the field of computer graphics. The roughly estimated dc spectrum in Eq. (43) for k ∈ L0 is modified by the Wiener filter shown in Eq. (45), where the noise is I 2b :  A0
k



ΨfjA†
kj2 ; I 2b gA†
k; k ∈ L0 ; 0; k ∉ L0 :

(49)

The filtered result is distributed around the origin as illustrated in Fig. 2(d). D.

1. Rough Estimation of Carrier Peak by Spectrum Shifting Technique The power distribution obtained by Eq. (50) is represented with the superposition of PC
k − s as P0twin
k ≃ PC
k − s  PC
−k − s;

where h ik denotes an average for the samples. When L
j b includes the points belonging to the dc peak, the average, I
j b , is larger than the background to be estimated. The points with large intensity, which belong to the dc peak, should be removed from Lb . The new set of samples can be updated by using
j I
j b and σ b :
j L
j1  fkj0 < jA†
kj < I
j b b  3σ b g:

where P0A
k  jA0
kj2 . The task in this stage is the isolation of the carrier peak, C
k, from the known functions of I 0twin
k and P0twin
k. Rough estimation of the carrier peak from the spectrum including twin peaks is similar to the case of dc peak isolation. However, the procedure to filter out the residual rough estimation is different from the case of the dc peak isolation process.

Isolation of Carrier Peak

Once the dc component, A0
k, is isolated, the power spectrum of twin peaks that is a sum of the carrier peak and the adjoint peak, I twin
k  C
k − s C
−k − s, can be easily evaluated, and the intensity of the twin peaks is evaluated by Wiener filter when A0
k is considered as noise, P0twin
k  P
k − P0A
k;

(50)

I 0twin
k  ΨfP
k; P0A
kgI
k;

(51)

(52)

which is illustrated in Fig. 3(a). Similar to the way to obtain the dc peak, the frequency shift to translate −2s is applied as P0twin
k  2s ≃ PC
k  s  PC
−k − 3s:

(53)

Subtraction of Eq. (53) from Eq. (52), and the use of properties of the power in Eqs. (36)–(39), gives us the following form: DsC
k ≜ P0twin
k − P0twin
k  2s ≃ PC
k − s  2PoC
−k − s − PC
−k − 3s:

(54)

Since the polarity of the second term cannot be determined, the positive component of this equation, Ps† C
k, includes an undesired spectrum as follows: s o Ps† C
k  PosfDC
kg ≃ PC
k − s  2PosfPC
−k − sg:

(55) Figures 3(b) and 3(c) illustrate these filtering steps. In the result shown in Fig. 3(c), the undesired part in kx < 0 remains, which corresponds to the second term on the right-hand side in Eq. (55). The estimated Fourier transform of the rough carrier, Cs†
k, is evaluated by applying the Wiener filter: Cs†
k 

PosfDsC
kg 0 I
k: P0twin
k twin

(56)

2. Symmetrical Point Comparison Filter In order to reduce the second term in Eq. (55) we propose a symmetrical point comparison filter. The adjoint peak spectrum cannot be completely removed by evaluation of Ps† C
k in Eq. (55); however, the adjoint peak spectrum is smaller than the original twin peak spectrum. This can be understood from the comparison between Figs. 3(a) and 3(c). Moreover, the carrier spectrum that should not be reduced by the spectrum shift filtering still remains. This means that we can determine which peak (carrier or adjoint) the arbitrary k belongs to by evaluation of the following quantity: 1 September 2014 / Vol. 53, No. 25 / APPLIED OPTICS

5625

interferogram obtained by an experiment in which the true phase distribution is unknown. A. Simulation

Fig. 3. Scheme of two-dimensional carrier peak isolation: (a) Ptwin
k, (b) DsC
k by spectrum shifting, (c) Ps† C
k, and 0 (d) jCs
kj2 by a symmetrical point comparison filter. The triplet of squares in kx > 0 and that in kx < 0 represent a carrier peak spectrum and an adjoint one, respectively. The carrier frequency, s, is positioned at the center of each triplet, and each spectrum is asymmetric with respect to the carrier frequency.

Δs
k ≜ PosfDsC
kg − PosfDsC
−kg:

(57)

The condition to determine the peak is Lc  fkjΔs
k > 0g;

(58)

La  fkjΔs
k < 0g;

(59)

where Lc and La are a set of points of the carrier peak and that of the adjoint peak, respectively. The above condition has an exception of the case in which Δs
k  0. In this case, the distances from the centers of two peaks, jk − sj and jk  sj, are evaluated, and the point with a shorter distance is selected. Every point belongs to either Lc or La. Furthermore, since the number of points belonging to each group is almost the same, the signal with high frequency can be represented by the sufficient number of points. The intensity of the carrier peak can be evaluated for the point set of the carrier peak as s0

C
k 



Cs†
k; k ∈ Lc ; 0; k ∉ Lc :

(60)

The spectrum after applying the symmetrical point comparison filter is illustrated in Fig. 3(d). 3. Numerical Results and Discussion

In order to demonstrate the validity of the proposed methods, three numerical results are shown in this section. Two of them are results from simulation data in which an interferogram is prepared from a given phase distribution; therefore, the unwrapped phase of filtering results can be compared with the given phase. The last one is a filtering for an 5626

APPLIED OPTICS / Vol. 53, No. 25 / 1 September 2014

1. Hash-Sign-Shaped Spectrum The first simulation is for the case in which the spectrum of the twin peaks surrounds the dc peak spectrum, but the twin peaks do not overlap with the dc peak. The simulation data shown in Figs. 4(a)–4(d) are prepared according to the definition shown in Appendix A. Since the assumed true phase distribution shown in Fig. 4(a) is given as a sum of two twodimensional elliptic Gaussian distributions with high aspect ratios, the Fourier transform shown in Fig. 4(b) spreads to higher frequencies as an oblique plus signed shape indicates. The image presented in Fig. 4(c) shows a fringe pattern that is a simulated result of the interferometer with a carrier modulation. The spectrum of Fig. 4(c) has a hash-signshaped distribution as shown in Fig. 4(d). The filtering process in the spectral domain by the spectrum shift method is shown in Figs. 4(e)–4(i). The first stage is the determination of a dc component by using two filters. The result shown in Fig. 4(e) is a roughly estimated dc spectrum by using the spectrum shift filter proposed in Subsection 2.C.1. We can see that the spectra around carrier peaks are not completely eliminated. Furthermore, the small spectrum spreads into a wide region like small islands. By applying the single connected support region filter described in Subsection 2.C.2, most of undesired spectra are filtrated out as shown in Fig. 4(f). The small tail is distributed around the origin. If we wish to remove the tail, we could control the filtering weight γ. However, the final results with γ  1; 2; …; 10 showed no significant differences for this simulated spectrum. The spectrum with twin peaks shown in Fig. 4(g) is computed by using the Wiener filter shown in Eq. (51). The second stage in the filtering process is the filtration for one of two carrier peaks. The result by using the spectrum shift filter proposed in Subsection 2.D.1 is shown in Fig. 4(h). In the result the adjoint carrier peak is reduced compared to the result shown in Fig. 4(g); however, undesired spectrum remains as speckles. 0 The isolated carrier peak, Cs
k, shown in Fig. 4(i) iϕ
r gj shown in Fig. 4(b), but the oriis similar to jF fe gin is shifted to k  s. By applying inverse Fourier 0 transform to Cs and taking the complex angle as Eqs. (7) and (8), we can obtain the wrapped phase with carrier shift as shown in Fig. 4(j). The image presented in Fig. 4(k) shows the wrapped phase after carrier removal by applying Eq. (9). The unwrapped phase, ϕ0
r, shown in Fig. 4(l) is the final result of the sequence in image processing. To compute the unwrapped phase, we applied a method that uses localized compensators [36]. The unwrapped phase agrees with the input true phase distribution, ϕ
r, given in Fig. 4(a). The images in Figs. 4(m)–4(o) are the results of applying the half-plane filter, which removes the

Fig. 4. Result of peak isolation and phase unwrapping for simulation data with a hash-sign-shaped spectrum. The input data for filtering are shown in (a)–(d); the images in (e)–(i) show the sequence of the filtering process; (j) and (k) show the wrapped phase with carrier modulation and that without modulation, respectively; and (l) is the unwrapped phase. The images shown in (m)–(o) are the results when using a half-plane filter, and those in (p)–(r) are the ones when using a circular bandpass filter. In the spectra, the displayed range for each axis is reduced to 50% of the range determined by the sampling theorem where the origin is the center of the image, and the brightness is converted to enhance the small spectra by an arc-tangent function. The filtering weight to obtain (e) is γ  1.

spectrum at the points where k · s < 0. Figure 4(m) shows the results by this filter, which is applied to the spectrum with twin peaks shown in Fig. 4(g). The unwrapped result shown in Fig. 4(o), ϕ0
rH , is different from the true phase shown in Fig. 4(a). In addition, the results by using a circular bandpass filter to filter out both the dc peak and the adjoint peak are shown in Figs. 4(p)–4(r). The unwrapped result shown in Fig. 4(r) is also different from the true phase. The root mean square errors between estimated ϕ0
r and ϕTrue
r by three filters were evaluated: 0.27, 7.6, and 4.3 rad for ϕ0
r, for ϕ0
rH, and for ϕ0
rC, respectively. From the above-mentioned discussion, it can be said that the spectrum, which cannot be filtered out by the region selecting method, can be filtered out successfully by the spectrum shift filter. 2. Partially Closed Fringe The second example is a case with a partially closed fringe. In its spectrum, the dc peak is included in the region where the carrier and the adjoint peaks are overlapping each other. Figure 5 demonstrates a result of filtering for both the spectrum shift filter and a half-plane filter. The fringe patterns, i
r, for the true phase with an isotropic Gaussian distribution were produced by simulations of an interferometer with various background fringes and noise. The definition is shown in Appendix A. The normalized

carrier frequency, kˆ bg , shows a carrier frequency normalized by the maximum gradient of the true phase: r


s e σ ϕ1  s: (61) kˆ bg ≜ 2 ϕ1 maxfj∇ϕTrue jg Closed fringe patterns are found with kˆ bg ≤ 1. When kˆ bg decreases, a closed pattern appears clearly in i
r, and overlaps of the twin peaks in jI
kj grow. In the case of the spectrum shift filter, the filtered 0 spectrum, jCs
kj, shows overfiltering of the carrier peak around the origin; however, the differences of phases (ϕw0
r and ϕ0
r) from the true phase are small. In contrast, the phases in the result by the half-plane filter have a large error at the portion where the fringe is closed. Figure 6 presents the error dependency on kˆ bg . The error for both filters increases as kˆ bg decreases. The errors by the spectrum shift filter are always smaller than those by the half-plane filter. Through this simulation, we can conclude that the proposed spectrum shift filter is superior to the halfplane filter for the partially closed fringe pattern. B. Experimental Fringe Pattern

In this subsection, we demonstrate a phase extraction for an experimental fringe pattern of which peaks are not clearly separated. Figure 7 shows the results in which the interferogram is measured 1 September 2014 / Vol. 53, No. 25 / APPLIED OPTICS

5627

Fig. 5. Result of peak isolation and phase unwrapping for simulation data with a partially closed interferogram. The assumed true , are shown at the top. The bottom tables show the simulated interferogram continuous phase, ϕTrue
r, and its wrapped phase, ϕTrue w as the input of filters, the results of the spectrum shift filter, and the results of the half-plane filter. Each row shows the dependency of the normalized carrier frequency, kˆ bg . The input image for each kˆ bg consists of three images: ibg
r (the background fringe), i
r (the 0 interferogram with noise), and jI
kj [the spectrum of i
r]. The results for both filters consist of three images: jCs
kj (the filtrated spectrum), ϕw0
r (the wrapped phase), and ϕ
r (the unwrapped phase). The scales of the brightness of each figure are shown at the top. In the images of the spectra, the displayed range for each axis is reduced to 25% of the range determined by the sampling theorem. The filtering weight in the spectrum shift filter is not considered, i.e., γ  1.

by using a Mach–Zehnder interferometer with a He– Ne laser (633 nm). The object is a horizontally placed burner flame. Since the burner flame fluctuates owing to the convection of heated gas, the exposure time is short as 0.2 ms. As a result, the obtained interferogram becomes noisy. Furthermore, since the dense background fringe is applied to avoid the closed fringe pattern, blurred portions are found in the interferogram shown in Fig. 7(b). Therefore, we only analyzed a portion with low fringe density determined by the box in the top figure in Fig. 7(b).

Fig. 6. Error of unwrapped phase for partially closed fringe simulation. The horizontal axis shows the normalized carrier frequency, and the vertical axis shows the root mean square of the error between the unwrapped phase and the true phase, i.e.,
hjϕ0
r − ϕ
rj2 ir 1∕2 . 5628

APPLIED OPTICS / Vol. 53, No. 25 / 1 September 2014

Figure 7(a) shows the schematic of the burner flame structure. The fuel gas is not burned at the exit of the burner nozzle, of which the unburned region is depicted by the tapered triangle in Fig. 7(a). When the fuel gas is burned, the flame can be observed; however, the boundary of the flame itself does not appear in the fringe pattern shown in Fig. 7(b). By the flame, temperatures of the fuel gas itself and those of surrounding gas rise, and the density of them decreases. The phase change depends on the density of the gases. The fuel gas density in the unburned region is higher than the density of unheated gas. As a result, the polarity of the phase gradient along the vertical direction is changed in the region around the unburned gas fuel, while it is not changed in the region apart from the nozzle. This means that choosing the half-plane filter to represent such a phase distribution is difficult. It can be found from the spectrum in Fig. 7(c) that the spectra corresponding to the carrier peak and the adjoint peak overlap. In the result by the spectrum shift filter shown in Fig. 7(e), the unwrapped phase, ϕ0
r, represents the higher gas density of the unburned gas fuel. In contrast, it does not show the increase of optical length owing to the unburned gas fuel in the result by the half-plane filter shown in Fig. 7(f). This is caused by the several small pieces of vertical patterns in the unburned region in ϕ0w
r. In the peak isolation for this experimental spectrum, of which the image size is 256 × 256, the computational time was almost 0.17 s with a PC equipped with an Intel Core 2 DUO CPU with 2.13 GHz clock in a single CPU operation mode.

Fig. 7. Phase extraction for an experimental fringe pattern. The fringe pattern is obtained by using an interferometer with an object that is flame of a horizontally placed gas burner. (a) Schematic view of image configuration. (b) Observed fringe pattern, i
r, with 512 × 512 pixels; the circle shows the laser beam profile with 50 mm diameter, and the bottom figure is a magnified image for the analyses with 256 × 256 pixels, whose region is determined by the box in the top image. (c) Spectrum, jI
kj, of (b). (d) Background fringe pattern without 0 object. (e) Filtered results by the spectrum shift filter; jCs
kj is the filtered spectrum, ϕ0w
r is the wrapped phase, and ϕ0
r is the unwrapped phase. (f) Filtered results by the half-plane filter. In the spectral images, the horizontal and vertical ranges are reduced to 25% of the range defined by the sampling theorem. The filtering weight in the spectrum shift filter is not considered, i.e., γ  1.

The computational time for different-sized images was also evaluated: 1.43 s for 254 × 254 size, 0.68 s for 512 × 512 size, and 10.5 s for 510 × 510 size. These differences are mainly caused by the difference of the algorithm of the Fourier transform to obtain three shifted spectra (two in the dc peak filtering, one in the carrier peak filtering), of which the method is shown in Eq. (21). The computational time for the Fourier transform is proportional to the number of times of multiplication. If the image is N by N square and if N  2p for an integer p, the fast Fourier transform is applicable in which the number of times is 2N 2 log2 N. In the case of N ≠ 2p, the time is proportional to N 3 . 4. Conclusion

We propose a new filtering method to obtain the wrapped phase from a single fringe pattern with spatial carrier frequency modulation. The Fourier spectrum of the fringe pattern consists of three peaks: a dc peak, a carrier peak, and an adjoint peak. The proposed filter can remove both the dc peak and the adjoint peak to obtain the wrapped phase, even if the three peaks are overlapping partially. The filtering procedure consists of two stages: dc peak filtering and adjoint peak filtering. A spectrum shift filter is employed in both stages, and it can remove most of undesired components. The remaining undesired components are removed by the single connected support region filter in the dc peak filtering stage, and by the symmetrical point comparison filter in the adjoint peak filtering stage. The validity of the proposed filtering method is demonstrated through analyses of the two sets of simulation data and a set of experimental data. In the simulations, the unwrapped phase computed

from the filtered spectrum by using a spectrum shift filter has smaller error than that by using a traditional half-plane filter. For the experimental data, the proposed filter can represent the nature of the optical length of heated gas around the gas burner. In the results shown in this article, we take no considerations to choose parameters to adjust the filter. Although results are not shown owing to the lack of pages, there are cases that require adjusting the weight of γ used in the rough estimation of the dc peak. These cases were found in other experimental data, of which the illuminated region is limited in a certain area. Even for these cases, the control parameter is the only one that is less than the other spectral domain filters. This means that the proposed method does not require most of the efforts made to find the right filtering parameters. Appendix A: Simulation of Interferogram

The two-dimensional interferograms to demonstrate the validity of the proposed filter shown in Figs. 4 and 5 are prepared by simulations of interferometer measurements. The phase and the amplitude of waves are represented by isotropic and elliptic Gaussian functions, which are defined as c 2 
y−yc 2 ∕σ 2

GI
r; rc ; σ ≜ e−

x−x  GE
r; rc ; σ; θ ≜ e−
X

;

(A1)

2
r;rc ;θ∕σ 2 Y 2
r;rc ;θ∕σ 2  X Y

;

c

(A2) c

c

where r 
x; y is the pixel position, r 
x ; y  is the center of the Gaussian, θ is the rotated angle, σ and σ 
σ X ; σ Y  are the widths of the Gaussian, and
X; Y represent the rotated position of
x; y; i.e., X 
x − xc  cos θ −
y − yc  sin θ and Y 
x − xc  sin θ 
y − yc  cos θ. 1 September 2014 / Vol. 53, No. 25 / APPLIED OPTICS

5629

The interferogram, i
r, is evaluated as the power of the sum of two waves [the object wave, f obj
r, and the reference wave, f ref
r]. After that, noise is added to simulate dark noise in a digital camera. The amplitudes of the two waves are isotropic Gaussians. The phase shift of the object wave, ϕ
r, is given as a sum of elliptic Gaussian functions, and is also added to the carrier shift, s · r. These are expressed by the following equations: ¯ σ n ; i
r  jf ref
r  f obj
rj2  Nn; c f ref
r  f ref max GI
r; rf ; σ f ; c j
s·rϕ
r ; f obj
r  f obj max GI
r; rf ; σ f e

ϕ
r 

Nϕ X n1

ϕn GE
r; rϕn ; σ ϕn ; θϕn ;

(A3)

¯ σ n  shows the noise with the normal diswhere Nn; tribution having the average n¯ and the standard deviation σ n . The carrier frequency, s, is specified by
N sx ; N sy , which expresses a one-dimensional number of cycles in the image size,
Lx ; Ly , as s∕2π 
N sx ∕Lx ; N sy ∕Lx . Both simulation data are evaluated in a square region with Lx  Ly  256 (r 
x; y fx; yg ∈ 0; ; 255). For the data in Fig. 4 two Gaussians are superposed (N ϕ  2). The other parameters are obj
f ref max ; f max  
1; 0.9;

¯ σ n  
0; 0;
n;
rcf ; σ f  

128; 128; 200;   1  4π
N sx ; N sy  
10.4; 10.4; ϕ 2   π c
rϕ1 ; σ ϕ1 ; θϕ1  
128; 108;
240; 10; ; 6
rcϕ2 ; σ ϕ2 ; θϕ2  

108; 128;
10; 240; 0: For Fig. 5 N ϕ  1, and the parameters are obj
f ref max ; f max  
1; 1;

¯ σ n  
0; 0.3;
n;
rcf ; σ f  

128; 128; 128; ϕ1  10π
rcϕ1 ; σ ϕ1 



128; 128;
100; 100;

and the three carrier frequencies are determined by N sx  N sy ∈ f4.0; 5.5; 7.0g. This research was supported by the Japan Society for the Promotion of Science (JSPS), Grant-in-Aid for Scientific Research (C), 24560216, 2012-2014. 5630

APPLIED OPTICS / Vol. 53, No. 25 / 1 September 2014

References 1. 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, 156–160 (1982). 2. W. W. Macy, Jr., “Two-dimensional fringe-pattern analysis,” Appl. Opt. 22, 3898–3901 (1983). 3. D. J. Bone, H.-A. Bachor, and R. J. Sandeman, “Fringe-pattern analysis using a 2-D Fourier transform,” Appl. Opt. 25, 1653–1660 (1986). 4. M. Takeda, “Spatial-carrier fringe-pattern analysis and its applications to precision interferometry and profilometry: an overview,” Ind. Metrol. 1, 79–99 (1990). 5. M. Takeda, “Fourier fringe analysis and its application to metrology of extreme physical phenomena: a review,” Appl. Opt. 52, 20–29 (2013). 6. Y. Arai, S. Yokozeki, K. Shiraki, and T. Yamada, “High precision two-dimensional spatial fringe analysis method,” J. Mod. Opt. 44, 739–751 (1997). 7. J. A. Ferrari and E. M. Frins, “Multiple phase-shifted interferograms obtained from a single interferogram with linear carrier,” Opt. Commun. 271, 59–64 (2007). 8. K. Creath, “Temporal phase measurement methods,” in Interferogram Analysis: Digital Fringe Pattern Measurement Techniques, D. W. Robinson and G. T. Reid, eds. (IOP, 1993), pp. 94–140. 9. J. C. Wyant, “Computerized interferometric surface measurements,” Appl. Opt. 52, 1–8 (2013). 10. J.-M. Desse, P. Picart, and P. Tankam, “Digital three-color holographic interferometry for flow analysis,” Opt. Express 16, 5471–5480 (2008). 11. K. Creath and G. Goldstein, “Dynamic quantitative phase imaging for biological objects using a pixelated phase mask,” Biomed. Opt. Express 3, 2866–2880 (2012). 12. D. I. Serrano-García, N.-I. Toto-Arellano, A. Martínez-García, and G. R. Zurita, “Radial slope measurement of dynamic transparent samples,” J. Opt. 14, 045706 (2012). 13. J. Heil, T. Bauer, T. Sure, and J. Wesner, “Iterative fullbandwidth wavefront reconstruction from a set of low-tilt Fizeau interferograms for high-numerical-aperture surface characterization,” Appl. Opt. 45, 4270–4283 (2006). 14. M. Kujawinska and J. Wójciak, “High accuracy Fourier transform fringe pattern analysis,” Opt. Lasers Eng. 14, 325–339 (1991). 15. C. Roddier and F. Roddier, “Interferogram analysis using Fourier transform techniques,” Appl. Opt. 26, 1668–1673 (1987). 16. R. W. Gerchberg, “Super-resolution through error energy reduction,” Opt. Acta 21, 709–720 (1974). 17. J. Yañez-Mendiola, M. Servín, and D. Malacara-Hernández, “Iterative method to obtain the wrapped phase in an interferogram with a linear carrier,” Opt. Commun. 178, 291–296 (2000). 18. S. Tomioka and S. Nishiyama, “Nondestructive threedimensional measurement of gas temperature distribution by phase tomography,” Proc. SPIE 8296, 829617 (2011). 19. D. Malacara, ed., “Interferometric optical profilers,” in Optical Shop Testing, 3rd ed. (Wiley, 2007), Chap. 15.4, pp. 695–702. 20. S. Kostianovski, S. G. Lipson, and E. N. Ribak, “Interference microscopy and Fourier fringe analysis applied to measuring the spatial refractive-index distribution,” Appl. Opt. 32, 4744–4750 (1993). 21. J.-F. Lin and X. Su, “Two-dimensional Fourier transform profilometry for the automatic measurement of threedimensional object shapes,” Opt. Eng. 34, 3297–3302 (1995). 22. J. H. Massig and J. Heppner, “Fringe-pattern analysis with high accuracy by use of the Fourier-transform method: theory and experimental tests,” Appl. Opt. 40, 2081–2088 (2001). 23. L.-C. Chen, H.-W. Ho, and X.-L. Nguyen, “Fourier transform profilometry (FTP) using an innovative band-pass filter for accurate 3-D surface reconstruction,” Opt. Lasers Eng. 48, 182–190 (2010). 24. M. A. Herráez, D. Burton, and M. Lalor, “Accelerating fast Fourier transform and filtering operations in Fourier fringe analysis for accurate measurement of three-dimensional surfaces,” Opt. Lasers Eng. 31, 135–145 (1999).

25. Z. Ge, F. Kobayashi, S. Matsuda, and M. Takeda, “Coordinatetransform technique for closed-fringe analysis by the Fouriertransform method,” Appl. Opt. 40, 1649–1657 (2001). 26. T. Kreis, “Digital holographic interference-phase measurement using the Fourier-transform method,” J. Opt. Soc. Am. A 3, 847–855 (1986). 27. Q. Kemao, “Two-dimensional windowed Fourier transform for fringe pattern analysis: principles, applications and implementations,” Opt. Lasers Eng. 45, 304–317 (2007). 28. M. Servin, J. L. Marroquin, and F. J. Cuevas, “Demodulation of a single interferogram by use of a two-dimensional regularized phase-tracking technique,” Appl. Opt. 36, 4540–4548 (1997). 29. M. Servin, J. L. Marroquin, and J. A. Quiroga, “Regularized quadrature and phase tracking from a single closed-fringe interferogram,” J. Opt. Soc. Am. A 21, 411–419 (2004). 30. J. C. Estrada, M. Servín, J. A. Quiroga, and J. L. Marroquín, “Path independent demodulation method for single image interferograms with closed fringes within the function space c2,” Opt. Express 14, 9687–9698 (2006).

31. C. Tian, Y. Yang, T. Wei, T. Ling, and Y. Zhuo, “Demodulation of a single-image interferogram using a Zernike-polynomialbased phase-fitting technique with a differential evolution algorithm,” Opt. Lett. 36, 2318–2320 (2011). 32. G. Rajshekhar and P. Rastogi, “Fringe demodulation using the two-dimensional phase differencing operator,” Opt. Lett. 37, 4278–4280 (2012). 33. G. Rajshekhar and P. Rastogi, “Multiple signal classification technique for phase estimation from a fringe pattern,” Appl. Opt. 51, 5869–5875 (2012). 34. I. Gurov and M. Volynsky, “Interference fringe analysis based on recurrence computational algorithms,” Opt. Lasers Eng. 50, 514–521 (2012). 35. W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, “Optimal (Wiener) filtering with the FFT,” in Numerical Recipes: the Art of Scientific Computing, 3rd ed. (Cambridge University, 2007), Chap. 13.3, pp. 649–652. 36. S. Tomioka and S. Nishiyama, “Phase unwrapping for noisy phase map using localized compensator,” Appl. Opt. 51, 4984–4994 (2012).

1 September 2014 / Vol. 53, No. 25 / APPLIED OPTICS

5631