1. No category

The Parabolic variance (PVAR), a wavelet
variance based on least-square fit
F. Vernotte, M. Lenczner, P.-Y. Bourgeois and E. Rubiola
arXiv:1506.00687v1 [physics.data-an] 1 Jun 2015
Final version – May 31, 2015
IEEE Transact. UFFC Special Issue to celebrate the
50th anniversary of the Allan Variance
Abstract
The Allan variance (AVAR) is one option among the wavelet variances. However a milestone in the analysis
of frequency fluctuations and in the long-term stability of clocks, and certainly the most widely used one, AVAR
is not suitable when fast noise processes show up, chiefly because of the poor rejection of white phase noise. The
modified Allan variance (MVAR) features high resolution in the presence of white PM noise, but it is poorer for slow
phenomena because the wavelet spans over 50% longer time.
This article introduces the Parabolic Variance (PVAR), a wavelet variance similar to the Allan variance, based
on the Linear Regression (LR) of phase data. The PVAR relates to the Ω frequency counter, which is the topics of a
companion article [the reference to the article, or to the ArXiv manuscript, will be provided later during the revision
process]. The PVAR wavelet spans over 2τ , the same of the AVAR wavelet.
After setting the theoretical framework, we analyze the degrees of freedom and the detection of weak noise
processes in a realistic case, where (at least) two types of noise are present simultaneously. This is for example
the case of the detection of white FM noise of an atomic standard in the presence of the white PM noise of the
instrument.
Comparing the three wavelet variances, PVAR features almost as well as AVAR in the detection of slow
phenomena, like random walk and drift. By contrast, PVAR exhibits the best ability to divide between fast noise
phenomena, from white PM to flicker FM, and it is in all cases superior to MVAR.
I. I NTRODUCTION
The Allan variance (AVAR) was the first of the wavelet-like variances used for the characterization of oscillators
and frequency standards [1], and despite 50 years of progress, is still the best at rendering the largest τ for a
given time series of measured data. This feature highly desired for monitoring the frequency standards used for
timekeeping. Unfortunately, this feature comes with poor resolution of white and flicker phase noise.
Michel Lenczner, Pierre-Yves Bourgeois and Enrico Rubiola are with the CNRS FEMTO-ST Institute, Department of Time and Frequency,
UBFC, UFC, UTBM, ENSMM, 26, chemin de l’´epitaphe, Besancon, France. E-mail {michel.lenczner|pyb2|rubiola}@femto-st.fr. Enrico’s home
page http://rubiola.org.
F. Vernotte is with UTINAM, Observatory THETA of Franche-Comt´e, University of Franche-Comt´e/UBFC/CNRS, 41 bis avenue de
l’observatoire - B.P. 1615, 25010 Besanc¸on Cedex - France (Email: francois.vernotte@obs-besancon.fr).
1
Conversely, high resolution in the presence of white and flicker phase noise is mandatory for the measurement
of short-term fluctuations (µs to s), and useful for medium-term fluctuations (up to days). This is the case of optics
and of the generation of pure microwaves from optics. The same features are of paramount importance for radars,
VLBI, geodesy, space communications, etc. As a fringe benefits, extending the time-domain measurements to lower
τ is useful to check on the consistency between variances and phase noise spectra. The modified Allan variance
[2]–[4] (MVAR) is suitable to the analysis of fast fluctuations, at a moderate cost in terms of computing power.
Frequency counters specialized for MVAR are available as a niche product, chiefly intended for research labs [5].
A sampling rate of 1/τ is sufficient for the measurement of AVAR, while a rate of 1/τ0 = m/τ is needed for
MVAR, where the rejection of white phase noise is proportional to m. MVAR is based on the simple averaging of
m almost-overlapped frequency data, before evaluating σ 2 (τ ).
The linear regression provides the lowest-energy (or lowest-power) fit of a data set, which is in most cases
considered the optimum approximation. For our purposes, the least-square fit finds an obvious application in the
estimation of the frequency from a time series of phase data, and opens the way to improvements in fluctuation
analysis. Besides, new digital hardware — like FPGAs and SoCs — provides bandwidth and computing power at
an acceptable complexity, and makes possible least-square fitting in real-time.
We apply least-square estimation of frequency to fast time stamping. The simplest estimator if this family is
the linear regression (LR) on phase data. The LR can be interpreted as a weight function applied to the measured
frequency fluctuations. The shape of such weight function is parabolic. We call the corresponding instrument1 ‘Ω
counter,’ for the graphical analogy of the parabola with the Greek letter, and in the continuity of the Π and Λ
counters [6], [7]. The Ω estimator is similar to the Λ estimator, but exhibits higher rejection of the instrument noise,
chiefly of white phase noise. This is important in the measurement of fast phenomena, where the cutoff frequency
fH is necessarily high, and the white phase noise is integrated over the wide analog bandwidth that follows.
In the same way as the Π and Λ estimators yield the AVAR and the MVAR, respectively, we define a variance
based on the Ω estimator. Like in the AVAR and MVAR, the weight functions are similar to wavelets, but for
the trivial difference that they are normalized for power-type signals. A similar use of the LR was proposed
independently by Benkler et al. [8] at the IFCS, where we gave our first presentation on the Ω counter and on our
standpoint about the PVAR. In a private discussion2 , we agreed on the name PVAR (Parabolic VARiance) for this
variance, superseding earlier terms.
After setting the theoretical framework of the PVAR, we provide the response to to noise described by the usual
polynomial spectrum. Then we calculate the degrees of freedom and confidence intervals, checking on the analytical
results against extensive simulations. Finally, we compare the performance of AVAR, MVAR and PVAR for the
detection of noise types, using the value of τ where σ 2 (τ ) changes law as an indicator. In most practical cases
PVAR turns out to be the fastest, to the extent that it enables such detection with the shortest data record.
1 Submitted
2 Private
for the IEEE T UFFC Special Issue on the 2015 IFCS.
discussion between E. Benkler, U. Sterr, F. Vernotte and E. Rubiola, occurred during the International Frequency Control Symposium,
Denver, CO, April 16, 2015.
2
II. S TATEMENT OF THE PROBLEM
The clock signal is usually written as
v(t) = V0 sin[2πν0 t + ϕ(t)]
where V0 is the amplitude, ν0 is the nominal frequency, and ϕ(t) is the random phase fluctuation. Notice that ϕ(t)
is allowed to exceed ±π. Alternatively, randomness is ascribed to the frequency fluctuation (∆ν)(t) = 2π ϕ(t).
˙
We introduce the normalized quantities
x(t) = t + x(t)
y(t) = 1 + y(t)
˙
where x(t) = ϕ(t)/2πν0 , and y(t) = x(t).
The quantity x(t) is the clock readout, which is equal to the time t plus
the random fluctuation x(t). Accordingly, the clock signal reads
v(t) = V0 sin[2πν0 x(t)]
= V0 sin[2πν0 t + 2πν0 x(t)]
For the layman, x is the time displayed by a watch, t is the ‘exact’ time from a radio broadcast, and x the watch
error. The error x is positive (negative) when the watch leads (lags). Similarly, y is the normalized frequency of
the watch’s quartz, and y its fractional error. For example, if y = +10 ppm (constant), the watch leads uniformly
by 1.15 s/day. For the scientist, x(t) is the random time fluctuation, often referred to as ‘phase time’ (fluctuation),
and y(t) is the random fractional-frequency fluctuation. The quantities x(t) and y(t) match exactly x(t) and y(t)
used in the general literature of time and frequency [citations].
The seed of this article is explained in Fig. 1. We use the linear regression of phase data to get a sequence {ˆ
y}
of data averaged on contiguous time intervals of duration τ , and in turn the sequence {ˆy} of fractional-frequency
fluctuation data. Two contiguous elements of {ˆ
y} and {ˆy} are shown in Fig. 1, from which we get one value of
1
2 (y2
− y1 )2 for the estimation of AVAR.
Most of the concepts underneath are expressed in both the continuous and the discrete settings with common
notations without risk of confusion. For example, the same expression x = t + x maps into xk = tk + xk in the
discrete case, and into x(t) = t + x(t) in the continuous case. The notations h . i, (., .) and || . || represent the
P
P
average, the scalar product and the norm. They are defined as hxi = n1 k xk , (x, y) = k xk yk , ||x|| = (x, x)1/2
R
R
where n is the number of terms of the sum in the discrete case, and as hxi = T1 x(t) dt, (x, y) = x(t)y(t) dt,
||x|| = (x, x)1/2 where T is the length of the interval of integration in the continuous case. The span of the sum
and the integral will be made precise in each case of application. The mathematical expectation and the variance
of random variables are denoted by E{ . } and V{ . }.
The linear regression problem consists in searching the optimum value yˆ of the slope η (dummy variable used
to avoid confusion) that minimizes the norm of the error x − ηt, i.e., yˆ = arg minη ||x − ηt||2 . Its solution is the
3
Fig. 1.
Principle of the measurements and notations
random variable
yˆ =
(x − hxi , t − hti)
.
||t − hti ||2
We recall some useful properties of yˆ as an estimator of the slope of x. Formulation is easier, with no loss of
generality, for the case where the time sequence is centered at zero, i.e., hti = 0. Proofs are omitted.
1) The estimator yˆ can be simplified as
yˆ =
(x, t)
.
||t||2
2) If in addition the component xk (or the values x(t)) are independent then the estimator variance
V{ˆ
y} =
σx2
.
||t||2
3) With a constant time step τ0 , the discrete time is tk = (k + 21 )τ0 for k ∈ {−p, ..., p}, m = 2p and τ = mτ0 .
For large m, we get
yˆ ≈ 1 +
12 (x, t)
mτ 2
and
V{ˆ
y} ≈
12σx2
mτ 2
4
4) With a signal that is continuous over a symmetric time interval (− τ2 , τ2 ), we get
yˆ = 1 +
12 (x, t)
.
τ3
(1)
The continuous form of the estimator yˆ can be expressed in the form of a weight average of x or y. To this end,
it is considered as a time dependent function defined over t ∈ (0, τ ),
Z
12 τ /2
yˆ(t) =
s x(t − τ /2 + s) ds
τ 3 −τ /2
Z
τ
12 0
(s + )x(t + s) ds
=
τ 3 −τ
2
Z τ
τ
12
( − s)x(t − s) ds.
=
τ3 0 2
(2)
III. T HE PARABOLIC VARIANCE : PVAR
A. Formulation in the time domain
Let us define T the duration of the data run, τ0 the sampling step, N the number of samples so T = N τ0 , and
n the ratio T /τ ; therefore N = m × n. We consider then, the series {ˆyi }i=1,..,n of frequency deviation estimates.
1) Mathematical expectation and estimators of PVAR:
a) Variances and wavelets: Discrete formulation for a generic variance: An unbiased estimator of the variance
2
Pn Pn
1
1
2
ˆ
ˆ
V {ˆy} is Sn−1
= n−1
y
−
y
, so
i
j
i=1
j=1
n
2 V {ˆy} = E Sn−1
.
2
Following the example of Allan [1], we replace the estimator Sn−1
by a two-sample variance by setting n = 2,
2
V {ˆy} = E S1 that we denote
o
1 n
2
(3)
σ 2 (τ ) = E (ˆy1 − ˆy2 ) .
2
and we define the estimator of σ 2 (τ )
E
1D
2
(ˆyi − ˆyi+1 )
2
σ
b2 (τ ) =
(4)
where the average is over the n − 1 differences ˆyi − ˆyi+1 .
Continuous formulation of PVAR: In the continuous setting, remarking
Z τ
Z τ
12
τ
τ
ˆy(t + τ ) − ˆy(t) =
(
(
−
s)x(t
+
τ
−
s)
ds
−
( − s)x(t − s) ds)
τ3 0 2
2
0
Z τ
12
τ
=
(|s| − )x(t − s) ds
τ 3 −τ
2
Z τ
12
τ
=
(|t − s| − )x(s) ds.
τ 3 −τ
2
the two-sample variance (3) is rewritten as σP2 (τ ) =
1
y(t
2 E{(ˆ
+ τ ) − ˆy(t))2 } independent of t can be expressed
under the form of a weighted running average
σP2 (τ )
(Z
+τ
=E
−τ
2 )
x(s)wx (s − t) ds
(5)
5
A
B
Fig. 2.
C
The wavelets associated to the different counters and variances: (A) Π counter and AVAR, (B) Λ counter and MVAR, (C) Ω counter
and PVAR.
with the even weight function wx (t) =
√
6 2
τ
τ 3 (|t| − 2 )χ(−τ,τ ) (t)
where χ(−τ,τ ) (t) is the characteristic function equals
to 1 for t ∈ (−τ, τ ) and equals to 0 otherwhere. Similarly, the estimator (4) is written
Z
1 T
σ
ˆP2 (τ ) =
(y(t) ∗ hy (t))2 dt.
T 0
(6)
√
Since y = x˙ , introducing the odd weight function wy (t) = − 3 τ 32t (|t| − τ ), σP2 (τ ) can also be expressed as a
weighted running average in y
σP2 (τ )
(Z
=E
2 )
y(s)wy (s − t) ds
.
R
For the purpose of operations with the Fourier transform, it is convenient to restate these expression in terms of
filters or convolution forms
σP2 (τ ) = E{(y ∗ hy )2 } = E{(x ∗ hx )2 }
(7)
with the even and odd kernels
hx (s)
=
hy (s)
=
√
6 2
τ
(|s| − )χ(−τ,τ ) (s)
3
τ
2
√
3 2s
(|s| − τ ) χ(−τ,τ ) (s).
τ3
(8)
(9)
Remark 1: The above formulations are identically stated for ˆy representing a more general estimator of the
frequency deviation obtained by either Π, Λ ou Ω counter. Depending on the type of counter which is used, i.e. Π,
Λ or Ω counter, (3) will involve a particular wavelet and thus a particular variance, i.e. AVAR, MVAR and PVAR.
The wavelets corresponding to these counters and variances are represented in Figure 2.
b) The approach of Lesage and Audoin: Following the approach of Lesage and Audoin [9], let us define the
√
variance elementary estimates αi = (ˆyi − ˆyi+1 )/ 2 as
σ
ˆ (τ ) =
M
1 X
2
(αi ) ,
M i=1
(10)
where the connection between αi and the N individual x(kτ0 ) time error measurements depends on the counter
which has been used for obtaining the ˆy estimates.
2) PVAR from phase data:
6
a) Expression of the variance elementary estimates: Let us denote xi = x(iτ0 ). The elementary estimates of
the variances are
AVAR:
1
αi = √ (−xi + 2xi+m − xi+2m )
2τ
(11)
and M = N − 2m.
MVAR:
αi = √
m−1
X
1
(−xi+k + 2xi+m+k − xi+2m+k )
2mτ k=0
(12)
and M = N − 3m.
PVAR: The discrete form of ˆy may be obtained from (2) by replacing the time integral by a summation with a τ0
step, τ by mτ0 , s by kτ0 , t by iτ0 , and x(t) by xi
m−1 m−1 12 X (m − 1)τ0
12 X m − 1
ˆyi = 3 3
− kτ0 xi−k τ0 = 2
− k xi−k
m τ0
2
m τ
2
k=0
k=0
Similarly,
ˆyi+1
m−1 12 X m − 1
− k xi+m−k
= 2
m τ
2
k=0
and
ˆyi − ˆyi+1
m−1 12 X m − 1
= 2
− k (xi−k − xi+m−k ) .
m τ
2
k=0
Since xi is defined for i ≥ 0, we ought to shift the {xi } sequence in such a way that ˆyi is defined even with i = 0.
We can then shift the {xi } sequence by m − 1
ˆyi − ˆyi+1
m−1 12 X m − 1
= 2
− k (xi+m−1−k − xi+2m−1−k ) .
m τ
2
k=0
Since the coefficient (m − 1)/2 − k is symmetrical for k and for m − 1 − k, we may exchange i + m − 1 − k by
i + m − 1 − (m − 1 − k) = i + k and i + 2m − 1 − k by i + 2m − 1 − (m − 1 − k) = i + m + k
m−1 12 X m − 1
ˆyi − ˆyi+1 = 2
− k (xi+k − xi+m+k ) .
m τ
2
k=0
Therefore, it comes
√ m−1 6 2 X m−1
− k (xi+k − xi+m+k )
αi = 2
m τ
2
(13)
k=0
and M = N − 2m.
b) Estimation of PVAR from phase data: Since we have N time error samples {xi } with a sampling period
τ , from (10) and (13), the estimate σ
ˆP2 (τ ) may be computed as
"
#2
M
−1 m−1
X
X m − 1
72
2
σ
ˆP (τ ) =
− k (xi+k − xi+m+k ) .
M m4 τ 2 i=0
2
k=0
Figure 3 (above) shows the shape of the phase weighting function of PVAR.
(14)
7
3) PVAR from frequency data: From (9) and (7), it comes
σ
ˆP2 (τ ) =
M −1
1 X
2
[y(ti ) ∗ hy (ti )] .
M i=0
(15)
Figure 3 (below) shows the shape of the frequency weighting function of PVAR.
B. Formulation in the frequency domain
1) Transfer function: The transfer function HP y (f ) of PVAR is the Fourier transform of its kernel hy (t). The
square of its modulus is given by
2
9 2 sin2 (πτ f ) − πτ f sin(2πτ f )
|HP y (f )| =
.
2(πτ f )6
2
(16)
2
Figure 4 shows that |HP y (f )| is an octave bandpass filter similar to that of AVAR, but with significantly smaller
side lobes.
2) Proof: The transfer function corresponding to hy is deduced by application of the Fourier transform. Since,
hy is an odd function,
Z
HP y (f ) =
−2iπf t
hy (t)e
R
But the primitive
Z
t (t − τ ) e−2iπf t dt
τ
t (t − τ ) e−2iπf t dt.
0
1
e−2iπtf −2iπ 2 t2 f 2 + 2iπ 2 τ tf 2 − 2πtf + πτ f + i
3
3
4π f
1
πτ f − ie−2iπτ f + πτ f e−2iπτ f + i
4π 3 f 3
τ
Z
Z
= −
=
so the integral
12i
Im
dt = √
2τ 3
t (t − τ ) e−2iπf t dt =
0
1
πτ f − ie−2iπτ f + πτ f e−2iπτ f + i .
3
3
4π f
Then,
HP y (f )
=
=
3i
Im πτ f − ie−2iπτ f + πτ f e−2iπτ f + i
3
3
3
2π τ f
3i
√
(1 − cos (2πτ f ) − πτ f sin (2πτ f )) .
2π 3 τ 3 f 3
√
Moreover 1 − cos (2πτ f ) = 1 − cos2 (πτ f ) + sin2 (πτ f ) = sin2 (πτ f ) + sin2 (πτ f ) = 2 sin2 (πτ f ), and we
conclude to
HP y (f ) = √
3i
2 sin2 (πτ f ) − πτ f sin (2πτ f ) ,
3
3
3
2π τ f
and to
|HP y (f )|2 =
9
2
(1 − cos (2πτ f ) − πτ f sin (2πτ f )) .
2π 6 τ 6 f 6
8
Fig. 3.
Time domain computation weight of PVAR from phase data (above) or frequency deviations (below) for τ = 8τ0 .
9
Fig. 4.
Transfer function of PVAR for 3 τ -values.
3) Convergence properties:
Remark 2: For small f ,
4
1
sin (πτ f ) ≈ πτ f − π 3 τ 3 f 3 + f 4 O(f ) and sin (2πτ f ) ≈ 2πτ f − π 3 τ 3 f 3 + f 4 O(f )
6
3
so
2
1 3 3 3
4 3 3 3
2 sin (πτ f ) − πτ f sin (2πτ f ) ≈ 2 πτ f − π τ f
− πτ f 2πτ f − π τ f
+ f 5 O(f )
6
3
1 4 4 4 2 2 2
≈
π τ f π τ f + 12 + f 5 O(f )
18
2
then
HP y (f ) ≈
We conclude that the module



