It is based on the development in Kerr Section 5 equations (6) and (8), using code from the godot implementation, specifically the function power_spectrum_fft which follows.
x = np.linspace(0, 4)
y = np.sin(x*2*np.pi)
dfp = PeakFinder(x,y)
plt.figure(figsize=(5,2))
plt.title('PeakFinder demo')
plt.plot(x,y, '.')
plt.plot(dfp.x,dfp.y, 'o');
dfp.find_closest(2.5)
dfp.query('2<x<3')
@ipynb_doc
def sinc_demo():
"""
### Sinc class demo
{fig}
"""
fig, ax = plt.subplots(figsize=(4,2))
ax.set(xlabel='Frequency ($s^{-1}$)', ylabel='power')
sinc = Sinc(1, 0.5, 0.1)
sinc.plot(ax=ax, pticks=[1,1.5,2,3,4], label='sinc', toplabel='Period (s)')
ax.legend()
return locals()
sinc_demo()
cd = None # the TimeSeries object for manipulation below
@ipynb_doc
def sim_test(amp=1, freq=5, T=50, srate=24, rng=None):
r"""
## Simulation Test
{setup_txt}
### Simulation/FFT details:
* Duration $T={T}\ d$,
* Sample rate ${srate}\ d^{-1}$
* Signal: $S(t) = 10^{-6}\ \big[1 + A \cos(2\pi f_S t)\big]$, where $f_S$ is the source
frequency {freq:.1f}, and the amplitude $A={amp}$.
* Background: fixed at $10^{-6}$.
* FFT: {nsamp} samples zero-padded to {nsize}, factor {resf:.1f}
### The periodogram
{fig1}
### Peak detection
Peak finding finds {nump} peaks.
The nearest peak to the expected signal:
<br>{signal_text}
{fig2}
"""
global cd
from wtlike.simulation import Simulation
def periodic(t):
return 1e-6 * ( 1+ amp*np.cos(2*np.pi*freq*t) )
tsamp = 1/srate
with capture_hide('Output during setup') as setup_txt:
sim = Simulation('periodic', src_flux=periodic, tstart=0, tstop=T, rng=rng)
sim.run()
cd = TimeSeries(sim, tsamp=tsamp, )
sim_info = cd.__repr__()
delf = 1/T # window resolution
fnyquist = 0.25/cd.tsamp
nsamp = int(T/cd.tsamp)
nsize = 4*len(cd.power_df)
resf = delf/cd.power_df.f[0]
def left(ax):
cd.power_plot(ax = ax, title='Kerr periodogram'); #fig, ax = plt.subplots(figsize=(8,4))
def right(ax):
ax.plot(cd.power_df.f, cd.power_df.p0, '.-');
ax.grid(True)
ax.axhline(0, color='grey')
xlim = (freq-5*delf,freq+5*delf)
ax.set(xlim =xlim,
xlabel='Frequency', ylabel='$P_0$',title='peak with sinc');
A, m,b = cd.power_df.p0.max(), freq, delf
sincsq = lambda x: A*np.sinc((x-m)/b)**2
x = np.linspace(*xlim)
ax.plot(x, sincsq(x), '--', lw=2 );
hh = np.array([freq-delf/2, freq+delf/2])
ax.plot(hh, sincsq(hh),'o-r', lw=2);
fig1 = figure(plt.figure(figsize=(12,3)),
caption=f"""Left plot: The full Kerr periodogram showing $P_B$ and $P_1$ up to the Nyquist frequency {fnyquist}.
Right plot: The $P_0$ spectral measure overlaid with the sinc function at the frequency $f_S$
adjusted to match the peak and with its expected width, $1/T={delf}$. (Note that the FFT processing uses
zero padding by at least a factor of 5 to improve the FFT resolution.)""",
width = 800)
plt.subplots_adjust(wspace=0)
gs = plt.GridSpec(1,10)
left( fig1.add_subplot(gs[:6]) )
right( fig1.add_subplot(gs[7:]) )
# find the peak
peak_df = cd.find_peaks()
nump = len(peak_df)
k = peak_df.find_closest(freq)
signal = peak_df.iloc[k] #peak_df.query(f'{freq-1/T}<x<{freq+1/T}')
signal_text = monospace(str(signal))
def peak_plot(fig):
y = cd.power_df.p0.values
x = cd.power_df.f.values
xp, yp = peak_df.iloc[:,0], peak_df.iloc[:,1]
ax = fig.add_subplot()
hkw = dict(bins=np.linspace(0,40,41), histtype='stepfilled', log=True)
ax.hist(yp, color='cornflowerblue', **hkw);
ax.hist(signal.p0, label='signal', color='orange', **hkw);
ax.set(xlabel='Power at peak')
ax.legend()
ax.grid(True)
fig2 = figure(plt.figure(figsize=(5,3)),
caption='Distribution of peak power values')
peak_plot(fig2)
return locals()
# need to bix
sim_test(amp=0.1, fignum=1, rng=24)
Simulation Test
Output during setup
generated 25902 photons
Simulation/FFT details:
- Duration $T=50\ d$,
- Sample rate $24\ d^{-1}$
- Signal: $S(t) = 10^{-6}\ \big[1 + A \cos(2\pi f_S t)\big]$, where $f_S$ is the source frequency 5.0, and the amplitude $A=0.1$.
- Background: fixed at $10^{-6}$.
- FFT: 1200 samples zero-padded to 8192, factor 6.8
The periodogram
Peak detection
Peak finding finds 189 peaks.
The nearest peak to the expected signal:
f 4.998702
p0 37.552967
Name: 158, dtype: float32
@ipynb_doc
def sim_phase():
"""
## Look at the phase
The plot below is of the differences of the phase measured in sequential frequencies, using the analysis above.
Subtracted from that is a constant offset, {offset:.2f} radians. This is $\pi$ divided by the FFT zero-padding factor {fftf:.2f}.
{fig1}
Here is a histogram of the phase differences.
{fig2}
Note that the periodic signal is not evident in the phase information--no visible change at its frequency.
"""
self = cd
fftf = cd.fft_factor
offset = np.pi/cd.fft_factor
dfx = cd.phase_df(offset=offset)
fig2, ax = plt.subplots(figsize=(4,3))
ax.hist(dfx.delphi, 100, log=False, histtype='step');
ax.grid();
fig1, ax1 = plt.subplots(figsize=(15,3))
ax1.plot(dfx.f, dfx.delphi, '.-');
ax1.set(xlabel='Frequency', ylabel='Phase change')
ax1.grid();#plt.xlim(4.8, 5.2);
return locals()
sim_phase()
Look at the phase
The plot below is of the differences of the phase measured in sequential frequencies, using the analysis above. Subtracted from that is a constant offset, 0.46 radians. This is $\pi$ divided by the FFT zero-padding factor 6.83.
Here is a histogram of the phase differences.
Note that the periodic signal is not evident in the phase information--no visible change at its frequency.
tms_list=[]
@ipynb_doc
def source_ffts(source_list, fs_list, xlim=None, xscale='log'):
"""
## Output spectra using `wtlike.TimeSeries`
{txt}
This analysis uses a time bin size of {tbin_h:.2f} h for a Nyquist frequency of {nyquist:.0f} $d^{-1}$
and a total interval of {tspan:.0f} d.
{fig1}
{txt2}
This set of plots show the results of analyses of the three sources that Matthew examined, plus Vela.
Like Kerr's Figure 6, the upper blue plot for each source is the profile likelihood estimator $P_1$,
and the bottom orange plot the background-fixed estimator $P_b$.
### Exposure power spectrum
It is also possible to extract only the exposure spectrum
{fig2x}
"""
global tms_list
plt.style.use('seaborn-ticks')
tbinsize=1/96; tbin_h = tbinsize*24
with capture_hide(f'Source setup output') as txt:
tms_list = [TimeSeries(name, tsamp=tbinsize) for name in source_list]
nyquist = tms_list[0].f_Nyquist
tspan = tms_list[0].tspan
N = len(source_list)
fig, axx = plt.subplots(N,1, figsize=(15,3*N+1), sharex=True,)
plt.subplots_adjust(hspace=0.05)
fig1 = figure(fig, width=600)
for name, ax, fs,tms in zip(source_list, axx.flatten(), fs_list, tms_list):
tms.power_spectrum(minf=1)
tms.power_plot(ax=ax, xlim=xlim or (2e-3,nyquist),
pmax=None, xscale=xscale, fs=fs, name=name)
fs = np.unique( np.array(sorted([1/53.05,1+1/365.25,15.1, 1/365.25, 4/365.25, 1/53.05, 4/365.25, 0.256,0.512, 5.01,15.1])))
df = pd.DataFrame.from_dict(dict(frequency=fs.round(3), period=(1/fs).round(2)))
df.loc[:,'comment']=['1/year','4/year', 'Fermi precession', 'LS 5039 orbital period',
'LS 5039 harmonic', 'sidereal day', 'Cyg X-3 orbital period', 'Fermi orbital frequency']
with capture_show('The tagged frequencies') as txt2:
print(df)
#--------------------
fig2, axx2 = plt.subplots(N,1, figsize=(15,2*N+1), sharex=True)
plt.subplots_adjust(hspace=0.05)
fig2x = figure(fig2, width=600)
for name, ax, fs, tms in zip(source_list, axx2.flatten(), fs_list, tms_list):
q = tms.power_spectrum(exp_only=True)
x,y = q.f,q.p_exp
ax.plot(x,y, '-')
pmax = y.max()*1.1
ap = dict(arrowstyle='->',color='k', lw=3)
for f in fs:
ax.annotate('', xy=(f, 0.85*pmax), xytext=(f, pmax),# transform=ax.transData,
arrowprops=ap);
if name is not None:
ax.text(0.02,0.96, name, va='top', transform=ax.transAxes)
ax.set(xscale='log', xlim=(2e-3,tms.f_Nyquist), ylim=(0,None),
xlabel=r'$\mathrm{Frequency\ (cycles\ d^{-1})}$',
ylabel='Exposure power')
ax.xaxis.set_major_formatter(
lambda val,pos: { 0.01:'0.01', 0.1:'0.1', 1.0:'1', 10.0:'10',}.get(val,''))
ax.grid(True)
return locals()
if Config().valid:
source_ffts(['Vela pulsar','Geminga','LS 5039', 'Cygnus X-3'],
[[1/53.05,1,2,15.1],[1/365.25, 4/365.25, 15.1],[1/53.05, 4/365.25, 0.256,0.512, 15.1],[5.01,15.1]])
Output spectra using wtlike.TimeSeries
Source setup output
SourceData: PSR J0835-4510: Restoring from cache with key "PSR J0835-4510_data"
SourceData: Source Vela pulsar with:
data: 3,587,233 photons from 2008-08-04 to 2022-05-26
exposure: 3,269,606 intervals, average effective area 3460 cm^2 for 97.7 Ms
rates: source 8.56e-06/s, background 2.05e-06/s, TS 7052015.5
CellData.rebin: Bin photon data into 484032 15-min bins from 54683.0 to 59725.0
SourceData: 4FGL J0633.9+1746: Restoring from cache with key "PSR J0633+1746_data"
SourceData: Source Geminga with:
data: 1,214,119 photons from 2008-08-04 to 2022-08-31
exposure: 3,412,848 intervals, average effective area 2691 cm^2 for 102.1 Ms
rates: source 3.31e-06/s, background 1.11e-06/s, TS 2290923.9
CellData.rebin: Bin photon data into 493344 15-min bins from 54683.0 to 59822.0
SourceData: LS 5039: Restoring from cache with key "P88Y5020_data"
SourceData: Source LS 5039 with:
data: 1,699,846 photons from 2008-08-04 to 2022-05-26
exposure: 3,186,410 intervals, average effective area 2940 cm^2 for 95.1 Ms
rates: source 2.90e-07/s, background 5.79e-06/s, TS 8946.5
CellData.rebin: Bin photon data into 484032 15-min bins from 54683.0 to 59725.0
SourceData: Cygnus X-3: Restoring from cache with key "504H-0317_data"
SourceData: Source Cygnus X-3 with:
data: 1,672,484 photons from 2008-08-04 to 2022-05-26
exposure: 3,584,619 intervals, average effective area 3330 cm^2 for 107.2 Ms
rates: source 3.72e-08/s, background 4.65e-06/s, TS 462.3
CellData.rebin: Bin photon data into 484032 15-min bins from 54683.0 to 59725.0
This analysis uses a time bin size of 0.25 h for a Nyquist frequency of 24 $d^{-1}$ and a total interval of 5042 d.
The tagged frequencies
frequency period comment
0 0.003 365.25 1/year
1 0.011 91.31 4/year
2 0.019 53.05 Fermi precession
3 0.256 3.91 LS 5039 orbital period
4 0.512 1.95 LS 5039 harmonic
5 1.003 1.00 sidereal day
6 5.010 0.20 Cyg X-3 orbital period
7 15.100 0.07 Fermi orbital frequency
Exposure power spectrum
It is also possible to extract only the exposure spectrum
