{'fraction_non_alphanumeric': 0.03517969958319134, 'fraction_numerical': 0.006343880252706635, 'mean_word_length': 4.702242152466368, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Nonlinear Digital Post-Processing to Mitigate Jitter in Sampling

This paper describes several new algorithms for estimating the parameters of a periodic bandlimited signal from samples corrupted by jitter (timing noise) and additive noise. Both classical (non-random) and Bayesian formulations are considered: an Expectation-Maximization (EM) algorithm is developed to compute the maximum likelihood (ML) estimator for the classical estimation framework, and two Gibbs samplers are proposed to approximate the Bayes least squares (BLS) estimate for parameters independently distributed according to a uniform prior. Simulations are performed to demonstrate the significant performance improvement achievable using these algorithms as compared to linear estimators. The ML estimator is also compared to the Cramer-Rao lower bound to determine the range of jitter for which the estimator is approximately efficient. These simulations provide evidence that the nonlinear algorithms derived here can tolerate 1.4-2 times more jitter than linear estimators, reducing on-chip ADC power consumption by 50-75 percent.


I. INTRODUCTION
Analog-to-digital converters (ADCs) are affected by different types of uncertainty that introduce noise into the samples of a signal. Two important causes of this noise are inaccurate sample times and variation in the measured signal amplitude. Many well-known approaches exist for separating noise due to amplitude error from the signal. Similarly, noise in the sample times, known as jitter, is also a well-studied phenomenon. However, the effect of jitter is assumed to be negligible in practice because of the low phase noise sampling circuitry employed in conventional ADCs.

A. Motivation
The increasing demand for ultra low power ADCs has renewed interest in mitigating jitter using digital post-processing techniques. Because the digital portion of mixed-signal systems like digital sensors and wireless receivers continues to shrink, the analog portion, including the clock generator for the ADC, dominates the size and power consumption of the chip. Such systems would benefit greatly from smaller, more power-efficient analog circuitry; however, the use of such circuitry, with lower voltages, significantly increases noise, such as jitter in the clock signal. The power consumed by an ADC is proportional to the desired accuracy and sampling rate [1]. In [2], the speed-power-accuracy tradeoff for high-rate ADCs is shown to satisfy Speed × (Accuracy (rms)) 2 Power ≈ constant.
Both [3] and [4] demonstrate that doubling the standard deviation of the jitter cuts down the effective number of bits (accuracy (rms) = 2 ENOB ) by one, so the power consumption would need to increase by a factor of four to achieve the same accuracy as before. Clearly, the impact of jitter is not negligible.
Incorporating a post-processing method that can recover some of this loss in accuracy should drastically reduce the power consumption of most high-speed ADCs. In addition, the ability to tolerate larger quantities of jitter allows the use of frequency-modulated clocks, which may yield additional benefits like lower electromagnetic interference and radiation [5].

B. Problem Formulation
In this article, nonlinear algorithms for improved classical and Bayesian estimation are developed to recover the parameters of a bandlimited signal. We restrict our attention to signals with a finite basis, leaving the development of streaming algorithms for the infinite-dimensional case a path for further investigation. In particular, we consider a finite segment of a periodic bandlimited signal with period K; the chosen basis is one of a periodic sinc function psinc K (t) ∆ = sin(πt) K sin(πt/K) and its integer shifts. One could also construct a suitable example by applying a smooth windowing function, rather than assuming periodicity. These algorithms are applicable for other signal classes, like splines. For this work, assume the Nyquist frequency is normalized to π, and the signal x(t) is sampled with an oversampling factor of M . The nth sample, y n , is described by the observation model: Let the jitter z n be zero-mean white Gaussian noise, with variance σ 2 z ; note that the jitter is normalized relative to the critical sampling period. The dominant source of additive noise is assumed to be a signalindependent external noise source (not quantization noise); this additive noise, denoted by w n , is also assumed to be zero-mean white Gaussian noise, with variance σ 2 w ; the additive noise is normalized to the scale of the parameters x k . All the noise sources are independent, and they are independent of the input signal parameters x k . Define [H(z)] n,k = psinc K (n/M + z n − k); then, To keep notation compact, denote by p(x), p(y; x), and p(y | x), the probability density function (pdf) of x, the pdf of y parameterized by the deterministic variable x, and the pdf of y conditioned on the random variable x, respectively. In this article, the random variable(s) associated with a pdf are given explicitly only when necessary to avoid ambiguity. Denote the d-dimensional multivariate uniform and normal distributions by and Then, the likelihood function corresponding to the observation model is classical: p(y; x) The classical estimation problem is to find the estimatorx(y) that minimizes the mean squared error (MSE), knowing only the samples y and the likelihood function: Computing this minimum MSE estimator (if such an estimator actually exists) is rarely straightforward; the maximum likelihood (ML) estimator is easier to approximate and will be used instead. The Bayesian estimation problem has the same objective, but with the additional knowledge of a prior on x: x(y) = arg min The Bayes least squares (BLS) estimator, E[x | y], is the solution to this minimization problem. However, neither the posterior density p(x | y) nor its mean have closed forms for the observation model in (3).

C. Related Work
The early literature focuses on linear reconstruction for signals affected by different types of random jitter: [6] develops a linear interpolation filter for jittered signals, and [7] explores the effects of sampling and reconstruction jitter. The sampling error due to jitter is examined in [8]: using low pass interpolation, the MSE is approximately linear in σ 2 z . In particular, when the jitter is Gaussian and small enough, and the input power spectral density S xx (jΩ) = 1 2ΩB is flat, the MSE is approximately 1 3 Ω 2 B σ 2 z . However, these works provide limited insight into the effectiveness of general nonlinear techniques.
More recently, [9] provides an iterative weighted least squares fitting algorithm for computing a cubic spline approximation to a jittered input signal. However, the weighted least squares fit uses only a secondorder approximation of the input signal. More along the lines of the contributions herein, [10] constructs an Expectation-Maximization (EM) algorithm to estimate a signal sampled by time-interleaved ADCs with timing offsets; however, the offset for each ADC stays constant over time, unlike jitter. Also, a Metropolis-Hastings algorithm is derived in [11] to estimate each sample's timing error and the jitter variance from a sequence of samples; this approach is similar to the Gibbs sampler used in this work.
Another related work, [12], implements several reconstruction methods for finite-length discrete-time signals with sparse coefficients; however, since the jitter is constrained to integer locations to preserve the structure of the DFT, combinatorial methods are employed in lieu of more general statistical approaches.
In the context of mitigating read-in and write-out jitter in data storage, a MAP-based estimator is proposed in [13] for jitter generated by a first-order Markov process and the bits corrupted by such jitter. December 14, 2009 DRAFT

D. Outline
In Section II, background on Gauss-Hermite quadrature and Monte Carlo methods is provided. For the classical estimation problem, a numerical method is developed in Section III to approximate the Cramér-Rao lower bound, and an Expectation-Maximization (EM) algorithm is presented to compute the ML estimate. Two Gibbs samplers are proposed in Section IV for the Bayesian case, one using rejection sampling, and the other using slice sampling. In Section V, simulation results are presented for all these algorithms, and conclusions and future directions for research are discussed in Section VI. The problem proposed here, along with further background material, is also discussed in detail in [14].

II. BACKGROUND
The likelihood function described in (6)  and (b) stochastic approximation, via rejection, Gibbs, and slice sampling.

A. Numerical Integration
Gauss quadrature approximates the univariate expectation f (x)w(x) dx by the sum n i=1 w i f (x i ), where x i and w i are fixed abscissas and weights based on the weighting function w(x), which is the pdf of X. The abscissas and weights can be precomputed for a choice of w(x) and n using the eigenvaluebased method in [15]. Quadrature generalizes to multiple input-multiple output functions f (x), but if the integrand is not separable, the computational complexity scales exponentially with the number of components in x.
One weighting function of particular interest is the pdf of the standard normal distribution. The related method of quadrature is known as Gauss-Hermite quadrature, since the abscissas and weights are derived from Hermite polynomials. Adapting the abscissas and weights for any univariate normal distribution is easy: where x i and w i are the abscissas and weights for the standard Normal weighting function. As discussed in [16], the error for Gauss-Hermite quadrature goes to zero exponentially fast as long as the nth derivative does not increase too fast with n; a sufficient condition for convergence is that f (x) does not increase more rapidly than e x 2 . in turns until the chain converges to the steady-state. Since the steady-state distribution is the joint pdf p(x), further samples from the transition probabilities will be generated as if they were generated by the joint distribution itself.

B. Stochastic Approximation
Stochastic approximation reduces the problem of calculating an expectation E[f (x)] to sampling from the distribution p(x) using methods like rejection sampling, Gibbs sampling, and slice sampling. If the normalizing constant for a target distribution is unknown, rejection sampling is still effective.
Suppose only the formp(x) is known, so P = p(x) dx is unknown (but assumed to be finite). Then, we choose c such that cq(x) >p(x), and rejection sampling works as before.
One disadvantage of rejection sampling is that if the proposal distribution is a poor approximation to the target distribution, the algorithm generally will reject a prohibitively large number of samples. The "envelope" requirement may make choosing a good proposal distribution very difficult. The expected number of iterations between accepted samples is c/P . The Gibbs and slice samplers overcome this issue by removing the enveloping requirement.
2) Gibbs Sampling: The Gibbs sampler is a Markov chain Monte Carlo method introduced in [17].
The desired sampling distribution p(x) is the stationary distribution of the constructed Markov chain.
Consider the joint probability density function p(x 1 , x 2 , . . . , x K ), and define a chain of random variables By augmenting these variables with the last K − 2 states to become x −2 , . . . , x −K , x −1 , the chain becomes a Markov chain; then, p(x) is the stationary distribution of this chain (see Figure 1).
Sufficient conditions for convergence of the underlying Markov chain are outlined in [18] and [19].
According to [20], separating correlated variables can slow convergence. The period until the chain converges is termed the "burn-in time." Once the chain has reached its steady-state, the generated samples can be treated as if they were generated from the joint distribution directly.
3) Slice Sampling: Slice sampling is a Markov chain Monte Carlo method described in [21] for sampling from p(x) by sampling uniformly from the region under the pdf. Considering the area under the pdf as a pair of random variables (X, Y ), where X is as before, and Y ∈ (0, p(X)), with joint distribution p(x, y). As in Figure 2, the joint pdf can be sampled by repeatedly drawing from the conditional distributions, like in Gibbs sampling. These two conditional densities correspond exactly to uniform distributions, the first p(y | x) being over a single interval, and the second p(x | y) over the "slice." Sampling from the slice {x : p(x) > y (i+1) } is generally difficult. If X is bounded, one could simply uniformly sample over the whole interval until the generated sample has sufficiently high probability p(x) to be in the slice. If the derivatives of p(x) are available, root-finding methods such as Newton's method or Halley's method may be used to locate boundaries of the slice. "Shrinkage," a simple accept-reject method to sample from the slice described in [21], is used in this paper.

