# Discrete Simulation of Power Law Noise

Neil Ashby
University of Colorado, Boulder CO
email<EMAIL_ADDRESS>
Affiliate, NIST, Boulder, CO

###### Abstract

A method for simulating power law noise in clocks and oscillators is presented
based on modification of the spectrum of white phase noise, then Fourier
transforming to the time domain. Symmetric real matrices are introduced whose
traces–the sums of their eigenvalues–are equal to the Allan variances, in
overlapping or non-overlapping forms, as well as for the corresponding forms
of the modified Allan variance. Diagonalization of these matrices leads to
expressions for the probability distributions for observing a variance at an
arbitrary value of the sampling or averaging interval $\tau$, and hence for
estimating confidence in the measurements. A number of applications are
presented for the common power-law noises.

##  1 Introduction

The characterization of clock performance by means of average measures such as
Allan variance, Hadamard variance, Theo, and modified forms of such variances,
is widely applied within the time and frequency community as well as by most
clock and oscillator fabricators. Such variances are measured by comparing the
times $t_{k}$ on a device under test, with the times at regular intervals
$k\tau_{0}$ on a perfect reference, or at least on a better reference.
Imperfections in performance of the clock under test are studied by analyzing
noise in the time deviation sequence $x_{k}=t_{k}-k\tau_{0}$, or the
fractional frequency difference during the sampling interval $\tau=s\tau_{0}$:

$\Delta_{k,s}^{(1)}=(x_{k+s}-x_{k})/(s\tau_{0}).$ (1)

The frequency spectrum of fractional frequency differences can usually be
adequately characterized by linear superposition of a small set of types of
power law noise. The frequency spectrum of the fractional frequency
differences of a particular noise type is given by a one-sided spectral
density

$S_{y}(f)=h_{\alpha}f^{\alpha},\quad f>0.$ (2)

(The units of $S_{y}(f)$ are Hz-1.) For the common power-law noise types,
$\alpha$ varies in integral steps from +2 down to -2 corresponding
respectively to white phase modulation, flicker phase modulation, white
frequency modulation, flicker frequency modulation, and random walk of
frequency.

Simulation of clock noise can be extremely useful in testing software
algorithms that use various windowing functions and Fourier transform
algorithms to extract spectral density and stability information from measured
time deviations, and especially in predicting the probability for observing a
particular value of some clock stability variance. This paper develops a
simple simulation method for a time difference sequence that guarantees the
spectral density will have some chosen average power law dependence.
Expressions for the common variances and their modified forms are derived here
that agree with expressions found in the literature, with some exceptions.
This approach also leads to predictions of probabilities for observing a
variance of a particular type at particular values of the sampling time. A
broad class of probability functions naturally arises. These only rarely
correspond to chi-squared distributions.

This paper is organized as follows. Sect. 2 introduces the basic simulation
method, and Sect 3 applies the method to the overlapping Allan variance. Sect.
4 shows how diagonalization of the averaged squared second-difference
operator, applied to the simulated time series, leads to expressions for the
probability of observing a value of the variance for some chosen value of the
sampling or averaging time. Expressions for the mean squared deviation of the
mean of the variance itself are derived in Sect. 6. The approach is used to
discuss the modified Allan variance in Sect. 7, and the non-overlapping form
of the Allan variance is treated in Sect. 8. Appendix 1 discusses evaluation
of a contour integral for the probability, first introduced in Sect. 4, in the
general case.

##  2 Discrete Time Series

We imagine the noise amplitudes at Fourier frequencies $f_{m}$ are generated
by a set of $N$ normally distributed random complex numbers $w_{n}$ having
mean zero and variance $\sigma$, that would by themselves generate a simulated
spectrum for white phase noise. These random numbers are divided by a function
of the frequency, $|f_{m}|^{\lambda}$, producing a spectral density that has
the desired frequency characteristics. For ordinary power law noise, the
exponent $\lambda$ is a multiple of $1/2$, but it could be anything. The
frequency noise is then transformed to the time domain, producing a time
series with the statistical properties of the selected power law noise. The
Allan variance, Modified Allan variance, Hadamard Variance, variances with
dead time, and other quantities of interest can be calculated using either the
frequency noise or the time series.

In the present paper we discuss applications to calculation of various
versions of the Allan variance. Of considerable interest are results for the
probability of observing a value of the Allan variance for particular values
of the sampling time $\tau$ and time series length $N$. The derivations in
this paper are theoretical predictions. A natural frequency cutoff occurs at
$f_{h}=1/(2\tau_{0})$, where $\tau_{0}$ is the time between successive time
deviations. This number is not necessarily related in an obvious way to some
hardware bandwidth. The measurements are assumed to be made at the times
$k\tau_{0}$, and the time errors or residuals relative to the reference clock
are denoted by $x_{k}$. The averaging or sampling time is denoted by
$\tau=s\tau_{0}$, where $s$ is an integer. The total length of time of the
entire measurement series is $T=N\tau_{0}$. The possible frequencies that
occur in the Fourier transform of the time residuals are

$f_{m}=\frac{m}{N\tau_{0}}\,.\quad\quad-\frac{N}{2}+1\leq m\leq\frac{N}{2}\,.$
(3)

Noise Sequences. In order that a set of noise amplitudes in the frequency
domain represent a real series in the time domain, the amplitudes must satisfy
the reality condition

$w_{-m}=(w_{m})^{\ast}.$ (4)

$N$ random numbers are placed in $N/2$ real and $N/2$ imaginary parts of the
positive and negative frequency spectrum. Thus if $w_{m}=u_{m}+iv_{m}$ where
$u_{m}$ and $v_{m}$ are independent uncorrelated random numbers, then
$(w_{m})^{\ast}=u_{m}-iv_{m}$. Since the frequencies $\pm 1/2\tau_{0}$
represent essentially the same contribution, $v_{N/2}$ will not appear. We
shall assume the variance of the noise amplitudes is such that

$\big{<}(w_{m})^{*}w_{n}\big{>}=\big{<}u^{2}+v^{2}\big{>}\delta_{mn}=2\sigma^{2}\delta_{mn};\quad
m\neq 0,N/2.$ (5)

Also, $<w_{m}^{2}>=<u^{2}-v^{2}+2iuv>=0{\rm\ for\ }m\neq N/2.$ The index $m$
runs from $-N/2+1$ to $N/2$. In order to avoid division by zero, we shall
always assume that the Fourier amplitude corresponding to zero frequency
vanishes. This only means that the average of the time residuals in the time
series will be zero, and has no effect on any variance that involves time
differences.

We perform a discrete Fourier transform of the frequency noise and obtain the
amplitude of the $k^{th}$ member of the time series for white PM:

$x_{k}=\frac{\tau_{0}^{2}}{\sqrt{N}}\sum_{m=-N/2+1}^{N/2}e^{-\frac{2\pi
imk}{N}}w_{m}.$ (6)

The factor $\tau_{0}^{2}$ is inserted so that the time series will have the
physical dimensions of time if $w_{m}$ has the dimensions of frequency. We
then multiply each frequency component by $|f_{0}/f_{m}|^{\lambda}$. This will
generate the desired power-law form of the spectral density. The time series
will be represented by

$X_{k}=\frac{\tau_{0}^{2}}{\sqrt{N}}\sum_{m=-N/2+1}^{N/2}\frac{|f_{0}|^{\lambda}}{|f_{m}|^{\lambda}}e^{-\frac{2\pi
imk}{N}}w_{m}.$ (7)

The constant factor $|f_{0}|^{\lambda}$ has been inserted to maintain the
physical units of the time series. The noise level is determined by $f_{0}$.
For this to correspond to commonly used expressions for the one-sided spectral
density, Eq. (2), we shall assume that

$\frac{\tau_{0}^{2}|f_{0}|^{\lambda}}{\sqrt{N}}=\sqrt{\frac{h_{\alpha}}{16\pi^{2}\sigma^{2}(N\tau_{0})}}\,.$
(8)

We shall show that if $2\lambda=2-\alpha$ the correct average spectral density
is obtained. The simulated time series is

$X_{k}=\sqrt{\frac{h_{\alpha}}{16\pi^{2}\sigma^{2}(N\tau_{0})}}\sum_{m}\frac{e^{-\frac{2\pi
imk}{N}}}{|f_{m}|^{\lambda}}w_{m}\,.$ (9)

