## Abstract

The conventional phase-shifting techniques commonly suffer from frequency aliasing because of their number of phase shifts below the critical sampling rate. As a result, fringe harmonics induce ripple-like artifacts in their reconstructed phase maps. For solving this issue, this paper presents an anti-aliasing phase-measuring technique. Theoretical analysis shows that, with phase-shifting, the harmonics aliased with the fundamental frequency component of a fringe signal depend on the greatest common divisor (GCD) of the used phase shifts. This fact implies a possibility of removing such aliasing effects by selecting non-uniform phase shifts that together with 2*π* have no common divisors. However, even if we do so, it remains challenging to separate harmonics from the fundamental fringe signals, because the systems of equations available from the captured fringe patterns are generally under-determined, especially when the number of phase shifts is very few. To overcome this difficulty, we practically presume that all the points over the fringe patterns have an identical characteristic of harmonics. Under this constraint, using an alternate iterative least-squares fitting procedure allows us to estimate the fringe phases and the harmonic coefficients accurately. Simulation and experimental results demonstrate that this proposed method enables separating high order harmonics from as few as 4 fringe patterns having non-uniform phase shifts, thus significantly suppressing the ripple-like phase errors caused by the frequency aliasing.

© 2022 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement

## 1. Introduction

Phase-shifting [1–6] is a high-resolution phase measuring technique widely used in laser interferometry [7–12], holography [13], moiré topography [14–16], fringe projection profilometry [17–20], and fringe deflectometry [21–24]. This technique assumes a fringe signal to be sinusoidal during phase-shifting, and therefore enables one to calculate the phase, amplitude, and background intensity of this signal from at least three samples having distinct phase shifts. In measurement practice, however, fringe signals are not always perfectly sinusoidal. For example, with interferometry, the nonlinear responses of the detectors induce harmonics in the recorded interferograms [4,5,7]; with fringe projection profilometry, the gamma nonlinearities of light sources severely distort the fringe profiles [25–33]; and in moiré topography, the moiré fringes are not sinusoidal but usually triangularly-shaped [34]. As a result, fringe harmonics induce ripple-like artifacts in the reconstructed phase maps. For this reason, the issue regarding fringe harmonics has been continuously motivating research on phase-shifting algorithms in the past decades.

To suppress the influence of harmonics, many phase-shifting algorithms that require capturing four or more fringe patterns have been developed. Among them, the synchronous detection algorithm invented by Bruning et al. [2] is recognized as the one having the best generality, and in other words, the popularly used three-step [3], four-step [4], and five-step [5] algorithms can be thought of as special cases of this method or their combination. The synchronous detection algorithm uses equal-spaced phase shifts over a full period of $2\pi $ radians and has a harmonic sensitivity depending on its number of phase shifts. Specifically, if the number of phase shifts is *N*, the captured fringes are aliased with the $({lN \pm 1} )\textrm{th}$ orders of harmonics with *l* being a positive integer [13,35]. This phenomenon means that, though the synchronous detection algorithm is known to be the most efficient in restraining fringe harmonics by use of uniform phase shifts, the aliasing effect makes it have an upper limitation of the $N - 2\textrm{th}$ order in restraining such harmonics.

With phase-shifting technique, using random or non-uniform phase shifts provides a possibility of avoiding aliasing in the captured fringe patterns. The reason is that the non-uniform phase shifts, corresponding to a varying sampling rate over time, make Fourier spectrum of the sampled signal lose its symmetry and periodicity. To retrieve the phases from fringe patterns having non-uniform phase shifts, Greivenkamp [6] proposed a least-squares phase-shifting algorithm. This algorithm is very contributive to enhancing the computational efficiency by taking advantage of a transformation for data linearization, but it has little effect in suppressing the influence of fringe harmonics. In fact, when the phase shifts are non-uniform, the available system of equations, when including the coefficients of harmonics as unknowns, is generally under-determined. This fact poses a challenge in separating harmonics from the captured fringe patterns, even if the aliasing effect has been removed by using non-uniform phase shifting.

To overcome the aforementioned issues, we propose an anti-aliasing phase measuring technique. First, we capture four or more fringe patterns having non-uniform phase shifts, with these phase shifts together with 2$\pi $ having no common divisors. As a result, the aliasing effect is avoided and the information of all the harmonics is preserved in the captured fringe patterns. Second, we practically presume that the points over the fringe patterns have an identical characteristic of harmonics, so that we achieve a determined system of equations. For solving this nonlinear system, we suggest an alternate iterative least-squares fitting algorithm, by which the fringe phases and the coefficients of fringe harmonics are estimated accurately. Simulation and experimental results demonstrate that this proposed method enables us to separate high order harmonics from as few as four fringe patterns, and further, to suppress the ripple-like phase errors caused by the frequency aliasing significantly.

## 2. Background

#### 2.1 Aliasing effect in the phase-shifting technique

Phase measuring techniques suffer from the issue of frequency aliasing because they generally involve sampling operation in spatial or temporal domain for getting discrete fringe signals. According to Shannon’s sampling theorem, to completely reconstruct a signal from a set of uniform samples, the sampling rate must be over twice the highest frequency in this signal, which is known as the Nyquist limit [36]. Somewhat differently, a sampling rate below the Nyquist limit is usually acceptable in phase measuring techniques. Its reason is that only the fundamental component of fringe signal needs to be measured and aliasing among high-order harmonics does not affect the measurement results. Even so, this weaker condition is not always satisfied in measurement practice because the fringe signals, as mentioned in the previous section, may have infinite harmonics caused by complicated factors.

