[time-nuts] Calculated ADEV does not match expectation

Wolfgang Wallner wolfgang-wallner at gmx.at
Sun Apr 24 14:01:21 UTC 2016


Hello Anders, hello timenuts,

This is in reply to the thread "[Announce] Simulation software for 
powerlaw noise and PTP clock synchronization" [1].
I have changed the mail subject as the discussion has moved quite far 
from the initial message.

[1] https://www.febo.com/pipermail/time-nuts/2016-April/097272.html


 > Do you get agreement between PLN time-series and calcluated ADEVs, as
 > per IEEE1139-2008 table B.2 ?
 >
 > I tried [3] using python libraries colorednoise [1] (also based on
 > Kasdin&Walter) and allantools [2], but it gives ADEVs that are
 > systematically higher than predicted by theory. OTOH the calculated
 >  MDEVs seem to fall on the line that theory predicts. Image here:
 > http://www.anderswallin.net/wp-content/uploads/2016/04
 > /colorednoise-1.png


I think it might be the best to discuss this with an actual example.

I had a look at the image that you provided and tried to compute a noise 
sample which is similar, and analyze it with the tools I use.

I will split the discussion into three parts, so that it is easier to 
follow the individual steps:

1) Calculating expected values
2) Generating a PLN sample
3) Analyzing the sample and compare with expectation


1) Calculating expected values
----------------------------------------------------

As already mentioned in my earlier reply, I had quite some troubles 
figuring out how all the formulas in IEEE 1139 and the Kasdin/Walter 
paper work, and I'm still not too confident. If any of the following 
steps raise doubts, please point it out.

I guessed the following values based on your image:

Sampling frequency f_s = 1Hz.
This implies that the highest measured frequency can be f_h = 0.5Hz
Also, it implies that the distance between two measurements Tau_0 = 1s.
For the example, I chose to use FFM noise.
This means alpha = -1 (Power S_y(f) ~ f^alpha), and mu = 0 (AVAR(Tau) ~ 
Tau^mu).
 From the top right figure I guess that S_y(10-3Hz) = 10^-19.5.

 From IEEE 1139 Table B.2 we know that: S_y(f) = h_alpha * f^alpha

In our example this means that S_y(10^-3) = h_-1 * (10^-3)^-1, which 
implies that h_-1 = 10^-22.5.

 From IEEE 1139 we also know that B = 2 ln(2) = 1.3863.

Thus, the expected AVAR(Tau) would be B * h_-1 * Tau^0 = 1.3863 * 
10^-22.5 = 4.3838 * 10^-23 for all Tau.
Using ADEV = Sqrt(AVAR) we get that the expected ADEV(Tau) = 9.9346 * 
10^-23 for all Tau.

The important values are thus:
h_-1 = 10^-22.5
ADEV(Tau) = 9.9346 * 10^-23

Judging from the lower left plot in your image, I guess that we agree here.


2) Generating a PLN sample
----------------------------------------------------

The Kasdin/Walter PLN-generation approach starts with white noise, and 
applies filters to it. The standard variance Qd is an important 
configuration parameter.

The original paper gives a formula as follows:

Qd = h_alpha/[ 2 * (2 pi)^alpha * Tau_0^(alpha-1) ]

This formula leads to unexpected values when I use it with my 
implementation. However, the following modified version works as expected:

Qd = h_alpha/[ 2 * (2 pi)^alpha * Tau_0^(alpha+1) ]

The difference is the sign of the 1 at the end. I don't know if this is 
an error in the original paper, or if I have an error in my 
implementation and both errors compensate each other.

Using my modified formula, and the values for h_alpha and Tau_0 as 
stated above, I get

Qd = 9.9345 * 10^-23

As stated in my last mail, LibPLN comes with a small console application 
in its 'Demos' folder, which is called PLN_Generator.
PLN_Generator provides a simple interface to the most important parts of 
LibPLN, and can be used to easily genrate PLN samples:

./PLN_Generator --sampling-frequency 1Hz --pln-filter-method 
kasdin-walter --qd 9.9346E-23 --alpha -1 --kw-filter-length 1000000 
--number-of-samples 1000000 --output-filename 'ffm.txt' --verbose

This generates a file called 'ffm.txt', which contains 10^6 samples of 
Time Derivative values (or phase observations as you call it).
The file is 35MB in size thus I haven't attached it to the mail.
But I have uploaded it into my Dropbox account, and you can find it here:

https://www.dropbox.com/s/k8tj1b0ywptat3z/ffm.txt?dl=0


3) Analyzing sample and compare with expectation
----------------------------------------------------

I have used my Matlab Scripts scripts to analyze this sample file.
Attached to the mail you can find the following images:

TimeDerivative.png: Plot of the time derivative over time (so basically 
the raw information that is contained in the sample file)