The average (two-sided) spectral density of the time residuals is obtained
from a single term in Eq. (9):

$s_{x}(f_{m})=\frac{h_{\alpha}}{16\pi^{2}\sigma^{2}(N\tau_{0})f_{m}^{2\lambda}}\bigg{<}\frac{w_{m}w_{m}^{*}}{\Delta
f}\bigg{>}=\frac{h_{\alpha}}{8\pi^{2}f_{m}^{2\lambda}}$ (10)

where $\Delta f=1/(N\tau_{0})$ is the spacing between successive allowed
frequencies. The average (two-sided) spectral density of fractional frequency
fluctuations is given by the well-known relation

$s_{y}(f)=(2\pi f)^{2}s_{x}(f)\,,$ (11)

and the one-sided spectral density is

$S_{y}(f)=\cases{0,&$f<0$;\cr 2s_{y}(f)=h_{\alpha}f^{\alpha},&$f>0$\,,\cr}$
(12)

where $2\lambda=2-\alpha$.

##  3 Overlapping Allan Variance

Consider the second-difference operator defined by

$\Delta_{j,s}^{(2)}=\frac{1}{\sqrt{2\tau^{2}}}(X_{j+2s}-2X_{j+s}+X_{j}).$ (13)

The fully overlapping Allan variance is formed by averaging the square of this
quantity over all possible values of $j$ from 1 to $N-2s$. Thus

$\sigma_{y}^{2}(\tau)=\bigg{<}\frac{1}{N-2s}\sum_{j=1}^{N-2s}(\Delta_{j,s}^{(2)})^{2}\bigg{>}.$
(14)

In terms of the time series, Eq. (9), the second difference can be reduced
using elementary trigonometric identities:

$\displaystyle\Delta_{j,s}^{(2)}=\sqrt{\frac{h_{\alpha}}{32\pi^{2}\tau^{2}\sigma^{2}(N\tau_{0})}}\hbox
to166.2212pt{}$
$\displaystyle\times\sum_{m}\frac{w_{m}}{|f_{m}|^{\lambda}}\bigg{(}e^{-\frac{2\pi
im(j+2s)}{N}}-2e^{-\frac{2\pi im(j+s)}{N}}+e^{-\frac{2\pi im(j)}{N}}\bigg{)}$
$\displaystyle=-\sqrt{\frac{h_{\alpha}}{2\pi^{2}\tau^{2}\sigma^{2}(N\tau_{0})}}\sum_{m}\frac{w_{m}}{|f_{m}|^{\lambda}}e^{-\frac{2\pi
im(j+s)}{N}}\bigg{(}\sin\frac{\pi ms}{N}\bigg{)}^{2}\,.$ (15)

We form the averaged square of $\Delta_{j,s}^{(2)}$ by multiplying the real
quantity times its complex conjugate, then averaging over all possible values
of $j$.

$\displaystyle\sigma_{y}^{2}(\tau)=\frac{h_{\alpha}}{2\pi^{2}\tau^{2}(N\tau_{0})(N-2s)}\sum_{m,n,j}\bigg{<}\frac{w_{m}w_{n}^{*}}{\sigma^{2}}\bigg{>}\frac{\bigg{(}\sin\big{(}\frac{\pi
ms}{N}\big{)}\sin\big{(}\frac{\pi
ns}{N}\big{)}\bigg{)}^{2}}{|f_{m}f_{n}|^{\lambda}}$ $\displaystyle\hbox
to-36.135pt{}\times e^{-\frac{2\pi i(m-n)(j+s)}{N}}\,.$ (16)

The average of the product of random variables only contributes $2\sigma^{2}$
when $m=n$ (see Eq. (5), except when $m=n=N/2$ where
$\big{<}|w_{N/2}|^{2}\big{>}=\sigma^{2}$. The Allan variance reduces to

$\sigma_{y}^{2}(\tau)=\frac{h_{\alpha}}{\pi^{2}\tau^{2}(N\tau_{0})}\bigg{(}\sum_{m}\frac{\bigg{(}\sin\big{(}\frac{\pi
ms}{N}\big{)}\bigg{)}^{4}}{|f_{m}|^{2\lambda}}+\frac{1}{2}\frac{\bigg{(}\sin\big{(}\frac{\pi
s}{2}\big{)}\bigg{)}^{4}}{\big{(}f_{N/2}\big{)}^{2\lambda}}\bigg{)}\,,$ (17)

since every term in the sum over $j$ contributes the same amount. The zero
frequency term is excluded from the sum. For convenience we introduce the
abbreviation

$K=\frac{2h_{\alpha}}{\pi^{2}\tau^{2}(N\tau_{0})}\,.$ (18)

If we sum over positive frequencies only, a factor of 2 comes in except for
the most positive frequency and so

$\sigma_{y}^{2}(\tau)=K\bigg{(}\sum_{m>0}^{N/2-1}\frac{\bigg{(}\sin\big{(}\frac{\pi
ms}{N}\big{)}\bigg{)}^{4}}{f_{m}^{2\lambda}}+\frac{1}{4}\frac{\bigg{(}\sin\frac{\pi
s}{2}\bigg{)}^{4}}{(f_{N/2})^{2\lambda}}\bigg{)}.$ (19)

If the frequencies are spaced densely enough to pass from the sum to an
integral, then $\Delta f=(N\tau_{0})^{-1}$ and

$\frac{1}{N\tau_{0}}\sum_{m}F(|f_{m}|)\rightarrow\int F(|f|)\,df$ (20)

and we obtain the well-known result[6]

$\sigma_{y}^{2}(\tau)=2\int_{0}^{f_{h}}\frac{S_{y}(f)\,df}{\pi^{2}\tau^{2}f^{2}}\bigg{(}\sin(\pi
f\tau)\bigg{)}^{4}\,.$ (21)

Similar arguments lead to known expressions for the non-overlapping version of
the Allan variance as well as for the modified Allan variance. These will be
discussed in later sections.

##  4 Confidence Estimates

In the present section we shall develop expressions for the probability of
observing a particular value $A_{o}$ for the overlapping Allan variance in a
single measurement, or in a single simulation run. $A_{o}$ is a random
variable representing a possible value of the overlapping variance. We use a
subscript “o” to denote the completely overlapping case. To save writing, we
introduce the following abbreviations:

$\displaystyle F_{m}^{j}=\frac{\bigg{(}\sin\bigg{(}\frac{\pi
ms}{N}\bigg{)}\bigg{)}^{2}}{|f_{m}|^{\lambda}}\cos\bigg{(}\frac{2\pi
m(j+s)}{N}{}\bigg{)}$ $\displaystyle
G_{m}^{j}=\frac{\bigg{(}\sin\bigg{(}\frac{\pi
ms}{N}\bigg{)}\bigg{)}^{2}}{|f_{m}|^{\lambda}}\sin\bigg{(}\frac{2\pi
m(j+s)}{N}{}\bigg{)}$ (22)

The dependence on $s$ is suppressed, but is to be understood. We write the
second difference in terms of a sum over positive frequencies only, keeping in
mind that the most positive and the most negative frequencies only contribute
a single term since $\sin(\pi(j+s))=0$. The imaginary contributions cancel,
and we obtain

$\Delta_{j,s}^{(2)}=\sqrt{K}\sum_{m>0}\big{(}F_{m}^{j}\frac{u_{m}}{\sigma}+G_{m}^{j}\frac{v_{m}}{\sigma}\big{)}\,.$
(23)

There is no term in $v_{N/2}$. It is easy to see that the overlapping Allan
variance is given by

$\sigma_{y}^{2}(\tau)=\frac{K}{N-2s}\sum_{j}\sum_{m>0}\bigg{(}(F_{m}^{j})^{2}+(G_{m}^{j})^{2}\bigg{)}\,.$
(24)

To compute the probability that a particular value $A_{o}$ is observed for the
Allan variance, given all the possible values that the random variables
$u_{1},v_{1},...u_{N/2}$ can have, we form the integral

$P(A_{o})=\int\delta\bigg{(}A_{o}-\frac{1}{N-2s}\sum_{j}\big{(}\Delta_{j,s}^{(2)}\big{)}^{2}\bigg{)}\prod_{m>0}\bigg{(}e^{-\frac{u_{m}^{2}+v_{m}^{2}}{2\sigma^{2}}}\frac{du_{m}dv_{m}}{2\pi\sigma^{2}}\bigg{)}\,.$
(25)

The delta function constrains the averaged second difference to the specific
value $A_{o}$ while the random variables
$u_{1},v_{1},...u_{m},v_{m},...u_{N/2}$ range over their (normally
distributed) values. There is no integral for $v_{N/2}$. Inspecting this
probability and the Eq. (23)for the second difference indicates that we can
dispense with the factors of $\sigma^{-1}$ and work with normally distributed
random variables having variance unity. Henceforth we set $\sigma=1$.

The exponent involving the random variables is a quadratic form that can be
written in matrix form by introducing the $N-1$ dimension column vector $U$
(the zero frequency component is excluded)

$U^{T}=[u_{1}\,v_{1}\,...u_{m}\,v_{m},...v_{N/2-1},u_{N/2}]\,.$ (26)

Then

$\frac{1}{2}\sum_{m>0}(u_{m}^{2}+v_{m}^{2})=\frac{1}{2}U^{T}U=\frac{1}{2}U^{T}{\bf
1}U,$ (27)

where ${\bf 1}$ represents the unit matrix. The delta-function in Eq. (25) can
be written in exponential form by introducing one of its well-known
representations, an integral over all angular frequencies $\omega$:[5]

$P(A_{o})=\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}e^{i\omega\big{(}A_{o}-\frac{1}{N-2s}\sum_{j}\big{(}\Delta_{j,s}^{(2)}\big{)}^{2}\big{)}}\prod_{m>0}\bigg{(}e^{-\frac{u_{m}^{2}+v_{m}^{2}}{2\sigma^{2}}}\frac{du_{m}dv_{m}}{2\pi\sigma^{2}}\bigg{)}\,.$
(28)