Phase-shifting is a temporally sampling technique. Assuming that the number of phase shifts is *N*, the $n\textrm{th}\; ({n = 1,2, \cdots N} )$ fringe pattern recorded by the camera is represented as

*k*th order harmonics, and $\phi ({x,y} )$ and ${\delta _n}$ $(0 \le {\delta _n} < 2\pi )$ represent the phase to be measured and a given phase shift, respectively.

To gain insight into aliasing in Eq. (1), we consider the harmonic of order *k*, i.e., ${a_k}({x,y}) \textrm{cos}[{k({\phi ({x,y} )+ {\delta_n}} )} ]$. If all of the phase shifts ${\delta _n}$ ($n = 1,2, \cdots ,N$) simultaneously satisfy

For examining the above analysis results, we take Bruning’s synchronous detection algorithm [2] as example to do a simple simulation. This algorithm uses equal-spaced phase shifts over a period, i.e., ${\delta _n} = 2({n - 1} )\pi /N$ for $n = 1,2, \cdots ,N$, and calculates phase values through

According to Eq. (6), the results of Eq. (7) are certainly affected by the harmonics of orders $lN \pm 1$ with $l = 1,2, \cdots , + \infty $. In Fig. 1, the top row illustrates its simulation result when $N = 4$. Figure 1(a) shows a sinusoidal fringe profile distorted by harmonics, in which the dot line contains only the second order harmonic and the solid curve is distorted by a nonlinear transformation of gamma 1.55 thus having infinite harmonics. Corresponding to the two intensity curves, Fig. 1(b) features their harmonics. Figure 1(c) gives the phases calculated by using Eq. (7). Figure 1(d) plots their phase errors, from which we observe that the influence of the second order harmonic is thoroughly removed by this 4-step algorithm, but the third, fifth and some higher orders of harmonics induce ripple-like artifacts in the calculated phases because aliasing exists.

The purpose of doing this analysis is to deduce a strategy of selecting phase shifts for avoiding aliasing effect. From Eq. (6), we know that decreasing the GCD in the denominator will increase the lower bound of order of the aliased harmonics. Especially, if we select non-uniform phase shifts and make these phase shifts together with $2\pi $ have no common divisors, aliasing effect caused by phase-shifting will be canceled.

#### 2.2 Under-determined system of equations

In Section 2.1, we deduced that using non-uniform phase-shifts provides a possibility of removing aliasing effect. The issue that arises here is how the phases are calculated from the fringe patterns having these non-uniform phase shifts.

Related to this issue, Greivenkamp’s algorithm [6] neglects the second and higher orders of harmonics from Eq. (1), and uses the transformation, ${b_1}({x,y} )= {a_1}({x,y} )\textrm{cos}[{\phi ({x,y} )} ]$ and ${b_2}({x,y} )={-} {a_1}({x,y} )\textrm{sin}[{\phi ({x,y} )} ]$, to linearize the equations. As a result, for each pixel $({x,y} )$, we have a system of linear equations as

Solving it in the least-squares sense for the three unknowns, i.e., ${a_0}({x,y} )$, ${b_1}({x,y} )$, and ${b_2}({x,y} )$, the phase at $({x,y} )$ is calculated by

Using this algorithm, we continue the simple simulation in Fig. 1 and show the result in its second row whose layout is parallel with the top row. The number of phase shifts keeps to be $N = 4$, but the phase shifts become ${\delta _n} = ({0,\; 1,\; 3,\; 5} )$, not uniformly distributed over a $2\pi $ period. Because these phase shifts together with $2\pi $ do not have common divisors, no aliasing occurs. The fringe profiles in Fig. 1(e) contain the same harmonics as in Fig. 1(a). Figure 1(f) shows these harmonics. Figure 1(g) plots the calculated phases. Figure 1(h) shows that the phase errors are not eliminated even if only the second order harmonic exists. These results demonstrate that Greivenkamp’s algorithm [6] is very contributive to enhancing the computational efficiency by taking advantage of linearization, but it has little effect in suppressing the influence of fringe harmonics when using non-uniform phase shifts.

A possible solution to the above issue is to include some coefficients of harmonics into the equation system, viz.,

*N*equations for each pixel $({x,y} )$, from which at most

*N*unknowns including ${a_0}\left( {x,y} \right)$, ${a_1}({x,y} )$, $\cdots $, ${a_{N - 2}}({x,y} )$, and $\phi ({x,y} )$ may be solved out. The system in Eq. (10) is nonlinear. Newton method [37] is a standard method for solving it and the appendix following this paper states the procedure. It is shown that, within a $2\pi $ phase period, there exist up to $\left( {N - 2} \right)\left( {N - 1} \right) - 2$ phase positions, at which the Jacobian matrices related to this system are singular. Near these phase positions, this system is ill-conditioned and does not have stable solutions. This fact means that the system in Eq. (10) is usually under-determined. Using the same parameters in Figs. 1(e) to 1(h), simulation of solving Eq. (10) is also implemented, and its result is given in Fig. 1 by the third row. The panels in this row show, from left to right, the fringe intensity profiles, the harmonics, the calculated phases, and the phase errors. From this result, we find that the method of solving Eq. (10) fails in recovering phase curves because of the singular points.

The analysis in this subsection reveals that, even if the aliasing effect has been removed by using non-uniform phase shifting, it remains challenging to separate harmonics from the captured fringe patterns and to solve the fringe phases accurately. We shall overcome this problem in the next section.

