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
_______________________________________________ time-nuts mailing list -- [email protected] To unsubscribe, go to https://www.febo.com/cgi-bin/mailman/listinfo/time-nuts and follow the instructions there.