The contour of integration goes along the real axis in the complex $\omega$
plane.

The squared second difference is a complicated quadratic form in the random
variables $u_{1},v_{1},...u_{m},v_{m},...u_{N/2}$. If this quadratic form
could be diagonalized without materially changing the other quadratic terms in
the exponent, then the integrals could be performed in spite of the imaginary
factor $i$ in the exponent. To accomplish this we introduce a column vector
$C^{j}$ that depends on $j,m,s,N$ and whose transpose is

$\displaystyle(C^{j})^{T}=[F_{1}^{j},G_{1}^{j},...F_{m}^{j},G_{m}^{j},...G_{N/2-1}^{j},F_{N/2}^{j}]\,.$
(29)

Dependence on $s$ is not written explicitly but is understood. The column
vector has $N-1$ real components. It contains all the dependence of the second
difference on frequency and on the particular power law noise. We use indices
$\\{m,n\\}$ as matrix (frequency) indices. The (scalar) second difference
operator can be written very compactly as a matrix product

$\Delta_{j,s}^{(2)}=\sqrt{K}(C^{j})^{T}U=\sqrt{K}U^{T}C^{j}.$ (30)

Then

$\frac{1}{N-2s}\sum_{j}\bigg{(}\Delta_{j,s}^{(2)}\bigg{)}^{2}=U^{T}\bigg{(}\frac{K}{N-2s}\sum_{j}C^{j}(C^{j})^{T}\bigg{)}U.$
(31)

The matrix

$H_{o}=\frac{K}{N-2s}\sum_{j}C^{j}(C^{j})^{T}$ (32)

is real and symmetric. $H_{o}$ is also Hermitian and therefore has real
eigenvalues. A real symmetric matrix can be diagonalized by an orthogonal
transformation,[1, 2] which we denote by $O$. Although we shall not need to
determine this orthogonal transformation explicitly, it could be found by
first finding the eigenvalues $\epsilon$ and eigenvectors $\psi_{\epsilon}$ of
$H_{o}$, by solving the equation

$H_{o}\psi_{\epsilon}=\epsilon\psi_{\epsilon}\,.$ (33)

The transformation $O$ is a matrix of dimension $(N-1)\times(N-1)$ consisting
of the components of the normalized eigenvectors placed in columns. Then

$H_{o}O=OE\,,$ (34)

where $E$ is a diagonal matrix with entries equal to the eigenvalues of the
matrix $H_{o}$. Then since the transpose of an orthogonal matrix is the
inverse of the matrix,

$O^{T}H_{o}O=E\,.$ (35)

The matrix $H_{o}$ is thus diagonalized, at the cost of introducing a linear
transformation of the random variables:

$\frac{K}{N-2s}\sum_{j}\big{(}\Delta_{j,s}^{(2)}\big{)}^{2}=U^{T}H_{o}U=U^{T}OO^{T}H_{o}OO^{T}U=(U^{T}O)E(O^{T}U)\,.$
(36)

We introduce $N-1$ new random variables by means of the transformation:

$V=O^{T}U\,.$ (37)

Then the term in the exponent representing the Gaussian distributions is

$-\frac{1}{2}U^{T}{\bf 1}U=-\frac{1}{2}U^{T}O1O^{T}U=-\frac{1}{2}V^{T}{\bf
1}V=-\frac{1}{2}\sum_{n=1}^{N-1}V_{n}^{2}\,.$ (38)

The Gaussian distributions remain basically unchanged.

Further, the determinant of an orthogonal matrix is $\pm 1$, because the
transpose of the matrix is also the inverse:

$\det(O^{-1}O)=1=\det(O^{T}O)=\big{(}\det(O)\big{)}^{2}.$ (39)

Therefore, changes in the volume element are simple since the volume element
for the new variables is

$\displaystyle dV_{1}dV_{2}...dV_{N-1}=\bigg{|}\det\bigg{(}\frac{\partial
V_{m}}{\partial U_{n}}\bigg{)}\bigg{|}dU_{1}dU_{2}...dU_{N-1}$
$\displaystyle=|\det(O)|dU_{1}dU_{2}...dU_{N-1}$
$\displaystyle=dU_{1}dU_{2}...dU_{N-1}.$ (40)

After completing the diagonalization,

$\frac{1}{N-2s}\sum_{j}\big{(}\Delta_{j,s}^{(2)}\big{)}^{2}=\sum_{i}\epsilon_{i}V_{i}^{2}\,.$
(41)

The probability is therefore

$P(A_{o})=\int\frac{d\omega}{2\pi}e^{i\omega\big{(}A_{o}-\sum_{k}\epsilon_{k}V_{k}^{2}\big{)}}\prod_{i}\bigg{(}e^{-\frac{V_{i}^{2}}{2}}\frac{dV_{i}}{\sqrt{2\pi}}\bigg{)}\,.$
(42)

An eigenvalue of zero will not contribute in any way to this probability since
the random variable corresponding to a zero eigenvalue just integrates out.

Let the eigenvalue $\epsilon_{i}$ have multiplicity $\mu_{i}$, by which is
meant that the eigenvalue $\epsilon_{i}$ is repeated $\mu_{i}$ times.
Integration over the random variables then gives a useful form for the
probability:

$P(A_{o})=\int_{-\infty}^{+\infty}\frac{d\omega}{2\pi}\frac{e^{i\omega
A_{o}}}{\prod_{i}(1+2i\epsilon_{i}\omega)^{\mu_{i}/2}}\,.$ (43)

Finally the contour integral may be deformed and closed in the upper half
complex plane where it encloses the singularities of the integrand. This is
discussed in Appendix 1. Knowing the probability, one may integrate with
respect to the variance to find the cumulative distribution, and then find the
limits corresponding to a 50% probability of observing the variance.

Properties of the eigenvalues. First, it is easily checked that the
probability is correctly normalized by integrating over all $A_{o}$ and using
properties of the delta-function:

$\displaystyle\int
P(A_{o})dA_{o}=\int_{-\infty}^{+\infty}\frac{d\omega}{2\pi}\frac{\int
e^{i\omega
A_{o}}dA_{o}}{\prod_{i}(1+2i\epsilon_{i}\omega)^{\mu_{i}/2}}=\int_{-\infty}^{+\infty}\frac{\delta(\omega)d\omega}{\prod_{i}(1+2i\epsilon_{i}\omega)^{\mu_{i}/2}}$
$\displaystyle=\int_{-\infty}^{+\infty}d\omega\delta(\omega)=1\,.$ (44)

Second, the eigenvalues are all either positive or zero. The eigenvalue
equation for the eigenvector labeled by $\epsilon$ is:

