Small-block frequency estimators to control step size in phase-shifting interferometry

Consider an interference signal

$$\begin{align} I_t = A \cos\left( \omega t + \phi \right) + B\end{align} \tag{ 1 }$$

oscillating with amplitude A and phase φ around a background or mean level B. The signal I_t is sampled at N different values of t with unit spacing, e.g. $t \in \{0, 1, \ldots, N-1\}$. In a PSI measurement, t labels the different samples and ω is the phase step from one sample to the next. More generally, we can view ω as the frequency of an arbitrary harmonic signal of which we know discrete values recorded with unit sampling period. If the instrument records data at several pixels, I_t , A, B, φ, and ω may all vary from one pixel to another. We make no assumptions about this variation, considering a single pixel at a time. Our goal is to estimate ω based on the set of samples $\{I_t\}$, without prior knowledge of A, B, or φ. Note that the Nyquist limit confines the estimate to the range $0 \lt \omega \lt {\pi}\,\mathrm{rad}$.

Equation (1) has four unknowns and so with N = 4 samples there is a unique solution for ω, given by Carré [2] as

$$\begin{align} \tan^2 \frac{\omega}{2} = \frac{3\left(I_1 - I_2\right) - \left(I_0 - I_3\right)}{\left(I_1 - I_2\right) + \left(I_0 - I_3\right)}.\end{align} \tag{ 2 }$$

Any expression for $\tan^2 \frac{\omega}{2}$ can be rewritten as an expression for $\cos\omega$, and in the latter form Carré’s solution becomes

$$\begin{align} \cos\omega = \frac{I_3- I_2 + I_1 - I_0}{2\left(I_2 - I_1\right)} = \frac{R_0}{S_0}\end{align} \tag{ 3 }$$

where, for later convenience, we define the shorthand

$$\begin{align} R_t = & I_{t+3} - I_{t+2} + I_{t+1} - I_{t}\end{align} \tag{ 4 }$$

$$\begin{align} S_t = & 2\left(I_{t+2} - I_{t+1}\right).\end{align} \tag{ 5 }$$

For certain phases of the signal, equations (2) and (3) yield an indeterminate result $(0/0)$. This was noted already by Hariharan et al [19], and as shown in figure 1 it is unavoidable: when the four samples are arranged symmetrically about an extremum of the oscillation, there are only two independent signal values ($I_0 = I_3$ and $I_1 = I_2$) but three unknowns (ω, A, and B; the phase being fixed by the symmetry), and oscillations of any frequency can interpolate the four samples exactly. In practice, this means that step-size estimates based on equation (2) or (3) become extremely sensitive to noise when the phase is near such a point of unfavourable symmetry. Thus, with evenly-spaced samples and an unknown backgound level B, at least five samples are needed for reliable step-size estimation. Since PSI systems using five samples [19, 20] are used in practical metrology [3, 10], we aim to solve the problem with N = 5.

Zoom In Zoom Out Reset image size

Estimating the frequency of a harmonic signal from discrete samples is a classic problem in signal processing. The maximum-likelihood frequency estimate for an isolated Fourier component can be found by locating a peak of the periodogram (discrete-time Fourier transform) [26], and computationally-efficient approximations for doing so based on the discrete Fourier transform (DFT) have been developed [27–31]. However the signal of equation (1) contains three Fourier components at $+\omega$, 0, and $-\omega$. Since N = 5 samples cover, at best, slightly over two periods of the oscillation, the DFT cannot cleanly resolve these three peaks: each is within two frequency bins (at most) of another. For the same reason, B cannot be eliminated simply by subtracting the average of the sample intensities: that average is the zero-frequency bin of the DFT, and is contaminated by spectral leakage from the signal peak. Spectral leakage also distorts the Fourier transform in the vicinity of the signal peak at $+\omega$, so frequency estimators designed for a single well-isolated Fourier component break down with a real signal and N = 5. A related problem afflicts frequency estimators based on fitting the phase of a complex signal [32–35]. The Hilbert transform used to obtain the complex signal from real data is distorted by spectral leakage from the $+\omega$ component to negative frequencies [36], corrupting the resulting frequency estimates.

There are multi-tone frequency estimation approaches suitable for poorly-resolved Fourier components [37–39], but many are too general for our purposes. The signal of equation (1) has a particular spectral structure, with one Fourier component known to be at zero frequency and the other two known to be a symmetric pair at ±ω. Estimators which assume nothing about the frequencies of the different Fourier components require additional data in order to rediscover these facts. By exploiting this structure in the design of the estimator we can achieve better performance for a given number of samples.