## 3. Anti-aliasing phase measuring technique

#### 3.1 Additional constraint to harmonics

In Section 2, we have confirmed that using non-uniform phase shifts, which together with 2$\pi $ have no common divisors, enables the captured fringe patterns to be free from aliasing. In fact, it is just a necessary condition for separating high-order harmonics from few fringe patterns having phase shifts. The difficulty is that the few fringe patterns cannot provide enough many equations to specify a unique solution. To make the system of equations determined, some additional constraints about the fringe harmonics have to be defined.

In different applications of optical measurement, the fringe harmonics may have different characteristics. These special characteristics can be exploited and used as the additional constraint. In this work, we practically presume that all the points over the fringe patterns have identical characteristic of harmonics, so that Eq. (1) is restated as

*N*fringe patterns having non-uniform phase shifts. The problem remains how we solve this system of equations which is not only nonlinear but also large-scaled.

#### 3.2 Alternate iterative least-squares algorithm

Following Section 3.1, we assume that a fringe pattern has *M* pixels, so we have $MN$ equations in total with *N* being the number of phase shifts. Given the rapid decline of the harmonics as their orders increase, we truncate the signals and keep only the first $K + 1$ frequency components, including the background intensities, the fundament fringe components, and the harmonics of orders from $2$ through *K*.

To solve such a large-scaled system of nonlinear equations, we suggest an alternate iterative least-squares algorithm. The procedure is as follows:

**Step 1.** Assuming the coefficients of harmonics are known, use Greivenkamp’s algorithm [6] to estimate ${a_0}({x,y} )$, ${a_1}({x,y} )$, and $\phi ({x,y} )$ for each pixel in the least-squares sense. In other words, solve the following normal equations

In the following iterations, we use

**Step 2.**As ${a_0}({x,y} )$, ${a_1}({x,y} )$, and $\phi ({x,y} )$ have been calculated by Step 1, Eq. (11) becomes a linear equation of ${r_k}\; (k = 2,3 \cdots K$). Here, ${r_1}$ keeps to be $1$. Estimate the coefficients of harmonics, ${r_k}\; (k = 2,3 \cdots K$), by fitting the pixels over the patterns in the least-squares sense. The normal equation system related to this over-determined linear system is where ${{\mathbf C}^{(i )}}$ is a $({K - 1} )\times ({K - 1} )$ symmetric matrix with its entry at the $p\textrm{th}$ ($p = 1,2, \cdots K - 1$) row and the $q\textrm{th}$ ($q = 1,2, \cdots K - 1$) column being

On the right-hand side of Eq. (19), ${{\mathbf H}^{(i )}}$ is a column vector with its $p\textrm{th}$ entry being

The method of solving Eq. (19) enables us to estimate the coefficients of harmonics, ${r_k}$, during iterations.

Solving Eq. (19) has a relatively high computational complexity. Note that, because of the orthogonality of sinusoidal functions, the matrix ${{\mathbf C}^{(i )}}$ is nearly a diagonal matrix as long as a sufficiently large number of points having different phases involved in its computation, and this condition is generally satisfied in practical measurement. Therefore, by using a diagonal matrix instead of the accurate ${{\mathbf C}^{(i )}}$, the coefficients, $r_{k}$ for $k = 2,3, \cdots K$, in the *i*th iteration can be directly estimated by

By calculating Eq. (23) instead of solving Eq. (19), the computational burden is significantly reduced at the expense of slightly decreased accuracy.

**Step 3.** Repeat implementing Steps 1 and 2 alternately until the algorithm converges.

In Fig. 1, the simple simulation result of implementing this anti-aliasing phase measuring technique is given by the bottom row. It has the same layout and the same predefined parameters as the two rows above it. The phase errors plotted in Fig. 1(p) demonstrated that this newly proposed technique enables effectively suppressing the influence of the aliasing on the phase measurement results.

## 4. Numerical simulations

In Sections 2 and 3, we have performed some simple numerical simulations as the aid of principle derivation. To verify the feasibility of the proposed method more completely, further simulations should be carried out.

Assuming that the number of phase shifts is four, three algorithms are used to measure the simulated fringe patterns. Their results are compared in Fig. 2 whose rows, from top to bottom, are with Bruning’s synchronous detection algorithm [2], Greivenkamp’s least-squares algorithm [6], and our newly proposed method, respectively. The phase-shifting patterns have harmonics because they are distorted by a nonlinearity of gamma 1.55. The first ones of the phase-shifting pattern sequences with the three methods are illustrated in Fig. 2 by the leftmost column. These three patterns all have the same size of $300 \times $300 pixels and the same appearances because they all have zero phase shifts. The background intensities and the modulations are not uniform across the patterns but Gaussian shaped having a 50% decrease at the image corners. Gaussian noise having a SD of 0.005 is added onto the patterns (The grey levels of the fringe patterns are normalized to have a range from 0 to1). When using Bruning’s algorithm, four phase-shifting fringe patterns have a fixed relative phase step of $\pi /2$ radians. Because the phase shifts are equal-spaced, there exists aliasing between the fundamental fringe component and the ($4l \pm 1)$th order harmonics with $l = 1,2, \cdots , + \infty $ according to Eq. (6). As a result, its calculated phases are shown in Fig. 2(b) whose errors are calculated by subtracting out the predefined phases and shown in Fig. 2(c). From this error map, we observe the ripple-like artifacts caused by the aliasing, which have a frequency four times higher than the fringe frequency. The second row gives the result of using Greivenkamp’s algorithm. The four phase shifts are defined as 0, 1, 3, 5 radians, which togather with $2\pi $ do not have common divisors, in order to avoid frequency aliasing in the captured patterns. In this case, the reconstructed phase map is shown in Fig. 2(e). As we see from the error map in Fig. 2(f), Greivenkamp’s algorithm does not suppress the effects of any harmonics. Using the same fringe patterns, we implement the proposed method. Here, we truncate the fringe signals and keep up to $7\textrm{th}$ order of harmonics. The reconstructed phase map is shown in Fig. 2(h), accompanyed by its error map in Fig. 2(i), from which it is evident that the phase errors induced by the harmonics are thoroughly removed.