$\frac{K}{N-2s}\sum_{j}C^{j}(C^{j})^{T}\psi_{\epsilon}=\epsilon\psi_{\epsilon}\,.$
(45)

Multiply on the left by $\psi_{\epsilon}^{T}$; assuming the vector has been
normalized, we obtain

$\epsilon=\frac{K}{N-2s}\sum_{j}\bigg{(}(C^{j})^{T}\psi_{\epsilon}\bigg{)}^{2}\geq
0\,.$ (46)

Thus every eigenvalue must be positive or zero.

Next let us calculate the trace of $H_{o}$. Since the trace is not changed by
an orthogonal transformation,

$\displaystyle{\rm Trace}(O^{T}H_{o}O)={\rm Trace}(H_{o}OO^{T})={\rm
Trace}(H_{o}OO^{-1})$ $\displaystyle={\rm
Trace}(H_{o})=\sum_{i}\epsilon_{i}\,.$ (47)

The sum of the diagonal elements of $H_{o}$ equals the sum of the eigenvalues
of $H_{o}$. If we then explicitly evaluate the sum of the diagonal elements of
$H_{o}$ we find

$\displaystyle\sum_{i}\epsilon_{i}=\frac{K}{N-2s}\sum_{j}{\rm
Trace}\big{(}C^{j}(C^{j})^{T})$
$\displaystyle=\frac{K}{N-2s}\sum_{j}\sum_{m>0}\bigg{(}(F_{m}^{j}\big{)}^{2}+\big{(}G_{m}^{j})^{2})\bigg{)}$
$\displaystyle=K\sum_{m>0}\frac{\bigg{(}\sin\frac{\pi
ms}{N}\bigg{)}^{4}}{|f_{m}|^{2-2\alpha}}=\sigma_{y}^{2}(\tau)\,.$ (48)

Every term labeled by $j$ contributes the same amount. We obtain the useful
result that the overlapping Allan variance is equal to the sum of the
eigenvalues of the matrix $H_{o}$. Similar results can be established for many
of the other types of variances.

Distribution of eigenvalues. The eigenvalue equation
$H_{o}\psi_{\mu}=\epsilon\psi_{\mu}$ produces many zero eigenvalues,
especially when $\tau$ is large. The dimension of the matrix $H_{o}$ is
therefore much larger than necessary. Numerical calculation indicates that for
the completely overlapping Allan variance, the eigenvalue equation has a total
of $N-1$ eigenvalues, but only $N-2s$ non-zero eigenvalues; the number of
significant eigenvalues is in fact equal to the number of terms in the sum
over $j$ in the equations:

$\frac{K}{N-2s}\sum_{n}\sum_{j}^{N-2s}(C_{m})^{j}(C_{n})^{j}\psi_{n}=\epsilon\psi_{m}\,.$
(49)

The factorized form of $H_{o}$, that arises on squaring a difference operator,
permits the reduction of the size of the matrix that is to be diagonalized. We
introduce the quantities

$\phi_{\mu}^{j}=\sum_{n}(C_{n})^{j}\psi_{n\mu}\,.$ (50)

We are using the Greek index $\mu$ to label a non-zero eigenvalue and the
index $\nu$ to label a zero eigenvalue. The eigenvalue equation becomes

$\frac{K}{N-2s}\sum_{j}(C_{m})^{j}\phi_{\mu}^{j}=\epsilon\psi_{m\mu}\,.$ (51)

Multiply by $(C_{m})^{l}$ and sum over the frequency index $m$. Then

$\frac{K}{N-2s}\sum_{m,j}(C_{m})^{l}(C_{m})^{j}\phi_{\mu}^{j}=\epsilon\phi_{\mu}^{l}\,.$
(52)

This is an eigenvalue equation with reduced dimension $N-2s$ rather than
$N-1$, since the number of possible values of $j$ is $N-2s$. The eigenvalue
equation can be written in terms of a reduced matrix $H_{red}$, given by

$(H_{red})^{lj}=\frac{K}{N-2s}\sum_{m}(C_{m})^{l}(C_{m})^{j}\,.$ (53)

The indices $l,j$ run from 1 to $N-2s$. Eigenvalues generated by Eq. (52) are
all non-zero. To prove this, multiply Eq. (52) by $\phi_{\mu}^{l}$ and sum
over $l$. We obtain

$\frac{K}{N-2s}\sum_{m}\bigg{(}\sum_{l}(C_{m})^{l}\phi_{\mu}^{l}\bigg{)}^{2}=\epsilon\sum_{l}\big{(}\phi_{\mu}^{l}\big{)}^{2}\,.$
(54)

The eigenvalue cannot be zero unless

$\sum_{l}(C_{m})^{l}\phi_{\mu}^{l}=0$ (55)

for every $m$. The number of such conditions however is larger than the number
$N-2s$ of variables, so the only way this can be satisfied is if
$\phi_{\mu}^{l}=0$, a trivial solution. Therefore to obtain normalizable
eigenvectors from Eq. (52), the corresponding eigenvalues must all be
positive. This is true even through some of these conditions may be trivially
satisfied if the factor $\sin(\pi ms/N)$ vanishes, which happens sometimes
when

$ms=MN$ (56)

where M is an integer. Every time a solution of Eq. (56) occurs, two equations
relating components of $\phi_{\mu}^{l}$ are lost. Suppose there were $n$
solutions to Eq. (56); then the number of conditions lost would be $2n$. The
number of variables is $N-2s$ and the number of conditions left in Eq. (55)
would be $N-1-2n$. The excess of conditions over variables is thus

$N-1-2n-(N-2s)=2(s-n)-1\,.$ (57)

In Appendix 2 we prove that under all circumstances $2(s-n)-1>0$.

We temporarily drop the subscript $o$ since the remainder of the results in
this section are valid for any of the variances. If the eigenvalues are found
and the appropriate matrix is diagonalized, we may compute the probability for
observing a value of the overlapping variance, denoted by the random variable
$A$, by

$\displaystyle
P(A)=\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}e^{i\omega\big{(}A-V^{T}EV\big{)}}\prod_{i}\bigg{(}\frac{e^{-V_{i}^{2}/2}dV_{i}}{\sqrt{2\pi}}\bigg{)}$
$\displaystyle=\int\frac{d\omega}{2\pi}\frac{e^{i\omega
A}}{\prod(1+2i\epsilon_{i}\omega)^{\mu_{i}/2}}\,.$ (58)

Case of a single eigenvalue. If a single eigenvalue occurs once only, the
general probability expression, Eq. (43), has a single factor in the
denominator:

$P(A)=\frac{1}{2\pi}\int_{-\infty}^{\infty}\frac{d\omega e^{i\omega
A}}{\sqrt{1+2i\omega\epsilon}}\,.$ (59)

The integral is performed by closing the contour in the upper half complex
$\omega$ plane. There is a branch point on the imaginary axis at
$\omega=i/(2\epsilon)$, and a branch line from that point to infinity.
Evaluation of the integral gives:

$P(A)=\frac{1}{\sqrt{2\pi\sigma_{y}^{2}(\tau)}}\frac{e^{-A/(2\sigma_{y}^{2}(\tau))}}{\sqrt{A}}$
(60)

This is a chi-squared distribution with exactly one degree of freedom. The
computation of the confidence interval for a given $s$ is simple. The
cumulative probability obtained from Eq. (60) is

$\Phi(A)=\int_{0}^{A}\frac{1}{\sqrt{2\pi\sigma_{y}^{2}(\tau)}}\frac{e^{-x/(2\sigma_{y}^{2}(\tau))}}{\sqrt{x}}dx={\rm
erf}\big{(}\frac{A}{\sqrt{2\sigma_{y}^{2}(\tau)}}\big{)}\,.$ (61)

The $\pm 25$% limits on the probability of observing a value $A$ are then
found to be $1.323\sigma_{y}^{2}$ and $0.1015\sigma_{y}^{2}$, respectively. An
example of this is plotted in Figure 2. ($A$ is a variance, not a deviation.)

Case of two distinct non-zero eigenvalues. For the overlapping variance, when
$s$ has its maximum value $N/2-1$ there are two unequal eigenvalues. The
probability integral can be performed by closing the contour in the upper half
plane and gives the expression

