Signal tracking approach for phase estimation in digital holographic interferometry Rahul G. Waghmare,1 Deepak Mishra,1 G. R. K. Sai Subrahmanyam,2 Earu Banoth,3 and Sai Siva Gorthi3,* 1
Department of Avionics, Indian Institute of Space Science and Technology, Trivandrum, Kerala 695547, India 2
Department of Earth & Space Science, Indian Institute of Space Science and Technology, Trivandrum, Kerala 695547, India
3
Department of Instrumentation & Applied Physics, Indian Institute of Science, Bangalore 560012, India *Corresponding author: [email protected] Received 26 March 2014; revised 21 May 2014; accepted 21 May 2014; posted 22 May 2014 (Doc. ID 204250); published 23 June 2014
In this research work, we introduce a novel approach for phase estimation from noisy reconstructed interference fields in digital holographic interferometry using an unscented Kalman filter. Unlike conventionally used unwrapping algorithms and piecewise polynomial approximation approaches, this paper proposes, for the first time to the best of our knowledge, a signal tracking approach for phase estimation. The state space model derived in this approach is inspired from the Taylor series expansion of the phase function as the process model, and polar to Cartesian conversion as the measurement model. We have characterized our approach by simulations and validated the performance on experimental data (holograms) recorded under various practical conditions. Our study reveals that the proposed approach, when compared with various phase estimation methods available in the literature, outperforms at lower SNR values (i.e., especially in the range 0–20 dB). It is demonstrated with experimental data as well that the proposed approach is a better choice for estimating rapidly varying phase with high dynamic range and noise. © 2014 Optical Society of America OCIS codes: (090.1995) Digital holography; (090.2880) Holographic interferometry; (100.2650) Fringe analysis; (100.5070) Phase retrieval; (120.5050) Phase measurement. http://dx.doi.org/10.1364/AO.53.004150
1. Introduction
Determination of the phase of a digital holographic interferogram plays a vital role in nondestructive testing with applications in displacement measurement, reliability assessment, safety monitoring, quality control, etc. [1]. In digital holographic interferometry (DHI), usually, two holograms are recorded corresponding to the object state before and after the deformation. Numerical reconstruction of these holograms provides their respective reconstructed object wave fields. Multiplication of one 1559-128X/14/194150-08$15.00/0 © 2014 Optical Society of America 4150
APPLIED OPTICS / Vol. 53, No. 19 / 1 July 2014
reconstructed object wave field with that of complex conjugate of the other, generates the reconstructed interference field [2]. The phase of this reconstructed interference field carries the information of the full-field deformation undergone by the object. The methods of phase estimation from a reconstructed interference field can be broadly categorized into two classes. One class of approaches uses the arctan function to get the phase, but the phase generated by this procedure is noisy and wrapped [3]. Hence, noise filtering [4,5] and unwrapping operations [6,7] become an essential part of phase estimation to get continuous phase. On the other hand, the second class of approaches deals with the piecewise polynomial approximation approach, where we get
unwrapped and continuous phase directly [2]. In these approaches, phase maps are modeled as piecewise polynomials of orders 2–4. Thus, the problem of phase estimation in such cases becomes that of parameter estimation and the coefficients of the phase polynomials becomes the parameters to be estimated. It utilizes various signal processing methods like maximum likelihood estimation (MLE) [2], discrete chirp-Fourier transform (DCFT) [8], highorder ambiguity function (HAF) [9], and extended Kalman filter (EKF) [10]. These methods are useful in DHI provided that the phase maps can be approximated with smaller order polynomials. If the order of the polynomial increases, the number of data points required per segment also increases. This makes the piecewise polynomial approximation approach unreliable for rapidly varying phase functions. In addition, as the phase map generated by these methods is discontinuous at the end points of each segment, phase stitching becomes an essential part of the algorithm. In this paper, we propose a novel approach based on signal tracking that provides more accurate, unwrapped, and continuous phase directly even if the phase is rapidly varying, has a larger dynamic range, or is corrupted by severe noise. In this approach, each column of the phase map is considered as an arbitrary function. The state space model is then derived based on the Taylor series expansion of this arbitrary function. We propose an unscented Kalman filter (UKF) as a signal tracking algorithm to process the derived state space model. The paper is organized as follows. The theoretical foundation for the proposed method is discussed in Section 2. Simulation and experimental results are presented in Sections 3 and 4, respectively, followed by the conclusions in Section 5. 2. Theoretical Foundation
ϕn 1 ϕn
1 0 1 ϕ n ϕ00 n w0 n: 1! 2!
(3)
Subsequently, we can calculate ϕ0 n 1 and ϕ00 n 1 as ϕ0 n 1 ϕ0 n
1 00 ϕ n w1 n; 1!
(4)
ϕ00 n 1 ϕ00 n w2 n:
(5)
Here, wi n; i 0; 1; 2 contains higher order terms of the Taylor series which are modeled as process noise. The amplitude is modeled as random walk so that we can write an 1 an wa n:
(6)
Now, the above equations can be written in matrix form as 2
an 1
3
2
1 0 0 0
32
an
3
2
wa n
3
6 ϕn 1 7 6 0 1 1 1 76 ϕn 7 6 w n 7 7 6 7 6 0 7 6 1! 2! 76 76 76 76 7; 6 0 4 ϕ n 1 5 4 0 0 1 1!1 54 ϕ0 n 5 4 w1 n 5 0 0 0 1 ϕ00 n 1 ϕ00 n w2 n (7) i.e., xn 1 Fxn wn:
(8)
Here, the vector xn denotes the state vector, the matrix F is the state transition matrix representing the relation between the present and the next state of the field, and wn denotes the process noise. Using the state vector we can generate the measurement signal with the relation Rfn an cosϕn;
(9)
The reconstructed interference field of DHI with variable amplitude embedded in noise can be expressed by the equation [2]
Ifn an sinϕn:
(10)
f m; n am; neiϕm;n ηm; n;
Hence, from above equations, the function h· used for the prediction of observation from state vector x can be written as
A.
State Space Model
(1)
where am; n is the amplitude, ϕm; n is the interference phase of the analytic signal, and ηm; n is the zero-mean complex additive white Gaussian noise (AWGN). Here, m and n are the rows and columns, respectively, of the N × N complex field. For any arbitrary column, Eq. (1) can be written as f n aneiϕn ηn:
(2)
In this proposed method, the phase is assumed to be a continuous and differentiable function. Hence, we can define the Taylor series expansion of the phase function ϕy in Eq. (2) as
x1 cosx2 hx ; x1 sinx2
(11)
where xk is the kth element in the state vector x. As the measurement signal is a complex field, the noise will also be in the complex form. We assume the measurement (observation) noise to be AWGN with zero mean and variance of σ 2ν . The noise νy is given by
Rηy νy : Iηy 1 July 2014 / Vol. 53, No. 19 / APPLIED OPTICS
(12) 4151
The noise covariance R can be estimated effectively from the measured signal as R kR σ 2v I;
The corresponding mean and the covariance of Z are approximated using transformed sigma points as
(13)
z¯
where kR is a scaling constant used for tuning the filter and I is 2 × 2 identity matrix. Hence, the observation model can be written as zn hxn νn: B.
Given the state model and the observation model, we propose the use of UKF for state estimation. The UKF uses a deterministic sampling technique, unscented transform (UT), to pick a minimal set of sample points (called sigma points) around the mean such that these points capture the mean and covariance of a prior random variable exactly, while approximating the mean and covariance of the transformed random variable up to the third order in the Taylor series [11]. Consider propagation of a random variable x (of dimension L), a state vector in our case, having mean x¯ and covariance Px , through a nonlinear function z hx. To calculate the first two moments of z, we form a matrix of sigma points as follows: χ 0 x¯ ;
(15)
p L λPx ; ∀ i 1; …; L; χ i x¯ i p χ i x¯ − L λPx ; ∀ i L 1; …; 2L; i
i0
(14)
Unscented Transformation
(16)
wμ 0
λ ; L λ
(17)
wc 0
λ 1 − α2 β; L λ
(18)
2L X
Pz
2L X i0
wμ i Zi;
(21)
¯ Z i − z¯ T : wc i Z i − z
(22)
These estimates in the mean are accurate up to the third order in the Taylor series, while the covariance estimates are accurate up to the fourth order in the Taylor series expansion [11,12]. C.
Algorithm
The UKF is a straightforward extension of the UT, presented in Section 2.B, to the recursive estimation of the state vector from Eqs. (8) and (14). The augmented state vector in UKF is redefined as the concatenation of the original state and noise variables. The UT sigma point selection scheme [Eqs. (15)–(19)] is applied to this new augmented state vector to calculate the prior sigma points. These sigma points are then passed through the state and observation functions to produce transformed sigma points. The statistics (mean and covariance) for the Kalman update equations are determined by using the transformed sigma points. Finally, the Kalman updates are employed. Let, xa xT uT vT be an augmented state vector, λ be the composite scaling parameter, L be dimension of augmented state, Pu be the process noise covariance, Pv be the measurement noise covariance, and wi be the weights as calculated in Eqs. (15)–(19). Then, the algorithm for the UKF as presented by Julier and Uhlmann [11] is used for state estimation. Initialize with: xˆ 0 Efx0 g; P0 Efx0 − xˆ 0 x0 − xˆ 0 T g;
c wμ i wi
1 ; ∀ i 1; …; 2L; 2L λ
(19)
and wc where χ i are the sigma points, wμ i i are the weights to calculate the mean and covariance, respectively. λ α2 L k − L, is a scaling factor, α determines the spread of the sigma points around x¯ and is usually set to a small positive value. β is used to incorporate the prior knowledge about distribup tion of x. L λPx i is the ith row of the matrix square root of L λPx. These sigma vectors are then passed through a transformation function to get Z i hχ i ; 4152
∀ i 0; …; 2L:
APPLIED OPTICS / Vol. 53, No. 19 / 1 July 2014
(20)
xˆ a0 Efxa g xˆ T0
0 0 T ;
Pa0 Efxa0 − xˆ a0 xa0 − xˆ a0 T g 3 2 0 P0 0 7 6 7 6 4 0 Pu 0 5; 0 0 Pv where Ef·g is an expectation operator. Now, ∀ k ∈ f1; …; Ng, where N is the number of pixels in a column. Calculate the Sigma Points: h p p i χ ak−1 xˆ ak−1 xˆ ak−1 L λPak−1 xˆ ak−1 − L λPak−1 :
Process Update Equations: χ xkjk−1 f χ xk−1 ; χ uk−1 ; xˆ −k P−k
2L X i0 2L X i0
x wμ i χ i;kjk−1 ; x ˆ −k χ xi;kjk−1 − xˆ −k T ; wc i χ i;kjk−1 − x
Z kjk−1 hχ xkjk−1 ; χ vk−1 ; zˆ −k
2L X i0
wμ i Z i;kjk−1 :
Measurement Update Equations: Pzz Pxz
2L X i0 2L X i0
ˆ −k Z i;kjk−1 − zˆ −k T ; wc i Z i;kjk−1 − z ˆ −k Z i;kjk−1 − zˆ −k T ; wc i χ i;kjk−1 − x
K Pxz P−1 zz ; xˆ k xˆ −k Kzk − zˆ −k ; Pk P−k − KPzz KT : We note that, unlike in EKF, no explicit calculation of the Jacobian or Hessian of the nonlinear function is necessary to implement this algorithm. Also, our algorithm has the capability of directly providing the unwrapped phase from the state vector, as the second element of the state vector is a phase function itself. In addition, the state vector [see Eq. (7)] includes first and second derivatives of the phase that provides direct access to phase derivatives without additional computational efforts. It is indeed a great advantage offered by this approach in nondestructive testing and evaluation applications, where calculation of displacement (phase) as well as strain/ curvature (its derivatives) are important. 3. Simulation Results
A complex interference field signal of size 512 × 512 was simulated in MATLAB. The zero-mean AWGN was then added at different SNRs using the AWGN function. The UKF algorithm was applied to each column with the initial guess taken as the arctan of the first element of each column. With this strategy of initialization, we generate the phase map which is unwrapped for a given column, but needs post processing to avoid wrapping along the rows. Instead of this, we have used the phase value, calculated for first element of the previous column as an initial guess of the phase for the current column so that neither unwrapping nor phase stitching is required at all. The process noise covariance matrix is taken to be Q diag 10−1 10−1 10−5 10−5 whereas
the measurement noise covariance matrix is taken to be R kR diag σ 2re σ 2im , where σ 2re and σ 2im are the variances estimated from the uniform region of the fringe pattern, and kR is a scaling factor generally taken as 2. To demonstrate the effectiveness of our approach, we have compared it with the state-of-the-art methods. The comparison is made at three different levels. At the first level, we modified our signal tracking approach so that we can generate the polynomial coefficients and then reconstruct the phase. This modification is just like the piecewise polynomial approximation approach, where we can fit secondorder polynomials for phase approximation in each segment. We call the estimated phase by this method as UKFPara . The comparison of UKFPara with different methods at slowly varying and lower dynamic range of phase maps at various noise levels are depicted in Section 3.A. Here, we have shown that the UKFPara performs very well as compared to different methods. At second level, in Section 3.B, we have compared our direct signal estimation approach with existing parameter estimation methods at larger dynamic range of the same phase map. Here, we have shown that the parameter estimation approach becomes unreliable for such phase maps. Finally, we have tested our algorithm for a rapidly varying phase map. Since no signal tracking algorithm has been proposed so far for the phase estimation, and to justify the choice of UKF versus EKF for realizing this approach, we have formulated EKF also as a signal tracking algorithm (EKFSig ). Comparison of the results of UKFSig and EKFSig is shown in Section 3.C. A. Parameter Estimation Approach
A phase map [Fig. 1(a)] is generated using the peaks function in MATLAB. The complex interference field is generated using Eq. (1) and the corresponding fringe pattern at a SNR of 20 dB is shown in Fig. 1(b). As shown in Section 2, the proposed approach is intended to estimate the phase signal directly as compared to the parameter estimation. But, for fair comparison with the existing parameter-estimationbased methods and the proposed approach, we have used the following relation [13,14] to get the coefficients of polynomials, which are then used to reconstruct the phase (UKFPara ):
Fig. 1. 3D mesh plot of the original phase and the corresponding fringe pattern at a SNR of 20 dB. 1 July 2014 / Vol. 53, No. 19 / APPLIED OPTICS
4153
θn : CF−n xn;
(23)
where θn An; a2 ; a1 ; a0 T is a parameter vector containing amplitude as its first element, and matrix C is a diagonal with elements 1, 1, 1, 0.5. We compared UKFPara with the following state-ofthe-art parameter estimation methods: the DCFTbased method [8], the IHAF-based method [9], and the EKF-based method (EKFPara ) [10]. While the number of segments is taken as 16 (Nw 16) for all the methods, the order of polynomial approximation for IHAF and EKFPara is taken as 4 (M 4), whereas it is 2 for DCFT, and UKFPara (M 2). Error in the estimated phase of the simulated pattern (Fig. 1) by different methods is shown in Fig. 2. Although all the methods perform equally well at a SNR of 20 dB, as can be seen from Figs. 2(a)–2(d), different methods start failing as the noise is increased. Figure 2(e) shows the root mean square error (RMSE) produced by different methods as SNR is varied from 0 to 30 dB. Since both EKF and UKF are the variants of Kalman filters, which have the tendency to diverge due to noisy situations, inaccurate state space model, and inaccurate initial conditions, we have compared the divergence rates of both the methods in Fig. 2(f). It can be seen from Fig. 2(f) that EKFPara quickly diverges below the cutoff SNR of 19 dB, while UKFPara does not diverge even at 0 dB. Divergence in the
50
1 DCFT IHAF EKF
0.9
UKF
0.7
Para
0.5 0.4 0.3 0.2
DCFT
EKFPara
IHAF
Proposed
Time (s)
204.3605
27.0874
9.1061
38.1811
EKFPara is due to the linearization in the measurement model and an inadequate process model. Whereas our process model is based on the Taylor series expansion, with an additional process noise parameter, thus the UKFPara works effectively even at very low SNR levels. Table 1 gives the comparison of different methods based on computational time. Each method is used to process the reconstructed interference field of size 512 × 512. These methods were evaluated on a Windows PC with Intel(R) Xeon(R) CPU E3-1225 V2 at 3.20 GHz, 16 GB RAM, and with MATLAB version 8.2.0.701 (R2013b). B. Comparison between Parameter and Signal Tracking Approaches
We have shown in Section 3.A that the UKF-based method applied for the parameter estimation approach (UKFPara ) outperforms state-of-the-art parameter estimation methods under noisy conditions. In this subsection, we compared the performance of the UKF-based method used as the proposed signal tracking approach (UKFSig ) with various parameter estimation methods. In order to quantify the performance variations among different methods in dealing with signals having higher dynamic range, we generated a phase map by using the peaks function but upscaled by a factor of 20, at the same resolution of 512 × 512. The 3D mesh plot of the simulated phase map is shown in Fig. 3(a), and the corresponding fringe pattern at a SNR of 20 dB is shown in Fig. 3(b). Comparison for the estimated phase with different methods along the middle column is shown in Fig. 4. As we increase the dynamic range (of values of the phase map), keeping the same resolution, parameterestimation-based methods underperform producing distortion in the estimated phase, even at higher SNR values. The reason is, at increased dynamic range of phase, even a small error in the estimation of parameters causes a larger change in the overall estimated phase signal.
Para
UKFPara
35 30 25 20 15 10
0.1
5 5
10
15 20 SNR (dB)
25
30
0
0
5
10
15 20 SNR (dB)
25
30
Fig. 2. (a)–(d) Error in phase estimation with different methods at a SNR of 20 dB. (e) RMSE versus SNR and (f) divergence rate versus SNR. 4154
Method
40
0.6
0 0
Comparison of Computational Time
EKF
45
Para
%age Divergence
RMSE (rad)
0.8
Table 1.
APPLIED OPTICS / Vol. 53, No. 19 / 1 July 2014
Fig. 3. 3D mesh plot of the original phase (20 × peaks) and the corresponding fringe pattern with a SNR of 20 dB.
True Phase DCFT IHAF EKFPara
Range of Phase (rad)
150
100
UKFSig 50
0
−50
−100 0
100
200
300
400
500
Pixels
Fig. 4. Estimated phase comparison along the middle column of the phase map at a SNR of 20 dB.
of the triangle and the slope of the side changes. This provides insights for the ability of the phase estimation technique to handle sudden change at the vertices of the triangles of different shapes. The fringe pattern shown in Fig. 6(b) is analyzed with the proposed (UKFSig ), and for the purpose of comparing its performance with any other possible signal-tracking-approach-based methods, we also implemented EKFSig and the error plots obtained with both the methods are shown in Figs. 6(c) and 6(d). It is found that both filters work almost identically above a SNR value of 12 dB. But in noisy conditions, EKFSig start to diverge owing to linearization of the measurement model. Figures 6(e) and 6(f) show the plots of RMSE and divergence rate against SNR for EKFSig and UKFSig , respectively. Figure 6 indicates that the UKF proves to be a better choice than EKF, in the successful realization of the proposed signal tracking approach for phase estimation in DHI. It is found that, though the process model follows the Taylor series expansion, which requires the function to be continuous and differentiable, the approach works well with nondifferentiable functions owing to signal tracking nature of the
Fig. 5. Phase estimation error for the phase map generated by 20 × peaks showing relatively poor performance of polynomialapproximation-based methods at a larger dynamic range of the phase map.
Dealing with Rapidly Varying Phase Signals
While deriving the state space model, we used the Taylor series expansion of the phase function which requires an assumption of continuity and differentiability at a given point. However, usually phase maps obtained in DHI may not always be well-behaved functions. To verify the ability of the proposed approach in handling such cases, we generated a phase map as shown in Fig. 6(a). This phase map consists of multiple triangles in each row. As we traverse through the columns, the angle between two sides
50
EKF Sig UKFSig
2.5
EKF Sig UKF
45
Sig
40 %age Divergence
C.
3
RMSE (rad)
Figure 5 shows the 3D plots of error in estimated phase by DCFT, IHAF, EKFPara , and UKFSig . It can be seen from Fig. 5 that if the phase has a high dynamic range, the parametric methods fail to estimate the overall phase map. It is observed that the computational time for both UKFPara and UKFSig are same.
2 1.5 1
35 30 25 20 15 10
0.5
5 0 0
5
10
15 20 SNR (dB)
25
30
0 0
5
10
15 20 SNR (dB)
25
30
Fig. 6. Comparison of EKFSig and UKFSig for signal-trackingbased phase estimation for rapidly varying phase maps at 20 dB. (a) and (b) The 3D mesh plot of the true phase map and corresponding 2D image of the fringe pattern. (c) and (d) The error in phase estimation using EKFSig and UKFSig . (e) and (f) Comparison of the performance of EKFSig and UKFSig for different values of SNR in terms of RMSE and divergence rate. 1 July 2014 / Vol. 53, No. 19 / APPLIED OPTICS
4155
algorithm and the use of recent observations in the Kalman update. 4. Experimental Results
The experimental setup used for recording the digital holographic interferograms is shown in Fig. 7. A laser (HNL150L-He–Ne laser, 632.8 nm, 15 mW, Thorlab) has been used as the light source. The beam is split into two equal halves (BE1, object beam; BE2, reference beam) with the help of a nonpolarizing beam splitter cube (BS010, Thorlab). Digital holograms are recorded with a monochrome CCD camera (DMK72BUC02 imaging source). The test specimen was made out of aluminum (0.5 mm thickness) and held on the optical table by clamping the plate on all four sides with a custom-designed mount. A threaded bore (M4 screw thread) was created at the center on the backside of the mount, such that driving an M4 screw through this bore would induce a deformation in the object (plate). Broadband dielectric mirrors (M1, M2, M3) (BB1-E02, Thorlab) mounted onto adjustable precision kinematic mirror mounts (KS1, Thorlab) were used to adjust the direction of the beams. The two beams reflected off the mirrors (M2, M3) were expanded using two concave lenses (f–80 mm). The beam of light scattered by the object interferes with the reference beam in a 3D volume of space. The interference patterns at the plane of the CCD sensor are recorded before and after inducing the deformation. These two digital holograms are then numerically reconstructed to obtain the complex amplitudes [A1 x; y; A2 x; y] by applying a discrete Fresnel transform. Figure 8 shows the intensity (jA1x; yj2 ) image of one of the reconstructed holograms. A portion of the numerical reconstruction (outlined by the bright rectangle in the figure) was used to perform digital holographic interference analysis. The complex amplitude of the predeformation of the object [A1 x; y] was multiplied with the complex conjugate of the postdeformation of the object [A2 x; y] to obtain the reconstructed interference field. The real part of the reconstructed interference field gives the fringe pattern, and the phase of it carries the information of the object deformation. The fringe patterns obtained under two such deformation experiments are shown in
Fig. 7. Experimental setup. 4156
APPLIED OPTICS / Vol. 53, No. 19 / 1 July 2014
Fig. 8. Intensity image of the reconstructed hologram.
Figs. 9(a) and 9(b), where Fig. 9(a) is the representative case for rapidly varying phase signals and Fig. 9(b) is that of the noisy data. EKFPara and UKFSig were applied to the reconstructed interference field to estimate the phase. It is observed that EKFPara , by being a piecewise polynomial approximation approach, performs poorly. Please note that although both methods directly provide the unwrapped phase distributions, fringe
Fig. 9. Comparison of EKFPara and UKFSig for experimental data. (a) and (b) The noisy fringe pattern of the reconstructed interference field. (c) and (d) The fringes corresponding to phase estimated by EKFPara whereas (e) and (f) are that of proposed UKFSig .
of process noise covariance into the state space model and use of UKF to estimate the state provides a significant improvement in the accuracy of phase estimation, especially in a noisy situation (SNR 0–20 dB). Also, since the derived state vector includes the phase and derivatives of the phase, it may give wide applicability of this algorithm in situations where simultaneous estimation of displacement (phase) as well as strain/curvature (derivatives of phase) are important. This work is funded by the IIST/ISRO Fellowship by the Department of Space, Govt. of India. We thank Dr. Rajaseker Gannavarpu for invaluable suggestions in the realization of [10]. We would also like to acknowledge the support extended by the Applied Photonics Initiative (API) of IISc in realizing the experimental setup. References
Fig. 10. 3D mesh plot of the estimated phase using the proposed method.
patterns are generated for the purpose of qualitative comparison. Figures 9(c) and 9(d) show fringe patterns generated using the estimated phase maps by EKFPara and Figs. 9(e) and 9(f) show that by UKFSig . The estimated phase by EKFPara shows distortion which is due to the divergence of EKF in a noisy situation. Figure 10 shows the 3D mesh plots of the estimated phase signals by the proposed UKFSig method. The proposed approach is able to estimate the phase even when it is nondifferentiable and rapidly varying [Fig. 10(a)] and in real noisy situations [Fig. 10(b)]. 5. Conclusions
This paper proposes a new approach, namely, the signal tracking approach, for phase estimation in DHI. Simulation and experimental results demonstrate that the proposed method yields more accurate phase estimates than the state-of-the-art approaches, especially at lower SNRs. It is observed that the method performs very well even for complicated phase signals where the polynomialapproximation-approach-based methods may have difficulty in the estimation of the parameters. The analysis results indicate that the introduction
1. U. Schnars and W. P. O. Jüptner, “Digital recording and numerical reconstruction of holograms,” Meas. Sci. Technol. 13, R85–R101 (2002). 2. S. S. Gorthi and P. Rastogi, “Piece-wise polynomial phase approximation approach for the analysis of reconstructed interference fields in digital holographic interferometry,” J. Opt. A Pure Appl. Opt. 11, 065405 (2009). 3. S. S. Gorthi and P. Rastogi, “Fringe projection techniques: whither we are?” Opt. Lasers Eng. 48, 133–140 (2010). 4. F. Palacios, E. Gonalves, J. Ricardo, and J. L. Valin, “Adaptive filter to improve the performance of phase-unwrapping in digital holography,” Opt. Commun. 238, 245–251 (2004). 5. M. J. Huang and W. Sheu, “Histogram-data-orientated filter for inconsistency removal of interferometric phase maps,” Opt. Eng. 44, 045602 (2005). 6. Y. Li, J. Zhu, and W. Shen, “Phase unwrapping algorithms, respectively, based on path-following and discrete cosine transform,” Optik 119, 545–547 (2008). 7. E. Zappa and G. Busca, “Comparison of eight unwrapping algorithms applied to Fourier-transform profilometry,” Opt. Lasers Eng. 46, 106–116 (2008). 8. S. S. Gorthi and P. Rastogi, “Estimation of phase derivatives using discrete chirp-Fourier-transform-based method,” Opt. Lett. 34, 2396–2398 (2009). 9. S. S. Gorthi and P. Rastogi, “Improved high-order ambiguityfunction method for the estimation of phase from interferometric fringes,” Opt. Lett. 34, 2575–2577 (2009). 10. G. Rajshekhar and P. Rastogi, “Phase estimation using a state-space approach based method,” Opt. Laser Eng. 51, 1004–1007 (2013). 11. S. J. Julier and J. K. Uhlmann, “Unscented filtering and nonlinear estimation,” Proc. IEEE 92, 401–422 (2004). 12. G. R. K. S. Subrahmanyam, “Recursive image estimation and inpainting in noise using non-Gaussian MRF prior,” Ph.D. thesis (Indian Institute of Technology, Madras, 2008). 13. J. Gal, A. Campeanu, and I. Nafornita, “Estimation of chirp signals in Gaussian noise by Kalman filtering,” in International Symposium on Signals, Circuits and Systems (IEEE, 2007), Vol. 1, pp. 1–4. 14. J. Gal, A. Campeanu, and I. Nafornita, “Identification of polynomial phase signals by extended Kalman filtering,” in Proceedings of the European Signal Processing Conference (EURASIP, 2008).
1 July 2014 / Vol. 53, No. 19 / APPLIED OPTICS
4157
Department of Avionics, Indian Institute of Space Science and Technology, Trivandrum, Kerala 695547, India 2
Department of Earth & Space Science, Indian Institute of Space Science and Technology, Trivandrum, Kerala 695547, India
3
Department of Instrumentation & Applied Physics, Indian Institute of Science, Bangalore 560012, India *Corresponding author: [email protected] Received 26 March 2014; revised 21 May 2014; accepted 21 May 2014; posted 22 May 2014 (Doc. ID 204250); published 23 June 2014
In this research work, we introduce a novel approach for phase estimation from noisy reconstructed interference fields in digital holographic interferometry using an unscented Kalman filter. Unlike conventionally used unwrapping algorithms and piecewise polynomial approximation approaches, this paper proposes, for the first time to the best of our knowledge, a signal tracking approach for phase estimation. The state space model derived in this approach is inspired from the Taylor series expansion of the phase function as the process model, and polar to Cartesian conversion as the measurement model. We have characterized our approach by simulations and validated the performance on experimental data (holograms) recorded under various practical conditions. Our study reveals that the proposed approach, when compared with various phase estimation methods available in the literature, outperforms at lower SNR values (i.e., especially in the range 0–20 dB). It is demonstrated with experimental data as well that the proposed approach is a better choice for estimating rapidly varying phase with high dynamic range and noise. © 2014 Optical Society of America OCIS codes: (090.1995) Digital holography; (090.2880) Holographic interferometry; (100.2650) Fringe analysis; (100.5070) Phase retrieval; (120.5050) Phase measurement. http://dx.doi.org/10.1364/AO.53.004150
1. Introduction
Determination of the phase of a digital holographic interferogram plays a vital role in nondestructive testing with applications in displacement measurement, reliability assessment, safety monitoring, quality control, etc. [1]. In digital holographic interferometry (DHI), usually, two holograms are recorded corresponding to the object state before and after the deformation. Numerical reconstruction of these holograms provides their respective reconstructed object wave fields. Multiplication of one 1559-128X/14/194150-08$15.00/0 © 2014 Optical Society of America 4150
APPLIED OPTICS / Vol. 53, No. 19 / 1 July 2014
reconstructed object wave field with that of complex conjugate of the other, generates the reconstructed interference field [2]. The phase of this reconstructed interference field carries the information of the full-field deformation undergone by the object. The methods of phase estimation from a reconstructed interference field can be broadly categorized into two classes. One class of approaches uses the arctan function to get the phase, but the phase generated by this procedure is noisy and wrapped [3]. Hence, noise filtering [4,5] and unwrapping operations [6,7] become an essential part of phase estimation to get continuous phase. On the other hand, the second class of approaches deals with the piecewise polynomial approximation approach, where we get
unwrapped and continuous phase directly [2]. In these approaches, phase maps are modeled as piecewise polynomials of orders 2–4. Thus, the problem of phase estimation in such cases becomes that of parameter estimation and the coefficients of the phase polynomials becomes the parameters to be estimated. It utilizes various signal processing methods like maximum likelihood estimation (MLE) [2], discrete chirp-Fourier transform (DCFT) [8], highorder ambiguity function (HAF) [9], and extended Kalman filter (EKF) [10]. These methods are useful in DHI provided that the phase maps can be approximated with smaller order polynomials. If the order of the polynomial increases, the number of data points required per segment also increases. This makes the piecewise polynomial approximation approach unreliable for rapidly varying phase functions. In addition, as the phase map generated by these methods is discontinuous at the end points of each segment, phase stitching becomes an essential part of the algorithm. In this paper, we propose a novel approach based on signal tracking that provides more accurate, unwrapped, and continuous phase directly even if the phase is rapidly varying, has a larger dynamic range, or is corrupted by severe noise. In this approach, each column of the phase map is considered as an arbitrary function. The state space model is then derived based on the Taylor series expansion of this arbitrary function. We propose an unscented Kalman filter (UKF) as a signal tracking algorithm to process the derived state space model. The paper is organized as follows. The theoretical foundation for the proposed method is discussed in Section 2. Simulation and experimental results are presented in Sections 3 and 4, respectively, followed by the conclusions in Section 5. 2. Theoretical Foundation
ϕn 1 ϕn
1 0 1 ϕ n ϕ00 n w0 n: 1! 2!
(3)
Subsequently, we can calculate ϕ0 n 1 and ϕ00 n 1 as ϕ0 n 1 ϕ0 n
1 00 ϕ n w1 n; 1!
(4)
ϕ00 n 1 ϕ00 n w2 n:
(5)
Here, wi n; i 0; 1; 2 contains higher order terms of the Taylor series which are modeled as process noise. The amplitude is modeled as random walk so that we can write an 1 an wa n:
(6)
Now, the above equations can be written in matrix form as 2
an 1
3
2
1 0 0 0
32
an
3
2
wa n
3
6 ϕn 1 7 6 0 1 1 1 76 ϕn 7 6 w n 7 7 6 7 6 0 7 6 1! 2! 76 76 76 76 7; 6 0 4 ϕ n 1 5 4 0 0 1 1!1 54 ϕ0 n 5 4 w1 n 5 0 0 0 1 ϕ00 n 1 ϕ00 n w2 n (7) i.e., xn 1 Fxn wn:
(8)
Here, the vector xn denotes the state vector, the matrix F is the state transition matrix representing the relation between the present and the next state of the field, and wn denotes the process noise. Using the state vector we can generate the measurement signal with the relation Rfn an cosϕn;
(9)
The reconstructed interference field of DHI with variable amplitude embedded in noise can be expressed by the equation [2]
Ifn an sinϕn:
(10)
f m; n am; neiϕm;n ηm; n;
Hence, from above equations, the function h· used for the prediction of observation from state vector x can be written as
A.
State Space Model
(1)
where am; n is the amplitude, ϕm; n is the interference phase of the analytic signal, and ηm; n is the zero-mean complex additive white Gaussian noise (AWGN). Here, m and n are the rows and columns, respectively, of the N × N complex field. For any arbitrary column, Eq. (1) can be written as f n aneiϕn ηn:
(2)
In this proposed method, the phase is assumed to be a continuous and differentiable function. Hence, we can define the Taylor series expansion of the phase function ϕy in Eq. (2) as
x1 cosx2 hx ; x1 sinx2
(11)
where xk is the kth element in the state vector x. As the measurement signal is a complex field, the noise will also be in the complex form. We assume the measurement (observation) noise to be AWGN with zero mean and variance of σ 2ν . The noise νy is given by
Rηy νy : Iηy 1 July 2014 / Vol. 53, No. 19 / APPLIED OPTICS
(12) 4151
The noise covariance R can be estimated effectively from the measured signal as R kR σ 2v I;
The corresponding mean and the covariance of Z are approximated using transformed sigma points as
(13)
z¯
where kR is a scaling constant used for tuning the filter and I is 2 × 2 identity matrix. Hence, the observation model can be written as zn hxn νn: B.
Given the state model and the observation model, we propose the use of UKF for state estimation. The UKF uses a deterministic sampling technique, unscented transform (UT), to pick a minimal set of sample points (called sigma points) around the mean such that these points capture the mean and covariance of a prior random variable exactly, while approximating the mean and covariance of the transformed random variable up to the third order in the Taylor series [11]. Consider propagation of a random variable x (of dimension L), a state vector in our case, having mean x¯ and covariance Px , through a nonlinear function z hx. To calculate the first two moments of z, we form a matrix of sigma points as follows: χ 0 x¯ ;
(15)
p L λPx ; ∀ i 1; …; L; χ i x¯ i p χ i x¯ − L λPx ; ∀ i L 1; …; 2L; i
i0
(14)
Unscented Transformation
(16)
wμ 0
λ ; L λ
(17)
wc 0
λ 1 − α2 β; L λ
(18)
2L X
Pz
2L X i0
wμ i Zi;
(21)
¯ Z i − z¯ T : wc i Z i − z
(22)
These estimates in the mean are accurate up to the third order in the Taylor series, while the covariance estimates are accurate up to the fourth order in the Taylor series expansion [11,12]. C.
Algorithm
The UKF is a straightforward extension of the UT, presented in Section 2.B, to the recursive estimation of the state vector from Eqs. (8) and (14). The augmented state vector in UKF is redefined as the concatenation of the original state and noise variables. The UT sigma point selection scheme [Eqs. (15)–(19)] is applied to this new augmented state vector to calculate the prior sigma points. These sigma points are then passed through the state and observation functions to produce transformed sigma points. The statistics (mean and covariance) for the Kalman update equations are determined by using the transformed sigma points. Finally, the Kalman updates are employed. Let, xa xT uT vT be an augmented state vector, λ be the composite scaling parameter, L be dimension of augmented state, Pu be the process noise covariance, Pv be the measurement noise covariance, and wi be the weights as calculated in Eqs. (15)–(19). Then, the algorithm for the UKF as presented by Julier and Uhlmann [11] is used for state estimation. Initialize with: xˆ 0 Efx0 g; P0 Efx0 − xˆ 0 x0 − xˆ 0 T g;
c wμ i wi
1 ; ∀ i 1; …; 2L; 2L λ
(19)
and wc where χ i are the sigma points, wμ i i are the weights to calculate the mean and covariance, respectively. λ α2 L k − L, is a scaling factor, α determines the spread of the sigma points around x¯ and is usually set to a small positive value. β is used to incorporate the prior knowledge about distribup tion of x. L λPx i is the ith row of the matrix square root of L λPx. These sigma vectors are then passed through a transformation function to get Z i hχ i ; 4152
∀ i 0; …; 2L:
APPLIED OPTICS / Vol. 53, No. 19 / 1 July 2014
(20)
xˆ a0 Efxa g xˆ T0
0 0 T ;
Pa0 Efxa0 − xˆ a0 xa0 − xˆ a0 T g 3 2 0 P0 0 7 6 7 6 4 0 Pu 0 5; 0 0 Pv where Ef·g is an expectation operator. Now, ∀ k ∈ f1; …; Ng, where N is the number of pixels in a column. Calculate the Sigma Points: h p p i χ ak−1 xˆ ak−1 xˆ ak−1 L λPak−1 xˆ ak−1 − L λPak−1 :
Process Update Equations: χ xkjk−1 f χ xk−1 ; χ uk−1 ; xˆ −k P−k
2L X i0 2L X i0
x wμ i χ i;kjk−1 ; x ˆ −k χ xi;kjk−1 − xˆ −k T ; wc i χ i;kjk−1 − x
Z kjk−1 hχ xkjk−1 ; χ vk−1 ; zˆ −k
2L X i0
wμ i Z i;kjk−1 :
Measurement Update Equations: Pzz Pxz
2L X i0 2L X i0
ˆ −k Z i;kjk−1 − zˆ −k T ; wc i Z i;kjk−1 − z ˆ −k Z i;kjk−1 − zˆ −k T ; wc i χ i;kjk−1 − x
K Pxz P−1 zz ; xˆ k xˆ −k Kzk − zˆ −k ; Pk P−k − KPzz KT : We note that, unlike in EKF, no explicit calculation of the Jacobian or Hessian of the nonlinear function is necessary to implement this algorithm. Also, our algorithm has the capability of directly providing the unwrapped phase from the state vector, as the second element of the state vector is a phase function itself. In addition, the state vector [see Eq. (7)] includes first and second derivatives of the phase that provides direct access to phase derivatives without additional computational efforts. It is indeed a great advantage offered by this approach in nondestructive testing and evaluation applications, where calculation of displacement (phase) as well as strain/ curvature (its derivatives) are important. 3. Simulation Results
A complex interference field signal of size 512 × 512 was simulated in MATLAB. The zero-mean AWGN was then added at different SNRs using the AWGN function. The UKF algorithm was applied to each column with the initial guess taken as the arctan of the first element of each column. With this strategy of initialization, we generate the phase map which is unwrapped for a given column, but needs post processing to avoid wrapping along the rows. Instead of this, we have used the phase value, calculated for first element of the previous column as an initial guess of the phase for the current column so that neither unwrapping nor phase stitching is required at all. The process noise covariance matrix is taken to be Q diag 10−1 10−1 10−5 10−5 whereas
the measurement noise covariance matrix is taken to be R kR diag σ 2re σ 2im , where σ 2re and σ 2im are the variances estimated from the uniform region of the fringe pattern, and kR is a scaling factor generally taken as 2. To demonstrate the effectiveness of our approach, we have compared it with the state-of-the-art methods. The comparison is made at three different levels. At the first level, we modified our signal tracking approach so that we can generate the polynomial coefficients and then reconstruct the phase. This modification is just like the piecewise polynomial approximation approach, where we can fit secondorder polynomials for phase approximation in each segment. We call the estimated phase by this method as UKFPara . The comparison of UKFPara with different methods at slowly varying and lower dynamic range of phase maps at various noise levels are depicted in Section 3.A. Here, we have shown that the UKFPara performs very well as compared to different methods. At second level, in Section 3.B, we have compared our direct signal estimation approach with existing parameter estimation methods at larger dynamic range of the same phase map. Here, we have shown that the parameter estimation approach becomes unreliable for such phase maps. Finally, we have tested our algorithm for a rapidly varying phase map. Since no signal tracking algorithm has been proposed so far for the phase estimation, and to justify the choice of UKF versus EKF for realizing this approach, we have formulated EKF also as a signal tracking algorithm (EKFSig ). Comparison of the results of UKFSig and EKFSig is shown in Section 3.C. A. Parameter Estimation Approach
A phase map [Fig. 1(a)] is generated using the peaks function in MATLAB. The complex interference field is generated using Eq. (1) and the corresponding fringe pattern at a SNR of 20 dB is shown in Fig. 1(b). As shown in Section 2, the proposed approach is intended to estimate the phase signal directly as compared to the parameter estimation. But, for fair comparison with the existing parameter-estimationbased methods and the proposed approach, we have used the following relation [13,14] to get the coefficients of polynomials, which are then used to reconstruct the phase (UKFPara ):
Fig. 1. 3D mesh plot of the original phase and the corresponding fringe pattern at a SNR of 20 dB. 1 July 2014 / Vol. 53, No. 19 / APPLIED OPTICS
4153
θn : CF−n xn;
(23)
where θn An; a2 ; a1 ; a0 T is a parameter vector containing amplitude as its first element, and matrix C is a diagonal with elements 1, 1, 1, 0.5. We compared UKFPara with the following state-ofthe-art parameter estimation methods: the DCFTbased method [8], the IHAF-based method [9], and the EKF-based method (EKFPara ) [10]. While the number of segments is taken as 16 (Nw 16) for all the methods, the order of polynomial approximation for IHAF and EKFPara is taken as 4 (M 4), whereas it is 2 for DCFT, and UKFPara (M 2). Error in the estimated phase of the simulated pattern (Fig. 1) by different methods is shown in Fig. 2. Although all the methods perform equally well at a SNR of 20 dB, as can be seen from Figs. 2(a)–2(d), different methods start failing as the noise is increased. Figure 2(e) shows the root mean square error (RMSE) produced by different methods as SNR is varied from 0 to 30 dB. Since both EKF and UKF are the variants of Kalman filters, which have the tendency to diverge due to noisy situations, inaccurate state space model, and inaccurate initial conditions, we have compared the divergence rates of both the methods in Fig. 2(f). It can be seen from Fig. 2(f) that EKFPara quickly diverges below the cutoff SNR of 19 dB, while UKFPara does not diverge even at 0 dB. Divergence in the
50
1 DCFT IHAF EKF
0.9
UKF
0.7
Para
0.5 0.4 0.3 0.2
DCFT
EKFPara
IHAF
Proposed
Time (s)
204.3605
27.0874
9.1061
38.1811
EKFPara is due to the linearization in the measurement model and an inadequate process model. Whereas our process model is based on the Taylor series expansion, with an additional process noise parameter, thus the UKFPara works effectively even at very low SNR levels. Table 1 gives the comparison of different methods based on computational time. Each method is used to process the reconstructed interference field of size 512 × 512. These methods were evaluated on a Windows PC with Intel(R) Xeon(R) CPU E3-1225 V2 at 3.20 GHz, 16 GB RAM, and with MATLAB version 8.2.0.701 (R2013b). B. Comparison between Parameter and Signal Tracking Approaches
We have shown in Section 3.A that the UKF-based method applied for the parameter estimation approach (UKFPara ) outperforms state-of-the-art parameter estimation methods under noisy conditions. In this subsection, we compared the performance of the UKF-based method used as the proposed signal tracking approach (UKFSig ) with various parameter estimation methods. In order to quantify the performance variations among different methods in dealing with signals having higher dynamic range, we generated a phase map by using the peaks function but upscaled by a factor of 20, at the same resolution of 512 × 512. The 3D mesh plot of the simulated phase map is shown in Fig. 3(a), and the corresponding fringe pattern at a SNR of 20 dB is shown in Fig. 3(b). Comparison for the estimated phase with different methods along the middle column is shown in Fig. 4. As we increase the dynamic range (of values of the phase map), keeping the same resolution, parameterestimation-based methods underperform producing distortion in the estimated phase, even at higher SNR values. The reason is, at increased dynamic range of phase, even a small error in the estimation of parameters causes a larger change in the overall estimated phase signal.
Para
UKFPara
35 30 25 20 15 10
0.1
5 5
10
15 20 SNR (dB)
25
30
0
0
5
10
15 20 SNR (dB)
25
30
Fig. 2. (a)–(d) Error in phase estimation with different methods at a SNR of 20 dB. (e) RMSE versus SNR and (f) divergence rate versus SNR. 4154
Method
40
0.6
0 0
Comparison of Computational Time
EKF
45
Para
%age Divergence
RMSE (rad)
0.8
Table 1.
APPLIED OPTICS / Vol. 53, No. 19 / 1 July 2014
Fig. 3. 3D mesh plot of the original phase (20 × peaks) and the corresponding fringe pattern with a SNR of 20 dB.
True Phase DCFT IHAF EKFPara
Range of Phase (rad)
150
100
UKFSig 50
0
−50
−100 0
100
200
300
400
500
Pixels
Fig. 4. Estimated phase comparison along the middle column of the phase map at a SNR of 20 dB.
of the triangle and the slope of the side changes. This provides insights for the ability of the phase estimation technique to handle sudden change at the vertices of the triangles of different shapes. The fringe pattern shown in Fig. 6(b) is analyzed with the proposed (UKFSig ), and for the purpose of comparing its performance with any other possible signal-tracking-approach-based methods, we also implemented EKFSig and the error plots obtained with both the methods are shown in Figs. 6(c) and 6(d). It is found that both filters work almost identically above a SNR value of 12 dB. But in noisy conditions, EKFSig start to diverge owing to linearization of the measurement model. Figures 6(e) and 6(f) show the plots of RMSE and divergence rate against SNR for EKFSig and UKFSig , respectively. Figure 6 indicates that the UKF proves to be a better choice than EKF, in the successful realization of the proposed signal tracking approach for phase estimation in DHI. It is found that, though the process model follows the Taylor series expansion, which requires the function to be continuous and differentiable, the approach works well with nondifferentiable functions owing to signal tracking nature of the
Fig. 5. Phase estimation error for the phase map generated by 20 × peaks showing relatively poor performance of polynomialapproximation-based methods at a larger dynamic range of the phase map.
Dealing with Rapidly Varying Phase Signals
While deriving the state space model, we used the Taylor series expansion of the phase function which requires an assumption of continuity and differentiability at a given point. However, usually phase maps obtained in DHI may not always be well-behaved functions. To verify the ability of the proposed approach in handling such cases, we generated a phase map as shown in Fig. 6(a). This phase map consists of multiple triangles in each row. As we traverse through the columns, the angle between two sides
50
EKF Sig UKFSig
2.5
EKF Sig UKF
45
Sig
40 %age Divergence
C.
3
RMSE (rad)
Figure 5 shows the 3D plots of error in estimated phase by DCFT, IHAF, EKFPara , and UKFSig . It can be seen from Fig. 5 that if the phase has a high dynamic range, the parametric methods fail to estimate the overall phase map. It is observed that the computational time for both UKFPara and UKFSig are same.
2 1.5 1
35 30 25 20 15 10
0.5
5 0 0
5
10
15 20 SNR (dB)
25
30
0 0
5
10
15 20 SNR (dB)
25
30
Fig. 6. Comparison of EKFSig and UKFSig for signal-trackingbased phase estimation for rapidly varying phase maps at 20 dB. (a) and (b) The 3D mesh plot of the true phase map and corresponding 2D image of the fringe pattern. (c) and (d) The error in phase estimation using EKFSig and UKFSig . (e) and (f) Comparison of the performance of EKFSig and UKFSig for different values of SNR in terms of RMSE and divergence rate. 1 July 2014 / Vol. 53, No. 19 / APPLIED OPTICS
4155
algorithm and the use of recent observations in the Kalman update. 4. Experimental Results
The experimental setup used for recording the digital holographic interferograms is shown in Fig. 7. A laser (HNL150L-He–Ne laser, 632.8 nm, 15 mW, Thorlab) has been used as the light source. The beam is split into two equal halves (BE1, object beam; BE2, reference beam) with the help of a nonpolarizing beam splitter cube (BS010, Thorlab). Digital holograms are recorded with a monochrome CCD camera (DMK72BUC02 imaging source). The test specimen was made out of aluminum (0.5 mm thickness) and held on the optical table by clamping the plate on all four sides with a custom-designed mount. A threaded bore (M4 screw thread) was created at the center on the backside of the mount, such that driving an M4 screw through this bore would induce a deformation in the object (plate). Broadband dielectric mirrors (M1, M2, M3) (BB1-E02, Thorlab) mounted onto adjustable precision kinematic mirror mounts (KS1, Thorlab) were used to adjust the direction of the beams. The two beams reflected off the mirrors (M2, M3) were expanded using two concave lenses (f–80 mm). The beam of light scattered by the object interferes with the reference beam in a 3D volume of space. The interference patterns at the plane of the CCD sensor are recorded before and after inducing the deformation. These two digital holograms are then numerically reconstructed to obtain the complex amplitudes [A1 x; y; A2 x; y] by applying a discrete Fresnel transform. Figure 8 shows the intensity (jA1x; yj2 ) image of one of the reconstructed holograms. A portion of the numerical reconstruction (outlined by the bright rectangle in the figure) was used to perform digital holographic interference analysis. The complex amplitude of the predeformation of the object [A1 x; y] was multiplied with the complex conjugate of the postdeformation of the object [A2 x; y] to obtain the reconstructed interference field. The real part of the reconstructed interference field gives the fringe pattern, and the phase of it carries the information of the object deformation. The fringe patterns obtained under two such deformation experiments are shown in
Fig. 7. Experimental setup. 4156
APPLIED OPTICS / Vol. 53, No. 19 / 1 July 2014
Fig. 8. Intensity image of the reconstructed hologram.
Figs. 9(a) and 9(b), where Fig. 9(a) is the representative case for rapidly varying phase signals and Fig. 9(b) is that of the noisy data. EKFPara and UKFSig were applied to the reconstructed interference field to estimate the phase. It is observed that EKFPara , by being a piecewise polynomial approximation approach, performs poorly. Please note that although both methods directly provide the unwrapped phase distributions, fringe
Fig. 9. Comparison of EKFPara and UKFSig for experimental data. (a) and (b) The noisy fringe pattern of the reconstructed interference field. (c) and (d) The fringes corresponding to phase estimated by EKFPara whereas (e) and (f) are that of proposed UKFSig .
of process noise covariance into the state space model and use of UKF to estimate the state provides a significant improvement in the accuracy of phase estimation, especially in a noisy situation (SNR 0–20 dB). Also, since the derived state vector includes the phase and derivatives of the phase, it may give wide applicability of this algorithm in situations where simultaneous estimation of displacement (phase) as well as strain/curvature (derivatives of phase) are important. This work is funded by the IIST/ISRO Fellowship by the Department of Space, Govt. of India. We thank Dr. Rajaseker Gannavarpu for invaluable suggestions in the realization of [10]. We would also like to acknowledge the support extended by the Applied Photonics Initiative (API) of IISc in realizing the experimental setup. References
Fig. 10. 3D mesh plot of the estimated phase using the proposed method.
patterns are generated for the purpose of qualitative comparison. Figures 9(c) and 9(d) show fringe patterns generated using the estimated phase maps by EKFPara and Figs. 9(e) and 9(f) show that by UKFSig . The estimated phase by EKFPara shows distortion which is due to the divergence of EKF in a noisy situation. Figure 10 shows the 3D mesh plots of the estimated phase signals by the proposed UKFSig method. The proposed approach is able to estimate the phase even when it is nondifferentiable and rapidly varying [Fig. 10(a)] and in real noisy situations [Fig. 10(b)]. 5. Conclusions
This paper proposes a new approach, namely, the signal tracking approach, for phase estimation in DHI. Simulation and experimental results demonstrate that the proposed method yields more accurate phase estimates than the state-of-the-art approaches, especially at lower SNRs. It is observed that the method performs very well even for complicated phase signals where the polynomialapproximation-approach-based methods may have difficulty in the estimation of the parameters. The analysis results indicate that the introduction
1. U. Schnars and W. P. O. Jüptner, “Digital recording and numerical reconstruction of holograms,” Meas. Sci. Technol. 13, R85–R101 (2002). 2. S. S. Gorthi and P. Rastogi, “Piece-wise polynomial phase approximation approach for the analysis of reconstructed interference fields in digital holographic interferometry,” J. Opt. A Pure Appl. Opt. 11, 065405 (2009). 3. S. S. Gorthi and P. Rastogi, “Fringe projection techniques: whither we are?” Opt. Lasers Eng. 48, 133–140 (2010). 4. F. Palacios, E. Gonalves, J. Ricardo, and J. L. Valin, “Adaptive filter to improve the performance of phase-unwrapping in digital holography,” Opt. Commun. 238, 245–251 (2004). 5. M. J. Huang and W. Sheu, “Histogram-data-orientated filter for inconsistency removal of interferometric phase maps,” Opt. Eng. 44, 045602 (2005). 6. Y. Li, J. Zhu, and W. Shen, “Phase unwrapping algorithms, respectively, based on path-following and discrete cosine transform,” Optik 119, 545–547 (2008). 7. E. Zappa and G. Busca, “Comparison of eight unwrapping algorithms applied to Fourier-transform profilometry,” Opt. Lasers Eng. 46, 106–116 (2008). 8. S. S. Gorthi and P. Rastogi, “Estimation of phase derivatives using discrete chirp-Fourier-transform-based method,” Opt. Lett. 34, 2396–2398 (2009). 9. S. S. Gorthi and P. Rastogi, “Improved high-order ambiguityfunction method for the estimation of phase from interferometric fringes,” Opt. Lett. 34, 2575–2577 (2009). 10. G. Rajshekhar and P. Rastogi, “Phase estimation using a state-space approach based method,” Opt. Laser Eng. 51, 1004–1007 (2013). 11. S. J. Julier and J. K. Uhlmann, “Unscented filtering and nonlinear estimation,” Proc. IEEE 92, 401–422 (2004). 12. G. R. K. S. Subrahmanyam, “Recursive image estimation and inpainting in noise using non-Gaussian MRF prior,” Ph.D. thesis (Indian Institute of Technology, Madras, 2008). 13. J. Gal, A. Campeanu, and I. Nafornita, “Estimation of chirp signals in Gaussian noise by Kalman filtering,” in International Symposium on Signals, Circuits and Systems (IEEE, 2007), Vol. 1, pp. 1–4. 14. J. Gal, A. Campeanu, and I. Nafornita, “Identification of polynomial phase signals by extended Kalman filtering,” in Proceedings of the European Signal Processing Conference (EURASIP, 2008).
1 July 2014 / Vol. 53, No. 19 / APPLIED OPTICS
4157
