Hello time-nuts community, I hope this is the right place for the following question :) I'm dealing with the simulation of powerlaw noise, and I stumbled upon something I cannot explain when I tried out some formulas of IEEE 1139 [1]: According to Table B.2 in [1] the one-sided power spectral density of fractional frequency data and the Allan variance of this data can be related as follows: PSD of FFD: S_y(f) = h_alpha * f1^alpha AVAR: Sigma_y(Tau) = K * h_alpha * f ^ x where K and x depend on the type of noise (alpha). For your convenience I made a screenshot of the formulas I'm referring to: http://postimg.org/image/6qcx3ggu9/ I have generated data sets with simulated powerlaw noise for different values of alpha, and did the following for each of these noise vectors: *) Calculated and plotted the Allan Variance *) Used to formulas of [1] to calculate h_alpha *) Calculated and plotted the FFD-PSD (the PSD plot is averaged the get better visual results) *) Added colored lines to both plots according to the calculated h_alpha values I would have expected that the colored lines would match each plot. However, this is only the case for White PM, Flicker PM, White FM and Flicker FM noise. To my surprise the calculated line for random walk noise does not match the PSD plot. For the random walk noise the expected line is off by a factor of exactly 2 from the calculated plot, and I don't know how to explain this behavior. I supposed that maybe the factor A is twice as large as it should be, and thus I searched in other powerlaw noise publications for different formulas. However, as far as I can see they agree with the definition given in [1]. I could only find one paper with another definition: [2] In that other paper A is defined as 2 pi^2/6 instead of 2 pi^2/3 (equation 24). Using this definition would result in a plot that matches what I would have expected. These are the graphs I'm referring to: White PM: http://postimg.org/image/fk059s243/full/ Flicker PM: http://postimg.org/image/6q71l03cp/full/ White FM: http://postimg.org/image/mxhxeszqx/full/ Flicker FM: http://postimg.org/image/3vzpj7jmf/full/ Random Walk: http://postimg.org/image/hxad6okwv/full/ <-- the bad guy Does anyone know the reason for the behavior I see? best regards, Wolfgang Wallner PS: I tried to keep this mail short. If I have left out any information that would be useful feel free to ask, please :) [1] 1139-2008 - IEEE Standard Definitions of Physical Quantities for Fundamental Frequency and Time Metrology [2] Gaderer, et al - Achieving a Realistic Notion of Time in Discrete Event Simulation

		Attila Kinali

Servus! On Thu, 05 Mar 2015 14:35:51 +0100 Wolfgang Wallner <wolfgang-wallner@gmx.at> wrote: > For the random walk noise the expected line is off by a factor of > exactly 2 from the calculated plot, and I don't know how to explain this > behavior. I'm probably the wrong one to answer, as I have never done any noise simulation or even read up the relevant papers, but... A factor of 2 sounds like the difference you would get between one sided and two sided noise PSD's. Attila Kinali -- < _av500_> phd is easy < _av500_> getting dsl is hard