III. CLASSICAL (NON-RANDOM) PARAMETER ESTIMATION
Minimizing the MSE without access to p(x) is difficult because the optimal estimator must effectively minimize the MSE for all possible values of x. Moreover, the MSE of a candidate estimator generally depends on x, so it is difficult even to evaluate the performance of any given estimator directly. As a basis of comparison, the Cramér-Rao bound (CRB) can provide a lower bound on the MSE performance of unbiased estimators, given that we can evaluate this bound for the noisy observation model.
The best linear unbiased estimator (BLUE) is explored as a first attempt to attain the Cramér-Rao bound.
To improve upon linear estimation, the ML estimator is derived, since it is guaranteed to asymptotically achieve the CRB as the oversampling factor, M , goes to infinity. While the likelihood function in (6) is separable, it has no closed form: where h T n (z n ) is the nth row vector of the matrix H(z). An EM algorithm is developed to simplify the problem of maximizing the likelihood function.

A. Cramér-Rao Lower Bound
The Cramér-Rao lower bound on the MSE for the class of unbiased estimators is equal to the trace of I y (x) −1 , where the Fisher information matrix I y (x) can be expressed as and l(x; y) ∆ = ln p(y; x) is the log-likelihood of y as a function x. An efficient estimator, which satisfies the CRB, can be found for the the no-jitter (z = 0) case using the general formula in [22]: However, for the random jitter case, the efficient estimator would depend on x, so it is not valid.
For the random jitter case, the lack of a closed-form for the likelihood function requires that the lower bound be approximated numerically. Because the likelihood function is separable, which means that the Fisher information in (11) can be re-written as The singleton likelihood function p(y n ; x) can be computed numerically using Gauss-Hermite quadrature: The derivative can also be approximated using quadrature: ∂p(y n ; x) ∂x Since the expression inside the expectation in (14) becomes a complicated non-linear function of y and x.
Since p(y n ; x) is approximately a Gaussian mixture, we can resort to simple stochastic approximation to approximate this expectation, and hence, the Fisher information matrix. For each data sample n, generate where The trace of the inverse of this matrix approximates Cramér-Rao lower bound.