|HP y (f )|2





 |H (f )|2
Py
≈
√
2iπτ f at low frequency.
2π 2 τ 2 f 2 at low frequency
and
decreases as (πτ f )−4 at high frequency.
10
Sy (f )
PVAR
h−2 f −2
26π 2 h−2 τ
35
h−1 f −1
2 [7 − ln(16)] h−1
5
h0 f 0
3h0
5τ
h+1 f +1
3 [ln(16) − 1] h+1
2π 2 τ 2
h+2 f +2
3h+2
2π 2 τ 3
y(t) = C1 · t
1 2 2
C τ
2 1
TABLE I
R ESPONSES OF PVAR FOR THE DIFFERENT TYPES OF NOISE AND FOR A LINEAR FREQUENCY DRIFT
Therefore
2
•
for small f : |HP y (f )| ≈ 2(πτ f )2 and then converges for f −2 FM
•
for large f : |HP y (f )| decreases as f −4 and then converges for f +2 FM.
2
C. Responses of σP2 (τ )
1) Calculation in the time domain: From (14) it comes
( M −1 "m−1 #
X X m−1
72
1
− k (xi+k − xi+m+k )
σP2 (τ ) = E σ
ˆP2 (τ ) = 4 2 E
m τ
M i=0
2
k=0
#)
"m−1 X m−1
− l (xi+l − xi+m+l )
×
2
l=0
m−1 m−1 72 X X m − 1
m−1
= 4 2
−k
−l
m τ
2
2
(17)
k=0 l=0
× [Rx ((k − l)τ0 ) − Rx ((k − m − l)τ0 ) − Rx ((m + k − l)τ0 ) + Rx ((k − l)τ0 )]
where Rx (θ) = E {x(t)x(t + θ)} is the autocorrelation function of x(t) (see table II).
2) Calculation in the frequency domain: The response of PVAR knowing the Power Spectral Density (PSD) of
frequency deviation Sy (f ) is:
σP2 (τ ) =
Z
∞
2
|HP y (f )| Sy (f )df.
(18)
0
Replacing Sy (f ) by the different power law noises, from h−2 f −2 (random walk FM) to h+2 f +2 (white PM), we
obtain the responses of PVAR for the different types of noises.
Table I gives the responses of PVAR for the different types of noise and for a linear frequency drift (y(t) = C1 ·t)
and Figure 5 shows their behavior versus the integration time τ .
The calculation can also be carried out in the time domain with (17). An example is given in Annex I for the
response of PVAR to white PM noise.
3) Monte-Carlo simulations: The response of PVAR may also be estimated by Monte-Carlo simulations. For
each type of noise, we realized 10 000 data runs of N = 2048 samples by using the noise simulator “bruiteur”3 .
The response of PVAR for each type of noise was estimated by averaging the 10 000 realizations from the data
runs.
3 This
software is available on request to the authors.
11
Fig. 5.
Responses of the P-variance for the different types of noise and for a linear frequency drift.
12
These 3 methods, from time domain, from frequency domain and by Monte-Carlo simulations give fully consistent
results.
IV. D EGREES OF FREEDOM AND CONFIDENCE INTERVALS
In this section, we compare the dispersion of the PVAR estimates with these of AVAR and MVAR.
A. Equivalent Degrees of Freedom (EDF)
We consider estimates of a generic variance σ 2 (τ ) that are assumed to follow a kχ2ν distribution, with k ∈ R,
whose Equivalent Degrees of Freedom (EDF) ν depends on the integration time τ . The mean and variance4 are
defined as:
  E σ 2 (τ )
 V σ 2 (τ )