The results of a second simulation are illustrated in Fig. 3 which has the same layout with Fig. 2. The difference is that the number of phase shifts becomes five. For the first row of Fig. 3, the relative phase step of the synchronous detection algorithm becomes $2\pi /5$ radians, leading to the aliasing of fundamental fringe component with the ($5l \pm 1)$th order harmonics for $l = 1,2, \cdots , + \infty $. In comparison with Fig. 2(c), the phase errors in Fig. 3(c) contain the artifacts having a slightly decreased amplitude and a frequency five times higher than the fringes. In the second and third rows, which are related with Greivenkamp’s algorithm and the proposed method, the phase shifts are defined to be 0, 1, 3, 4.5, 6 radians thus no aliasing existing in the fringe patterns. Comparing Figs. 3(f) and 3(i), we know that, the fringe harmonics induce large errors in the calculated phase map by Greivenkamp’s algorithm, but they do not affect the results of our proposed method.

It is well known that, when using equal-spaced phase shifts, Greivenkamp’s algorithm gives results the same as those with Bruning’s, because in this case its coefficient matrix of normal system is a diagonal one. It is interesting that, using equal-spaced phase shifts, our proposed method still works, but it delivers results almost the same as those with Bruning’s. This phenomenon means that selecting non-uniform is a necessary condition for avoiding possible aliasing.

The effectiveness of the proposed method is possibly affected by random noise. We implement the above simulations repeatedly under different noise conditions. The root-mean-square (RMS) values of their phase errors are calculated and listed in Table 1. From this table, we find that Bruning’s algorithm enables restraining phase errors to a certain extent because it is immune to up to the $({N - 2} )\textrm{th}\; $order harmonic with *N* being the number of phase shifts, but it is affected by aliasing of some higher order harmonics. Under the same condition, Greivenkamp’s least-squares algorithm, when using non-uniform phase shifts, has little effect in suppressing the influence of fringe harmonics. The proposed technique which uses non-uniform phase shifts is much more effective than others in canceling the errors caused by aliasing, especially when the fringes have low level noise. In the high noise situations, the phase errors of all these three methods increase, but the proposed method is still valid whose residual errors are mainly caused by noise rather than by the fringe harmonics.

Another issue is regarding convergence of this proposed technique, which is directly related to the computational efficiency of the method. Figure 4 plots the curves of the RMS phase errors decline with the number of iterations, under different noise conditions and using different numbers of phase shifts. Under the noise free condition, the curves rapidly converge to 0 after about 20 iterations. In the presence of noise, the curves converge to certain values depending on the noise level. These values are residual errors caused by noise rather than by fringe harmonics. Comparing Figs. 4(a) and 4(b), we find that using more phase shifts makes the algorithm converge faster, and simultaneously, to smaller residual errors. This phenomenon occurs because the phases estimated from more samples generally have smaller variances according to the theory of probability and statistics.

## 5. Experiments and discussions

#### 5.1 Phase measuring in fringe projection profilometry

Besides numerical simulations, it is necessary to confirm the performance of this anti-aliasing phase measuring technique experimentally. As we mentioned in the previous sections, our proposed method is suitable for use in different optical measurement techniques like interferometry, fringe projection, and moiré topography.

Firstly, we use this proposed method in fringe projection profilometry, with which the captured fringe patterns typically contain harmonics because the nonlinearities of the projector and the camera exist. The measurement system consists of a camera (AVT G192B, 1600 × 1200) and a projector (LiyingLDLP-500). The projector casts sinusoidal fringes onto the measured object, and the camera captures the deformed fringe patterns. We calculate the fringe phases and then convert them into the object depths through a mapping relationship obtained by the system calibration [38,39]. Because both the projector and the camera have nonlinearities, harmonics are introduced into the captured fringe patterns. Here, we use three different algorithms, including Bruning’s synchronous detection algorithm [2], Greivenkamp’s least squares algorithm [6], and our proposed method, to reterive the fringe phases.

