/
deltanu_estimators.py
215 lines (170 loc) · 9.72 KB
/
deltanu_estimators.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
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
"""Helper functions for estimating deltanu from periodograms."""
from __future__ import division, print_function
import numpy as np
from scipy.signal import find_peaks
from matplotlib import pyplot as plt
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
from astropy import units as u
from .utils import SeismologyQuantity
from . import utils
from .. import MPLSTYLE
__all__ = ['estimate_deltanu_acf2d', 'diagnose_deltanu_acf2d']
def estimate_deltanu_acf2d(periodogram, numax):
"""Returns the average value of the large frequency spacing, DeltaNu,
of the seismic oscillations of the target, using an autocorrelation
function.
There are many papers on the topic of autocorrelation functions for
estimating seismic parameters, including but not limited to:
Roxburgh & Vorontsov (2006), Roxburgh (2009), Mosser & Appourchaux (2009),
Huber et al. (2009), Verner & Roxburgh (2011) & Viani et al. (2019).
We base this approach first and foremost off the approach taken in
Mosser & Appourchaux (2009). Given a known numax, a window around this
numax is taken of one estimated full-width-half-maximum (FWHM) of the
seismic mode envelope either side of numax. This width is chosen so that
the autocorrelation includes all of the visible mode peaks.
The autocorrelation (numpy.correlate) is given as::
C = sum(s * s)
where s is a window of the signal-to-noise spectrum. When shifting
the spectrum over itself, C will increase when two mode peaks are
overlapping. Because of the method of this calculation, we need to first
rescale the power by subtracting its mean, placing its mean around 0. This
decreases the noise levels in the ACF, as the autocorrelation of the noise
with itself will be close to zero.
As is done in Mosser & Appourchaux, we rescale the value of C in terms
of the noise level in the ACF spectrum as::
A = (|C^2| / |C[0]^2|) * (2 * len(C) / 3) .
The method will autocorrelate the region around the estimated numax
expected to contain seismic oscillation modes. Repeating peaks in the
autocorrelation implies an evenly spaced structure of modes.
The peak closest to an empirical estimate of deltanu is taken as the true
value. The peak finding algorithm is limited by a minimum spacing
between peaks of 0.5 times the empirical value for deltanu.
Our empirical estimate for numax is taken from Stello et al. (2009) as::
deltanu = 0.294 * numax^0.772
If `numax` is None, a numax is calculated using the estimate_numax()
function with default settings.
NOTE: This function is intended for use with solar like Main Sequence
and Red Giant Branch oscillators only.
Parameters
----------
numax : float
An estimated numax value of the mode envelope in the periodogram. If
not given units it is assumed to be in units of the periodogram
frequency attribute.
Returns
-------
deltanu : `.SeismologyQuantity`
The average large frequency spacing of the seismic oscillation modes.
In units of the periodogram frequency attribute.
"""
# Run some checks on the passed in numaxs
# Ensure input numax is in the correct units
numax = u.Quantity(numax, periodogram.frequency.unit)
fs = np.median(np.diff(periodogram.frequency.value))
if numax.value < fs:
raise ValueError("The input numax can not be lower than"
" a single frequency bin.")
if numax.value > np.nanmax(periodogram.frequency.value):
raise ValueError("The input numax can not be higher than"
"the highest frequency value in the periodogram.")
# Calculate deltanu using the method by Stello et al. 2009
# Make sure that this relation only ever happens in microhertz space
deltanu_emp = u.Quantity((0.294 * u.Quantity(numax, u.microhertz).value ** 0.772)*u.microhertz,
periodogram.frequency.unit).value
window_width = 2*int(np.floor(utils.get_fwhm(periodogram, numax.value)))
aacf = utils.autocorrelate(periodogram, numax=numax.value, window_width=window_width)
acf = (np.abs(aacf**2)/np.abs(aacf[0]**2)) / (3/(2*len(aacf)))
fs = np.median(np.diff(periodogram.frequency.value))
lags = np.linspace(0., len(acf)*fs, len(acf))
#Select a 25% region region around the empirical deltanu
sel = (lags > deltanu_emp - .25*deltanu_emp) & (lags < deltanu_emp + .25*deltanu_emp)
#Run a peak finder on this region
peaks, _ = find_peaks(acf[sel], distance=np.floor(deltanu_emp/2. / fs))
#Select the peak closest to the empirical value
best_deltanu_value = lags[sel][peaks][np.argmin(np.abs(lags[sel][peaks] - deltanu_emp))]
best_deltanu = u.Quantity(best_deltanu_value, periodogram.frequency.unit)
diagnostics = {'lags':lags, 'acf':acf, 'peaks':peaks, 'sel':sel,
'numax':numax, 'deltanu_emp':deltanu_emp}
result = SeismologyQuantity(best_deltanu,
name="deltanu",
method="ACF2D",
diagnostics=diagnostics,
diagnostics_plot_method=diagnose_deltanu_acf2d)
return result
def diagnose_deltanu_acf2d(deltanu, periodogram):
"""Returns a diagnostic plot which elucidates how deltanu was estimated.
[1] Scaled correlation metric vs frequecy lag of the autocorrelation
window, with inset close up on the determined deltanu and a line
indicating the determined deltanu.
For details on the deltanu estimation, see the `estimate_deltanu()`
function. The calculation performed is identical.
NOTE: When plotting, we exclude the first two frequency lag bins, to
make the relevant features on the plot clearer, as these bins are close to
the spectrum correlated with itself and therefore much higher than the rest
of the bins.
Parameters
----------
deltanu : `.SeismologyResult` object
The object returned by `estimate_deltanu_acf2d()`.
Returns
-------
ax : `~matplotlib.axes.Axes`
The matplotlib axes object.
"""
with plt.style.context(MPLSTYLE):
fig, axs = plt.subplots(2, figsize=(8.485, 8))
ax = axs[0]
periodogram.plot(ax=ax, label='')
ax.axvline(deltanu.diagnostics['numax'].value, c='r', linewidth=1,
alpha=.4, ls=':')
ax.text(deltanu.diagnostics['numax'].value, periodogram.power.value.max()*0.45,
'{} ({:.1f} {})'.format(r'$\nu_{\rm max}$', deltanu.diagnostics['numax'].value, deltanu.diagnostics['numax'].unit.to_string('latex')),
rotation=90, ha='right', color='r', alpha=0.5, fontsize=8)
ax.text(.025, .9, 'Input Power Spectrum', horizontalalignment='left',
transform=ax.transAxes, fontsize=11)
window_width = 2*int(np.floor(utils.get_fwhm(periodogram, deltanu.diagnostics['numax'].value)))
frequency_spacing = np.median(np.diff(periodogram.frequency.value))
spread = int(window_width/2/frequency_spacing) # spread in indices
a = np.argmin(np.abs(periodogram.frequency.value - deltanu.diagnostics['numax'].value)) + spread
b = np.argmin(np.abs(periodogram.frequency.value - deltanu.diagnostics['numax'].value)) - spread
a = [periodogram.frequency.value[a] if a < len(periodogram.frequency)
else periodogram.frequency.value[-1]][0]
b = [periodogram.frequency.value[b] if b > 0
else periodogram.frequency.value[0]][0]
ax.axvline(a, c='r', linewidth=2, alpha=.4, ls='--')
ax.axvline(b, c='r', linewidth=2, alpha=.4, ls='--')
h = periodogram.power.value.max() * 0.9
ax.annotate("", xy=(a, h), xytext=(a - (a-b), h),
arrowprops=dict(arrowstyle="<->", color='r', alpha=0.5), va='bottom')
ax.text(a - (a-b)/2, h, r"2 $\times$ FWHM", color='r', alpha=0.7, fontsize=10, va='bottom', ha='center')
ax.set_xlim(b - ((a-b)*0.2), a + ((a-b)*0.2))
ax = axs[1]
ax.plot(deltanu.diagnostics['lags'][2:], deltanu.diagnostics['acf'][2:])
ax.set_xlabel("Frequency Lag [{}]".format(periodogram.frequency.unit.to_string('latex')))
ax.set_ylabel(r'Scaled Auto Correlation', fontsize=11)
axin = inset_axes(ax, width="50%",height="50%", loc="upper right")
axin.set_yticks([])
axin.plot(deltanu.diagnostics['lags'][deltanu.diagnostics['sel']],
deltanu.diagnostics['acf'][deltanu.diagnostics['sel']])
axin.scatter(deltanu.diagnostics['lags'][deltanu.diagnostics['sel']][deltanu.diagnostics['peaks']],
deltanu.diagnostics['acf'][deltanu.diagnostics['sel']][deltanu.diagnostics['peaks']],
c='r', s=5)
mea_label = r'Measured {} {:.1f} {}'.format(
r'$\Delta\nu$',
deltanu.value,
periodogram.frequency.unit.to_string('latex'))
ax.axvline(deltanu.value,c='r', linewidth=2, alpha=.4, label=mea_label)
emp_label = r'Empirical {} {:.1f} {}'.format(
r'$\Delta\nu$',
deltanu.diagnostics['deltanu_emp'],
periodogram.frequency.unit.to_string('latex'))
ax.axvline(deltanu.diagnostics['deltanu_emp'], c='b', linewidth=2,
alpha=.4, ls='--', label=emp_label)
axin.axvline(deltanu.value,c='r', linewidth=2, alpha=.4)
axin.axvline(deltanu.diagnostics['deltanu_emp'], c='b', linewidth=2,
alpha=.4, ls='--')
ax.text(.025, .9, 'Scaled Auto Correlation Within 2 FWHM',
horizontalalignment='left',
transform=ax.transAxes, fontsize=11)
ax.legend(loc='lower right', fontsize=10)
return ax