Specialized frequency estimators for real harmonic signals using four or five samples have been developed for the background-free case B = 0, including specializations of Prony’s method, discrete energy separation algorithms, the modified covariance method, and reformed Pisarenko harmonic decomposition [40–44]. One fruitful approach has been to describe the samples with a second-order auto-regressive or auto-regressive-moving-average model [45] and derive a closed-form estimator from the associated linear predictor [41–44]. A second-order model cannot describe the signal of equation (1) with an unknown B ≠ 0, but a third-order model can.

We therefore derive an estimator similar to the modified covariance estimator [41, 42, 44, 46], adapting it to a third-order model. Note that, for a given frequency ω, the signal of equation (1) obeys the linear recurrence

$$\begin{align} I_t = \left(1 + 2\cos\omega\right)\left(I_{t-1} - I_{t-2}\right) + I_{t-3}.\end{align} \tag{ 6 }$$

This relation holds for arbitrary A, B, and φ, and can be used to predict I_t based on the preceding three samples and an estimate of the frequency ${\widehat\omega}$. The prediction error is

$$\begin{align} e_t = I_t - \left[\left(1+2\cos{\widehat\omega}\right)\left(I_{t-1} - I_{t-2}\right) + I_{t-3}\right].\end{align} \tag{ 7 }$$

With N samples we can make and test N − 3 predictions. The estimate ${\widehat\omega}$ which minimizes the sum of squared prediction errors satisfies

$$\begin{align} \cos{\widehat\omega} &= \frac{ \sum_{t = 0}^{N - 4}\left(I_{t+2} - I_{t+1}\right) \left(I_{t+3} - I_{t+2} + I_{t+1} - I_{t}\right)} {2\sum_{t = 0}^{N - 4} \left(I_{t+2} - I_{t+1}\right)^2}\end{align} \tag{ 8 }$$

$$\begin{align} &= \frac{\sum_{t = 0}^{N - 4} S_t R_t}{\sum_{t = 0}^{N - 4} S_t^2}.\end{align} \tag{ 9 }$$

For N = 4 this reduces to the Carré solution of equation (3), as it must. For $N\unicode{x2A7E}5$, however, equation (9) yields well-defined estimates for any phase φ. It is easy to check that, for a signal obeying equation (1), these estimates are independent of A, B, and φ.

A different approach to the problem takes the well-known Carré solution as a starting point. Given N = 5 samples, we have two overlapping four-sample subsequences, namely the first four and last four samples. We can construct weighted combinations of the Carré solutions for these two subsequences as

$$\begin{align} \cos{\widehat\omega} = \frac{\alpha R_0 + \left(1-\alpha\right) R_1}{\alpha S_0 + \left(1-\alpha\right) S_1} = w \frac{R_0}{S_0} + \left(1-w\right) \frac{R_1}{S_1}\end{align} \tag{ 10 }$$

with α or w specifying the relative weighting of the two solutions in producing the combined estimate. Either weighting coefficient can be derived from the other:

$$\begin{align} w = & \frac{\alpha S_0}{\alpha S_0 + \left(1-\alpha\right) S_1}\end{align} \tag{ 11 }$$

$$\begin{align} \alpha = & \frac{w S_1}{w S_1 + \left(1-w\right) S_0}.\end{align} \tag{ 12 }$$

If we choose $\alpha = w = 0$ or $\alpha = w = 1$, we obtain the usual Carré solution based respectively on the first or last four samples respectively. Choosing $\alpha = 1/2$ or $w = S_0 / (S_0 + S_1)$ yields another four-sample estimator, suggested by Hariharan et al [19], which uses the first two and last two samples:

$$\begin{align} \cos{\widehat\omega} = \frac{R_0 + R_1}{S_0 + S_1} = \frac{I_4 - I_0}{2\left(I_3 - I_1\right)}.\end{align} \tag{ 13 }$$

As they noted, this estimator also sometimes gives indeterminate results. In fact, for any fixed choice of the coefficient α, the denominator $\alpha S_0 + (1-\alpha) S_1$ will vanish for some phase of the oscillating signal, yielding ill-defined and noise-sensitive frequency estimates. To avoid this problem we can choose $\alpha = S_0 / (S_0 + S_1)$ or $w = S_0^2 / (S_0^2 + S_1^2)$. This ensures that, whenever one of the two subsequences approaches a degenerate point as shown in figure 1, the corresponding estimate has vanishing weight in the combination. The result is

$$\begin{align} \cos{\widehat\omega} = \frac{S_0 R_0 + S_1 R_1}{S_0^2 + S_1^2},\end{align} \tag{ 14 }$$

which is just the N = 5 case of equation (9). Thus the estimator we derived from linear prediction can also be understood as a weighted average of Carré estimates, with the weights chosen to ignore indeterminate or noise-sensitive estimates.