$P(A)=\frac{1}{2\sqrt{\epsilon_{1}\epsilon_{2}}}e^{-\frac{A}{4}\big{(}\frac{1}{\epsilon_{1}}+\frac{1}{\epsilon_{2}}\big{)}}I_{0}\bigg{(}\frac{A}{4}\big{(}\frac{1}{\epsilon_{2}}-\frac{1}{\epsilon_{1}}\big{)}\bigg{)}\,.$
(62)

where $I_{0}$ is the modified Bessel function of order zero. The probability
is correctly normalized. It differs from a chi-squared distribution in that
the density does not have a singularity at $A=0$. This is illustrated in
Figure 1.

Evaluation of the contour integral when there are more than two distinct
eigenvalues is discussed in Appendix 1. If an eigenvalue occurs an even number
$2n$ times, the corresponding singularity becomes a pole of order $n$ and a
chi-squared probability distribution may result; this has only been observed
to occur for white PM.

##  5 Variance of the Variance

One of the goals of this investigation is to estimate the uncertainty in the
variance when some value of the variance is measured. This can be rigorously
defined if the probability $P(A_{o})$ is available, but it may be difficult to
evaluate the integral in Eq. (43). In this case, one might be interested in a
measure such as the rms value of the variance, measured with respect to the
mean variance $\sigma_{y}^{2}(\tau)$. One method for obtaining this measure is
obtained from the probability without introducing a Fourier integral
representation for the delta function:

$\displaystyle\int A_{o}^{2}P(A_{o})dA_{o}=\int
A_{o}^{2}\delta(A_{o}-\sum_{i}\epsilon_{i}V_{i}^{2})\prod_{m}\bigg{(}e^{-V_{m}^{2}/2}\frac{dV_{m}}{\sqrt{2\pi}}\bigg{)}dA_{o}$
$\displaystyle=\int\bigg{(}\sum_{i}\epsilon_{i}V_{i}^{2}\bigg{)}^{2}\delta(A_{o}-\sum_{i}\epsilon_{i}V_{i}^{2})\prod_{m}\bigg{(}e^{-V_{m}^{2}/2}\frac{dV_{m}}{\sqrt{2\pi}}\bigg{)}dA_{o}$
$\displaystyle=\int\bigg{(}\sum_{i}\epsilon_{i}V_{i}^{2}\bigg{)}^{2}\prod_{m}\bigg{(}e^{-V_{m}^{2}/2}\frac{dV_{m}}{\sqrt{2\pi}}\bigg{)}\,.$
(63)

Expanding the square in the integrand, there are $N$ fourth-order terms and
$N(N-1)$ cross terms so we get

$\displaystyle\int
A_{o}^{2}P(A_{o})dA_{o}=\sum_{i}\epsilon_{i}^{2}\big{<}V^{4}\big{>}+2\sum_{i\neq
j}\epsilon_{i}\epsilon_{j}\big{<}V^{2}\big{>}^{2}$
$\displaystyle=3\sum_{i}\epsilon^{2}+2\sum_{i\neq
j}\epsilon_{i}\epsilon_{j}\,.$ (64)

To obtain the variance of the variance, we must subtract

$\big{<}A_{o}\big{>}^{2}=\sum_{i}\epsilon_{i}^{2}+2\sum_{i\neq
j}\epsilon_{i}\epsilon_{j}\,.$ (65)

The result is

$\big{<}A_{o}^{2}\big{>}-\big{<}A_{o}\big{>}^{2}=2\sum_{i}\epsilon_{i}^{2}\,.$
(66)

Therefore the rms deviation of the variance from the mean variance equals
$\sqrt{2}$ times the square root of the sum of the squared eigenvalues. This
result is not very useful if there are only a small number of eigenvalues,
since the factor $\sqrt{2}$ can make the rms deviation from the mean of the
variance larger than the variance. This is because the probability
distribution of the variance is not a normal distribution, and can lead to
unreasonably large estimates of confidence intervals.

If the eigenvalues or the probabilities are not readily available, a similar
confidence estimate can be obtained from an alternative form of Eq. (66) by
considering the trace of the square of $H_{o}$:

$2\sum_{i}\epsilon_{i}^{2}=2\sum_{i,j}\big{(}\sum_{m>0}(F_{m}^{i}F_{m}^{j}+G_{m}^{i}G_{m}^{j})\big{)}^{2}\,.$
(67)

To prove this, consider diagonalization of $H_{o}$ by the orthogonal
transformation $O$, and compute

$\displaystyle 2\sum_{i}\epsilon_{i}^{2}=2{\rm Trace}(E^{2})=2{\rm
Trace}(O^{T}H_{o}OO^{T}H_{o}O)=2{\rm Trace}H_{o}^{2}$
$\displaystyle=\frac{2K^{2}}{(N-2s)^{2}}\sum_{m,n}\bigg{(}\sum_{i}\big{(}(C_{m})^{i}(C_{m})^{i}\big{)}\sum_{j}\big{(}(C_{n})^{j}(C_{n})^{j}\big{)}\bigg{)}$
$\displaystyle=\frac{2K^{2}}{(N-2s)^{2}}\sum_{i,j}\bigg{(}\sum_{m}\big{(}(C_{m})^{i}(C_{m})^{j}\big{)}\sum_{n}\big{(}(C_{n})^{j}(C_{n})^{i}\big{)}\bigg{)}$
$\displaystyle=\frac{2K^{2}}{(N-2s)^{2}}\sum_{i,j}\bigg{(}\sum_{m}\big{(}F_{m}^{i}F_{m}^{j}+G_{m}^{i}G_{m}^{j}\big{)}\bigg{)}^{2}\,.$
(68)

If there are not too many terms in the sums over $i$ and $j$ then the sum over
the frequency index in Eq. (67) becomes

$\sum_{m}\big{(}F_{m}^{i}F_{m}^{j}+G_{m}^{i}G_{m}^{j}\big{)}=\sum_{m>0}\frac{\big{(}\sin\frac{\pi
ms}{N}\big{)}^{4}}{f_{m}^{2\lambda}}\cos\frac{2\pi m(i-j)}{N}\,.$ (69)

This leads to some useful approximation schemes, but these will not be
discussed further in this paper.

Usually the variance of the variance is larger than the range of possible
values of the variance computed from the $\pm 25$% limits obtained from
calculated probabilities. An example is given in Figure 1, where the
confidence intervals are plotted for the non-overlapping variance for
$N=1024$. For $s$ ranging between 342 and 511 there is only one eigenvalue, so
the rms deviation of the variance from the variance is $\sqrt{2}$ times the
variance. This gives an upper limit on the confidence interval that is 1.414
times the variance, larger than that obtained from the actual probability
distribution. The medium heavy lines in Figure 1 show the true $\pm 25$%
probability limits obtained from calculations such as Eq. (10)

Figure 1: Comparisons of average non-overlapping Allan variance (light line)
with rms deviations from the mean of the variance (heavy lines) and true $\pm
25$% limits (medium heavy lines) for $N=1024$ data items in the time series.
Results for the Allan variance for two independent simulation runs are plotted
for comparison.

##  6 Modified Allan Variance

The modified Allan variance is defined so that averages over time are
performed before squaring and averaging. We use a subscript $m$ to distinguish
this form from the overlapping form of the variance. The definition is

$\sigma_{m}^{2}(\tau)=\bigg{<}\bigg{(}\frac{1}{s}\sum_{j=1}^{s}\Delta_{j,s}^{(2)}\bigg{)}^{2}\bigg{>}\,.$
(70)

Summing the expression for the second-difference given in Eq. (3),

$\displaystyle\frac{1}{s}\sum_{j}^{s}\Delta_{j,s}^{(2)}=-\sqrt{\frac{h_{\alpha}}{2s^{2}\pi^{2}\tau^{2}\sigma^{2}(N\tau_{0})}}\sum_{m}\frac{w_{m}}{|f_{m}|^{\lambda}}\big{(}\sin\frac{\pi
ms}{N}\big{)}^{2}$ $\displaystyle\times\bigg{(}\frac{\sin\frac{\pi
ms}{N}}{\sin\frac{\pi m}{N}}\bigg{)}e^{-\frac{\pi im(3s+1)}{N}}\,.$ (71)

Squaring this, we write the complex conjugate of one factor and obtain the
ensemble average

