はじめに
天文学における「位相折り畳み (epoch folding)」は、観測された時系列データに対して、ある周期 $P$ を仮定し、その周期でデータを「折り畳んで」一周期のプロファイルを再構築する方法です。
HEASOFT の XRONOS efold【公式マニュアル】はこの処理を行う標準ツールで、これについては「すざく」の XRONOSの使い方などが参考になるかと思います。
efsearch は周期探索をするツールです。この記事では、python で簡単な周期サーチの実装も紹介します。
本記事ではその簡易版を Python で実装し、さらに「周期を10%わざとずらした場合にどう見えるか」を比較できるように工夫した例と、簡易的な周期探索の例を紹介します。
宇宙の時系列データ解析に特化分析ライブラリstingray を使っても良いでしょうが、宇宙のデータは観測的な事情により、自分で簡易的にコーディングできた方が便利なことが多いでしょう。
宇宙の実際のデータを使いたい方は、MAXIデータなど、
こちらの記事で紹介してますので、手に入れてください。この記事では擬似だけで完結する簡単な例を紹介します。
数式でみる「folding」の仕組み
位相折り畳みの考え方はシンプルです。
観測された各時刻を $t_i$、仮定する周期を $P$、基準時刻 (epoch) を $T_0$ とすると、位相 (phase) $\phi_i$ は次式で定義されます：
\phi_i = \frac{(t_i - T_0) \bmod P}{P}, \quad 0 \leq \phi_i < 1
ここで：
- $(t_i - T_0) \bmod P$ は「時刻差を周期 $P$ で割った余り」
- $\phi_i$ は「その時刻が周期のどこに位置するか」を示す無次元量
この $\phi_i$ を横軸に、対応する観測値（カウント数や光度）$x_i$ を縦軸にしてプロットすると、全ての周期が 0–1 の範囲に折り畳まれることになります。
さらに、位相を $N$ 個のビンに分け、それぞれのビンに入ったデータ点の平均を取ると：
\bar{x}_k = \frac{1}{N_k} \sum_{i \in \text{bin }k} x_i
という 折り畳みプロファイル（平均光度曲線）が得られます。
ここで $N_k$ はビン $k$ に含まれるデータ点の数です。
この手法は特に パルサーのパルスプロファイルや連星系の軌道変調の解析で頻繁に使われます。
実装と「ずらし検証」
今回の実装では、以下を行います：
- 正しい周期 $P$ で折り畳んでプロファイルを表示
- 周期を わざと 10% ずらす（$P' = 1.1P$）
- 折り畳み結果を比較して、信号が崩れる様子を確認
コード例
import numpy as np
import matplotlib.pyplot as plt
from astropy.io import fits
# 擬似データの生成
np.random.seed(42)
period_true = 10.0 # 真の周期 [s]
time = np.linspace(0, 300, 1000) # 0～300秒まで1000点
rate = 5.0 + 2.0 * np.sin(2 * np.pi * time / period_true) + np.random.normal(0, 0.5, len(time))
# プロット
plt.figure(figsize=(8,4))
plt.plot(time, rate, '.')
plt.xlabel("Time [s]")
plt.ylabel("Count rate")
plt.title("Simulated Light Curve")
plt.show()
# ASCIIに保存
np.savetxt("sim_lightcurve.txt", np.column_stack([time, rate]))
def efold(filename, period, epoch=None, pfrac=0.1, n_phase_bins=20,
time_col=0, rate_col=1, rate_col_error=2, skiprows=0, flag_use_rate_error=False, sn_limit=0.05):
# ファイル形式の判定
file_is_fits = filename.lower().endswith('.fits')
if file_is_fits:
with fits.open(filename) as hdul:
data = hdul[1].data
times = data['TIME']
if 'RATE' in data.columns.names:
intensities = data['RATE']
elif 'COUNT' in data.columns.names:
intensities = data['COUNT']
else:
intensities = data.field(1)
else:
raw = np.loadtxt(filename, skiprows=skiprows)
times = raw[:, time_col]
intensities = raw[:, rate_col]
if flag_use_rate_error:
intensities_error = raw[:, rate_col_error]
cutid = np.where(intensities_error/intensities < sn_limit)
print(f"event cut : {len(intensities)} --> {len(cutid[0])}")
times = times[cutid]
intensities = intensities[cutid]
intensities_error = intensities_error[cutid]
if epoch is None:
epoch = times[0]
def fold_and_bin(times, intensities, period, epoch, n_phase_bins):
# 位相を 0–1 に正規化
phases = ((times - epoch) % period) / period
# 2周期分に拡張（0–2）
phases_extended = np.concatenate([phases, phases + 1.0])
intensities_extended = np.concatenate([intensities, intensities])
# 並べ替え
idx = np.argsort(phases_extended)
phases_sorted = phases_extended[idx]
intensities_sorted = intensities_extended[idx]
# bin の設定: 0.0, 0.5, 1.0, 1.5, 2.0 に点が来るように
bin_edges = np.linspace(0, 2, n_phase_bins + 1)
bin_centers = 0.5 * (bin_edges[:-1] + bin_edges[1:])
# 平均値を for ループで計算
bin_means = []
for i in range(n_phase_bins):
if i < n_phase_bins - 1:
mask = (phases_sorted >= bin_edges[i]) & (phases_sorted < bin_edges[i+1])
else:
# 最後の bin は右端も含める
mask = (phases_sorted >= bin_edges[i]) & (phases_sorted <= bin_edges[i+1])
if np.any(mask):
mean_val = intensities_sorted[mask].mean()
else:
mean_val = np.nan # データが無い場合
bin_means.append(mean_val)
bin_means = np.array(bin_means)
return phases_sorted, intensities_sorted, bin_centers, bin_means
# 通常の折り畳み
phases, ints, centers, means = fold_and_bin(times, intensities, period, epoch, n_phase_bins)
# 軌道周期を pfrac だけわざとずらす(検証用)。
period_pfrac = pfrac * period
phases_off, ints_off, centers_off, means_off = fold_and_bin(times, intensities, period_pfrac, epoch, n_phase_bins)
# プロット
fig, axes = plt.subplots(1, 2, figsize=(12,4), sharey=True)
# 正しい位相
axes[0].scatter(phases, ints, s=1, alpha=0.3, label='Folded data')
axes[0].plot(centers, means, '-o', color='orange', label='Binned profile')
axes[0].set_title("Correct Epoch")
axes[0].set_xlabel("Normalized Orbital Phase")
axes[0].set_ylabel("Intensity")
axes[0].legend()
# ずらした位相
axes[1].scatter(phases_off, ints_off, s=1, alpha=0.3, label='Offset-Period data')
axes[1].plot(centers_off, means_off, '-o', color='orange', label='Binned profile')
axes[1].set_title(f"Period offset by {pfrac*100:.1f}% of P")
axes[1].set_xlabel("Normalized Orbital Phase")
axes[1].legend()
plt.tight_layout()
plt.show()
# efold 実行
efold("sim_lightcurve.txt", period=10.0, epoch=None, pfrac=0.1, n_phase_bins=40)
実行結果
擬似的に生成したライトカーブ
folding したライトカーブ
黄色は位相幅に含まれるデータ内で平均した値で、青は単にライトカーブを位相(time%period)で表示しただけである。このように、平均化するときは、平均化する前の生データも表示して、全データが正しそうであることを可視化して確認しましょう。
-
正しい周期 $P=10.0$ s
→ プロファイルがきれいな正弦波状に揃う。 -
ずれた周期 $P'=11.0$ s
→ データ点がバラけ、平均プロファイルも潰れてしまう。
つまり、「folding」が鋭く揃っているかどうかが、周期が正しいかどうかの判断材料になります。
この例はあまりにわかりやすいので、有り難みがわからないと思いますが、食のある連星系について、MAXIのライトカーブの folding などを行う時には、統計誤差によるデータが揺らぐので、このような方法で真の変動なのか、偽物か、判断材料として使います。
まとめ (foldingに関して)
- folding は「$t \mapsto (t-T_0) \bmod P$」という単純な操作。
- 正しい周期で折り畳むと信号が際立ち、ずれると崩れる。
- Python 実装により、XRONOS
efoldの基本を手軽に理解できる。
【周期サーチを行いたい場合】
可視光観測などで、周期が厳密に決まってる場合は周期を固定値にする場合もあるが、周期にまだ曖昧さがある場合は、時間変化してる可能性がある場合は、周期サーチも重要です。その場合についても考えてみましょう。
まずは周期探索アルゴリズムを数式で説明します。ここでは、与えられた観測データ列を 仮定周期ごとに位相折り畳みし、そのプロファイルの「鋭さ」 を χ² 統計量で評価する手法で、周期サーチ方法の中でも、最も基本的となるカイ2乗検定を用いた方法を紹介します。
folding の部分は上記とやや重なりますが、それも含めて数式で概要を整理してみましょう。
1. データと仮定周期
観測データが
\{\, (t_i, x_i, \sigma_i) \,\}_{i=1}^N
で与えられるとします。
- $t_i$ : 観測時刻
- $x_i$ : 観測された強度（count rate など）
- $\sigma_i$ : 強度の誤差（無ければ $\sigma_i = 1$ と置く）
仮定する周期を $P$、基準時刻（epoch）を $T_0$ とします。
2. 位相の計算
各観測点に対して「位相」$\phi_i$ を計算します：
\phi_i = \frac{(t_i - T_0) \bmod P}{P}, \quad 0 \leq \phi_i < 1.
3. 位相ビン分け
位相範囲 [0,1) を $M$ 個のビンに等分します：
\Delta\phi = \frac{1}{M}, \quad \text{bin } k: \; [k\Delta\phi, (k+1)\Delta\phi).
各観測点をその位相に応じてビンに振り分けます。
4. ビン平均の計算
ビン $k$ に属するデータ点の集合を $\mathcal{I}_k$ とすると、ビン平均値 $\bar{x}_k$ を
- 誤差がある場合：重み付き平均
\bar{x}_k = \frac{\sum\limits_{i \in \mathcal{I}_k} \dfrac{x_i}{\sigma_i^2}}
{\sum\limits_{i \in \mathcal{I}_k} \dfrac{1}{\sigma_i^2}}
- 誤差が無い場合：単純平均
\bar{x}_k = \frac{1}{|\mathcal{I}_k|} \sum_{i \in \mathcal{I}_k} x_i
と定義します。
注意点としては、誤差がついていれば、常に誤差で重みづけすべきとは限らない点です。誤差が正しさと不定性も含めて考える必要があります。もし誤差が異常に小さいデータが含まれていたら、その誤差を使って重みづけするほうが間違えてしまうこともあります。誤差の扱いの基本は、なるべく単純簡潔に記述して、誤差の不定性によらないデータ解析を心がけると良いでしょう。
5. χ² の計算
Stingray の “Pulsation search with epoch folding” と同様に，フラットなモデル（＝無変調）に対する χ² を使う形を紹介します。
位相で M 分割した折り畳みプロファイル ($x_k$) を作り，全位相での（重み付き）平均$ (\langle x\rangle)$ を平坦モデルとして
\chi^2(P)=\sum_{k=0}^{M-1}\frac{\bigl(x_k-\langle x\rangle\bigr)^2}{\sigma_k^2}
を各試行周期 (P) について計算し，最小値（≒最も「平坦」に見える周期）/ あるいは最大値（≒最も「非平坦」＝有意なパルス）を指標にします。無変調仮説の下では ($\chi^2$) は自由度 (M-1) の ($\chi^2$) 分布に従う（($\langle x\rangle$) をデータから推定しているため 1 自由度減る）ことに対応します。
コードは def efold_search()（flat-model χ²版） です。
-
位相 0–1 で M ビンに分けて ビン平均 ($x_k$) と ビン平均の標準誤差 ($\sigma_k$) を作ります。
-
誤差列がある場合：各点の重み $w_i=1/\sigma_i^2$。
- ビン平均 $x_k=\sum w_i x_i / \sum w_i$
- ビン平均誤差 $\sigma_k=\sqrt{1/\sum w_i}$
-
誤差列が無い場合：そのビン内の標準偏差から 標準誤差 ($\sigma_k=\mathrm{std}/\sqrt{n_k}$) を推定（データが少ないときは ($\sigma_k$) を小さく見積もりすぎないよう最小値クリップ）。
-
全ビンの（重み付き）平均 ($\langle x\rangle=\sum (x_k/\sigma_k^2) / \sum(1/\sigma_k^2)$) を求め，式のとおり χ² を計算します。
-
周期は P×(1±pfrac) を n_sweep_period ステップで掃引します。
使い分けの目安
- 検出目的（「どこで最も非平坦か？」）なら χ² 最大の周期を採用。
- 平坦性確認（折り畳みの位相合わせのデバッグ等）なら χ² 最小の周期も参照。
- 無変調仮説の下での有意性は，($\chi^2$) と DOF (4$\approx M-1$) から p 値を評価できます（必要なら scipy.stats.chi2.sf を使ってください）。
6. 周期の探索
探索範囲は
P_{\min} = P \,(1 - p_{\mathrm{frac}}), \quad
P_{\max} = P \,(1 + p_{\mathrm{frac}})
とし、これを $n_{\mathrm{sweep}}$ ステップで分割して各候補周期 $P_j$ について χ² を計算します：
P_j = P_{\min} + \frac{j}{n_{\mathrm{sweep}} - 1}\,(P_{\max} - P_{\min}), \quad j=0,\dots,n_{\mathrm{sweep}}-1.
7. 最良周期の決定
得られた χ² の中で最小値を与える周期を最良周期とみなします：
P_{\mathrm{best}} = \arg\min_{P_j} \chi^2(P_j).
コードの紹介
import numpy as np
import matplotlib.pyplot as plt
from astropy.io import fits
def efold_search(
filename,
period, # 中心周期 P（探索の基準となる周期）
epoch=None, # 基準時刻（与えられない場合はデータ先頭を使用）
pfrac=0.1, # 探索幅（中心周期の ±pfrac 倍）
n_sweep_period=101, # 周期ステップ数（探索の分解能。奇数推奨）
n_phase_bins=20, # 位相ビン数 (M)。1周期を M 個に分割
time_col=0,
rate_col=1,
rate_col_error=2,
skiprows=0,
flag_use_rate_error=False,
sn_limit=None, # ASCIIのみ: S/Nカットのしきい値（例: 0.05）
debug=False # デバッグ情報を表示するかどうか
):
"""
【目的】
Stingray の epoch folding 法に基づき、
無変調（フラット）なモデルに対する χ² を計算して周期探索を行う関数。
【手法】
与えられた周期 P のまわりに ±pfrac の範囲で周期を振って試行し、
各周期に対してデータを位相で折り畳む。
折り畳み後のプロファイルを「全体平均（フラットモデル）」と比較し、
χ²(P) = Σ_k ( (x_k - <x>)² / σ_k² ) を計算する。
M = 位相ビン数、x_k = ビン平均、σ_k = その誤差。
【戻り値】
periods: 試行周期の配列
chi2s: 各周期に対する χ² 値
best_period: χ² 最小値・最大値を与える周期（min, max）
best_index: そのインデックス（min, max）
"""
# -------------------------------
# 1. データの読み込み
# -------------------------------
# 入力が FITS か ASCII かを判定し、それぞれ対応する方法で読み込む。
file_is_fits = filename.lower().endswith('.fits')
if file_is_fits:
with fits.open(filename) as hdul:
data = hdul[1].data
times = np.array(data['TIME'], dtype=float)
# 強度（カウント率やカウント数）のカラム名はデータに依存するので順に探索
if 'RATE' in data.columns.names:
intensities = np.array(data['RATE'], dtype=float)
elif 'COUNT' in data.columns.names:
intensities = np.array(data['COUNT'], dtype=float)
elif 'COUNTS' in data.columns.names:
intensities = np.array(data['COUNTS'], dtype=float)
else:
intensities = np.array(data.field(1), dtype=float)
# 誤差カラム（もしあれば使用）
sigmas = None
for cand in ['ERROR','RATE_ERR','ERR','RATEERROR','COUNTS_ERR','COUNT_ERR']:
if cand in data.columns.names:
sigmas = np.array(data[cand], dtype=float)
break
else:
# ASCII の場合
raw = np.loadtxt(filename, skiprows=skiprows)
times = raw[:, time_col].astype(float)
intensities = raw[:, rate_col].astype(float)
sigmas = None
if flag_use_rate_error and raw.shape[1] > rate_col_error:
sigmas = raw[:, rate_col_error].astype(float)
# S/N カット：誤差/強度 がしきい値以上の点を除外する
if sn_limit is not None:
n_before = len(times)
with np.errstate(divide='ignore', invalid='ignore'):
mask_sn = (sigmas / intensities) < sn_limit
times = times[mask_sn]
intensities = intensities[mask_sn]
sigmas = sigmas[mask_sn]
n_after = len(times)
n_cut = n_before - n_after
if debug:
print(f"S/N cut applied (limit={sn_limit}): {n_before} -> {n_after} (cut {n_cut})")
# -------------------------------
# 2. epoch の設定
# -------------------------------
# 基準時刻 epoch が与えられていない場合は、データの最初を採用。
if epoch is None:
epoch = times[0]
# -------------------------------
# 3. データのクリーニング
# -------------------------------
# NaN や Inf を含む点を除外して解析の安定性を確保。
n_total = len(times)
if sigmas is None:
m = np.isfinite(times) & np.isfinite(intensities)
else:
m = np.isfinite(times) & np.isfinite(intensities) & np.isfinite(sigmas)
n_valid = np.sum(m)
n_invalid = n_total - n_valid
times = times[m]
intensities = intensities[m]
if sigmas is not None:
sigmas = sigmas[m]
if debug:
print("=== efold_search DEBUG INFO ===")
print(f"Total data points : {n_total}")
print(f"Valid data points : {n_valid}")
print(f"Excluded (NaN/Inf) : {n_invalid}")
print(f"Time range : {times.min():.6g} - {times.max():.6g}")
print(f"Intensity mean ± std : {np.mean(intensities):.6g} ± {np.std(intensities):.6g}")
if sigmas is not None:
print(f"Error mean ± std : {np.mean(sigmas):.6g} ± {np.std(sigmas):.6g}")
print(f"Min error : {np.min(sigmas):.6g}")
print(f"Max error : {np.max(sigmas):.6g}")
else:
print("No error column used.")
print("===============================")
# -------------------------------
# 4. χ² を計算する内部関数
# -------------------------------
def chi2_flat_for_period(candidate_period):
# 各データ点の位相 (0–1) を計算
phases = ((times - epoch) % candidate_period) / candidate_period
# 位相を M ビンに分割
bin_edges = np.linspace(0.0, 1.0, n_phase_bins + 1)
bin_index = np.searchsorted(bin_edges, phases, side='right') - 1
bin_index[bin_index < 0] = 0
bin_index[bin_index >= n_phase_bins] = n_phase_bins - 1
# 各ビンの平均強度 x_k とその誤差 σ_k を計算
xk = np.full(n_phase_bins, np.nan)
sk = np.full(n_phase_bins, np.nan)
for k in range(n_phase_bins):
mk = (bin_index == k)
if not np.any(mk):
continue
x_in = intensities[mk]
if sigmas is None:
# 誤差列が無い場合: 標準誤差 = 標準偏差/√n
xk[k] = np.mean(x_in)
n_k = x_in.size
if n_k > 1:
sk[k] = np.std(x_in, ddof=1) / np.sqrt(n_k)
else:
sk[k] = max(1e-6, 0.1 * abs(xk[k]) + 1e-6)
else:
# 誤差列がある場合: 重み付き平均とその誤差
s_in = sigmas[mk]
good = np.isfinite(s_in) & (s_in > 0) & np.isfinite(x_in)
if not np.any(good):
continue
w = 1.0 / (s_in[good] ** 2)
xg = x_in[good]
xk[k] = np.sum(w * xg) / np.sum(w)
sk[k] = np.sqrt(1.0 / np.sum(w))
if sk[k] < 1e-12:
sk[k] = 1e-12
# データが存在するビンのみ採用
ok = np.isfinite(xk) & np.isfinite(sk) & (sk > 0)
xk = xk[ok]
sk = sk[ok]
if xk.size < 2:
return np.inf # ビンが少なすぎる場合は χ² 無限大
# 全体平均 <x> を計算（重み = 1/σ²）
wk = 1.0 / (sk ** 2)
x_mean = np.sum(wk * xk) / np.sum(wk)
# χ² 計算：フラットモデルに対する残差の二乗和
chi2 = np.sum(((xk - x_mean) ** 2) / (sk ** 2))
return chi2
# -------------------------------
# 5. 周期を掃引して χ² を計算
# -------------------------------
pmin = period * (1.0 - pfrac) # 最小試行周期
pmax = period * (1.0 + pfrac) # 最大試行周期
periods = np.linspace(pmin, pmax, n_sweep_period) # 試行周期リスト
chi2s = np.array([chi2_flat_for_period(p) for p in periods])
# χ² が最小・最大となる周期を抽出
best_index_min = int(np.argmin(chi2s))
best_period_min = periods[best_index_min]
best_index_max = int(np.argmax(chi2s))
best_period_max = periods[best_index_max]
# -------------------------------
# 6. プロット
# -------------------------------
plt.figure(figsize=(7,4))
plt.plot(periods, chi2s, '-', lw=1.5)
plt.axvline(best_period_max, ls='--', alpha=0.5,
label=f'max χ² @ {best_period_max:.6g}')
plt.axvline(best_period_min, ls='--', alpha=0.5,
label=f'min χ² @ {best_period_min:.6g}')
plt.xlabel("Trial Period")
plt.ylabel(r"$\chi^2$ vs flat profile")
plt.title(f"Epoch-folding χ² (flat model), M={n_phase_bins} bins")
plt.legend()
plt.tight_layout()
plt.show()
return periods, chi2s, (best_period_min, best_period_max), (best_index_min, best_index_max)
# 例：中心周期10秒、±10%の範囲を 201 ステップで走査、位相ビンは40
periods, chi2s, bestP, idx = efold_search(
"sim_lightcurve.txt",
period=10.0,
epoch=None,
pfrac=0.1,
n_sweep_period=201,
n_phase_bins=40,
time_col=0,
rate_col=1,
debug=True
# 誤差列がある場合のみ
# flag_use_rate_error=True, rate_col_error=2, sn_limit=0.05
)
print("Best period (best_period) =", bestP[1])
実行結果
周期10.0秒のsin波を入れてるので、10.0秒のところのカイ2乗値が大きい(non pulse と大きく異なる)はずだが、、10.01秒がミニマムになっている。統計の揺らぎですね。
まとめ (周期サーチ)
- 観測データを位相折り畳みして 「各ビン内の平均プロファイル」 を作成。
- そのプロファイルからのズレの大きさを χ² で定量化。
- 周期を少しずつ変えながら χ² を比較し、最大値を与える周期を採用。
これにより、正しい周期ではデータが揃って χ² が直線から大きく外れて大きくなる一方、間違った周期では揃わず直線的になるので χ² が小さくなる、という性質を利用して周期探索を行います。
上級者向け：Z² 検定から H-test, Lomb–Scargle 法へ
Z² 検定の定義
χ² 折り畳み法は「位相ビン平均とのズレ」を見る方法でしたが、フーリエ成分を直接評価する方法として Z² 検定 がよく使われます。
観測位相を $\phi_i$ ($i=1, \dots, N$) とすると、次数 $m$ の Z² 統計量は次式で定義されます：
Z^2_m = \frac{2}{N} \sum_{k=1}^m \left[ \left( \sum_{i=1}^N \cos(2\pi k \phi_i) \right)^2 + \left( \sum_{i=1}^N \sin(2\pi k \phi_i) \right)^2 \right].
これは「位相分布が一様であれば小さく、特定の位相に集中していれば大きくなる」統計量で、特にパルサーのような周期現象検出に広く使われています。次数 $m$ を変えることで、単純な正弦波から複雑なパルス波形まで柔軟に対応できます。
stingray では下記に解説があります。
これをマグネター、研究に応用した例もあります(牧島先生の天文月報)。
「マグネターの自由歳差運動の観測的研究」牧島 一夫
H-test
Z² 検定では次数 $m$ をあらかじめ決める必要がありますが、H-test は $m$ を 1〜20 程度で自動的に走査し、最も有意なものを選びます。定義は以下の通りです：
H = \max_{1 \leq m \leq m_{\max}} \left( Z^2_m - 4(m-1) \right).
この「罰則項 $4(m-1)$」により、高調波を不必要に増やして χ² を稼ぐことを防ぎます。H-test は波形の形状に依存せず最適な次数を選んでくれるため、パルサー解析で標準的に利用されます。
Python ではパルサータイミングの解析ツールで出来るかもしれませんが、使ったことはないです。。
Lomb–Scargle Periodogram
一方、位相折り畳みではなく、観測データを直接フーリエ的に走査する方法が Lomb–Scargle periodogram (LS法) です。不均一サンプリングの時系列 $(t_i, x_i)$ に対して、各試行角周波数 $\omega$ ごとに
P_{\mathrm{LS}}(\omega) = \frac{1}{2\sigma^2} \left\{
\frac{\left[ \sum_i (x_i - \bar{x}) \cos \omega (t_i - \tau) \right]^2}{\sum_i \cos^2 \omega (t_i - \tau)} +
\frac{\left[ \sum_i (x_i - \bar{x}) \sin \omega (t_i - \tau) \right]^2}{\sum_i \sin^2 \omega (t_i - \tau)}
\right\},
というパワースペクトルを計算します。ここで $\tau$ は位相シフト項です。Astropy のastropy.timeseries.LombScargle や SciPy の scipy.signal.lombscargle で簡単に利用でき、特に観測間隔が不規則な光度曲線の周期探索に有効です。エイリアスには十分な注意が必要な解析です。
あるいは、下記の記事を参考ください。
補足: 時系列データ解析ライブラリ群と中性子星タイミング解析
出来合いのコード群としては Stingray の他にも、Stingrayをベースに使用しつつ NuSTAR など現在のX線望遠鏡のデータ解析に特化した HENDRICS があります。
さらに、これらのライブラリは中性子星の界隈で広く使われている
といったパルサータイミング解析用のライブラリの上に構築されています。
なぜ専用ライブラリが必要なのか？
中性子星の周期探査は、最短で ミリ秒周期 のシグナルを対象にする必要があります。単純な畳み込み解析に加えて、非常に精密な 時刻補正 が必須です。
代表的な補正は以下のようなものです：
-
連星運動による相対論効果
連星系にあるパルサーでは、観測されるパルスの到来時間が重力ポテンシャルや運動によって変動します。例えば、シャピロ遅延 (Shapiro delay) は以下のように近似されます：
\Delta t_{\text{Shapiro}} \simeq -\frac{2 G M_c}{c^3}
\ln \bigl[ 1 - e \cos E \bigr]
ここで $M_c$ は伴星質量、$E$ は離心近点角です。マイクロ秒精度での補正が必要になります。
-
星間物質による分散 (Dispersion Measure, DM)
パルスは星間媒質を通過する際に周波数依存の遅延を受けます：
\Delta t_{\text{DM}} \propto \frac{\text{DM}}{\nu^2}
($\nu$ は観測周波数)。X線帯は影響は小さいですが、電波と組み合わせる際は重要になります。
このように、中性子星の周期探査では
- 畳み込み解析（Stingrayなど）
- 高精度時刻補正（PINT, TEMPO2など）
を組み合わせる必要があります。結果として、Stingray → HENDRICS → PINT/TEMPO2 といった「積み重なった」ライブラリ構造になっています。こうした背景を理解しておくと、なぜ中性子星分野の解析環境が他のX線源より複雑に見えるのかが分かると思います。
まとめ (上級者向け)
- Z² 検定: フーリエ成分を用いた位相一様性検定
- H-test: Z² の高調波次数を自動で最適化
- Lomb–Scargle: 不均一サンプリングに強いフーリエ法
- パルサーの時系列解析はもっと高精度な補正が必要なこともある
χ² 折り畳み法を理解したうえで、これらの手法を使い分ければ、周期信号探索の精度と堅牢性は飛躍的に向上します。ただし、どの手法もデータの統計誤差と系統誤差によってはエイリアス(擬似信号)との戦いになるので、感度ギリギリの周期性探査をやる場合には装置の深い理解も必要であることを記しておきたいと思います。