3738

OPTICS LETTERS / Vol. 36, No. 19 / October 1, 2011

Application of complex-lag distributions for estimation of arbitrary order phase derivatives in digital holographic interferometry Gannavarpu Rajshekhar and Pramod Rastogi* Applied Computing and Mechanics Laboratory, Ecole Polytechnique Fédérale de Lausanne, 1015-Lausanne, Switzerland *Corresponding author: [email protected] Received July 26, 2011; accepted August 26, 2011; posted August 30, 2011 (Doc. ID 151811); published September 19, 2011 This Letter proposes a method to estimate phase derivatives of arbitrary order in digital holographic interferometry. Based on the desired order, the generalized complex-lag distribution is computed from the reconstructed interference field. Subsequently, the phase derivative is estimated by tracing the peak of the distribution. Simulation and experimental results are presented to validate the method’s potential. © 2011 Optical Society of America OCIS codes: 120.2880, 090.1995, 090.2880.

For many nondestructive testing applications using digital holographic interferometry (DHI), estimation of firstand higher-order derivatives of interference phase is desired since these are related to the physical quantities, such as strain and curvature. Recently, ridge-detection methods based on space-frequency distributions, such as windowed Fourier transform (FT) [1], wavelet transform [2], and Wigner–Ville distributions [3–5], have been developed to measure the first-order phase derivative. In these methods, the ridge or peak location of the distribution provides information about the phase derivative. The main advantage of using space-frequency distributions is direct estimation of the phase derivative without requiring any phase unwrapping or numerical differentiation operations. However, these methods cannot be applied to directly measure higher-order derivatives. In this Letter, a generalized complex-lag-distribution (GCD)-based method is presented that permits the direct estimation of phase derivatives of any order. The complex reconstructed interference field in DHI can be expressed as Γðx; yÞ ¼ aðx; yÞ exp½ jϕðx; yÞ
þ ηðx; yÞ;

ð1Þ

where aðx; yÞ is the amplitude, assumed to be slowly varying, ϕðx; yÞ is the interference phase, and ηðx; yÞ is the noise. For a given row y, the above equation could be written as ΓðxÞ ¼ aðxÞ exp½ jϕðxÞ
þ ηðxÞ:

ϕ

∂m ϕðxÞ ðxÞ ¼ : ∂xm

ð3Þ

To estimate ϕðmÞ ðxÞ, the generalized complex-lag moment (GCM) is used, which is given as [6] GCMm N ½Γ
ðx; τÞ

¼ hðτÞ

N−1 Y k¼0

ωN−m N;k

Γ




GCMm N ½Γ
ðx; τÞ "

# rffiffiffiffiffiffiffiffiffi! m m! N−m ϕ x þ ωN;k ¼ hðτÞ exp j τ ωN;k : N k¼0 N −1 X

ð5Þ

Using Taylor series expansion, we have rffiffiffiffiffiffiffiffiffi! X rffiffiffiffiffiffiffiffiffi!p p ∞ ω N;k m m! m m! ϕðpÞ ðxÞ ϕ x þ ωN;k τ ¼ τ ; ð6Þ p! N N p¼0 and Eq. (5) becomes GCMm N ½Γ
ðx; τÞ " ¼ hðτÞ exp j

ωpþN−m ϕðpÞ ðxÞ N;k p! p¼0

N −1 X ∞ X k¼0

rffiffiffiffiffiffiffiffiffi!p # m m! : τ N

ð7Þ

ð2Þ

The mth-order phase derivative ϕðmÞ ðxÞ with respect to x is given as ðmÞ

where ωN;k ¼ exp½j2πk=N
, N is a parameter such that N ≥ m, and h is a real symmetric window having a lowpass behavior. For simplicity of analysis, the amplitude term within the window region is assumed constant and neglected along with the noise term. So we have

rffiffiffiffiffiffiffiffiffi m m! τ ; ð4Þ x þ ωN;k N

0146-9592/11/193738-03$15.00/0

pþN−m ¼ exp½j2πkðp þ N − mÞ=N
, we have Noting ωN;k

N −1 X k¼0

pþN−m ωN;k