$\
\sigma_{m}^{2}(\tau)=\frac{h_{\alpha}}{2s^{2}\pi^{2}\tau^{2}\sigma^{2}(N\tau_{0})}\sum_{m,n}\frac{\big{<}w_{m}w_{n}^{*}\big{>}}{|f_{m}f_{n}|^{\lambda}}\frac{\bigg{(}\sin\frac{\pi
ms}{N}\bigg{)}^{6}}{\big{(}\sin\frac{\pi m}{N}\big{)}^{2}}\,.$ (72)

Using Eq. (5)and writing the result in terms of a sum over positive
frequencies only,

$\sigma_{m}^{2}(\tau)=2\sum_{m>0}\frac{S_{y}(f_{m})}{\pi^{2}s^{2}\tau^{2}f_{m}^{2}(N\tau_{0})}\frac{\bigg{(}\sin\frac{\pi
ms}{N}\bigg{)}^{6}}{\bigg{(}\sin\frac{\pi m}{N}\bigg{)}^{2}}\,.$ (73)

Passing to an integral when the frequency spacing is sufficiently small,

$\sigma_{m}^{2}(\tau)=2\int_{0}^{f_{h}}\frac{dfS_{y}(f)}{\pi^{2}s^{2}\tau^{2}f^{2}}\frac{\big{(}\sin(\pi
f\tau)\big{)}^{6}}{\big{(}\sin(\pi f\tau_{0})\big{)}^{2}}\,.$ (74)

This agrees with previously derived expressions[12]. The appearance of the
sine function in the denominator of this expression makes the explicit
evaluation of the integral difficult. In Table 1 we give the leading
contributions to the modified Allan variance for very large values of the
averaging time $\tau$. In general there are additional oscillatory
contributions with small amplitudes that are customarily neglected. These
results do not agree however with those published in[7]; the leading terms do
agree with those published in[8].

Noise Type | $S_{y}(f)$ | Mod $\sigma_{y}^{2}(\tau)$ | $\lambda$ | $\alpha$
---|---|---|---|---
White PM | $h_{2}f^{2}$ | $\frac{3h_{2}}{\pi^{2}\tau^{3}}\big{(}1+\frac{5}{18s^{3}}\big{)}$ | 0 | 2
Flicker PM | $h_{1}f$ | $\frac{3h_{1}\ln(256/27)}{8\pi^{2}\tau^{2}}$ | $\frac{1}{2}$ | 1
White FM | $h_{0}$ | $\frac{h_{0}}{4\tau}\big{(}1+\frac{1}{2s^{2}}\big{)}$ | 1 | 0
Flicker FM | $h_{-1}f^{-1}$ | $2h_{-1}\ln\big{(}\frac{3^{27/16}}{4}\big{)}$ | $\frac{3}{2}$ | -1
Random Walk | $h_{-2}f^{-2}$ | $h_{-2}\pi^{2}\tau\big{(}\frac{11}{20}+\frac{1}{12s^{2}}+\frac{1}{40s^{4}}\big{)}$ | 2 | -2
Table 1: Asymptotic expressions for the Modified Allan Variance, in the limit
of large sampling times $\tau=s\tau_{0}$.

Eigenvalue structure for modified Allan variance. For the modified case,
instead of Eq. (32) there are two summations; we use a subscript $m$ on $H$ to
denote the modified form:

$H_{m}=\frac{K}{s^{2}}\sum_{j}C^{j}\sum_{k}(C^{k})^{T}\,.$ (75)

We then seek solutions of the eigenvalue equation

$H_{m}\psi_{\epsilon}=\frac{K}{s^{2}}\sum_{j}C^{j}\sum_{k}(C^{k})^{T}\psi_{\epsilon}=\epsilon\psi_{\epsilon}\,.$
(76)

Multiply by $\sum_{l}(C^{l})^{T}$ and sum over the frequency index.

$\frac{K}{s^{2}}\sum_{l,j,m}\big{(}(C_{m})^{l})^{T}(C_{m})^{j}\big{(}\sum_{k}(C^{k})^{T}\psi_{\epsilon}\big{)}=\epsilon\big{(}\sum_{l}(C^{l})^{T}\psi_{\epsilon}\big{)}\,.$
(77)

The matrix equation has been reduced to a scalar equation for the quantity

$\phi=\sum_{k}(C^{k})^{T}\psi_{\epsilon}\,.$ (78)

Here the matrix product of $(C^{k})^{T}$ with $\psi_{\epsilon}$ entails a sum
over all frequencies. Since we have a scalar eigenvalue equation for
$\epsilon$, there can be one and only one eigenvalue, which is easily seen to
be the same as given by Eq. (73) for each sampling time $\tau$:

$\epsilon=\frac{K}{s^{2}}\sum_{l,k,m}\big{(}(C_{m})^{l}\big{)}^{T}(C_{m})^{k}=\sigma_{m}^{2}(\tau)\,.$
(79)

The probability distribution will be of the form of Eq. (60), a chi-squared
distribution with one degree of freedom.

##  7 Non-overlapping Allan variance

In the non-overlapping form of the Allan variance, the only change with
respect to Eq. (14) is that the sum over $j$ goes from 1 in steps of $s$ up to
$j_{max}\leq N-2s$. We denote the number of values of $j$ by $n_{max}$. The
average non-overlapping variance is the same as that for the overlapping case,
but the probability distributions are different. The matrix $H_{no}$ takes the
form

$(H_{no})_{mn}=\frac{K}{n_{max}}\sum_{j=1,1+2s..}^{j_{max}}(C_{m})^{j}(C_{n})^{j}\,.$
(80)

The subscript “no” labels a non-overlapping variance, and the indices $m$ and
$n$ label frequencies. Here there are only $n_{max}$ terms in the sum since
the values of $j$ skip by $s$. $n_{max}$ is given by:

$n_{max}=\lfloor\frac{N-1}{s}\rfloor-1\,,$ (81)

where $\lfloor Q\rfloor$ denotes the largest integer less than or equal to
$Q$. In order to diagonalize $H_{no}$ and compute probabilities, we look for
eigenvalues by seeking solutions of:

$\sum_{n}(H_{no})_{mn}\psi_{n}=\epsilon\psi_{m}.$ (82)

We define a reduced eigenvector by

$\phi_{j}=\sum_{n}(C_{n})^{j}\psi_{n}.$ (83)

The eigenvalue equation reduces to

$\frac{K}{n_{max}}\sum_{j}(C_{m})^{j}\phi_{j}=\epsilon\psi_{m}.$ (84)

Multiply this equation by $(C_{m})^{l}$ and sum over the frequency labels $m$.
Then

$\frac{K}{n_{max}}\sum_{j}\sum_{m}(\big{(}C_{m})^{l}(C_{m})^{j}\big{)}\phi_{j}=\epsilon\phi_{l}\,.$
(85)

Let

$J_{lj}=\frac{K}{n_{max}}\sum_{m}(C_{m})^{l}(C_{m})^{j}=\frac{K}{n_{max}}(C^{l})^{T}C^{j};$
(86)

the sum over frequencies is accomplished by the matrix multiplication. Just as
for the overlapping case, if we calculate the trace of the matrix $J_{lj}$, we
find that since each term in the sum over $j$ contributes the same amount,

${\rm
Trace}(J_{lj})=\frac{K}{n_{max}}\sum_{j}C_{j}^{T}C_{j}=\sigma_{no}^{2}(\tau)\,.$
(87)

Because the trace remains unchanged under an orthogonal transformation, the
non-overlapping Allan variance will be equal to the sum of the eigenvalues of
Eq. (85). The probability functions for a given $\tau$ can differ from those
for the overlapping or modified cases because the number of eigenvalues and
their multiplicities may be different. The eigenvalues are still all greater
than or equal to zero.

If the eigenvalues are found and the matrix $J_{lj}$ is diagonalized, we may
compute the probability for observing a value of the overlapping variance,
denoted by the random variable $A_{no}$, by

$\displaystyle
P(A_{no})=\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}e^{i\omega\big{(}A_{no}-U^{T}JU\big{)}}\prod_{i}\frac{e^{-U_{i}^{2}/2}}{\sqrt{2\pi}}$
$\displaystyle=\int\frac{d\omega}{2\pi}\frac{e^{i\omega
A_{no}}}{\prod(1+2i\epsilon_{i}\omega)^{\mu_{i}/2}}\,.$ (88)

For example, when there are only two distinct non-zero eigenvalues, The matrix
$J_{lk}$ will have elements