B. Linear Estimation
It is well known that for any linear observation model y = Hx + w with zero-mean additive noise  (3) is the pseudoinversê The question naturally arises as to how to find the best linear unbiased estimator (BLUE), in terms of minimizing the MSE. In [22], the BLUE is derived for a fixed matrix H. As derived in [14], the BLUE for a linear observation model with random matrix H iŝ where the covariance matrix of the data Λ y depends on the value of the parameters: A sufficient condition for (21) to be valid is for Λ y to be a scalar matrix (a scalar times the identity matrix). In general, and for the observation model and prior on the jitter used here, the diagonal elements of Λ y are not equal, and Λ y is not a scalar matrix. For the no-jitter case, H(z) is deterministic, and the BLUE is equal to the efficient linear estimator in (12).
This efficient linear estimator for the no-jitter case is used as a baseline to compare the performance of the ML estimator proposed in the next section. This estimator would be optimal in the MSE sense if the jitter were removed from the observation model.

C. ML Estimation using the EM Algorithm
Operating under the assumption that a linear estimator is too restrictive to capture the complicated behavior of the random jitter observation model, but unable to formulate an efficient estimator directly, maximum likelihood estimation appears to be a promising alternative. However, maximizing the likelihood function is a difficult task, because the likelihood function is not convex in x, does not have a closed form, and is not separable. One approach would be to utilize a gradient method, or a hill-climbing method, to find a local maximum iterating over one variable at a time. One related method, investigated in [23], involves maximizing the joint likelihood p(y, z; x) by alternating maximizing p(z | y; x) and p(y | z; x).
Maximum likelihood estimation is also a classic application of the Expectation-Maximization (EM) algorithm, as described in [24]. For each iteration i, the EM algorithm consists of maximizing the expectation Q(x;x (i−1) ) = E log p(y, z; x) | y;x (i−1) with respect to the unknown parameters x; i.e.x (i) = arg max x Q(x;x (i−1) ). Here, z is the missing data. Repeated iterations of this step results in the maximizing x converging to a local maximum of the likelihood function.
Maximizing (23) will result in finding a local maximum of the likelihood function; however, nothing has been said about how fast the algorithm converges to a stationary point, or if this local maximum coincides with the true ML estimate. The rate of convergence depends on the choice of complete data [24].
Proven in [25] is the fact that the rate of convergence is worse if the CRB for the incomplete data set is much greater than the CRB for the augmented data set. The SEM algorithm in [26] also approximates the observed Fisher information, which can be used to establish the quality of the estimate. Since the likelihood function is not concave in x, the estimate depends on initial conditions, the effects of which are studied in [14].
If the values of the jitter were known, the problem would decompose into a linear one, and the solution is readily available. Then, selecting the jitter to be the missing data, Expanding and substituting into the expectation in (23) yields Since this equation is quadratic in x, the maximum valuex (i) satisfies the linear system of equations The Hessian matrix is negative-definite, so the equation is strictly concave, and this maximum is unique.
To solve this linear system of equations, the expectations in (26) need to be approximated. Using Bayes rule, If the jitter or additive noise were not independent (white), the expectation would not be separable. Using Gauss-Hermite quadrature (with I terms) on the left and right expectations, n,: and p(y n ; x (i−1) ) is computed again using (15). The complexity of performing these approximations is linear in the number of samples, so for a large number of parameters, the limiting step is actually solving the system once the expectations are computed.