The results of measuring the first object are given in Fig. 5, with its first and second rows corresponding to the cases that the numbers of phase shifts are four and five, respectively. In each row, the first panel is shared by the three methods as the first frames of their phase-shifting fringe patterns. When using Bruning’s algorithm, the phase shifts are equal-spaced over the $2\pi $ period, and in other words it has a relative phase step of $\pi /2\; $ radians in the first row and $2\pi /5$ radians in the second row, respectively. According to Section 2.1, frequency aliasing occurs between the fundamental fringe component and the ($Nl \pm 1)$th order harmonics with *N* being the number of phase shifts and $l = 1,2, \cdots , + \infty $. As a result, their reconstructed depth maps, shown in Figs. 5(b) and 5(f), contain visible ripple-like artifacts having frequencies four- and five- times, respectively, higher than the fringes. To avoid aliasing effect, we use non-uniform phase shifts 0, 1, 3, 5 radians in the first row and 0, 1, 3, 4.5, 6 radians in the second row, respectively, which together with $2\pi $ do not have any common divisors. Using these two sets of phase shifts, the measurement results of Greivenkamp’s algorithm are illustrated in Figs. 5(c) and 5(g), respectively. The large artifacts in them show that Greivenkamp’s algorithm, when using non-uniform phase shifts, has little effect in suppressing the influence of fringe harmonics. Using the same sets of non-uniform phase shifts, we implement the proposed method. Here, we truncate the fringe signals and keep up to $7\textrm{th}$ order of harmonics. Its measurement results are shown in Figs. 5(d) and 5(h), where the ripple-like artifacts caused by harmonics have almost disappeared. These results demonstrate that the additional constraint in Eq. (11) properly models the characteristic of harmonics in this optical measurement technique, and that, with the suggested non-uniform phase shifts and the alternate iterative algorithm, the effects of aliasing have been removed.

The measurement results of measuring the second object are given in Fig. 6. Its layout and the used parameters are the same as those in Fig. 5, but the object has a larger depth and a more complex shape. The phenomena observed from Fig. 6 are similar to those from Fig. 5, confirming the validity of the proposed technique in restraining effects of harmonics once again.

Through Figs. 5 and 6, we have checked the performance of our proposed technique in suppressing the effects of harmonics on phase measuring in fringe projection technique, by comparing the appearances of its reconstructed depth maps with those of other two methods. To quantitively test how accurate this technique achieves, we do an additional experiment that measures a flat plate. Figure 7 shows the reconstructed profiles along the same cross-section perpendicular to the fringes on the measured plate. The three columns of this figure, from left to right, are the results of using Bruning’s, Greivenkamp’s, and our proposed methods, respectively. The parameters used for measurement, e.g., phase shifts, are the same as those used in Figs. 5 and 6, and the first and second rows correspond to the results of using four and five phase shifts, respectively. The measured profiles of using Bruning’s algorithm shown in Figs. 7(a) and 7(d) contain periodic errors having frequencies *N* times higher than fringes with *N* being the number of phase shifts. These errors are typically induced by frequency aliasing because of uniform sampling. When using Greivenkamp’s algorithm, large errors exist as we see from Figs. 7(b) and 7(e). It is known that Greivenkamp’s algorithm is contributive to enhancing the computational efficiency, but it has little effect in suppressing the influence of fringe harmonics, even if aliasing in the captured fringes has been avoided by using non-uniform phase shifts. The rightmost column shows the measurement results of using our proposed technique, where the periodic errors caused by frequency aliasing have been removed and the residual errors are induced by noise rather than by harmonics. By using the measurement results of this flat plate, we calculate the RMS values of the depth profiles deviating from a fitted plane and quantitatively list them in Table 2. Comparisons between the data in Table 2 demonstrate that this proposed method is superior in anti-aliasing over the existing techniques.

#### 5.2 Phase measuring in interferometry

As the aforementioned, the constraint in Section 3.1 is also suitable for modelling the effect of the nonlinear detector on the interferograms, provided the illumination is enough uniform. For examining this model, we experimentally measure the phase of practical interferograms using the proposed algorithm. The interferograms are recorded by using a microscopic interferometer which uses a light source having a wavelength of 530 nm and has a phase shifter driven by piezoelectric ceramic. The number of phase shifts is four and the nominal phase shifts are set to be $0$, $\pi /2$, $\pi $, and $3\pi /2$ radians. Each captured interferogram has a size of $512 \times 512$ pixels. When measuring a flat optical surface, the first interferogram is shown in Fig. 8(a).

By using the nominal phase shifts which are uniform over a $2\pi $ period, we employ Bruning’s algorithm to measure the phase map of the interferograms. The profile of the tested piece recovered from the phase map is given in Fig. 8(b). From this result, we observe that the ripple-like artifacts have a more complex appearance in comparison with those in Figs. 5(b) and 6(b). The reason for this phenomenon is that, in addition to the nonlinearity of detector, miscalibration of the phase shifter also induces phase errors. It is known that linear or some low-order nonlinear errors in phase shifts can be effectively coped with by using phase-shifting algorithms based on averaging strategy [5,40,41]. In this experiment, however, the phase shifter has large random errors and therefore their effect cannot be suppressed by using these algorithms. To solve this problem, we employ the self-calibrating phase-shifting algorithm [12] to estimate the true phase shifts as 0, 1.58, 3.65, and 5.17 radians for the four patterns in turn, and then calculate the phases using Greivenkamp’s least-squares algorithm. Its result is shown in Fig. 8(c), from which the periodic artifacts are most likely induced by the nonlinearity of detector. Note that the true phase shifts just estimated are non-uniform over a $2\pi $ period, so our proposed technique is used to eliminate the effects of harmonics. Our result is given in Fig. 8(d), from which it is evident the ripple-like artifacts have been suppressed significantly. Figure 9 compares the cross-sections of the measured profiles along the image diagonal roughly perpendicular to the fringes. The periodic errors with Bruning’s and Greivenkamp’s algorithms, shown in Figs. 9(a) and 9(b), respectively, are smoothed in Fig. 9(c) by using the proposed technique. Table 3 lists the RMS values of the profiles with different algorithms, where the proposed method gives the smallest value. These results experimentally demonstrate that our newly proposed method performs very well in interferometry to restrain the phase errors caused by the detector nonlinearities, especially when only very few interferograms are recoded.