The fact that the two approaches yield the same estimator is not a coincidence. In section 2 we estimated $\cos{\widehat\omega}$ by minimizing the squared prediction errors of equation (7) across all four-sample subsequences, where each subsequence contains three samples used to make the prediction and one to test it. The solution to such a single-parameter linear-least-squares optimization is generically a weighted average of the solutions that would be computed by considering each prediction individually, with weights proportional to the square of the sensitivity of the error to the fit parameter. At points of unfavourable symmetry the prediction error is independent of $\cos{\widehat\omega}$ and this sensitivity (which is proportional to S_t ) vanishes. Thus the linear-prediction formalism imposes automatically the constraint we imposed deliberately in this section, of giving no weight to estimates which provide no information on the step size.

For ideal, noiseless data obeying equation (1), the estimator of equation (9) is exact. For real, imperfect data, the uncertainty of the estimate will depend on the nature and scale of the imperfection. As a first test case, we simulated independent normally-distributed additive noise on each sample. Figure 2 shows the resulting root-mean-square estimation error (i.e. standard uncertainty) as a function of ω. The simulations are for a signal-to-noise ratio of 40 dB, using N = 5 samples for the estimate. Each point is an average of 6284 simulated estimates with uniformly distributed φ and independent realisations of the noise.

Zoom In Zoom Out Reset image size

If B is known and can be subtracted from the signal, then the modified covariance estimator [41, 42, 44, 46] is directly applicable (brown empty squares), and yields results close to the Cramér-Rao bound for this scenario [47] (black dashed line). Applying the modified-covariance estimator to data with B ≠ 0 yields large errors, even if the average of the samples is subtracted in a preliminary step (red filled squares), showing that the assumption B = 0 was crucial to its good performance. The estimator of equation (9) was derived without this assumption, and its performance (black crosses) is independent of B. Ignorance of B makes this a harder estimation problem, and the corresponding Cramér-Rao bound is higher (black solid line), particularly for small step sizes. Our proposed estimator’s performance is generally close to this benchmark.

The estimator proposed here, like the modified-covariance estimator, is statistically consistent but slightly biased [43, 44, 48], particularly near ω = 0 and $\omega = \pi$ where the non-linear behaviour of the cosine function distorts the error distribution. It is this bias which allows the estimators to outperform the Cramér-Rao bounds in some cases, since the bounds are strictly valid only for unbiased estimators. The bias could be removed by a non-linear transformation [43, 44], at the cost of increasing the variance of the estimate. In this work we omit this extra step and accept the bias in exchange for reduced variance and lower overall RMS error.

Figure 3 compares the effect of different noise types on the proposed estimator, at a consistent signal-to-noise ratio of 40 dB. Additive noise (black crosses, same data as in figure 2) corresponds to random fluctuation of the nominally-constant parameter B, here with a standard deviation of $10^{-2}\cdot A/\sqrt{2}$. Multiplicative noise (blue Xs) corresponds to fluctuations of A with a standard deviation of $10^{-2}\cdot A$, and phase noise (green dots) to fluctuations of φ with a standard deviation of ${10^{-2}}\,\mathrm{rad}$. From a signal-processing perspective, multiplicative noise is in-phase with the interference signal, phase noise is in quadrature with the interference signal, and additive noise is a mix of both. We thus expect the estimator’s response to additive noise to be a good guide to its overall performance. It is interesting to note, however, that for small step sizes $\omega\lesssim{0.4\pi}\,\mathrm{rad}$ it is in-phase noise which dominates the uncertainty of the estimate, while for large step sizes $\omega\gtrsim{0.8\pi}\,\mathrm{rad}$ it is quadrature noise which has the most effect. In general, the estimation uncertainty is lowest near $\omega\approx{0.6\pi}\,\mathrm{rad}$ and grows substantially worse for small step sizes (ω → 0), where the signal is almost constant across the samples and the differences between them are dominated by noise. This is not a particular feature of this estimator, but a general consequence of the amount of frequency information encoded in the N = 5 samples, as seen in the Cramér-Rao bound (figure 2, black line). The inset to figure 3 plots the RMS error as a function of signal-to-noise ratio around $\omega = {(\pi/2)}\,\mathrm{rad}$, showing that the estimation uncertainty is inversely proportional to the signal-to-noise ratio and converges to zero in the limit of noise-free signals.

Zoom In Zoom Out Reset image size

With very noisy data, equation (9) can yield meaningless estimates with $\lvert \cos{\widehat\omega} \rvert \gt 1$, especially at small step sizes (ω → 0). Under these conditions step-size estimation performance will be poor no matter how the data are handled. This is different from the situation with the Carré estimator, which can give indeterminate results even with a perfect noiseless signal. The method used to handle meaningless estimates has little impact on practical performance. Here we simply clamp the estimate of $\cos{\widehat\omega}$ to the range $[-1, 1]$ before inverting the cosine function, but we have checked that other procedures, including simply discarding the out-of-range estimates, have no noticeable effect on the plots of figures 2 and 3.

