[time-nuts] AVAR <-> S_Y conversion

[Wolfgang Wallner]

Hello time-nuts community,

thanks to your feedback and that of others, I could make additional
progress on my journey to understand powerlaw noise :)

I would like to reiterate what I have learned so far. Please comment on
anything that you think is done wrong.

-----------------------------------------------
0) Conventions:
-----------------------------------------------

My main goal is to simulate powerlaw noise and analyze it as described
in IEEE1139 [1]. I will use the following conventions:

f_s     sampling frequency of the simulated noise
tau_0   time interval between two samples (directly related as 1/f_s)

y(t)    Fractional frequency deviation at time t
S_y(f)  one-sided PSD of y, has a shape of as h_alpha * f^alpha

alpha has one of the following values:

2   White PM noise
1   Flicker PM noise
0   White FM noise
-1  Flicker FM noise
-2  Random Walk noise

-----------------------------------------------
1) Using the formulas in IEEE 1139:
-----------------------------------------------

As I said, I mainly use the formulas given in IEEE1139, especially those
in Table B.2, which define the relationship between Allan Variance and
the S_y. I have attached a screenshot of those formulas to this mail.

I'm not sure what f_h should be in the calculation of the terms D and E.
IEEE 1139 defines it as the "high-frequency cutoff of an infinitely
sharp low-pass filter".

I don't sample my noise from a real device, but simulate it.
Am I right that I can use f_h = f_s/2 (Nyquist theorem)? I haven't found
a publication that explicitly states this, but in my experiments this
assumptions works well.

-----------------------------------------------
2) Noise generation, calculation of h_alpha:
-----------------------------------------------

I generate powerlaw noise according to the publication by Kasdin and
Walter [2]. As you might have noticed in my initial mail, I estimated
h_alpha from an Allan Variance plot. This is not how it should be done.
The better way would be to estimate it from my noise configuration.

The reason why I went the other way is that I had trouble to estimate
h_alpha from my noise configuration. The approach described in [2]
generates white noise, and filters it to get the required PSD shape. The
relationship between the standard variance Qd of the white input noise
and the scaling factor h_alpha of the powerlaw noise is given in [2] as
follows (equ 39):

Qd = h_alpha / (2 * (2*pi)^alpha * tau_0^(alpha-1))

However, this definition never worked for me to predict the relationship
between h_alpha and Qd. I think the formula should be modified as

Qd = h_alpha / (2 * (2*pi)^alpha * tau_0^(alpha+1))

I changed the sign of the 1 in the exponent of tau_0.
Using this change, the formula now works for calculating h_alpha and Qd
from each other, and the results match if I do a counter-check and
estimage h_alpha from the AVAR or the PSD.

This change also makes the formula more consistent (e.g. the AVAR is
defined so that the standard variance of White FM noise should match the
AVAR for tau_0, this holds with the modified formula).

-----------------------------------------------
3) PSD estimation
-----------------------------------------------

I tried to implement the PSD estimation as a mixture of the information
found in [3] and [4]. However, I'm a novice when it comes to PSDs, and
my approach had some error (I still don't know exactly what was wrong).
As I know now, the 3dB difference that I saw for RW noise in my initial
mail is a bug in my clumsy implementation.

I tried to learn more about PSDs, and [5] proved to be very useful.
I know now that the PSD estimation approach that I tried to use is known
as the 'Welch's method' and supported in Matlab as 'pwelch'.

--> Using this tool the PSD estimate converges to the expected value for
all 5 types of noise! :)

As PSD estimation configuration I use non-overlapping segments with a
Hann window. There is no deeper reason for this choice (as I said, I'm
new to these topics, so I would'nt know any better), it's just what is
used in the example in [5] and it provides a PSD estimate as I would
expect it.

-----------------------------------------------
4) All summed up
-----------------------------------------------

With the assumptions and concepts of 1-3 I have finally been able to
generate powerlaw noise in a way that the results match what I had
configured :)

I tried once more to generate the 5 noise types, and compare them with
my expectations. I have included the resulting plots in a PDF file,
which is available here:

https://www.dropbox.com/s/lrdbpxrghkca0y8/Relationship_PSD_AVAR.pdf?dl=0

However, I also found new things that confuse me :P

When I try to estimate the PSD with the 'pwelch' function in Matlab, I
can select the number of non-overlapping segments my raw data should be
divided into. Using a larger number of segments leads to a nicer plot,
which converges to clean lines at some point.

However, I also see that the lower part of the frequency spectrum
increases with the number of configured segments. I attached a plot
(PSD_NumSegments.png) to this mail to visualize what I mean.

I don't know why this is the case. Maybe someone of you can point me to
an explanation. I have already been told that the power in the frequency
domain and in the time domain has to be the same (i.e.
(mean(abs(ffd)^2)) == (sum(psd*fs/n_dft))). This condition holds for any
configured number of segments, and would be violated if the lower part
of the spectrum would not increase, so I guess these things have
something to do with each other.

If I would state it with my novice words, I would say "the left out
spectrum parts due to the segmentation process leak into the lower
frequency bins". However, I would be glad to here any real explanation :D

Thanks so far,
Wolfgang

[1] IEEE 1139
[2] Kasdin and Walter, Discrete Simulation of Power Law noise, 1992
[3] http://www.dspguide.com/ch9/1.htm
[4] http://de.mathworks.com/help/signal/ug/psd-estimate-using-fft.html
[5] Spectrum and spectral density estimation by the Discrete Fourier
transform (DFT), G. Heinzel et al

-------------- next part --------------
A non-text attachment was scrubbed...
Name: IEEE1139.png
Type: image/png
Size: 56227 bytes
Desc: not available
URL: <http://www.febo.com/pipermail/time-nuts/attachments/20150322/935be2e6/attachment-0002.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: PSD_NumSegments.png
Type: image/png
Size: 20121 bytes
Desc: not available
URL: <http://www.febo.com/pipermail/time-nuts/attachments/20150322/935be2e6/attachment-0003.png>