#### 5.3 Advantages and limitations

This proposed method has a superior performance to others in restraining effects of harmonics because it avoids frequency aliasing in the early stage of measurement when recording fringe patterns. As comparisons, there have been many techniques using uniform phase shifts developed, which suppress nonlinearity-induced phase errors through data postprocessing [25–33]. Since frequency aliasing has occurred in their fringe patterns, these techniques generally require additional conditions or knowledge about fringes or phase errors besides the same constraint in Section 3.1regarding the identical characteristic of harmonics over fringe patterns. For example, they can depend on a one-parameter Gamma function for system nonlinearity [25,26], on a pre-calibrated look-up table for phase errors [27,28], or on an assumption of uniform distribution for phase probability [29]. Compared with them, this proposed method offers some new performances by taking advantage of anti-aliasing sampling. First, the additional conditions and knowledge just mentioned about fringes or phases are not required. Second, for each pixel, not only the phase of fundamental signal but also amplitude of each frequency component can be retrieved accurately, which is useful in some applications but not available from the uniform samples. Third, the underlying principle of this method is easy to extend to form new algorithms to adapt to special situations. As the aforementioned, with different measurement techniques, the fringe harmonics may have different characteristics. These characteristics can serve as the additional constraint instead of the one assumed in Section 3.1 thus forming new algorithms.

On the other hand, with this method, there exist some limitations. Firstly, this method focuses on dealing with the aliasing issue caused by sampling, and therefore it cannot remove the errors existing in the original fringe signals. For example, some physical factors, such as the illumination fluctuation of light source [42] or the miscalibration of phase shifter [12], will not only induce harmonics in a fringe signal but also distort the fundamental component of this signal. In this case, the phase error of the fundamental component cannot be eliminated by increasing the sampling rate or using non-uniform sampling. Secondly, this method requires an additional constraint about fringe harmonics to make the available equation system determined. Typically, the non-sinusoidal fringes caused by system nonlinearities or the defocused binary fringes can be modeled using the constraint defined in Section 3.1. However, this constraint cannot express the characteristics of the harmonics caused by fringe saturation [19,43] because they are highly localized.

With this method, there exist other issues worth discussing. For example, random noise in fringe patterns is always a main factor affecting performance of a phase measuring technique. Noise induced by complicated physical factors is generally modeled with a Gaussian distribution [44] according to the central limit theorem from the probability theory. Gaussian noise affects the proposed method mainly by degrading its capability of recognizing high order harmonics. It is known that Gaussian white noise has a flat power spectrum, but the intensities of harmonics rapidly decline with their orders. This fact means that the fringe patterns have significantly decreased signal-noise-ratios (SNRs) in high frequency band. For removing noise, a low-pass smoothing filter is usually used, but it is not suitable for use with our method because the low-pass filter simultaneously removes the high order harmonics. An alternative solution is to capture more images and average them. This method decreases the variance of noise on the fringe patterns at the expense of time duration for image capturing.

Phase shift selection is an issue related to noise. According to the theory from numerical analysis, when the phase shifts are uniform over a 2$\pi $ period, the resulting system of linear equations like Eq. (8) is the best-conditioned. In other words, the solution of this system will have the least sensitivity to noise. In contrast, extremely uneven phase shifts may lead to biased values and uneven variance of the reconstructed phases. With this proposed method, however, a set of non-uniform phase shifts is required. Therefore, we have to trade off these two aspects, making the selected phase shifts non-uniform and simultaneously as dispersed as possible within a 2$\pi $ period. This requirement is easy to meet in practical measurement. When using very few (e.g., 4 or 5) phase shifts, simply selecting integers below 2$\pi $ as the phase shifts enables delivering satisfied measurement results. If more phase shifts are required, we can draw some random values, from the standard uniform distribution on the interval from 0 to 2$\pi $, as phase shifts instead.

This proposed method involves an alternate iterative least-squares fitting procedure. A concern about it is its computational efficiency. Using the experimental data related to Fig. 7 and Table 2 in Section 5.1, Fig. 10(a) investigates the convergence of the proposed method by plotting the curves of RMS errors against the number of iterations. It shows that, in comparison with using four phase shifts, using five phase shifts makes the algorithm converges slightly faster, and simultaneously, to smaller residual errors. Through only 20 iterations, the both curves sufficiently converge. Using a laptop Desktop-MAPV02 K with Intel Core i7-10710U CPU and 16 GB memory, the overall computational time for processing fringe patterns having $500 \times 500$ pixels, when the number of phase shifts are four and five, are 4.3 seconds and 4.9 seconds, respectively. Figure 10(b) plots the curves of profile RMS values in Section 5.2 against the number of iterations, from which we observe the similar convergence trend. Optimizing the algorithm is helpful for further improve the computational efficiency, for example, using Eq. (23) to estimate the harmonic coefficients or involving few pixels in computation in the intermediate iterations. In addition, using more advanced hardware also enables shortening the computational time.

## 6. Conclusion