8 > > < N; if p ¼ 0 and N ¼ m ¼ N; if p ¼ Nr þ m ; > > : 0; otherwise

where r is an integer such that r ¼ 0; 1; 2; …. Consequently, the GCM in Eq. (7) vanishes unless p ¼ Nr þ m or p ¼ 0 and N ¼ m. © 2011 Optical Society of America

October 1, 2011 / Vol. 36, No. 19 / OPTICS LETTERS

GCD since the scaling term is independent of τ and j exp½jNϕðxÞ
j ¼ 1. Hence, Eq. (13) remains valid. So, the proposed method can be summarized in the following steps.

Considering the case when N ≠ m, we have GCMm N ½Γ
ðx; τÞ "

rffiffiffiffiffiffiffiffiffi!Nrþm # ∞ X ϕðNrþmÞ ðxÞ m m! ¼ hðτÞ exp jN τ ðNr þ mÞ! N r¼0

ð8Þ

¼ hðτÞ exp½ jϕðmÞ ðxÞτ
exp½ jUðx; τÞ
; where ϕðmÞ ðxÞ is the desired phase derivative and rffiffiffiffiffiffiffiffiffi!Nrþm ∞ X ϕðNrþmÞ ðxÞ m m! Uðx; τÞ ¼ N : τ ðNr þ mÞ! N r¼1

ð9Þ

Within the given window region, the phase usually varies slowly and its (Nr þ m)th-order derivatives, i.e., ϕðNrþmÞ ðxÞ, are negligible for r ≥ 1. Hence, Uðx; τÞ ≈ 0 and Eq. (8) can be written as ðmÞ GCMm ðxÞτ
: N ½Γ
ðx; τÞ ¼ hðτÞ exp½jϕ

3739

ð10Þ

The FT of the GCM provides the GCD: Z ∞ ½Γ
ðx; ΩÞ ¼ hðτÞ exp½jϕðmÞ ðxÞτ
exp½−jΩτ
∂τ; GCDm N −∞

ð11Þ

1. The complex reconstructed interference field of an arbitrary row or column is considered as in Eq. (2). 2. The GCM for the desired phase derivative order m is computed using Eq. (4). 3. The GCD is obtained through the FT of the GCM as in Eq. (11). 4. The mth-order phase derivative is calculated by tracing the peak of the GCD using Eq. (13). 5. The above steps are repeated for all rows or columns to obtain the overall phase derivative. To test its effectiveness, the proposed method was analyzed for the estimation of the first- and second-order phase derivatives. A complex reconstructed interference field signal (size 256 × 256 pixels) of the following form was simulated in MATLAB: Γðx; yÞ ¼ exp½jϕðx; yÞ
;

ð17Þ

and noise was added at a signal-to-noise ratio of 25 dB using MATLAB’s “awgn” function. The original phase ϕðx; yÞ in radians is shown in Fig. 1(a). The

which is essentially the FT of hðτÞ modulated by the term exp½jϕðmÞ ðxÞτ
. For the frequency domain Ω, the FT of the ^ window hðτÞ is denoted as HðΩÞ and the frequency shifting property of the FT yields ðmÞ ^ ðxÞÞ: GCDm N ½Γ
ðx; ΩÞ ¼ HðΩ − ϕ

ð12Þ

Since the window function is chosen to have a low-pass ^ behavior, GCD is maximum for Hð0Þ, i.e., when Ω ¼ ðmÞ ϕ ðxÞ. Hence, the mth-order phase derivative ϕðmÞ ðxÞ can be obtained by finding the frequency Ω at which the GCD spectrum attains its maximum. So we have ϕðmÞ ðxÞ ¼ arg maxjGCDm N ½Γ
ðx; ΩÞj: Ω

ð13Þ

These equations were obtained for the case N ≠ m. When N ¼ m, an extra phase term corresponding to p ¼ 0 in Eq. (7) would arise and consequently, Eqs. (8), (10), and (12) are modified as ðmÞ GCMm ðxÞτ
N ½Γ
ðx; τÞ ¼ hðτÞ exp½jϕ

× exp½jNϕðxÞ
exp½jUðx; τÞ
;

ð14Þ

