/usr/lib/python3/dist-packages/astroML/time_series/_periodogram.py is in python3-astroml 0.3-6.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
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 | import numpy as np
def lomb_scargle(t, y, dy, omega, generalized=True,
subtract_mean=True, significance=None):
"""
(Generalized) Lomb-Scargle Periodogram with Floating Mean
Parameters
----------
t : array_like
sequence of times
y : array_like
sequence of observations
dy : array_like
sequence of observational errors
omega : array_like
frequencies at which to evaluate p(omega)
generalized : bool
if True (default) use generalized lomb-scargle method
otherwise, use classic lomb-scargle.
subtract_mean : bool
if True (default) subtract the sample mean from the data before
computing the periodogram. Only referenced if generalized is False
significance : None or float or ndarray
if specified, then this is a list of significances to compute
for the results.
Returns
-------
p : array_like
Lomb-Scargle power associated with each frequency omega
z : array_like
if significance is specified, this gives the levels corresponding
to the desired significance (using the Scargle 1982 formalism)
Notes
-----
The algorithm is based on reference [1]_. The result for generalized=False
is given by equation 4 of this work, while the result for generalized=True
is given by equation 20.
Note that the normalization used in this reference is different from that
used in other places in the literature (e.g. [2]_). For a discussion of
normalization and false-alarm probability, see [1]_.
To recover the normalization used in Scargle [3]_, the results should
be multiplied by (N - 1) / 2 where N is the number of data points.
References
----------
.. [1] M. Zechmeister and M. Kurster, A&A 496, 577-584 (2009)
.. [2] W. Press et al, Numerical Recipies in C (2002)
.. [3] Scargle, J.D. 1982, ApJ 263:835-853
"""
t = np.asarray(t)
y = np.asarray(y)
dy = np.asarray(dy) * np.ones_like(y)
assert t.ndim == 1
assert y.ndim == 1
assert dy.ndim == 1
assert t.shape == y.shape
assert y.shape == dy.shape
w = 1. / dy / dy
w /= w.sum()
# the generalized method takes care of offset automatically,
# while the classic method requires centered data.
if (not generalized) and subtract_mean:
# subtract MLE for mean in the presence of noise.
y = y - np.dot(w, y)
omega = np.asarray(omega)
shape = omega.shape
omega = omega.ravel()[np.newaxis, :]
t = t[:, np.newaxis]
y = y[:, np.newaxis]
dy = dy[:, np.newaxis]
w = w[:, np.newaxis]
sin_omega_t = np.sin(omega * t)
cos_omega_t = np.cos(omega * t)
# compute time-shift tau
# S2 = np.dot(w.T, np.sin(2 * omega * t)
S2 = 2 * np.dot(w.T, sin_omega_t * cos_omega_t)
# C2 = np.dot(w.T, np.cos(2 * omega * t)
C2 = 2 * np.dot(w.T, 0.5 - sin_omega_t ** 2)
if generalized:
S = np.dot(w.T, sin_omega_t)
C = np.dot(w.T, cos_omega_t)
S2 -= (2 * S * C)
C2 -= (C * C - S * S)
tan_2omega_tau = S2 / C2
tau = np.arctan(tan_2omega_tau)
tau *= 0.5
tau /= omega
# compute components needed for the fit
omega_t_tau = omega * (t - tau)
sin_omega_t_tau = np.sin(omega_t_tau)
cos_omega_t_tau = np.cos(omega_t_tau)
Y = np.dot(w.T, y)
YY = np.dot(w.T, y * y) - Y * Y
wy = w * y
YCtau = np.dot(wy.T, cos_omega_t_tau)
YStau = np.dot(wy.T, sin_omega_t_tau)
CCtau = np.dot(w.T, cos_omega_t_tau * cos_omega_t_tau)
SStau = np.dot(w.T, sin_omega_t_tau * sin_omega_t_tau)
if generalized:
Ctau = np.dot(w.T, cos_omega_t_tau)
Stau = np.dot(w.T, sin_omega_t_tau)
YCtau -= Y * Ctau
YStau -= Y * Stau
CCtau -= Ctau * Ctau
SStau -= Stau * Stau
p_omega = (YCtau * YCtau / CCtau + YStau * YStau / SStau) / YY
p_omega = p_omega.reshape(shape)
if significance is not None:
N = t.size
M = 2 * N
z = (-2.0 / (N - 1.)
* np.log(1 - (1 - np.asarray(significance)) ** (1. / M)))
return p_omega, z
else:
return p_omega
|