kE χ2ν
= k 2 V χ2ν
=
= kν
=
2k 2 ν.
Thus, the degrees of freedom ν may be estimated from
2
2E σ 2 (τ )
.
(19)
ν=
V {σ 2 (τ )}
Knowing ν, it is easy to define a confidence interval around E σ 2 (τ ) with a p % confidence. As an example,
Figure 6 shows the histogram of the estimates of PVAR(512τ0 ) obtained from 10 000 realizations of a white FM
sequence.
We got the mean µ = 2.3500 · 10−3 and the variance σ 2 = 3.3264 · 10−6 . From (19), σ
bP2 (512τ0 ) is a random
bP2 (512τ0 ) = kχ2ν with ν = 2µ2 /σ 2 = 3.3204 and k = σ 2 /(2µ) =
variable proportional to a χ2ν random variable σ
7.0774 · 10−4 . It results the quantiles B2.5% = 0.29191 and B97.5% = 9.9381. The confidence interval over
σ
bP2 (512τ0 ) is then bounded by kB2.5% and kB97.5% . Thus, the most probable value of σ
bP2 (512τ0 ) is µ and 2.0659 ·
10−4 < σ
bP2 (512τ0 ) < 7.0336 · 10−3 with 95% or confidence. This is the classical way for defining the error bars
over the variance estimates.
For applying this result to PVAR, we have then to calculate the variance of PVAR.
B. Variance of variance
The variance of the estimate σ
ˆ 2 (τ ) is given by
"
( M −1 )#2 
−1

 1 M
