forked from jleute/colorednoise
-
Notifications
You must be signed in to change notification settings - Fork 1
/
colorednoise.py
140 lines (113 loc) · 4.33 KB
/
colorednoise.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
'''
Generate Dicrete Colored Noise
python / numpy implementation of
Kasdin, N.J., Walter, T., "Discrete simulation of power law noise
[for oscillator stability evaluation],"
Frequency Control Symposium, 1992. 46th., Proceedings of the 1992 IEEE,
pp.274,283, 27-29 May 1992
http://dx.doi.org/10.1109/FREQ.1992.270003
'''
# author: Julia Leute
import numpy as np
def noiseGen(nr, Qd, b):
''' Generates discrete colored noise
required inputs:
nr (number of points, must be power of two)
Qd (discrete variance)
b (slope of the noise) '''
mhb = -b/2.0
Qd = np.sqrt(Qd)
# fill sequence wfb with white noise
wfb = np.zeros(nr*2)
wfb[:nr] = np.random.normal(0, Qd, nr)
# generate the coefficients hfb
hfb = np.zeros(nr*2)
hfb[0] = 1.0
indices = np.arange(nr-1)
hfb[1:nr] = (mhb+indices)/(indices+1.0)
hfb[:nr] = np.multiply.accumulate(hfb[:nr])
# discrete Fourier transforms of white noise and coefficients,
# multiplication of resulting complex vectors,
# inverse Fourier transform
colorednoise = np.fft.irfft(np.fft.rfft(wfb)*np.fft.rfft(hfb))[:nr]
return colorednoise
def phase_psd_from_qd(qd, b, tau0):
""" return phase power spectral density coefficient g_b from QD
Colored noise generated with (qd, b, tau0) parameters will
show a phase power spectral density of
S_x(f) = Phase_PSD(f) = g_b * f^b
Kasdin & Walter eqn (39)
"""
g_b = qd*2.0*pow(2.0*np.pi, b)*pow(tau0, b+1.0)
return g_b
def frequency_psd_from_qd(qd, b, tau0):
""" return frequency power spectral density coefficient h_a from QD
Colored noise generated with (qd, b, tau0) parameters will
show a frequency power spectral density of
S_y(f) = Frequency_PSD(f) = h_a * f^a
where the slope a comes from the phase PSD slope b:
a = b + 2
Kasdin & Walter eqn (39)
"""
a = b + 2.0
h_a = qd*2.0*pow(2.0*np.pi, a)*pow(tau0, a-1.0)
return h_a
def adev_from_qd(qd, b, tau0, tau):
""" prefactor for Allan deviation from QD and slope
Colored noise generated with (qd, b, tau0) parameters will
show an Allan variance of:
AVAR = prefactor * h_a * tau^c
where a = b+2 is the slope of the frequency PSD.
and h_a is the frequency PSD prefactor S_y(f) = h_a * f^a
The relation between a, b, c is:
a b c(AVAR) c(MVAR)
-----------------------
-2 -4 1 1
-1 -3 0 0
0 -2 -1 -1
+1 -1 -2 -2
+2 -2 -2 -3
"""
g_b = phase_psd_from_qd(qd,b,tau0)
f_h = 0.5/tau0
gamma = 0.577
if b == 0:
coeff = 3.0*f_h / (4.0*pow(np.pi,2)) # E, White PM, tau^-1
elif b == -1:
coeff = (np.log(2.0*np.pi*f_h*tau)+gamma)/(2.0*pow(np.pi,2))# D, Flicker PM, tau^-1
elif b == -2:
coeff = 0.5 # C, white FM, 1/sqrt(tau)
elif b == -3:
coeff = 2*np.log(2) # B, flicker FM, constant ADEV
elif b == -4:
coeff = 2.0*pow(np.pi,2)/3.0 # A, RW FM, sqrt(tau)
return np.sqrt(coeff*g_b*pow(2.0*np.pi,2))
def mdev_from_qd(qd, b, tau0, tau):
""" prefactor for Modified Allan deviation from QD and slope
Colored noise generated with (qd, b, tau0) parameters will
show an Modified Allan variance of:
MVAR = prefactor * h_a * tau^c
where a = b+2 is the slope of the frequency PSD.
and h_a is the frequency PSD prefactor S_y(f) = h_a * f^a
The relation between a, b, c is:
a b c(AVAR) c(MVAR)
-----------------------
-2 -4 1 1
-1 -3 0 0
0 -2 -1 -1
+1 -1 -2 -2
+2 -2 -2 -3
"""
g_b = phase_psd_from_qd(qd,b,tau0)
f_h = 0.5/tau0
if b == 0:
coeff = 3.0/(8.0*pow(np.pi,2)) # E, White PM, tau^-{3/2}
elif b == -1:
coeff = (24.0*np.log(2)-9.0*np.log(3))/8.0/pow(np.pi,2) # D, Flicker PM, tau^-1
elif b == -2:
coeff = 0.25 # C, white FM, 1/sqrt(tau)
elif b == -3:
coeff = 27.0/20.0*np.log(2) # B, flicker FM, constant MDEV
elif b == -4:
coeff = 11.0/20.0*pow(np.pi,2) # A, RW FM, sqrt(tau)
return np.sqrt(coeff*g_b*pow(2.0*np.pi,2))