Many other imperfections are possible besides uncorrelated random noise. While a full study of the impact of arbitrary signal distortions is beyond the scope of this work, we show several illustrative examples in figure 4, using the same axes as figures 2 and 3 for easier comparison. The plot shows the RMS estimation error, averaging over signal phases, for three noiseless but distorted signals: one where the background level B drifts linearly at a rate of $0.04\cdot A$ from one frame to the next (black solid line), one where the amplitude A drifts linearly at the same rate (blue dashed line), and one with higher harmonics where the harmonic amplitudes fall off by a factor of 0.04 at each order (green chain-dotted line). This last case approximates the distortion expected due to multi-wave interference in a Fizeau interferometer with an uncoated glass beam splitter and a highly-reflective end mirror. Using equation (9) with a higher N improves the robustness to distortions, but where significant signal distortion is expected it may be useful to develop a specialized estimator for the expected signal form [49].

Zoom In Zoom Out Reset image size

The NRC Gauge-Block Interferometer is used for the calibration of gauge blocks and end bars as part of the dissemination of the SI metre. In 2016–2017 the instrument was upgraded to use PSI for the readout of interferometric phase, with the phase steps generated by a piezo-driven stage that displaces the interferometer’s reference mirror. Because the reference mirror tilts somewhat as it travels, the effective phase step varies by as much as $\pm{25}{\%}$ across the beam. To reduce our sensitivity to this inhomogeneous step size we use Servin’s broadband 9-frame PSI algorithm [15], which can accomodate wide variations from the nominal step $\omega_0 = {(\pi/2)}\,\mathrm{rad}$. In our setup there is also a noticeable non-linearity in the response of the piezo to its control voltage.

The upper-left image in figure 5 shows the typical fringe pattern seen in a single frame (camera image) of the interferometer. Servin’s PSI algorithm provides estimates of the interference amplitude (upper right) and phase (not shown). We can apply the estimator of equation (9) to any five-frame subsequence to obtain a map of step size vs. position (lower left, using the first five frames). This map clearly shows the overall gradient of step size due to the tilt of the reference mirror, as well as a pattern of lines parallel to the interference fringes. The latter correspond to phase-dependent errors associated with uneven step sizes. By comparing the estimates ${\widehat\omega}$ for different five-frame subsequences, we can model and correct the non-linearity in the piezo response to generate more uniform steps. Since the estimates are computed from the same set of control voltages as would be used in a phase measurement, hysteresis or path-dependent dynamics in the actuator are included in the estimates ${\widehat\omega}$ and can be corrected along with other repeatable non-linearities. Alternatively, we can use all nine frames to estimate the average step-size for the whole scan (lower right), where equation (9) is used directly with N = 9, reducing the sensitivity to noise.

Zoom In Zoom Out Reset image size

Prior to this work, the response of the piezo-driven mount to the control voltage v was calibrated using non-linear least-squares fitting to the model $I = B + A\cos\left( \psi\left( v \right) + \phi \right)$, where we modelled A, B, and φ as per-pixel parameters and $\psi\left( v \right)$ as a single cubic polynomial for the entire camera image. Assuming uniform phase shifts across the $612\times514$ pixel camera frames amounts to neglecting the tilt of the mirror as it travels, and even with this simplification the fits would take over five minutes of computation time to converge, making it impractical to monitor the step size during routine measurements. Image-based estimators such as those of [22–24] could have significantly reduced this comuptation time, but would still have required assuming a spatially uniform step size. After implementing the estimator of equation (9) the instrument can now report independent step-size estimates at every pixel and for every contiguous five-frame subset of the nine-frame PSI scan, thus allowing an assessment of the mirror tilt and as well as the shifter non-linearity all in less than two seconds of computation time. This check is now performed routinely for every measurement performed on the instrument, providing additional confidence in its repeatable and reliable operation.

When the background level B is unknown, a minimum of five samples are required for reliable estimation of the step-size in single-pixel PSI, or equivalently for the estimation of frequency from regularly sampled data. We have presented a reliable five-sample estimator, relating it both to the modified-covariance estimator and the broader tradition of auto-regressive frequency estimators on the one hand, and to the Carré step-size estimator familiar in the PSI community on the other. The new estimator has a modest and predictable computational cost while saturating the relevant Cramér-Rao bound in the presence of additive Gaussian noise. It can be used to measure and correct step-size errors in practical phase-shifting systems, as well as for ongoing performance monitoring and process control. The NRC Gauge Block Interferometer takes advantage of these possibilities to maintain confidence in the reliable performance of its PSI implementation.

The estimator presented here is implicitly optimized for additive Gaussian noise. It should be possible to design other estimators optimized for different signal imperfections such as higher harmonics. We are not aware of a systematic theory of step-size estimation as complete as the existing frameworks for robust phase estimation [14–18], and developing such a theory would be an interesting direction for future work.