ADEV.png: Plot of the calculated ADEV and the expected values as 
calculated in Section 1

PSD.png: Plot of the calculated PSD and the expected values as given in 
section 1.
The PSD is plotted in dB, calculated as PSD_db = 10 log10( PSD )

I plan to publish my Matlab scripts as GPL, so that they can be used 
together with LibPLN. But I haven't done so yet, mostly because they are 
just ugly prototypes right now.

For the calculation of the Allan Deviation I use the script 'allan' from 
Matlab's file exchange:
http://www.mathworks.com/matlabcentral/fileexchange/13246-allan

In my scripts I use it as follows:

   % Calculation of actual ADEV
   tau = generate_tau( length(td) / 1000, f_s );
   s   = struct ('phase', td, 'rate', f_s );

   [AdevAct, ~, ~, Tau] = allan(s, tau, '', 0);

generate_tau is a simple script that generates logarithmic values with 
the common money-step values (1, 2, 5, 10, 20, 50, ...).

For the calculation of the PSD I use the Matlab function pwelch with a 
hanning window. For the generated PSD plot I split the sample into 100 
segments for the Welch's method (e.g. the Hanning window length was 
floor(NumSamples/100)).


As you can see in the attached plots, the expected and actual PSD and 
ADEV values agree. What would be the results if you apply your analyzing 
methods to the same sample file?


4) Additional comments
----------------------------------------------------

In the ADEV plot in your image, the lines for the expected ADEV of FPM 
and WPM noise seem to be exactly parallel. Is this the case? It is hard 
to judge from the figure.

Because I think it should not be the case:
E only depends on f_h, and is thus constant for all Tau.
But D depends on both f_h and Tau, and thus it changes for different Tau 
values.


regards, Wolfgang


On 04/19/2016 07:46 AM, Wolfgang Wallner wrote:
> Hello Anders,
>
>> Do you get agreement between PLN time-series and calcluated ADEVs, as per
>> IEEE1139-2008 table B.2 ?
>
> Yes, I explicitly cite IEEE1139-2008 table B.2 in my thesis. I will
> write a more detailed reply with example data later when I get home.
>
> If you would like to try out the noise generations of libPLN: I also
> wrote a small command line utility for LibPLN, called PLN_Generator. You
> can find it in the Demo/ directoy [1].
>
>> I tried [3] using python libraries colorednoise [1] (also based on
>> Kasdin&Walter) and allantools [2], but it gives ADEVs that are
>> systematically higher than predicted by theory. OTOH the calculated MDEVs
>> seem to fall on the line that theory predicts. Image here:
>> http://www.anderswallin.net/wp-content/uploads/2016/04/colorednoise-1.png
>
> I have only used ADEV and AVAR so far, so I don't know about the
> MDEV-behavior of my simulated noise. I will try out your allan tools
> when I get home.
>
> But the conversion between AVAR and power spectral density has been
> quite confusing for me, see [2] for an example. I thought I had figured
> it out at the end, but maybe I'm still wrong :)
>
>> On your libPLN page the time-series graph is called 'time deviation'
>> which
>> could be a bit confusing as there is also an ADEV-like statistic called
>> time deviation. Perhaps time-series or 'phase observations' or similar is
>> better?
>
> Yes, I agree that the wording is very unfortunate. IEEE 1139 first calls
> x(t) time fluctuations and time derivative, but later refers to it as
> time deviation in the text.
>
> I only realized quite late that there is a name conflict with another
> statistic measure.
>
> I guess the easiest solution would be to replace all occurrences of
> 'time deviation' with e.g. 'time derivative'.
>
> regards, Wolfgang
>
> [1] https://github.com/w-wallner/libPLN/tree/master/Demos/PLN_Generator
> [2] https://www.febo.com/pipermail/time-nuts/2015-March/091098.html
> _______________________________________________
> time-nuts mailing list -- time-nuts at febo.com
> To unsubscribe, go to
> https://www.febo.com/cgi-bin/mailman/listinfo/time-nuts
> and follow the instructions there.
-------------- next part --------------
A non-text attachment was scrubbed...
Name: ADEV.png
Type: image/png
Size: 8782 bytes
Desc: not available
URL: <http://febo.com/pipermail/time-nuts_lists.febo.com/attachments/20160424/1ab888c8/attachment.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: PSD.png
Type: image/png
Size: 7306 bytes
Desc: not available
URL: <http://febo.com/pipermail/time-nuts_lists.febo.com/attachments/20160424/1ab888c8/attachment-0001.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: TimeDerivative.png
Type: image/png
Size: 8635 bytes
Desc: not available
URL: <http://febo.com/pipermail/time-nuts_lists.febo.com/attachments/20160424/1ab888c8/attachment-0002.png>


More information about the Time-nuts_lists.febo.com mailing list