In this paper, we have proposed an anti-aliasing phase-measuring technique by use of non-uniform phase-shifting. It is demonstrated that selecting non-uniform phase shifts, which together with $2\pi $ have no common divisors, enables avoiding aliasing phenomenon in the captured fringe patterns, and that the constraint about identical characteristic of harmonics, which is practical in optical measurement applications, allows us to achieve a determined equation system having a stable solution. By employing an alternate iterative least-squares fitting procedure, we solve this system and separate the harmonics from fundamental fringes. Simulation and experimental results demonstrate that this proposed method allows us to significantly suppress the ripple-like phase errors caused by aliasing and to estimate the fringe phases accurately from as few as 4 fringe patterns. This proposed method is suitable for measuring fringe phases with various optical measurement techniques including interferometry, moiré topography, fringe projection, and fringe deflectometry.

## Appendix

We use Newton method [37] to solve the nonlinear system in Eq. (10). Defining the left side of Eq. (10) to be ${F_n}$ with $n = 1,2, \cdots N$ and ${\mathbf w} = {\left[ {\begin{array}{ccccc} {{a_0}}&{{a_1}}& \cdots &{{a_{N - 2}}}&\phi \end{array}\; } \right]^\textrm{T}}$ to be the unknown vector, we solve the linear system

Note that ${\mathbf J}({\mathbf w} )$ may become singular depending on the phase $\phi $. This phenomenon implies that the nonlinear system ${\mathbf F}({\mathbf w} )= 0$ does not always have stable solutions. For revealing this fact, we consider the determinant of this Jacobian matrix as

## Funding

National Natural Science Foundation of China (51975345).

## Disclosures

The authors declare no conflicts of interest.

## Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

## References

**1. **P. Carré, “Installation et utilization du comparateure photoélectrique et intertérentiel du Bureau International des Poids et Mesures,” Metrologia **2**(1), 13–23 (1966). [CrossRef]

**2. **J. H. Bruning, D. R. Herriott, J. E. Gallagher, D. P. Rosenfeld, A. D. White, and D. J. Brangccio, “Digital wavefront measuring interferometer for testing optical surfaces and lenses,” Appl. Opt. **13**(11), 2693–2703 (1974). [CrossRef]

**3. **J. C. Wyant, C. L. Koliopoulos, B. Bhushan, and O. E. George, “An optical profilometer for surface characterization of magnetic media,” ASLE Trans. **27**(2), 101–113 (1984). [CrossRef]

**4. **J. C. Wyant, “Interferometric optical metrology: basic principles and new systems,” Laser Focus **18**(5), 65–71 (1982).

**5. **P. Hariharan, B. F. Oreb, and T. Eiju, “Digital phase-shifting interferometry: a simple error-compensating phase calculation algorithm,” Appl. Opt. **26**(13), 2504–2506 (1987). [CrossRef]

**6. **J. E. Greivenkamp, “Generalized data reduction for heterodyne interferometry,” Opt. Eng. **23**(4), 350–352 (1984). [CrossRef]

**7. **K. Creath, “Temporal phase measurement method,” in * Interferogram Analysis: Digital Fringe Pattern Measurement*, D. W. Robinson and G. Reid, eds. (IOP, Bristol, UK, 1993), pp. 94–140.

**8. **M. Chen, H. Guo, and C. Wei, “Algorithm immune to tilt phase-shifting error for phase-shifting interferometers,” Appl. Opt. **39**(22), 3894–3898 (2000). [CrossRef]

**9. **H. Guo and M. Chen, “Least-squares algorithm for phase-stepping interferometry with an unknown relative step,” Appl. Opt. **44**(23), 4854–4859 (2005). [CrossRef]

**10. **H. Guo, Z. Zhao, and M. Chen, “Efficient iterative algorithm for phase-shifting interferometry,” Opt. Lasers Eng. **45**(2), 281–292 (2007). [CrossRef]

**11. **P. Gao, B. Yao, N. Lindlein, K. Mantel, I. Harder, and E. Geist, “Phase-shift extraction for generalized phase-shifting interferometry,” Opt. Lett. **34**(22), 3553–3555 (2009). [CrossRef]

**12. **H. Guo and Z. Zhang, “Phase shift estimation from variances of fringe pattern differences,” Appl. Opt. **52**(26), 6572–6578 (2013). [CrossRef]

**13. **K. A. Stetson and W. R. Brohinsky, “Electro-optic holography and its application to hologram interferometry,” Appl. Opt. **24**(21), 3631–3637 (1985). [CrossRef]

**14. **J. J. Dirckx, W. F. Decraemer, and G. Dielis, “Phase shift method based on object translation for full field automatic 3-D surface reconstruction from Moiré topograms,” Appl. Opt. **27**(6), 1164–1169 (1988). [CrossRef]

**15. **T. Yoshizawa and T. Tomisawa, “Shadow moiré topography by means of the phase-shift method,” Opt. Eng. **32**(7), 1668–1674 (1993). [CrossRef]

**16. **Y. B. Choi and S. W. Kim, “Phase-shifting grating projection moiré topography,” Opt. Eng. **37**(3), 1005–1010 (1998). [CrossRef]

**17. **V. Srinivasan, H. C. Liu, and M. Halioua, “Automated phase-measuring profilometry of 3-D diffuse objects,” Appl. Opt. **23**(18), 3105–3108 (1984). [CrossRef]

**18. **S. Zhang, “Recent progresses on real-time 3d shape measurement using digital fringe projection techniques,” Opt. Lasers Eng. **48**(2), 149–158 (2010). [CrossRef]

**19. **H. Guo and B. Lü, “Phase-shifting algorithm by use of Hough transform,” Opt. Express **20**(23), 26037–26049 (2012). [CrossRef]