IV. BAYESIAN ESTIMATION
The classical estimation framework is suitable for designing a general post-processor for a sampling system when nothing is known about the process generating the parameters of the input signal. However, if the generating process is known, Bayesian inference leverages that information to optimize the performance of the post-processor. The Bayes least squares (BLS) estimator minimizes the MSE, and approximating that estimator is the objective of this section.
How should the prior distribution be chosen? A white Gaussian process is attractive for its simplicity.
The conjugate prior is another prior employed for convenience. However, since the likelihood function is already intractable, a conjugate prior is unhelpful for this problem. The "least-informative" prior is another suitable choice, but in general, it is hard to find. The approach used here is the maximum entropy prior. The maximum entropy model is justified in that it imposes minimal structure on the parameters; in this way, the maximum entropy model is the most "challenging" prior to use. Therefore, the parameters should be independent, and the entropy-maximizing distribution over a finite interval is the uniform prior [27]. For simplicity, we will assume that the parameters lie between −1 and 1.
Like the likelihood function in the previous section, the resulting posterior density has no closed form.
It can be expressed in terms of the singleton likelihoods in (10):

A. LLS Estimation
The linear least squares (LLS) estimator is defined to be the linear estimator with minimum MSE. For the random jitter observation model in (3), the LLS estimator iŝ The expectations in (31) can be computed using Gauss-Hermite quadrature or some other numerical method. However, the optimal linear estimator has already been studied and does not show much improvement over the LLS estimator assuming no jitter is present. In [6], the optimal filter coefficients for a low pass interpolator are derived for the random jitter problem. As mentioned in the introduction, the MSE of a low pass interpolator is linear in the variance of the jitter [8].
When no jitter is assumed, the above estimator simplifies tô This linear estimator is the best linear operation that can be performed in the absence of jitter. Hence, the no-jitter LLS estimator is the baseline estimator against which the BLS estimators derived later will be measured. The error covariance of this estimator is