List -- I apologize for letting Wolfgang's posting earlier today slip though the filter. Normally we pre-screen all URL's to prevent exposure to questionable "free image hosting" sites. Wolfgang -- I have taken your 6 images and attached them to this posting. This way members can view your equations and plots without being exposed to postimg.org and its third party affiliates. The original text of your excellent posting is included below, with attachments instead of spam URL's. Thanks, /tvb ----- Original Message ----- From: "Wolfgang Wallner" <wolfgang-wallner@gmx.at> To: "Time-Nuts" <time-nuts@febo.com> Sent: Thursday, March 05, 2015 5:35 AM Subject: [time-nuts] AVAR <-> S_Y conversion > > Hello time-nuts community, > > I hope this is the right place for the following question :) > > I'm dealing with the simulation of powerlaw noise, and I stumbled upon > something I cannot explain when I tried out some formulas of IEEE 1139 [1]: > > According to Table B.2 in [1] the one-sided power spectral density of > fractional frequency data and the Allan variance of this data can be > related as follows: > > PSD of FFD: S_y(f) = h_alpha * f1^alpha > > AVAR: Sigma_y(Tau) = K * h_alpha * f ^ x > > where K and x depend on the type of noise (alpha). > > For your convenience I made a screenshot of the formulas I'm referring > to: [IEEE1139_Table_B_2.png] > > I have generated data sets with simulated powerlaw noise for different > values of alpha, and did the following for each of these noise vectors: > > *) Calculated and plotted the Allan Variance > *) Used to formulas of [1] to calculate h_alpha > *) Calculated and plotted the FFD-PSD (the PSD plot is averaged the get > better visual results) > *) Added colored lines to both plots according to the calculated h_alpha > values > > I would have expected that the colored lines would match each plot. > However, this is only the case for White PM, Flicker PM, White FM and > Flicker FM noise. To my surprise the calculated line for random walk > noise does not match the PSD plot. > > For the random walk noise the expected line is off by a factor of > exactly 2 from the calculated plot, and I don't know how to explain this > behavior. > > I supposed that maybe the factor A is twice as large as it should be, > and thus I searched in other powerlaw noise publications for different > formulas. However, as far as I can see they agree with the definition > given in [1]. I could only find one paper with another definition: [2] > In that other paper A is defined as 2 pi^2/6 instead of 2 pi^2/3 > (equation 24). Using this definition would result in a plot that matches > what I would have expected. > > These are the graphs I'm referring to: > > White PM: [White_PM.png] > Flicker PM: [Flicker_PM.png] > White FM: [White_FM.png] > Flicker FM: [Flicker_FM.png] > > Random Walk: [Random_Walk.png] <-- the bad guy > > Does anyone know the reason for the behavior I see? > > best regards, Wolfgang Wallner > > PS: I tried to keep this mail short. If I have left out any information > that would be useful feel free to ask, please :) > > [1] 1139-2008 - IEEE Standard Definitions of Physical Quantities for > Fundamental Frequency and Time Metrology > [2] Gaderer, et al - Achieving a Realistic Notion of Time in Discrete > Event Simulation > _______________________________________________ > time-nuts mailing list -- time-nuts@febo.com > To unsubscribe, go to https://www.febo.com/cgi-bin/mailman/listinfo/time-nuts > and follow the instructions there.

On 03/05/2015 07:23 PM, Attila Kinali wrote: > Servus! Servus :) > On Thu, 05 Mar 2015 14:35:51 +0100 > Wolfgang Wallner <wolfgang-wallner@gmx.at> wrote: > >> For the random walk noise the expected line is off by a factor of >> exactly 2 from the calculated plot, and I don't know how to explain this >> behavior. > > I'm probably the wrong one to answer, as I have never done any noise > simulation or even read up the relevant papers, but... > A factor of 2 sounds like the difference you would get between one sided > and two sided noise PSD's. > I calculate the one-sided PSD of the FFD data as described in [1] (first paragraph), so the code looks like this: xdft = fft(x); xdft = xdft(1:N/2+1); psdx = (1/(Fs*N)) * abs(xdft).^2; psdx(2:end-1) = 2*psdx(2:end-1); Remark: Before calculating the PSD, I split the data into parts of equal size, calculate the PSD for each one, and average over the set of PSDs. This improves the graphical visualization a lot. As the result matches my expectation exactly for 4 different kinds of noise, I would have assumed that this PSD calculation approach is quite reasonable. As I see the unexpected behavior only with random walk noise, and the main difference in the calculation is the term A, I would suspect that it has something to do with it. However, I'm a novice in this field, so any hint is very appreciated. regards, Wolfgang [1] http://de.mathworks.com/help/signal/ug/psd-estimate-using-fft.html

xdft = fft(x);
xdft = xdft(1:N/2+1);
psdx = (1/(Fs*N)) * abs(xdft).^2;
psdx(2:end-1) = 2*psdx(2:end-1);

Wolfgang, I have checked several sources, and they match up with the IEEE 1139 in this regard. I have also evaluated the equation for Allan variance for the random walk noise, and it matches up with the references and what I put here: https://en.wikipedia.org/wiki/Allan_variance#Power-law_noise So, the A formula you have matches up. You will need to find another source of the mismatch. Cheers, Magnus On 03/06/2015 11:04 AM, Wolfgang Wallner wrote: > > > On 03/05/2015 07:23 PM, Attila Kinali wrote: >> Servus! > > Servus :) > >> On Thu, 05 Mar 2015 14:35:51 +0100 >> Wolfgang Wallner <wolfgang-wallner@gmx.at> wrote: >> >>> For the random walk noise the expected line is off by a factor of >>> exactly 2 from the calculated plot, and I don't know how to explain this >>> behavior. >> >> I'm probably the wrong one to answer, as I have never done any noise >> simulation or even read up the relevant papers, but... >> A factor of 2 sounds like the difference you would get between one sided >> and two sided noise PSD's. >> > > I calculate the one-sided PSD of the FFD data as described in [1] (first > paragraph), so the code looks like this: > > xdft = fft(x); > xdft = xdft(1:N/2+1); > psdx = (1/(Fs*N)) * abs(xdft).^2; > psdx(2:end-1) = 2*psdx(2:end-1); > > Remark: Before calculating the PSD, I split the data into parts of equal > size, calculate the PSD for each one, and average over the set of PSDs. > This improves the graphical visualization a lot. > > As the result matches my expectation exactly for 4 different kinds of > noise, I would have assumed that this PSD calculation approach is quite > reasonable. > > As I see the unexpected behavior only with random walk noise, and the > main difference in the calculation is the term A, I would suspect that > it has something to do with it. > > However, I'm a novice in this field, so any hint is very appreciated. > > regards, Wolfgang > > > [1] http://de.mathworks.com/help/signal/ug/psd-estimate-using-fft.html > _______________________________________________ > time-nuts mailing list -- time-nuts@febo.com > To unsubscribe, go to https://www.febo.com/cgi-bin/mailman/listinfo/time-nuts > and follow the instructions there. >