n
o
X
1 X 2
2
.
V σ
ˆ 2 (τ ) = E σ
ˆ 2 (τ ) − E σ
ˆ 2 (τ )
=E
αi2 − E
αi

 M
M
i=0
i=0
Developing (20) yields
M −1 M −1
2 1 X X 2 2
V σ
ˆ (τ ) = 2
E αi αj − E αi2 E αj2 .
M i=0 j=0
4 Here,
it means the variance of the P-variance!
(20)
13
Fig. 6.
Histogram of 10 000 realizations of a variance estimator compared to the theoretical χ2 probability density function. The estimator is
2 (τ ) for τ = 512τ and was applied to 10 000 white FM noise sequences with h = 1 Hz−1 , τ = 1 s and 2048 samples. In this example,
σP
0
0
0
the number of Equivalent Degrees of Freedom are ν = 3.3. Therefore, the 2.5 % and 97.5 % bounds are respectively 2.1 · 10−4 and 7.0 · 10−3 .
From Isserlis’s theorem [10]–[12], if z and w are two centered, jointly Gaussian random variables, we know that
2
E z 2 w2 − E z 2 − E w2 = 2 [E {zw}] .
Assuming that x is a Gaussian process, αi and αj are two centered jointly Gaussian random variables, it comes
M −1 M −1
2 2 X X
2
[E {αi αj }] .
V σ
ˆ (τ ) = 2
M i=0 j=0
2 Thus, in order to calculate V σ
ˆ (τ ) , we need to calculate E {αi αj }.
(21)
14
C. Equivalent Degrees of Freedom of PVAR
From (13), we can calculate the following mathematical expectation
("m−1 #
X m−1
72
E
E {αi αj } =
− k (xi+k − xi+m+k )
m4 τ 2
2
k=0
#)
"m−1 X m−1
×
− l (xj+l − xj+m+l )
2
l=0
which yields after development
E {αi αj }
m−1 m−1 72 X X m − 1
m−1
−k
− l {2Rx [(i + k − j − l)τ0 ]
m4 τ 2
2
2
=
k=0 l=0
−Rx [(i + k − j − m − l)τ0 ] − Rx [(i + m + k − j − l)τ0 ]} .
(22)
Thanks to the relationships (21) and (22), it is possible to calculate the variance of the variance for each type of
noise from the autocorrelation function of this type of noise. For example, for a white PM noise, we found
m 2
2
9h2+2
m
m
V σ
ˆP (τ ) =
− 12
23
− 175 2 .
70π 4 τ 6
M
M
M
(23)
The calculation of the EDF yields
ν=
35
.
23m/M − 12(m/M )2 − 175m/M 2
(24)
Annex II gives the details of this formal calculation.
D. Numerical computation of EDF
The EDF may be computed from (21) and (22). The main issue of this later relationship is the calculation
of the autocorrelation function Rx (τ ) of the time error data, which must be defined for all types of noise. The
autocorrelation function may be defined as the Fourier transform of the Power Spectral Density. Since the PSD is
real an even [13], [14]:
Z
Rx (τ ) =
+∞
Sx (f ) cos(2πτ f )df.
(25)
0
Replacing Sx (f ) by the classical power laws ranging from white PM to f −4 PM, Mathematica software gives the
results displayed in Table II5 .
1
(Nyquist
We considered the following case N = 2048 data, τ0 = 1 s, the high cut-off frequency fh =
2τ0
1
frequency), the low cut-off frequency fl =
(for details about the physical meaning of fl see [14]) and
256N τ0
τ ∈ {τ0 , 2τ0 , 4τ0 , . . . , 1024τ0 }. Thus, we computed (22) for all types of noise and then to deduce V σP2 (τ ) from
(21). The corresponding degrees of freedom were obtained from (38). Figure 7 shows the EDF for all types of
noises.
5 Since
the functions Ci(x) and Si(x) are generally badly tabulated, it is strongly recommended to replace them by the following approximations
limx→0 Ci(x) = C + ln(x) −
limx→∞ Si(x) =
π
.
2
x2
96
where C ≈ 0.5772 is the Euler-Mascheroni constant, limx→∞ Ci(x) = 0, limx→0 Si(x) = x −
x3
9
and
15
Sx (f )
Rx (0)
k0
k0 fh
k−1 f −1
k−1
Rx (τ )
(for τ 6= 0)
0
1
+ ln(fh /fl )
2
k−1
cos(2πfl τ ) − 1 + 2πfl τ sin(2πfl τ )
(2πfl τ )2
+ Ci(2πτ fh ) − Ci(2πτ fl )
k−2 f −2
k−2
1
1
−
fl
fh
k−2
cos(2πfl τ )
cos(2πfh τ )
−
fl
fh
+2πτ [Si(2πfl τ ) − Si(2πfh τ )]
k−3 f −3
k−3
2
"
1
1
− 2
fl2
fh
(
#
k−3
cos(2πfl τ )
cos(2πfh τ )
−
2fl2
2fh2
+2π 2 τ 2 [Ci(2πfl τ ) − Ci(2πfh τ )]
sin(2πfl τ )
sin(2πfh τ )
−
+πτ
fh
fl
k−4
f −4
k−4
3
"
1
1
− 3
fl3
fh
#
(
k−4
−
(2π 2 fh2 τ 2 − 1) cos(2πfh τ ) + πfh τ sin(2πfh τ )
3fh3
(2π 2 fl2 τ 2 − 1) cos(2πfl τ ) + πfl τ sin(2πfl τ )
3fl3
4π 3 τ 3
+
[Si(2πfh τ ) − Si(2πfl τ )]
3
TABLE II
AUTOCORRELATION FUNCTION OF THE TIME ERROR DATA FOR THE DIFFERENT TYPES OF NOISES . Ci(x) AND Si(x) REPRESENT
RESPECTIVELY THE
C OSINE I NTEGRAL FUNCTION AND THE S INE I NTEGRAL FUNCTION .
E. Monte-Carlo simulations
Another possibility to assess the EDF of PVAR estimates consists in simulating series of 10 000 sequences of
N = 2048 samples for each types of noise by using the noise simulator “bruiteur”. It is then easy to obtain the
EDF by using (19). The results we obtain were very close to the ones obtained by numerical computation: generally
less than a few percents, except for τ = τ0 and 2τ0 because of the calculation of PVAR which cannot follow the
parabolic shape of the weighting function with only one or two samples (unless the samples were obtained by an
Ω-counter). Figure 8 compares the EDF given by the formal expression (24), by numerical simulations and the
experimental EDF obtained by Monte-Carlo simulations over 10 000 realizations of white phase noise. We observe
a good agreement between the three methods.
Table III and Figure 9 compare the EDF of PVAR to the ones of AVAR and MVAR.
F. Comparison of AVAR, MVAR and PVAR
PVAR exhibits slightly higher EDF than MVAR. On the other hand, the EDF of PVAR are generally lower than
those of AVAR, except for the random walk noise.
Table IV presents the response of the normalized variances. PVAR seems to be less efficient than MVAR since
16
Fig. 7.
Equivalent Degrees of Freedom of PVAR for all types of noise obtained from numerical computations. A zoom around τ = 256τ0
and 512τ0 allows us to separate the EDF for the different noises.
its noise response is always higher. However, we must keep in mind that PVAR has lower EDF and then tighter
confidence interval than MVAR. Moreover, PVAR may be calculated over τ -values 1.5 higer than MVAR. The
combination of all these factors must be carefully analyzed. The next section yields an efficiency way for comparing
the performances of the variances.
V. C ORNER DETECTION
We address the question of which variance is the most suitable to the detection of the noise processes of Table IV.
The criterion we choose is the lowest level of a process ‘X’ we can detect in the presence of a process ‘Y,’ with a
given data series. In a realistic example, we may search for the white FM noise (X) in the presence of the white
PM noise (Y) of the time-stamp counter. This is typical of passive atomic standards, as Cs beam or Rb cell.
For the comparison to be fair, all the variances are re-normalized for the same noise response. For example, it
17
Fig. 8.
Comparison of the EDF calculated by the formal expression (24) in Annex II, by numerical computation (see §IV-D) and assessed by
Monte-Carlo simulations over 10 000 realizations of white phase noise (see §IV-E).
is seen on Table IV that the response to white frequency noise Sy = h0 is h0 /2τ for the AVAR, h0 /4τ for the
MVAR, and 3h0 /5τ for the PVAR. Accordingly, a coefficient of 2, 4, or 5/3 is applied. The coefficients for the
other noise types are obtained from Table IV in the same way. Of course this normalization is useful at comparing
the variances, not for reporting on data.
Our comparison is based on a simulation with N = 2048 samples uniformly spaced by τ0 = 1 s. So, the lowest
value of τ is equal to 1 s, and the largest is equal to N τ0 /2 = 1024 s for AVAR and PVAR, and N τ0 /3 = 682 s
for MVAR. For each case, we averaged the results on 104 runs.
The detection criterion the is based on the uncertainty bars at 97.5% confidence level. We explain the process
referring to Fig. 10 A, which shows the detection of white FM noise σ 2 (τ ) ∝ 1/τ in the presence of white PM
noise. We use PVAR (blue red plot) as the reference, and we set the white FM noise (grey slanted line) at the
lowest level that PVAR can detect with 97.5% probability, i.e. at the upper point of the two-sigma uncertainty bar
at τ = 1014 s. This is highlighted by a grey circle at τ = 1014 s. We notice that the lowest value of MVAR (green
plot) at 97.5% confidence (grey circle at τ = 682 s) exceeds the reference white FM noise. Thus MVAR cannot
18
τ /τ0
1
2
4
8
16
32
64
128
256
512
682
1024
White PM (f +2 FM)
AVAR
892
1060
1020
1010
955
953
922
896
811
652
MVAR
891
970
685
355
173
82.5
38.9
17.3
7.48
2.88
PVAR
892
1150
824
419
202
99.1
46.9
22.0
10.0
4.13
1.03
0.930
Flicker PM
(f +1
0.981
1.02
FM)
AVAR
1090
1140
984
728
523
340
209
127
69.5
33.8
MVAR
1090
1020
544
258
126
62.1
29.3
13.9
5.73
2.09
PVAR
1090
1300
701
329
165
79.4
38.2
18.4
8.42
3.36
1.04
1.05
White FM (f 0 FM)
AVAR
1380
1200
716
372
186
91.7
45.3
21.8
10.2
4.07
MVAR
1380
1060
505
247
119
58.4
28.6
13.2
5.71
1.87
1.01
PVAR
1380
1390
680
319
157
76.7
37.5
18.2
8.43
3.32
1.01
1.02
1.04
Flicker FM (f −1 FM)
AVAR
1780
1200
595
299
150
72.8
36.1
17.1
7.58
3.05
MVAR
1780
1030
484
241
120
57.9
28.5
12.9
5.32
1.58
PVAR
1780
1470
648
319
159
77.8
38.2
18.2
8.01
3.16
Random walk FM
(f −2
AVAR
1990
1020
480
238
117
57.9
28.1
13.3
5.93
2.29
MVAR
1990
861
398
197
96.5
47.1
22.6
10.3
4.26
1.31
PVAR
1990
1290
548
266
131
64.3
31.2
14.8
6.53
2.49
1.02
1.02
FM)
1.01
1.02
1.02
TABLE III
C OMPARISON OF THE EDF OF AVAR, MVAR AND PVAR VERSUS τ AND THE NOISE TYPE . A LL NOISE SEQUENCES WERE SIMULATED
WITH A UNITY COEFFICIENT NOISE ,
CALCULATION SEQUENCE IS
2048 SAMPLES AND A SAMPLING FREQUENCY OF 1 H Z . S INCE THE EXTENSION OF THE MVAR
3τ , THE LARGER τ - VALUE IS THE LENGTH OF THE SEQUENCE DIVIDED BY 3 (682 HERE ) AND NOT BY 2
(1024 HERE ) AS FOR AVAR AND PVAR.
detect the such noise. As expected, AVAR shows poor behavior in the presence of white PM noise.
The same analysis is done for other noise types, as shown in Fig. 10.
A. Detection of White FM in the presence of white PM noise
White FM noise is present in all atomic standards and its fast detection is an important issue. Figure 10-A, alredy
discussed, shows that PVAR exhibits the highest sensitivity in the detection of white FM in the presence of white
PM. The classical AVAR is a poor option because of its τ −2 response law, versus the τ −3 of the other variances.
B. Detection of Linear Frequency Drift in the presence of white PM noise
We see on Fig. 10-B that PVAR exhibit superior sensitivity, while AVAR suffers from the same limitation as
above. Frequency drift is a severe limitation in cavity stabilized lasers, and in other precision oscillators based on
the mechanical properties of an artifact.
19
Fig. 9.
Comparison of the Equivalent Degrees of Freedom of AVAR, MVAR and PVAR for the different types of noise. All noise sequences
were simulated with a unity coefficient noise, 2048 samples and a sampling frequency of 1 Hz.
C. Detection of Flicker FM noise in the presence of white FM noise
Figure 10-C shows that the 3 variances detects the flicker FM almost at the same level, with a slight superiority
of AVAR and PVAR ex aequo. Since Cs clocks have no random walk and drift, flicker of frequency is the ultimate
limitation to long-term stability, thus to time keeping accuracy. In commercial standards, flicker FM exceeds the
white FM at approximately 1 day integration time. Thus, early detection of flicker enables early estimation of the
long term behavior, and a useful diagnostic. In this case however, since white PM is rejected, AVAR is the favorite
tool of time keepers.
20
Short term detection:
A. White FM in white PM
B. Drift in white PM
Long term detection:
C. Flicker FM in white FM
E. Drift in flicker FM
Fig. 10.
D. Random walk in flicker FM
F. Drift in random walk
Corner detection for the most common noise types (grey circles). The 97.5 % upper bounds of the confidence intervals over the
variance estimates are figured by dashed lines (blue for AVAR, green for MVAR and red for PVAR). The lowest detected noise or drift by
PVAR is represented as a solid black line.
21
AVar(τ )
MVar(τ )
2 (τ )
σP
2 (τ )
σP
MVar(τ )
h−2 f −2
2π 2 h−2 τ
3
11π 2 h−2 τ
20
26π 2 h−2 τ
35
1.4
h−1 f −1
2 ln(2)h−1
[27 ln(3) − 32 ln(2)] h−1
8
2 [7 − ln(16)] h−1
5
1.8
h0
2τ
h0
4τ
3h0
5τ
2.4
h+1 f +1
[1.038 + 3 ln(πτ /τ0 )]
4π 2 τ 2
[24 ln(2) − 9 ln(3)] h+1
8π 2 τ 2
3 [ln(16) − 1] h+1
2π 2 τ 2
3.2
h+2 f +2
3h+2
8π 2 τ0 τ 2
3h+2
8π 2 τ 3
3h+2
2π 2 τ 3
4
Drift C1 t
1 2 2
C τ
2 1
1 2 2
C τ
2 1
1 2 2
C τ
2 1
1
Sy (f )
(s)
h0 f 0
TABLE IV
C OMPARISON OF THE RESPONSES OF AVAR, MVAR AND PVAR FOR THE DIFFERENT TYPES OF NOISE AND A LINEAR FREQUENCY DRIFT.
D. Detection of slow phenomena
Figure 10-D-E-F show the detection of random walk and drift in the presence of flicker M, and the detection of
drift in the presence of random walk. This is typical of Rb clocks and H masers. In all cases AVAR is superior,
and yet close to PVAR, while MVAR is the poorest choice.
VI. C ONCLUSION
PVAR is a powerful tool for short term analysis because of its optimal rejection of white PM. But this study
proves that it is also almost as efficient as AVAR for long term analysis. It is even better than AVAR for studying
the long term stability of Cs standards.
VII. ACKNOWLEDGEMENTS
This work is supported by the ANR Programme d’Investissement d’Avenir in progress at the TF Departments of
FEMTO-ST Institute and UTINAM (Oscillator IMP, First-TF and Refimeve+), and by the R´egion de Franche-Comt´e.
We wish to thank Charles Greenhall for his valuable help concerning Isserlis’theorem.
22
R EFERENCES
[1] D. W. Allan, “Statistics of atomic frequency standards,” Proceedings of the IEEE, vol. 54, no. 2, pp. 222–231, 1966.
[2] J. J. Snyder, “Algorithm for fast digital analysis of interference fringes,” Applied Optics, vol. 19, no. 4, pp. 1223–1225, Apr. 1980.
[3] D. W. Allan and J. A. Barnes, “A modified “Allan variance” with increased oscillator characterization ability,” in Proc. 35 IFCS,
Ft. Monmouth, NJ, May 1981, pp. 470–474.
[4] J. J. Snyder, “An ultra-high resolution frequency meter,” in Proc. 35th Annual Frequency Control Symposium, May 27–29 1981, pp.
464–489.
[5] G. Kramer and K. Klitsche, “Multi-channel synchronous digital phase recorder,” in Proceedings of the 2001 IEEE International Frequency
Control Symposium, Seattle (WA), USA, June 2001, pp. 144–151.
[6] E. Rubiola, “On the measurement of frequency and of its sample variance with high-resolution counters,” Review of Scientific Instruments,
vol. 76, no. 054703, May 2005, article no. 054703.
[7] S. T. Dawkins, J. J. McFerran, and A. N. Luiten, “Considerations on the measurement of the stability of oscillators with frequency counter,”
IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 54, no. 5, pp. 918–925, May 2007.
[8] E. Benkler, C. Lisdat, and U. Sterr, “On the relation between uncertainties of weighted frequency averages and the various types of Allan
deviations,” arXiv:1504.00466, April 2015. [Online]. Available: http://arxiv.org/pdf/1504.00466v3.pdf
[9] P. Lesage and C. Audoin, “Characterization of frequency stability: uncertainty due to the finite number of measurements,” IEEE Transactions
on Instrumentation and Measurement, vol. 22, no. 2, pp. 157–161, June 1973.
[10] L. Isserlis, “On certain probable errors and correlation coefficients of multiple frequency distributions with skew regression,” Biometrika,
vol. 11, pp. 185–190, 1916.
[11] ——, “On a formula for the product-moment coefficient of any order of a normal frequency distribution in any number of variables,”
Biometrika, vol. 12, pp. 134–139, 1918.
[12] C. Greenhall, “‘Greenhall’s proof’ in wikipedia entry ‘Isserlis’ theorem’,” March 2015. [Online]. Available: http://en.wikipedia.org/wiki/
Isserlis%27 theorem
[13] F. Vernotte, J. Delporte, M. Brunet, and T. Tournier, “Uncertainties of drift coefficients and extrapolation errors: Application to clock error
prediction,” Metrologia, vol. 38, no. 4, pp. 325–342, December 2001.
[14] F. Vernotte and E. Lantz, “Metrology and 1/ f noise: linear regressions and confidence intervals in flicker noise context,” Metrologia,
vol. 52, no. 2, pp. 222–237, 2015. [Online]. Available: http://stacks.iop.org/0026-1394/52/i=2/a=222
A NNEX I: R ESPONSE OF THE PVAR FOR A WHITE NOISE
In the discrete case, we calculate the variance σP2 (τ ) in the case of a white noise from the expression (17). For
a white noise, Rx (t) = 0 if t 6= 0 and Rx (0) = k0 fh where k0 is the noise level and fh the high cut-off frequency.
Since l > k − m and l < k + m for all k, we obtain
σP2 (τ )
=
2
m−1 72 X m − 1
− k 2Rx (0)
m4 τ 2
2
=
144k0 fh m(m2 − 1)
12k0 fh
≈
.
4
2
m τ
12
mτ 2
k=0
Because of sampling, fh = 1/(2τ0 ). On the other hand, the phase noise level k0 may be written as the corresponding
frequency noise level: h+2 = 4π 2 k0 . Finally, for large τ -values (m 1), it comes:
σP2 (τ ) ≈
3h+2
.
2π 2 τ 3
In the continuous setting, for a noise with power spectral density Sy (f ) = h2+2 |f |2 we can show that σP2 (τ ) =
Indeed,
σP2 (τ ) =
1
2
Z
R
F −1 (Sy )(t)
Z
hy (t − s)hy (−s) ds dt
R
(26)
3h2+2
2π 2 τ 3 .
23
where F(u)(f ) =
R
R
u(t)e−2iπtf dt is the Fourier transform. Moreover, in this case F −1 (Sy )(t) = − 4π1 2 δ 00 (t) so
Z
Z
h2+2
2
00
σP (τ ) = − 2
δ (t) hy (t − s)hy (−s) ds dt
8π R
R
Z
Z
h2+2
δ(t) h0y (t − s)h0y (−s) ds dt,
=
8π 2 R
R
2 Z
h+2
=
(h0 )2 (s) ds.
(27)
8π 2 R y
Notice that the second step has been obtained by converting the derivative h0y (t − s) into the partial deriva√
tive ∂hy (t − s)/∂s and applying an integration by parts. Finally, using h0y (s) = − 3τ 32 (2s + τ ) χ(−τ,0) (s) −
√
3 2
τ3
(τ − 2s) χ(0,τ ) (s) because hy vanishes at the discontinuity points gives σP2 (τ ) =
3h2+2
2π 2 τ 3
as expected.
A NNEX II: C ALCULATION OF THE EDF OF PVAR IN THE PRESENCE OF WHITE NOISE
A. Autocorrelation of the elementary estimates in the presence of white noise
We restart from (22) for a white noise. Since E {αi αj } does not depend on the sign of i − j, let us define
d = |i − j|. We must distinguish three cases,
•
0 ≤ d < m:
E {αi αj }
"m−d−1 #
X
m−1
m−1
72k0 fh
−k
−k−d
=
m4 τ 2
2
2
k=0
" m−1 #
X
m−1
m−1
×
−k
−k−d+m
2
2
k=m−d
=
•
6k0 fh
× 6d3 − 6md2 − 3m2 d + 3d + 2m3 − 2m
4
2
m τ
(28)
m ≤ d < 2m:
E {αi αj }
=
=
72k0 fh
m4 τ 2
"2m−d−1 X
m−1
k=0
2
−k
#
m−1
−k−d+m
2
6k0 fh
× (2m − d)(1 + 2d2 − 2md − m2 )
m4 τ 2
(29)
m ≥ 2m:
E {αi αj } = 0.
2
We can now calculate V σ
ˆP (τ ) from (21).
B. Calculation of the variance of PVAR in the presence of white noise
To explain the principle of calculation, let us begin with the simple case m = N/4.
(30)
24
1) Case m = N/4: Since E {αi αj } = E {αM −i αM −j }, we may calculate only the first half of the autocorrelations:
M/2−1 M −1
2
4 X X
2
V σ
ˆP (τ ) = 2
[E {αi αj }] .
M
i=0
j=0
(31)
Since m = N/4, the number of estimates is M = N − 2m = 2m. Hence, αi should vary from α0 to αm−1 , i.e.
m possible values. We have then to distinguish two regions for αj :
•
from j = 0 to m + i − 1, |i − j| < m and (28) must be used
•
from j = m + i to j = M − 1 = 2m − 1, m ≤ |i − j| < 2m and we ought to use (29).
Let us denote:
em
=
6d3 − 6md2 − 3m2 d + 3d + 2m3 − 2m
(32)
e2m
=
(2m − d)(1 + 2d2 − 2md − m2 )
(33)
the numerical parts of (28) and (29).
The variance of PVAR may be calculated by:


2 m−1
2m−1
X m+i−1
X
X
2
6k0 fh
4

e2m +
e22m  .
V σ
ˆP (τ ) = 2
M
m4 τ 2
i=0
j=0
j=m+i
(34)
The computation of (34) with Mathematica yields:
2
2
4
6k0 fh
65m2 − 98m4 − 35m6 + 68m8
V σ
ˆP (τ ) = 2
M
m4 τ 2
35
with M = 2N :
2
36k 2 f 2
65m2 − 98m4 − 35m6 + 68m8
V σ
ˆP (τ ) = 100 4h ×
m τ
35
(35)
Since m 1, we may neglect the terms in m2 , m4 and m6 regarding m8 and (35) becomes:
2
2448k02 fh2
V σ
ˆP (τ ) =
.
35m2 τ 4
(36)
By using k0 = h+2 /(4π 2 ) and fh = 1/(2τ0 ), it comes:
2
153h2+2
V σ
ˆP (τ ) =
.
140π 4 τ 6
(37)
Now, we can deduce the Equivalent Degrees of Freedom (EDF) of PVAR estimator. The EDF ν may be calculated
from:
ν=
2
2
2 E σ
ˆP (τ )
.
V {ˆ
σP2 (τ )}
(38)
From (26), (37) and (38), we obtain:
ν=
70
≈ 4.1176.
17
A numerical computation of (21) and (22) gives the following numerical value: νnum ≈ 4.1203, i.e. a difference of
0.06 %.
By Monte-Carlo simulations over 10 000 realizations of white phase noise, we got an experimental value νexp ≈
4.1290, i.e. a difference of 0.28 %.
25
2) Case m ≤ N/6: Following the previous example, we see that the domain of i should be divided in three
regions:
1) 0 ≤ i < m in which the domain of j has to be separated in
•
0 ≤ j < i + m where em must be used
•
i + m ≤ j < i + 2m where e2m must be used
2) m ≤ i < 2m − 1 in which the domain of j has to be separated in
•
0 ≤ j < i − m where e2m must be used
•
i − m ≤ j < i + m where em must be used
•
i + m ≤ j < i + 2m − 1 where e2m must be used
3) 2m ≤ i < M/2 − 1 in which the domain of j has to be separated in
•
i − 2m ≤ j < i − m where e2m must be used
•
i − m ≤ j < i + m where em must be used
•
i + m ≤ j < i + 2m − 1 where e2m must be used.
Therefore, the variance of PVAR may be calculated by:



2 m−1
i+2m
X
X
X i+m−1
2
4
6k0 fh

e22m 
e2m +
V σ
ˆP (τ )
=

M 2 m4 τ 2
j=i+m
j=0
i=0


2m−1
i+m−1
i+2m
X i−m
X
X
X

+
e22m +
e2m +
e22m 
+
i=m
j=0
M/2−1

X
i−m
X

i=2m
j=i−2m
j=i−m
e22m +
i+m−1
X
j=i−m
j=i+m
i+2m
X
e2m +
j=i+m


e22m  .

(39)
Mathematica gives the following result:
2
2
6k0 fh
1 4
V σ
ˆP (τ )
=
−m 165m + 966m3 + 1365m5 + 24m7
×
2
4
2
M
m τ
35
+m6 (350 − 46M ) − 45M + 35m4 (52 + M ) + 14m2 (25 + 4M )
Keeping only the terms of higher degrees, it comes:
m 2
2
9h2+2
m
m
V σ
ˆP (τ ) =
−
12
−
175
23
.
70π 4 τ 6
M
M
M2
The calculation of the EDF yields:
ν=
35
.
23m/M − 12(m/M )2 − 175m/M 2