B. BLS Estimation using Gibbs Sampling
The BLS estimator is well-known to have the minimum MSE of any estimator for the Bayesian framework. However, optimal performance comes at a price: the BLS estimator involves taking the expected value of the posterior distribution, which is generally difficult to compute, as it is for this problem. However, stochastic approximation can be helpful here. Consider the BLS estimator for (x, z).
The expectation E[x, z | y] is an obvious application of Gibbs sampling. Sampling from p(x, z | y) directly is a terrible idea due to the high dimensionality of the sample space, even if it were trivial to generate samples. Instead, Gibbs sampling can be used to generate samples of one variable at a time, conditioned on the rest.
Once enough samples have been taken so that the Markov chain is sufficiently close to the stationary distribution, additional samples are averaged to approximate the BLS estimator. In the description below, I b represents the "burn-in time", the number of iterations until the Markov chain has reached its steady state, and I represents the number of samples to generate after convergence has been achieved.
Require: y, I, I b Consider generating samples z n from the distribution p(· | z −n , x, y). Using Bayes rule, and the independence of z n and w n , The full-conditional distribution of z n is independent of z −n , so sampling the jitter values can be easily grouped together. Because this functional form is enveloped by the prior on z n , rejection sampling is an obvious choice to produce samples. The proposal density is q(z n ) = N zn (0, σ 2 z ), and the scaling factor and if y 2 n > 2y n x 2 , c can be multiplied by exp − y 2 n −2yn x 2 2σ 2 w , and cq(z n ) would more tightly envelope (34).
Unlike the independent z n 's, the conditional distribution on x k does depend on the other parameters Denote by H k (z) and H −k (z) the columns of H(z) that multiply x k and x −k in (36), respectively. The truncated normal distribution in (36) can be modified to be a distribution explicitly on x k : where Sampling from a truncated normal distribution is a typical application of cdf inversion. When the mean µ k lies between −1 and 1, the cdf is almost linear, and the cdf inversion technique works well. However, when the mean lies far outside this range, the cdf is practically flat, and the technique is limited by the precision of the computer. A rejection sampling-based approach for producing samples of a truncated standard normal distribution is developed in [28]. For a standard normal distribution truncated to the with variance σ 2 z ) is also shown for reference. Note that even for values of z0 likely according to the proposal density q(z0), the target density is many orders of magnitude smaller, so rejection sampling will reject an extraordinarily large number of samples before accepting a sample.
interval [a, b], where a > 0, the proposal distribution should be an exponential distribution shifted to cover x ≥ a, whose optimal scale factor The case where [a, b] lies to the left of zero, generating samples is the mirror-image of this problem, so the optimal exponential distribution is flipped and shifted to cover x ≤ b, and α * = − 1 2 (b − √ b 2 + 4).
According to [28], a uniform distribution U x (a, b) is the preferred proposal distribution when a and b are too close together, so that a disproportionately large number of samples are not rejected due to the The parameters x k are all highly correlated, so one might wonder if it is possible to use a similar technique to generate realizations of the multivariate truncated normal distribution p(x | y, z) directly.
However, as [28] laments, the rejection-sampling method has no natural extension to multivariate truncated normal distributions when the variables are correlated. Gibbs sampling is recommended to reduce the problem down to univariate sampling, which is what is proposed here anyway. Sampling x k one at a time has the added benefit of extending naturally for a less trivial generating process for x.