xdft = fft(x);
xdft = xdft(1:N/2+1);
psdx = (1/(Fs*N)) * abs(xdft).^2;
psdx(2:end-1) = 2*psdx(2:end-1);

On 03/06/2015 10:29 PM, Magnus Danielson wrote: > I have checked several sources, and they match up with the IEEE 1139 in > this regard. > > I have also evaluated the equation for Allan variance for the random > walk noise, and it matches up with the references and what I put here: > https://en.wikipedia.org/wiki/Allan_variance#Power-law_noise Thanks a lot for your effort! > So, the A formula you have matches up. > > You will need to find another source of the mismatch. I will take a step back and describe the overall picture of what I'm doing. Maybe someone can help me spot where I do something wrong. (As stated later, the part where I'm quite unsure what I'm doing is the PSD estimation part.) My main goal is to simulate powerlaw noise. I then analyze the generated noise to check if my simulation is reasonable. So the basic workflow would be the following: 1) Generate noise 2) Analyze the noise in the time and frequency domain 3) See that everything agrees and be happy :) Step 1: Noise generation ----------------------------------------------------------- I generate powerlaw noise with the method described by Kasdin and Walter in [1]. So basically I generate white noise and apply a filter as described in [1] to get a PSD shape corresponding to the different values of alpha. The part of the PSD that will have the correct shape depends on the filter length and the simulated sampling frequency. Basically: the length of a simulation I would like to carry out specifies a lower bound on the filter length to get correct results. For WPM, WFM and RW noise I can use a shortcut: for these types of noise the filter coefficients are basically a discrete derivative, an identity filter and a cumulative sum. This is expected, as it agrees with [2], which states that integration of powerlaw noise decreases alpha by 2 (chapter 3.4 in [2]). Thus for even values of alpha I can even skip the expensive convolution to apply the filter and implement the filters directly. As input white noise I use a Gaussian distribution, mainly because that is what is used in the original paper. (I have also found another implementation [3] that optionally provides a uniform distribution). I'm quite confident that the noise generation part works as expected. However, even if I do something wrong here, it should not influence the analyzing part. Step 2: Analyzing noise ----------------------------------------------------------- 2.1 Time domain To analyze powerlaw noise in the time domain, I use a Matlab script called 'allan' [4], which calculates the Allan Deviation. I also found another Matlab tool called 'Stability Analyzer' [5], which can also calculate ADEV values. These two tools are developed by different authors and expect different input formats, but their results agree for any noise example I have tried so far. Thus I would say both of them can be trusted to work as expected. 2.2 Frequency domain IEEE 1139[6] defines S_y as: "frequency spectrum Sy(f): One-sided spectral density of the normalized frequency fluctuations, as defined in normalized frequency fluctuations y(t)." However, I'm not sure how to calculate this measure for a given noise sample. Anything I describe below is just based on 'I think this might work'. If anyone knows a better way of calculating S_y, or tools that can be used for this task, I would be glad to hear about it :) As already stated in the earlier mail I use the method described in [7] to estimate the one-sided PSD of my noise data in FFD format. These plots are quite noisy, and to improve the graphical presentation I use the averaging method described in [8]. I split the noise vector in parts of equal length, calculate the individual PSDs and average over them. Using this averaging method, the PSD plots converge to lines on a log-log plot with the expected slopes. I have an example figure attached to the mail that shows the effect of the averaging (PSD_Average.png). Step 3: Comparing time and frequency domain results ----------------------------------------------------------- At this point I have plots for both the Allan Deviation and the FFD-PSD, and would like to compare them. As first step I estimate h_alpha from the Allan Deviation plot (I'm aware that I need to take care for the Allan Deviation <-> Allan Variance conversion). Then I try to estimate the expected PSD values and compare them with my actual plot using the formulas from IEEE 1139. However, at this point a see that RW noise behaves unexpected :( Numerical Example: ----------------------------------------------------------- Suppose the figure attached as 'Numeric_example.png': At Tau = 0.1s the ADEV plot has a value of 0.005849, so de AVAR would be 3.4211e-05 at this point. The constant A is 2 * pi^2/3 = 6.5797. Thus the value of h_-2 could be roughly estimated as AVAR / (Tau * A) = ~5.2e-05. This would lead to an expected S_y value at a frequency f = 10Hz of h_-2 * f = 5.2000e-07, or -62.84dB The actual plot value is at -59.83, so its ~3dB too larger than expected. regards, Wolfgang [1] Kasdin and Walter, Discrete Simulation of Power Law noise, 1992 [2] Riley, NIST SP 1065: Handbook of Frequency Stability Analysis, 2008 [3] http://people.sc.fsu.edu/~jburkardt%20/c_src/cnoise/cnoise.html [4] http://de.mathworks.com/matlabcentral/fileexchange/13246-allan [5] http://de.mathworks.com/matlabcentral/fileexchange/31319-stability-analyzer-53230a [6] IEEE 1139 [7] http://de.mathworks.com/help/signal/ug/psd-estimate-using-fft.html [8] http://www.dspguide.com/ch9/1.htm > > Cheers, > Magnus > > On 03/06/2015 11:04 AM, Wolfgang Wallner wrote: >> >> >> On 03/05/2015 07:23 PM, Attila Kinali wrote: >>> Servus! >> >> Servus :) >> >>> On Thu, 05 Mar 2015 14:35:51 +0100 >>> Wolfgang Wallner <wolfgang-wallner@gmx.at> wrote: >>> >>>> For the random walk noise the expected line is off by a factor of >>>> exactly 2 from the calculated plot, and I don't know how to explain >>>> this >>>> behavior. >>> >>> I'm probably the wrong one to answer, as I have never done any noise >>> simulation or even read up the relevant papers, but... >>> A factor of 2 sounds like the difference you would get between one sided >>> and two sided noise PSD's. >>> >> >> I calculate the one-sided PSD of the FFD data as described in [1] (first >> paragraph), so the code looks like this: >> >> xdft = fft(x); >> xdft = xdft(1:N/2+1); >> psdx = (1/(Fs*N)) * abs(xdft).^2; >> psdx(2:end-1) = 2*psdx(2:end-1); >> >> Remark: Before calculating the PSD, I split the data into parts of equal >> size, calculate the PSD for each one, and average over the set of PSDs. >> This improves the graphical visualization a lot. >> >> As the result matches my expectation exactly for 4 different kinds of >> noise, I would have assumed that this PSD calculation approach is quite >> reasonable. >> >> As I see the unexpected behavior only with random walk noise, and the >> main difference in the calculation is the term A, I would suspect that >> it has something to do with it. >> >> However, I'm a novice in this field, so any hint is very appreciated. >> >> regards, Wolfgang >> >> >> [1] http://de.mathworks.com/help/signal/ug/psd-estimate-using-fft.html >> _______________________________________________ >> time-nuts mailing list -- time-nuts@febo.com >> To unsubscribe, go to >> https://www.febo.com/cgi-bin/mailman/listinfo/time-nuts >> and follow the instructions there. >> > _______________________________________________ > time-nuts mailing list -- time-nuts@febo.com > To unsubscribe, go to > https://www.febo.com/cgi-bin/mailman/listinfo/time-nuts > and follow the instructions there.