$\displaystyle J_{1,1}=J_{1+s,1+s}=\frac{1}{2}\sigma_{no}^{2};\hbox
to72.26999pt{}$ $\displaystyle
J_{1,1+s}=J_{1+s,1}=\int_{0}^{f_{h}}\frac{S_{y}(f)df}{(\pi\tau
f)^{2}}\big{(}\sin(\pi\tau f)\big{)}^{4}\cos(2\pi\tau f)\,.$ (89)

The eigenvalues are then

$\displaystyle\epsilon_{1}=\frac{1}{2}\sigma_{no}^{2}+|J_{1,1+s}|$
$\displaystyle\epsilon_{2}=\frac{1}{2}\sigma_{no}^{2}-|J_{1,1+s}|\,.$ (90)

The probability is of the form of Eq. (62).

Figure 2: Simulation of Allan variance with $N=64$, $s=19$. For this case
there are two distinct eigenvalues in the matrix for the non-overlapping case.
The variance for $s=19$ was extracted from each of 4000 independent runs and a
histogram of the values obtained was constructed for comparison with the
probability, Eq. (62). Chi-squared distributions with 1, 2, and 3 degrees of
freedom are plotted for comparison.

Figure 1 shows an example of this for flicker FM noise, for $N=64$ items in
the time series. The variance corresponding to $\tau=19\tau_{0}$ was extracted
from each of 4000 independent simulation runs and a histogram of the resulting
values was plotted. Good agreement with the predicted probability distribution
can be seen.

##  8 Other variances

The analysis methods developed in this paper can be extended to other
variances. For example, the Theo variance[13] is defined by an overlapping
average,

$\displaystyle\sigma_{Th}^{2}(s,\tau_{0},N)=\frac{1}{N-s}\sum_{i=1}^{N-m}\frac{4}{3\tau^{2}}\times\hbox
to108.405pt{}$
$\displaystyle\sum_{\delta=0}^{s/2-1}\frac{1}{s/2-\delta}\bigg{(}(X_{i}-X_{i-\delta+s/2})-(X_{i+s}-X_{i+\delta+s/2})\bigg{)}^{2}\,.$
(91)

where $s$ is assumed to be even, and $\tau=3s\tau_{0}/4$. For any of the power
law noises, and after passing to an integral, this expression can be
transformed with the aid of elementary trigonometric identities to

$\sigma_{Th}^{2}(\tau,\tau_{0},N)=\frac{8}{3}\int_{0}^{f_{h}}\frac{S_{y}(f)df}{(\pi\tau
f)^{2}}\sum_{\kappa=1}^{2\tau/3\tau_{0}}\frac{1}{\kappa}\bigg{(}\sin(\pi
f\kappa\tau_{0})\sin\big{(}\pi
f(\frac{4}{3}\tau-\kappa\tau_{0})\big{)}\bigg{)}^{2}\,.$ (92)

The integer $\kappa$ in the denominator makes this variance more difficult to
evaluate, but a factorized quadratic form can still be constructed and the
eigenvalue structure and resulting probabilities can be analyzed.

The Hadamard variance is defined in terms of a third difference, and is widely
used to characterize clock stability when clock drift is a significant
issue[14]. A third difference operator may be defined as

$\Delta_{j,s}^{(3)}=\frac{1}{\sqrt{6\tau^{2}}}\big{(}X_{j+3s}-3X_{j+2s}+X_{j+s}-X_{j}\big{)}\,.$
(93)

The completely overlapping Hadamard variance is the average of the square of
this third difference:

$\displaystyle\sigma_{H}^{2}(\tau)=\frac{1}{N-3s}\sum_{j=1}^{N-3s}\big{(}\Delta_{j,s}^{(3)}\big{)}^{2}$
$\displaystyle\rightarrow\frac{8}{3}\int_{0}^{f_{h}}\frac{S_{y}(f)df}{(\pi\tau
f)^{2}}\big{(}\sin(\pi\tau f)\big{)}^{6}\,.$ (94)

These methods can be applied to cases in which there is dead time between
measurements of average frequency during the sampling intervals. Suppose for
example that the measurements consist of intervals of length $\tau=s\tau_{0}$
during which an average frequency is measured, separated by dead time
intervals of length $D-\tau$ during which no measurements are available. Let
the index $j$ label the measurement intervals with $j=1,2,...N$. An
appropriate variance can be defined in terms of the difference between the
average frequency in the $j^{th}$ interval and that in the interval labeled by
$j+r$:

$\Delta_{j,r,s}^{(2)}=\frac{1}{\sqrt{2}}\big{(}\overline{y}_{j+r,s}-\overline{y}_{j,s}\big{)}\,,$
(95)

where $\overline{y}_{j,s}$ is the average frequency in the interval $j$ of
length $s\tau_{0}$. Then an appropriate variance can be defined as

$\Psi(\tau,D)=\bigg{<}\big{(}\Delta_{j,r,s}^{(2)}\big{)}^{2}\bigg{>}\,.$ (96)

If the measurements are sufficiently densely spaced that it is possible to
pass to an integral, this can be shown to reduce to

$\Psi(\tau,D)=\frac{2}{D\tau}\int_{0}^{f_{h}}df\frac{S_{y}(f)}{(\pi
f)^{2}}\big{(}\sin(\pi frD)\big{)}^{2}\big{(}\sin(\pi f\tau)\big{)}^{2}\,.$
(97)

When $D=\tau$ and $r=1$ there is no real dead time and this variance reduces
to the ordinary Allan variance.

##  9 Summary and Conclusion

In this paper a method of simulating time series for the common power-law
noises has been developed and applied to several variances used to
characterize clock stability. These include overlapping and non-overlapping
forms of the Allan variance, and the Modified Allan variance. Diagonalization
of quadratic forms for the average variances leads to expressions for the
probabilities of observing particular values of the variance for a given
sampling time $\tau=s\tau_{0}$. The probabilities are expressed in terms of
integrals depending on the eigenvalues of matrices formed from squares of the
second differences that are used to define the variance. Generally speaking,
the number of eigenvalues is equal to the number of terms occurring in the sum
used to define averages of the second-difference operator, and this number
gets smaller as the sampling time $\tau$ gets larger. The probability
distribution $P(A)$ for some variance $A$ is useful in estimating the $\pm
25$% confidence interval about the average variance. The eigenvalues are
usually distinct; only for white PM have eigenvalues been observed to occur
with multiplicities other than unity. This must happen in order that chi-
square probability distributions result. Methods for computing the
probabilities have been presented in a few useful cases. Results presented in
this paper are confined to various forms of the Allan variance.

Other methods of simulating power law noise have been published; the present
approach differs from that of[3] in that no causality condition is imposed.
The application of the methods developed here will be applied to other
variances in future papers.

##  10 Appendix 1. Evaluation of contour integrals for probability

In almost all cases except for white PM, the eigenvalues are all distinct. We
show here how then the probability function, Eq. (43), can be reduced to a sum
of real integrals. For each of the eigenvalues, we introduce the quantity

$r_{k}=1/(2\epsilon_{k})\,.$ (98)

The integral becomes

$P(A)=\prod_{k}\big{(}\frac{\sqrt{r_{k}}}{e^{\pi
i/4}}\big{)}\int\frac{d\omega}{2\pi}\frac{e^{i\omega A}}{\prod_{k}(\omega-
ir_{k})^{1/2}}\,.$ (99)

By Jordan’s lemma[11], the contour of integration can be deformed into the
upper half plane. The addition of a circular arcs at radius $|\omega|=\infty$
contributes nothing. Each of the square root factors in Eq. (99) has a branch
point at $\omega=ir_{k}$ with a branch line extending to $+i\infty$ along the
imaginary axis. We define the complex argument of each such factor by

$-\frac{3\pi}{2}<\arg(\omega-ir_{k})<\frac{\pi}{2}\,.$ (100)

Figure 3: Contour deformed to run along branch line on the $y-$axis. A pair of
segments is shown with an odd number of singularities below the segments. The
branch line required by square root factors of the form $\sqrt{\omega-r_{k}}$
is defined by the angles–either $\pi/2$ or $-3\pi/2$–of the directed segments
relative to the Re$(\omega)$ axis. Contributions from pairs of segments that
are above an even number of branch points cancel out.