ðmÞ ðxÞτ
exp½jNϕðxÞ
; ð15Þ GCMm N ½Γ
ðx; τÞ ¼ hðτÞ exp½jϕ ðmÞ ^ GCDm ðxÞÞ: N ½Γ
ðx; ΩÞ ¼ exp½jNϕðxÞ
HðΩ − ϕ

ð16Þ

The above equations show that the additional scaling term exp½jNϕðxÞ
appears when N ¼ m. However, it has no effect on the peak-detection strategy of the

Fig. 1. (Color online) (a) Original phase ϕðx; yÞ in radians. (b) Fringe pattern. (c) Estimated first-order phase derivative ϕð1Þ ðx; yÞ in radians/pixel. (d) First-order phase derivative estimation error. (e) Estimated second-order phase derivative ϕð2Þ ðx; yÞ in radians=pixel2 . (f) Second-order phase derivative estimation error.

3740

OPTICS LETTERS / Vol. 36, No. 19 / October 1, 2011

Fig. 2. (Color online) (a) Experimental fringe pattern. (b) Estimated first-order phase derivative in radians/pixel. (c) Estimated second-order phase derivative in radians=pixel2 .

corresponding fringe pattern constituted by the real part of the complex reconstructed interference field is shown in Fig. 1(b). For the estimation of first-order phase derivative ϕð1Þ ðx; yÞ, we used m ¼ 1 and N ¼ 2 to compute the GCM in Eq. (4). So, the corresponding GCM is given as GCM12 ½Γ
ðx; τÞ ¼ hðτÞΓðx þ τ=2ÞΓ−1 ðx − τ=2Þ:

ð18Þ

This expression leads to a pseudo Wigner–Ville-like distribution [4]. Hence, a pseudo-Wigner–Ville distribution could be considered a special case of the GCD. Similarly, for the estimation of the second-order phase derivative ϕð2Þ ðx; yÞ, we used m ¼ 2 and N ¼ 2 and the corresponding GCM is given as GCM22 ½Γ
ðx; τÞ ¼ hðτÞΓðx þ

pffiffiffi pffiffiffi τÞΓðx − τÞ:

ð19Þ

For the analysis, a Gaussian window of a length of 65 samples was generated using MATLAB’s “gausswin” function. The estimated first-order phase derivative and the corresponding absolute estimation errors in radians/ pixel are shown in Figs. 1(c) and 1(d). Similarly, the estimated second-order phase derivative and the corresponding absolute estimation errors in radians=pixel2 are shown in Figs. 1(e) and 1(f). The root mean square errors for the estimation of the first- and second-order phase derivatives using the proposed method were 0.008 radians/pixel and 2:443 × 10−4 radians=pixel2 . The pixels near the borders were neglected for the analysis. To test its practical applicability in DHI, the proposed method was applied on an experimental fringe pattern, as shown in Fig. 2(a). The estimated first-order phase derivative in radians/pixel after median filtering is shown in Fig. 2(b). Finally, the estimated second-order phase derivative in radians=pixel2 after median filtering is shown in Fig. 2(c). To conclude, this Letter proposes a GCD-based method for the estimation of phase derivatives of arbitrary order in DHI. The method provides a generalized approach to designing distributions based on the desired phase derivative order. This enables direct retrieval of the required phase derivative. In addition, the proposed method does not require phase unwrapping or numerical differentiation operations. Simulation and experimental results validate the strength of the proposed method. References 1. Q. Kemao, Opt. Lasers Eng. 45, 304 (2007). 2. L. R. Watkins, S. M. Tan, and T. H. Barnes, Opt. Lett. 24, 905 (1999). 3. C. A. Sciammarella and T. Kim, Opt. Eng. 42, 3182 (2003). 4. G. Rajshekhar, S. S. Gorthi, and P. Rastogi, Rev. Sci. Instrum. 80, 093107 (2009). 5. G. Rajshekhar, S. S. Gorthi, and P. Rastogi, Opt. Express 18, 18041 (2010). 6. C. Cornu, S. Stanković, C. Ioana, A. Quinquis, and L. Stanković, IEEE Trans. Signal Process. 55, 4831 (2007).