Hi Wolfgang, Have a look at the 20 plots in: http://leapsecond.com/pages/allan/Exploring_Allan_Deviation_v2.pdf This shows phase/frequency/ADEV+MDEV and PSD for 5 noise types. Zoom the PDF 400x if necessary. This was generated with Stable32 and should be 100% correct. See if your results agree. The raw data is at: http://leapsecond.com/pages/allan/ /tvb ----- Original Message ----- From: "Wolfgang Wallner" <wolfgang-wallner@gmx.at> To: <time-nuts@febo.com> Sent: Monday, March 09, 2015 8:57 AM Subject: Re: [time-nuts] AVAR <-> S_Y conversion On 03/06/2015 10:29 PM, Magnus Danielson wrote: > I have checked several sources, and they match up with the IEEE 1139 in > this regard. > > I have also evaluated the equation for Allan variance for the random > walk noise, and it matches up with the references and what I put here: > https://en.wikipedia.org/wiki/Allan_variance#Power-law_noise Thanks a lot for your effort!

 xdft = fft(x);
 xdft = xdft(1:N/2+1);
 psdx = (1/(Fs*N)) * abs(xdft).^2;
 psdx(2:end-1) = 2*psdx(2:end-1);

Wolfgang, Remember to scale the double-integral right: d(t) = D y(t) = integrate(d(t),t) = y_0 + Dt x(t) = integrate(y(t),t) = x_0 + y_0*t + D/2*t^2 Did you miss the 1/2 factor somewhere? That would make sense for the Random Walk phase noise. Cheers, Magnus On 03/09/2015 04:57 PM, Wolfgang Wallner wrote: > > > On 03/06/2015 10:29 PM, Magnus Danielson wrote: >> I have checked several sources, and they match up with the IEEE 1139 in >> this regard. >> >> I have also evaluated the equation for Allan variance for the random >> walk noise, and it matches up with the references and what I put here: >> https://en.wikipedia.org/wiki/Allan_variance#Power-law_noise > > Thanks a lot for your effort! > >> So, the A formula you have matches up. >> >> You will need to find another source of the mismatch. > > I will take a step back and describe the overall picture of what I'm doing. > Maybe someone can help me spot where I do something wrong. > (As stated later, the part where I'm quite unsure what I'm doing is the PSD estimation part.) > > My main goal is to simulate powerlaw noise. I then analyze the generated noise to check if my simulation is reasonable. > > So the basic workflow would be the following: > > 1) Generate noise > 2) Analyze the noise in the time and frequency domain > 3) See that everything agrees and be happy :) > > Step 1: Noise generation > ----------------------------------------------------------- > > I generate powerlaw noise with the method described by Kasdin and Walter in [1]. > So basically I generate white noise and apply a filter as described in [1] to get a PSD shape corresponding to the different values of alpha. > The part of the PSD that will have the correct shape depends on the filter length and the simulated sampling frequency. > Basically: the length of a simulation I would like to carry out specifies a lower bound on the filter length to get correct results. > > For WPM, WFM and RW noise I can use a shortcut: for these types of noise the filter coefficients are basically a discrete derivative, an identity filter and a cumulative sum. > This is expected, as it agrees with [2], which states that integration of powerlaw noise decreases alpha by 2 (chapter 3.4 in [2]). > Thus for even values of alpha I can even skip the expensive convolution to apply the filter and implement the filters directly. > > As input white noise I use a Gaussian distribution, mainly because that is what is used in the original paper. > (I have also found another implementation [3] that optionally provides a uniform distribution). > > I'm quite confident that the noise generation part works as expected. > However, even if I do something wrong here, it should not influence the analyzing part. > > Step 2: Analyzing noise > ----------------------------------------------------------- > > 2.1 Time domain > > To analyze powerlaw noise in the time domain, I use a Matlab script called 'allan' [4], which calculates the Allan Deviation. > I also found another Matlab tool called 'Stability Analyzer' [5], which can also calculate ADEV values. > These two tools are developed by different authors and expect different input formats, but their results agree for any noise example I have tried so far. > Thus I would say both of them can be trusted to work as expected. > > 2.2 Frequency domain > > IEEE 1139[6] defines S_y as: "frequency spectrum Sy(f): One-sided spectral density of the normalized frequency fluctuations, as defined in normalized frequency fluctuations y(t)." > > However, I'm not sure how to calculate this measure for a given noise sample. > Anything I describe below is just based on 'I think this might work'. > > If anyone knows a better way of calculating S_y, or tools that can be used for this task, I would be glad to hear about it :) > > As already stated in the earlier mail I use the method described in [7] to estimate the one-sided PSD of my noise data in FFD format. > These plots are quite noisy, and to improve the graphical presentation I use the averaging method described in [8]. > I split the noise vector in parts of equal length, calculate the individual PSDs and average over them. > Using this averaging method, the PSD plots converge to lines on a log-log plot with the expected slopes. > I have an example figure attached to the mail that shows the effect of the averaging (PSD_Average.png). > > Step 3: Comparing time and frequency domain results > ----------------------------------------------------------- > > At this point I have plots for both the Allan Deviation and the FFD-PSD, and would like to compare them. > As first step I estimate h_alpha from the Allan Deviation plot (I'm aware that I need to take care for the Allan Deviation <-> Allan Variance conversion). > Then I try to estimate the expected PSD values and compare them with my actual plot using the formulas from IEEE 1139. > > However, at this point a see that RW noise behaves unexpected :( > > Numerical Example: > ----------------------------------------------------------- > > Suppose the figure attached as 'Numeric_example.png': > > At Tau = 0.1s the ADEV plot has a value of 0.005849, so de AVAR would be 3.4211e-05 at this point. > The constant A is 2 * pi^2/3 = 6.5797. > > Thus the value of h_-2 could be roughly estimated as AVAR / (Tau * A) = ~5.2e-05. > This would lead to an expected S_y value at a frequency f = 10Hz of > > h_-2 * f = 5.2000e-07, or -62.84dB > > The actual plot value is at -59.83, so its ~3dB too larger than expected. > > regards, Wolfgang > > > [1] Kasdin and Walter, Discrete Simulation of Power Law noise, 1992 > [2] Riley, NIST SP 1065: Handbook of Frequency Stability Analysis, 2008 > [3] http://people.sc.fsu.edu/~jburkardt%20/c_src/cnoise/cnoise.html > [4] http://de.mathworks.com/matlabcentral/fileexchange/13246-allan > [5] http://de.mathworks.com/matlabcentral/fileexchange/31319-stability-analyzer-53230a > [6] IEEE 1139 > [7] http://de.mathworks.com/help/signal/ug/psd-estimate-using-fft.html > [8] http://www.dspguide.com/ch9/1.htm > >> >> Cheers, >> Magnus >> >> On 03/06/2015 11:04 AM, Wolfgang Wallner wrote: >>> >>> >>> On 03/05/2015 07:23 PM, Attila Kinali wrote: >>>> Servus! >>> >>> Servus :) >>> >>>> On Thu, 05 Mar 2015 14:35:51 +0100 >>>> Wolfgang Wallner <wolfgang-wallner@gmx.at> wrote: >>>> >>>>> For the random walk noise the expected line is off by a factor of >>>>> exactly 2 from the calculated plot, and I don't know how to explain >>>>> this >>>>> behavior. >>>> >>>> I'm probably the wrong one to answer, as I have never done any noise >>>> simulation or even read up the relevant papers, but... >>>> A factor of 2 sounds like the difference you would get between one sided >>>> and two sided noise PSD's. >>>> >>> >>> I calculate the one-sided PSD of the FFD data as described in [1] (first >>> paragraph), so the code looks like this: >>> >>> xdft = fft(x); >>> xdft = xdft(1:N/2+1); >>> psdx = (1/(Fs*N)) * abs(xdft).^2; >>> psdx(2:end-1) = 2*psdx(2:end-1); >>> >>> Remark: Before calculating the PSD, I split the data into parts of equal >>> size, calculate the PSD for each one, and average over the set of PSDs. >>> This improves the graphical visualization a lot. >>> >>> As the result matches my expectation exactly for 4 different kinds of >>> noise, I would have assumed that this PSD calculation approach is quite >>> reasonable. >>> >>> As I see the unexpected behavior only with random walk noise, and the >>> main difference in the calculation is the term A, I would suspect that >>> it has something to do with it. >>> >>> However, I'm a novice in this field, so any hint is very appreciated. >>> >>> regards, Wolfgang >>> >>> >>> [1] http://de.mathworks.com/help/signal/ug/psd-estimate-using-fft.html >>> _______________________________________________ >>> time-nuts mailing list -- time-nuts@febo.com >>> To unsubscribe, go to >>> https://www.febo.com/cgi-bin/mailman/listinfo/time-nuts >>> and follow the instructions there. >>> >> _______________________________________________ >> time-nuts mailing list -- time-nuts@febo.com >> To unsubscribe, go to >> https://www.febo.com/cgi-bin/mailman/listinfo/time-nuts >> and follow the instructions there. >> >> >> _______________________________________________ >> time-nuts mailing list -- time-nuts@febo.com >> To unsubscribe, go to https://www.febo.com/cgi-bin/mailman/listinfo/time-nuts >> and follow the instructions there.

 xdft = fft(x);
 xdft = xdft(1:N/2+1);
 psdx = (1/(Fs*N)) * abs(xdft).^2;
 psdx(2:end-1) = 2*psdx(2:end-1);