All branch lines extend along the positive $y-$axis to infinity. The largest
eigenvalues give singularities closest to the real axis. Fig. 3 illustrates
the resulting contour. Around each branch point is a circular portion which
contributes nothing because as the radius $\delta$ of each circle approaches
zero, the contribution to the integral approaches zero as $\sqrt{\delta}$. The
integral then consists of straight segments where $\omega=iy\pm\delta$, where
$\delta$ approaches zero. Two such straight segments are illustrated in Fig.
3, one on each side of the branch line along the $y-$axis. Suppose the
interval of interest is placed so that out of a total of $M$ eigenvalues, $n$
of them are below the interval (three are shown in the figure). For the
contribution to the integral on the left, $y$ is decreasing and we can account
for this by reversing the limits and introducing a minus sign. The phase
factor contributed by $n$ factors with branch points below the interval is

$\frac{1}{\big{(}e^{-\frac{3\pi i}{4}}\big{)}^{n}}=\bigg{(}e^{\frac{3\pi
i}{4}}\bigg{)}^{n}\,.$ (101)

The contribution from branch points above the interval of interest is the same
on both sides of the $y-$axis, and is

$\bigg{(}e^{\frac{i\pi}{4}}\bigg{)}^{M-n}\,.$ (102)

The factor in front of the integral sign in Eq. (99) is the same on both sides
of the $y-$axis and is

$\bigg{(}e^{-\frac{i\pi}{4}}\bigg{)}^{M}\,.$ (103)

The total phase factor of this contribution, including a factor $i$ that comes
from setting $\omega=iy$ is thus

$-i\bigg{(}e^{-\frac{i\pi}{4}}\bigg{)}^{M}\bigg{(}e^{\frac{i\pi}{4}}\bigg{)}^{M-n}\bigg{(}e^{\frac{3\pi
i}{4}}\bigg{)}^{n}=-i\bigg{(}e^{\frac{\pi i}{2}}\bigg{)}^{n}\,.$ (104)

For the integral along the segment on the positive side of the $y-$axis, the
only difference is that the phase of each of the $n$ contributions from branch
points below the interval changes from $3\pi/4$ to $-\pi/4$. The phase factor
for this part of the contour is thus

$+i\bigg{(}e^{-\frac{i\pi}{4}}\bigg{)}^{M}\bigg{(}e^{\frac{i\pi}{4}}\bigg{)}^{M-n}\bigg{(}e^{\frac{-\pi
i}{4}}\bigg{)}^{n}=+i\bigg{(}e^{\frac{-\pi i}{2}}\bigg{)}^{n}\,.$ (105)

If $n$ is even, the two contributions cancel. If $n=2m+1$ is odd, then the
contributions add up with a factor $2(-1)^{m}$. The probability is thus always
real and consists of contributions with alternating signs, with every other
interval left out.

In summary, the contour integral contributions from portions of the imaginary
axis in the complex $\omega$ plane that have an even number of branch points
below the interval will not contribute to the integral. For example, if there
are four distinct eigenvalues the probability will reduce to

$\displaystyle
P(A)=\frac{\sqrt{r_{1}r_{2}r_{3}r_{4}}}{\pi}\int_{r_{1}}^{r_{2}}\frac{e^{-yA}dy}{\sqrt{(y-r_{1})(r_{2}-y)(r_{3}-y)(r_{4}-y)}}$
$\displaystyle-\frac{\sqrt{r_{1}r_{2}r_{3}r_{4}}}{\pi}\int_{r_{3}}^{r_{4}}\frac{e^{-yA}dy}{\sqrt{(y-r_{1})(y-r_{2})(y-r_{3})(r_{4}-y)}}\,.$
(106)

Such results have been used to evaluate the probabilities for certain sampling
intervals for flicker fm noise in Sect. 4.

##  11 Appendix 2. Proof that Eq. (52) generates positive eigenvalues

In this Appendix we show that under all circumstances Eq. (56) gives rise to
positive eigenvalues.

Obviously if $N$ is a prime number Eq. (56) can never be satisfied ($n=0$).
Consider the case $s=N-1$. Then $m=MN$ and there are no solutions since $m\leq
N/2$. In general both $m$ and $s$ may contain factors that divide $M$ or $N$.
Suppose $s=ab$ where at least one of the factors $a,b$ is greater than 1,
where $a$ divides $M$ and $b$ divides $N$. Then let

$M=am_{a},\quad N=bn_{b}\,.$ (107)

Then $m=m_{a}n_{b}$. $m=m_{a}n_{b}>N/2=bn_{b}/2$ cannot happen because $m\leq
N/2$. If $m_{a}n_{b}=N/2=bn_{b}$, then the number of solutions is $n=1$ and
$2(s-n)-1=2(ab-1)-1\geq 2(2-1)-1>0$. If $m=m_{a}n_{b}<N/2=bn_{b}/2$, then
there may be solutions $m=m_{a}n_{b},m=2m_{a}n_{b},m=3m_{a}n_{b},...$up to
$bn_{b}/2$. The number of such solutions is

$n=\lfloor\frac{bn_{b}}{m_{a}n_{b}}\rfloor=\lfloor\frac{b}{m_{a}}\rfloor$
(108)

where $\lfloor x\rfloor$ means the largest integer less than or equal to $x$.
The excess of conditions is then

$2(s-n)-1=2(ab-\lfloor\frac{b}{m_{a}}\rfloor)>1\,,$ (109)

which proves the assertion.

## References

  * [1] O. Strang, _Linear Algebra and its Applications_ , $2^{nd}$ ed., Academic Press (New York) 1980, Sect. 5.5.
  * [2] R. Stoll,_Linear Algebra and Matrix Theory_ , McGraw-Hill (New York) 1980.
  * [3] N. J. Kasdin and T. Walter, _Discrete Simulation of Power Law Noise_ , Proc. 1992 IEEE Frequency Control Symposium, pp 274-283.
  * [4] _Characterization of Frequency and Phase Noise_ , Report 580 of the International Radio Consultative Committee (C.C.I.R.), 1986, pp 142-150, Reprinted in NIST Technical Note 1337, pp TN-162-170. (see [12]
  * [5] M. J. Lighthill, _Introduction to Fourier Analysis and Generalized Functions_ , Cambridge University Press (1958).
  * [6] J. A. Barnes, A. R. Chi, L. S. Cutler, D. J. Healey, D. B. Leeson, T. E. McGunigal, J. A. Mullen, Jr., W. L. Smith, R. L. Sydnor, R. F. C. Vessot, and G. M. R. Winkler, _Characterization of Frequency Stability_ , IEEE Trans. Inst. Meas. 1M-20(2), May 1971; Sect. III. discusses the one-sided spectral density.
  * [7] P. Lesage and T. Ayi, _Characterization of Frequency Stability: analysis of the Modified Allan Variance and Properties of Its Estimate_ , IEEE Trans. Inst. Meas. Vol. IM-33 (4), i pp 332-336, December 1984; (See TN-1337, pp TN-259-TN263).
  * [8] D. A. Howe and F. Vernotte, _Generalization of the Total variance approach to the modified Allan variance_ , Proc. 1999 PTTI Mtg., pp. 267-276.
  * [9] S. Bregni and L. Jmoda, _Improved Estimation of the Hurst Parameter of Long-Range Dependent Traffic Using the Modified Hadamard Variance_ , Proceedings of the IEEE ICC, June 2006.
  * [10] M. Abramowitz and I. Stegun, _Handbook of Mathematical Functions_ , NBS
  * [11] V. I. Smirnov,_Complex Variables-Special Functions_ , III.2, “A Course of Higher Mathematics,” Pergamon Press 1964 Sect. 60, pp 232-233.
  * [12] D. B. Sullivan, D. W. Allan, D. A. Howe and F. L. Walls, eds., _Characterization of Clocks and Oscillators_ NIST Technical Note 1337 March 1990 U. S. Government Printing Office, Washington, D.C. p. TN-9
  * [13] D. A. Howe and T. Peppler, “Very long-term frequency stability: Estimation using a special-purpose statistic,” in _Proc. 2003 Joint Meeting IEEE Int. Frequency Control Symp. and European Frequency and Time Forum Conf.,_ April 2004, pp. 233-238.
  * [14] S. A. Hutsell, J. A. Buisson, J. D. Crum, H. S. Mobbs and W. G. Reid, _Proc. 27 th Annual PTTI Systems and Applications Meeting_, San Diego, CA. Nov. 29-Dec. 1 1995 pp. 291-301