C. Improvement via Slice Sampling
While the rejection sampler provides an exact method for producing realizations of a target distribution, a significant disparity between the shape of the proposal and target distributions causes a large number of samples to be rejected. In the case of rejection sampling employed in the Gibbs sampler to generate jitter values z n from the density (34), finding the smallest value for c is difficult, and the shape of the target distribution can vary depending on the parameters. Empirical evidence, shown in Figure 3, portrays the extent of the problem with rejection sampling, when σ z ≫ σ w . High oversampling only compounds this phenomenon.
Slice sampling is not susceptible to such problems since no tightly enveloping proposal density is necessary; the ability to evaluate an unnormalized form of the target distribution is sufficient. Thus, the expression in (34) can still be used. Each iteration of slice sampling consists of two uniform sampling problems: 1) Choose a slice u uniformly from [0,p(z 2) Sample z (i+1) n uniformly from the slice S ∆ = {z n :p(z n | y, x) ≥ u}.
The first step is trivial, since we are sampling from a single interval. The second step is more difficult.
However, since u ≤p(z n | y, x) for all z n in the slice, Solving for z n , the range of possible z n is bounded: Using these extreme points for the initial interval containing the slice, and the "shrinkage" method specified in [21] to sample from the slice by repeatedly shrinking the interval, slice sampling becomes a relatively efficient alternative to rejection sampling. The "shrinkage" method decreases the size of the This expectation is quadratic in x 0 , so the maximum occurs at the extreme point x 0 = (R + L)/2. The maximum value is Concavity implies that the minima are at the two endpoints x 0 = L and x 0 = R. In both cases, the expected size of the interval is (R − L)/2. Therefore, which implies that at worst, the size of the interval shrinks to 3/4 its previous size per iteration, on average. Then, given the initial interval [L 0 , R 0 ] and previous point x 0 , the expected size of the interval [L I , R I ] after I iterations of the shrinkage algorithm is If the target distribution p(x) is continuous, the algorithm is guaranteed to terminate once the search interval is small enough. Since the interval size shrinks exponentially fast, on average, the number of "shrinkage" iterations is approximately proportional to the log of the fraction of the initial interval contained in the slice. In a discussion included in [21], a binary search-like shrinkage algorithm is proposed that can converge faster on the slice than the algorithm used here. Incorporating such an approach to accelerate the slice sampler merits future investigation.
Using slice sampling to generate the jitter values when σ z is large relative to σ w improves the speed of the Gibbs sampler. However, the addition of new auxiliary variables through slice sampling can be expected to slow the Gibbs sampler's overall rate of convergence. Thus, both algorithms are included for simulation.

V. SIMULATION RESULTS
The objectives of this section are: (a) to determine when jitter significantly affects the samples of a signal and (b) to demonstrate that nonlinear post-processing can improve upon linear estimation of the signal parameters. The convergence behavior of all the included algorithms, as well as their sensitivity to initial conditions, is presented in [14]. This article focuses on the performance analysis of these algorithms, using Monte Carlo simulation to determine the MSE of each algorithm for the different parameter choices (specified in each figure). As in the rest of the paper, σ z is relative to the critical sampling period T = 1, and σ w is relative to the scale of the signal parameters x k ∼ U x (−1, 1).
A. Evaluating the EM Algorithm 1) Convergence: The convergence properties of the EM algorithm were discussed briefly in the background section. Intuitively, increasing M , increasing σ z , or decreasing σ w increases the additional information from knowing the missing data z, slowing the rate of convergence. This behavior is observed in the first simulations in [14]. In addition, these simulations suggest that the proposed algorithm converges exponentially fast.
The other major sticking point for using an EM algorithm is the local nature of the solution. If the likelihood function has many local maxima, the EM algorithm may not necessarily converge to the global maximum. Instead, the local maximum is determined by the choice of initial conditions; i.e. the initial value of the parameters x. Simulations in [14] suggest the sensitivity to initial conditions increases with the non-concavity of the likelihood function. Larger M , larger σ z , or smaller σ w increase the rate and severity of oscillations in the likelihood function, precipitating more local maxima. In cases when the sensitivity is an issue, using multiple random starting points or approaches like the deterministic annealing EM algorithm [29] can help alleviate the problem.