**20. **C. Zuo, S. Feng, L. Huang, T. Tao, W. Yin, and Q. Chen, “Phase shifting algorithms for fringe projection profilometry: a review,” Opt. Lasers Eng. **109**, 23–59 (2018). [CrossRef]

**21. **M. C. Knauer, J. Kaminski, and G. Hausler, “Phase measuring Deflectometry: a new approach to measuring specular free-form surfaces,” Proc. SPIE **5457**, 366–376 (2004). [CrossRef]

**22. **H. Guo, P. Feng, and T. Tao, “Specular surface measurement by using least squares light tracking technique,” Opt. Lasers Eng. **48**(2), 166–171 (2010). [CrossRef]

**23. **L. Huang, M. Idir, C. Zuo, and A. Asundi, “Review of phase measuring deflectometry,” Opt. Lasers Eng. **107**, 247–257 (2018). [CrossRef]

**24. **X. Liu, Z. Zhang, N. Gao, and Z. Meng, “3D shape measurement of diffused/specular surface by combining fringe projection and direct phase measuring deflectometry,” Opt. Express **28**(19), 27561–27574 (2020). [CrossRef]

**25. **H. Guo, H. He, and M. Chen, “Gamma correction for digital fringe projection profilometry,” Appl. Opt. **43**(14), 2906–2914 (2004). [CrossRef]

**26. **K. Liu, Y. Wang, D. L. Lau, Q. Hao, and L. G. Hassebrook, “Gamma model and its analysis for phase measuring profilometry,” J. Opt. Soc. Am. A **27**(3), 553–562 (2010). [CrossRef]

**27. **S. Zhang and S.-T. Yau, “Generic nonsinusoidal phase error correction for three dimensional shape measurement using a digital video projector,” Appl. Opt. **46**(1), 36–43 (2007). [CrossRef]

**28. **B. Pan, Q. Kemao, L. Huang, and A. Asundi, “Phase error analysis and compensation for nonsinusoidal waveforms in phase-shifting digital fringe projection profilometry,” Opt. Lett. **34**(4), 416–418 (2009). [CrossRef]

**29. **X. Yu, Y. Liu, N. Liu, M. Fan, and X. Su, “Flexible gamma calculation algorithm based on probability distribution function in digital fringe projection system,” Opt. Express **27**(22), 32047 (2019). [CrossRef]

**30. **S. Ma, C. Quan, R. Zhu, L. Chen, B. Li, and C. J. Tay, “A fast and accurate gamma correction based on Fourier spectrum analysis for digital fringe projection profilometry,” Opt. Commun. **285**(5), 533–538 (2012). [CrossRef]

**31. **F. Lü, S. Xing, and H. Guo, “Self-correction of projector nonlinearity in phase-shifting fringe projection profilometry,” Appl. Opt. **56**(25), 7204–7216 (2017). [CrossRef]

**32. **S. Xing and H. Guo, “Directly recognizing and removing the projector nonliearity errors from a phase map in phase-shifting fringe projection profilometry,” Opt. Commun. **435**(15), 212–220 (2019). [CrossRef]

**33. **C. Jiang, S. Xing, and H. Guo, “Fringe harmonics elimination in multi-frequency phase-shifting fringe projection profilometry,” Opt. Express **28**(3), 2838–2856 (2020). [CrossRef]

**34. **J. Jiang and H. Guo, “One-shot dual-projection topography enhanced by phase-shifting logical moiré,” Appl. Opt. **60**(19), 5507–5516 (2021). [CrossRef]

**35. **H. Guo and M. Chen, “Fourier analysis of the sampling characteristics of the phase-shifting algorithm,” Proc. SPIE **5180**, 437–444 (2003). [CrossRef]

**36. **A. J. Jerri, “The Shannon sampling theorem—Its various extensions and applications: A tutorial review,” Proc. IEEE **65**(11), 1565–1596 (1977). [CrossRef]

**37. **A. Quarteroni, R. Sacco, and F. Saleri, * Numerical mathematics* (Springer Science, 2007), pp. 285–331.

**38. **H. Guo, H. He, Y. Yu, and M. Chen, “Least-squares calibration method for fringe projection profilometry,” Opt. Eng. **44**(3), 033603 (2005). [CrossRef]

**39. **S. Xing and H. Guo, “Iterative calibration method for measurement system having lens distortions in fringe projection profilometry,” Opt. Express **28**(2), 1177–1195 (2020). [CrossRef]

**40. **J. Schwider, O. R. Falkenstoerfer, H. Schreiber, A. Zoeller, and N. Streibl, “New compensating four-phase algorithm for phase-shift interferometry,” Opt. Eng. **32**(8), 1883–1885 (1993). [CrossRef]

**41. **J. Schmit and K. Creath, “Extended averaging technique for derivation of error-compensating algorithms in phase-shifting interferometry,” Appl. Opt. **34**(19), 3610–3619 (1995). [CrossRef]

**42. **Y. Lu, R. Zhang, and H. Guo, “Correction of illumination fluctuations in phase-shifting technique by use of fringe histograms,” Appl. Opt. **55**(1), 184–197 (2016). [CrossRef]

**43. **S. Feng, C. Zuo, L. Zhang, W. Yin, and Q. Chen, “Generalized framework for non-sinusoidal fringe analysis using deep learning,” Photonics Res. **9**(6), 1084–1098 (2021). [CrossRef]

**44. **H. Guo, “A simple algorithm for fitting a Gaussian function,” IEEE Signal Process. Mag. **28**(5), 134–137 (2011). [CrossRef]