Hello Magnus, I'm sorry, but I can't follow you here. I know that time deviation values and fractional frequency values are related via integration, so I think I understand the third line. But I don't know what is meant especially with the first one: What is d(t)? What is D? Is this the frequency drift? - Wolfgang On 03/14/2015 11:12 AM, Magnus Danielson wrote: > Wolfgang, > > Remember to scale the double-integral right: > > d(t) = D > y(t) = integrate(d(t),t) = y_0 + Dt > x(t) = integrate(y(t),t) = x_0 + y_0*t + D/2*t^2 > > Did you miss the 1/2 factor somewhere? > That would make sense for the Random Walk phase noise. > > Cheers, > Magnus > > On 03/09/2015 04:57 PM, Wolfgang Wallner wrote: >> >> >> On 03/06/2015 10:29 PM, Magnus Danielson wrote: >>> I have checked several sources, and they match up with the IEEE 1139 in >>> this regard. >>> >>> I have also evaluated the equation for Allan variance for the random >>> walk noise, and it matches up with the references and what I put here: >>> https://en.wikipedia.org/wiki/Allan_variance#Power-law_noise >> >> Thanks a lot for your effort! >> >>> So, the A formula you have matches up. >>> >>> You will need to find another source of the mismatch. >> >> I will take a step back and describe the overall picture of what I'm >> doing. >> Maybe someone can help me spot where I do something wrong. >> (As stated later, the part where I'm quite unsure what I'm doing is >> the PSD estimation part.) >> >> My main goal is to simulate powerlaw noise. I then analyze the >> generated noise to check if my simulation is reasonable. >> >> So the basic workflow would be the following: >> >> 1) Generate noise >> 2) Analyze the noise in the time and frequency domain >> 3) See that everything agrees and be happy :) >> >> Step 1: Noise generation >> ----------------------------------------------------------- >> >> I generate powerlaw noise with the method described by Kasdin and >> Walter in [1]. >> So basically I generate white noise and apply a filter as described in >> [1] to get a PSD shape corresponding to the different values of alpha. >> The part of the PSD that will have the correct shape depends on the >> filter length and the simulated sampling frequency. >> Basically: the length of a simulation I would like to carry out >> specifies a lower bound on the filter length to get correct results. >> >> For WPM, WFM and RW noise I can use a shortcut: for these types of >> noise the filter coefficients are basically a discrete derivative, an >> identity filter and a cumulative sum. >> This is expected, as it agrees with [2], which states that integration >> of powerlaw noise decreases alpha by 2 (chapter 3.4 in [2]). >> Thus for even values of alpha I can even skip the expensive >> convolution to apply the filter and implement the filters directly. >> >> As input white noise I use a Gaussian distribution, mainly because >> that is what is used in the original paper. >> (I have also found another implementation [3] that optionally provides >> a uniform distribution). >> >> I'm quite confident that the noise generation part works as expected. >> However, even if I do something wrong here, it should not influence >> the analyzing part. >> >> Step 2: Analyzing noise >> ----------------------------------------------------------- >> >> 2.1 Time domain >> >> To analyze powerlaw noise in the time domain, I use a Matlab script >> called 'allan' [4], which calculates the Allan Deviation. >> I also found another Matlab tool called 'Stability Analyzer' [5], >> which can also calculate ADEV values. >> These two tools are developed by different authors and expect >> different input formats, but their results agree for any noise example >> I have tried so far. >> Thus I would say both of them can be trusted to work as expected. >> >> 2.2 Frequency domain >> >> IEEE 1139[6] defines S_y as: "frequency spectrum Sy(f): One-sided >> spectral density of the normalized frequency fluctuations, as defined >> in normalized frequency fluctuations y(t)." >> >> However, I'm not sure how to calculate this measure for a given noise >> sample. >> Anything I describe below is just based on 'I think this might work'. >> >> If anyone knows a better way of calculating S_y, or tools that can be >> used for this task, I would be glad to hear about it :) >> >> As already stated in the earlier mail I use the method described in >> [7] to estimate the one-sided PSD of my noise data in FFD format. >> These plots are quite noisy, and to improve the graphical presentation >> I use the averaging method described in [8]. >> I split the noise vector in parts of equal length, calculate the >> individual PSDs and average over them. >> Using this averaging method, the PSD plots converge to lines on a >> log-log plot with the expected slopes. >> I have an example figure attached to the mail that shows the effect of >> the averaging (PSD_Average.png). >> >> Step 3: Comparing time and frequency domain results >> ----------------------------------------------------------- >> >> At this point I have plots for both the Allan Deviation and the >> FFD-PSD, and would like to compare them. >> As first step I estimate h_alpha from the Allan Deviation plot (I'm >> aware that I need to take care for the Allan Deviation <-> Allan >> Variance conversion). >> Then I try to estimate the expected PSD values and compare them with >> my actual plot using the formulas from IEEE 1139. >> >> However, at this point a see that RW noise behaves unexpected :( >> >> Numerical Example: >> ----------------------------------------------------------- >> >> Suppose the figure attached as 'Numeric_example.png': >> >> At Tau = 0.1s the ADEV plot has a value of 0.005849, so de AVAR would >> be 3.4211e-05 at this point. >> The constant A is 2 * pi^2/3 = 6.5797. >> >> Thus the value of h_-2 could be roughly estimated as AVAR / (Tau * A) >> = ~5.2e-05. >> This would lead to an expected S_y value at a frequency f = 10Hz of >> >> h_-2 * f = 5.2000e-07, or -62.84dB >> >> The actual plot value is at -59.83, so its ~3dB too larger than expected. >> >> regards, Wolfgang >> >> >> [1] Kasdin and Walter, Discrete Simulation of Power Law noise, 1992 >> [2] Riley, NIST SP 1065: Handbook of Frequency Stability Analysis, 2008 >> [3] http://people.sc.fsu.edu/~jburkardt%20/c_src/cnoise/cnoise.html >> [4] http://de.mathworks.com/matlabcentral/fileexchange/13246-allan >> [5] >> http://de.mathworks.com/matlabcentral/fileexchange/31319-stability-analyzer-53230a >> >> [6] IEEE 1139 >> [7] http://de.mathworks.com/help/signal/ug/psd-estimate-using-fft.html >> [8] http://www.dspguide.com/ch9/1.htm >> >>> >>> Cheers, >>> Magnus >>> >>> On 03/06/2015 11:04 AM, Wolfgang Wallner wrote: >>>> >>>> >>>> On 03/05/2015 07:23 PM, Attila Kinali wrote: >>>>> Servus! >>>> >>>> Servus :) >>>> >>>>> On Thu, 05 Mar 2015 14:35:51 +0100 >>>>> Wolfgang Wallner <wolfgang-wallner@gmx.at> wrote: >>>>> >>>>>> For the random walk noise the expected line is off by a factor of >>>>>> exactly 2 from the calculated plot, and I don't know how to explain >>>>>> this >>>>>> behavior. >>>>> >>>>> I'm probably the wrong one to answer, as I have never done any noise >>>>> simulation or even read up the relevant papers, but... >>>>> A factor of 2 sounds like the difference you would get between one >>>>> sided >>>>> and two sided noise PSD's. >>>>> >>>> >>>> I calculate the one-sided PSD of the FFD data as described in [1] >>>> (first >>>> paragraph), so the code looks like this: >>>> >>>> xdft = fft(x); >>>> xdft = xdft(1:N/2+1); >>>> psdx = (1/(Fs*N)) * abs(xdft).^2; >>>> psdx(2:end-1) = 2*psdx(2:end-1); >>>> >>>> Remark: Before calculating the PSD, I split the data into parts of >>>> equal >>>> size, calculate the PSD for each one, and average over the set of PSDs. >>>> This improves the graphical visualization a lot. >>>> >>>> As the result matches my expectation exactly for 4 different kinds of >>>> noise, I would have assumed that this PSD calculation approach is quite >>>> reasonable. >>>> >>>> As I see the unexpected behavior only with random walk noise, and the >>>> main difference in the calculation is the term A, I would suspect that >>>> it has something to do with it. >>>> >>>> However, I'm a novice in this field, so any hint is very appreciated. >>>> >>>> regards, Wolfgang >>>> >>>> >>>> [1] http://de.mathworks.com/help/signal/ug/psd-estimate-using-fft.html >>>> _______________________________________________ >>>> time-nuts mailing list -- time-nuts@febo.com >>>> To unsubscribe, go to >>>> https://www.febo.com/cgi-bin/mailman/listinfo/time-nuts >>>> and follow the instructions there. >>>> >>> _______________________________________________ >>> time-nuts mailing list -- time-nuts@febo.com >>> To unsubscribe, go to >>> https://www.febo.com/cgi-bin/mailman/listinfo/time-nuts >>> and follow the instructions there. >>> >>> >>> _______________________________________________ >>> time-nuts mailing list -- time-nuts@febo.com >>> To unsubscribe, go to >>> https://www.febo.com/cgi-bin/mailman/listinfo/time-nuts >>> and follow the instructions there. > _______________________________________________ > time-nuts mailing list -- time-nuts@febo.com > To unsubscribe, go to > https://www.febo.com/cgi-bin/mailman/listinfo/time-nuts > and follow the instructions there.