2) MSE Performance:
The EM algorithm is compared against two linear estimators. First, to demonstrate the MSE improvement attainable by not restricting ourselves to linear estimators, the EM algorithm is pitted against the linear unbiased estimator in (20). However, a major motivating factor for developing these algorithms is to reduce the power consumption due to clock accuracy. By comparing the EM algorithm against the optimal linear estimator for the no-jitter observation model, the EM algorithm can be shown to achieve the same MSE as the linear estimator for a substantially larger jitter variance, reducing the clock's power consumption.
When the additive noise dominates the jitter (σ z ≪ σ w ), the improvement can be expected to be minimal, since the system is nearly linear, and the jitter is statistically insignificant. As the amount of jitter increases, the density function p(z | y; x) used in each iteration of the EM algorithm becomes more non-linear in z, and the quadrature becomes less accurate for a given number of terms. Therefore, the EM algorithm generally takes longer to converge, and the result should be a less accurate approximation to the true ML estimator. This behavior can be observed in Figure 4, where the EM algorithm is compared More extensive plots can be found in [14].
against both the linear unbiased estimator in (20) and the efficient no-jitter linear estimator in (12). The linear unbiased estimator has lower MSE for higher M , and the EM algorithm generally has lower MSE than either linear estimator.
To answer the question of how much more jitter can be tolerated for the same desired MSE using the EM algorithm, the maximum proportional increase is plotted as a function of M and σ w in Figure 5. The maximum proportional increase for a choice of M and σ w is computed by approximating log-log domain MSE curves, like those in Figure 4, with piece-wise linear curves and interpolating the maximum distance between them over the range of σ z up to 0.5. The proportion of improvement increases logarithmically as M increases, and the improvement stays approximately the same for different values of σ w . However, the maximum improvement shown for σ w = 0.5 is low because no σ z > 0.5 is tested. The maximum improvement factor shown corresponds to the power consumption dropping by 50 to 75 percent.
As a final experiment, the Cramér-Rao lower bound is compared to the unbiased linear estimator (20) and EM algorithm to measure the efficiency of the algorithms. Although computational difficulties prevent a complete comparison for every possible value of x, carrying out a comparison for a few randomly chosen values of x provide a measure of the quality of the algorithms. As the curves in Figure 6 demonstrate for one such random choice of x, both algorithms are approximately efficient for small σ z , but the EM algorithm continues to be efficient for larger values of σ z than the linear estimator. Since the ML estimator is asymptotically efficient, as M → ∞, the EM algorithm should be efficient for all σ z with large enough oversampling.   Additional data can be found in [14].

B. The Gibbs Samplers 1) Convergence:
The Gibbs sampler is known to converge if the Markov chain is irreducible and aperiodic; however, the rate of convergence needs to be established. These results can be used to tune the "burn-in" time I b and the number of intervals I after the chain has sufficiently converged. In [14], the convergence rates of the Gibbs sampler and Gibbs sampler using slice sampling are simulated to determine trends in M , σ z , and σ w . The essential result shown is that the convergence rate increases when jitter or additive noise power decrease, or when the oversampling factor increases.
2) MSE Performance: The MSE performance of both Gibbs samplers is determined for a wide range of oversampling factors M , jitter variance σ 2 z , and additive noise variance σ 2 w . Like the EM algorithm, the Gibbs samplers are expected to outperform both linear estimators when the jitter dominates the additive noise (σ z ≫ σ w ) and when the jitter is not so large that the Gibbs sampling methods become inaccurate for the given number of samples. The Gibbs sampler using slice sampling is also expected to be more accurate for large σ z , since rejection sampling fails (the acceptance probability becomes too small). To summarize the performance of the Gibbs samplers, the maximum factor of improvement in σ z attainable by using these sampling algorithms is plotted against increasing M and σ w in Figure 8. These plots portray the same trends as observed for the EM algorithm; however, these algorithms appear to handle the high σ w case better than the EM algorithm does.

VI. CONCLUSION
The results of the previous section are very encouraging from a power-consumption standpoint. A maximum improvement of between 1.4 to 2 times the jitter translates to a two-to-fourfold decrease in power consumption, according to (1). To put the magnitude of such an improvement in context, consider the digital baseband processor for ultra-wideband communication in [30]. This processor incorporates an ADC and a PLL, which consume 86 mW and 45 mW, respectively, out of a 271 mW budget for the chip.
Reducing by a factor of two the power consumed by the ADC alone would decrease the total power consumption of the chip by almost sixteen percent.
While effective, both the EM algorithm and the Gibbs samplers are computationally expensive. One benefit of digital post-processing is that these algorithms can be performed off-chip, on a computer or other system with less limited computational resources. For real-time on-chip applications, Kalman-filterlike versions of the EM algorithm or Gibbs sampler would be more practical; this extension is a topic for further investigation. Related to real-time processing is developing streaming algorithms for the infinite-dimensional case, extending this work for general real-time sampling systems. Another future direction involves modifying these algorithms for correlated or periodic jitter.
Sampling jitter mitigation is actually just one application of these new algorithms. In the frequency domain, jitter maps to uncertainty in frequency; using algorithms such as these should produce more reliable Fourier transforms for systems like spectrum analyzers. In higher dimensions, timing noise becomes location jitter in images or video. Greater tolerance of the locations of pixels in images would allow camera manufacturers to place smaller pixels closer together, enabling higher quality images to be acquired using conventional photosensor technology. This paper shows that significant improvements over the best linear post-processing are possible; thus, further work may impact these and other applications.