１ はじめに
科学や工学のさまざまな分野において、実験データの取得と測定信号の処理は欠かせません。その過程でノイズの影響を避けることは難しく、信号を正確に抽出するためには適切なフィルタ処理が不可欠です。特に、信号雑音比（Signal-to-Noise Ratio; SNR）が低い場合、単純なローパスフィルタでは十分な効果が得られず、ノイズのばらつきが誤差に直結してしまいます。このような状況では、 最適フィルタ処理（Optimal Filtering） を用いることでSNRを最大化し、ノイズの影響を最小限に抑えることができます。
実用例として、2023年に打ち上げられたX線分光撮像衛星 XRISM に搭載されているX線マイクロカロリメータ Resolve では、軌道上のオンボードCPUで最適フィルタ処理が行われています。1
本記事では、パルス信号の 波高値（Pulse Height; PH） について、テンプレート波形を用いた最適フィルタで補正する基本的な考え方を解説します。さらに、Pythonを用いた実装例と、実際のデータに適用した際の効果についても紹介します。
２ サンプル
コードを見れば分かる方は、Google Colabで作成した以下の例を参照ください。
３ 最適フィルタ処理とは
原理的には非常に高いエネルギー分解能を実現できるX線マイクロカロリメータですが、実際にはパルス波形がノイズによって変形するため、単純に波高値をエネルギーに較正するだけでは十分な分解能が得られません。そこで、期待される理想的な波形に最もよく一致する信号部分を強調して波高値を取り出すため、 最適フィルタ （ Matched filter とも）というフィルタを適用し、テンプレート波形を作ります。これにより、テンプレート波形と一致する部分が強く現れ、ノイズの影響が減少します。
ここでは、まず最適フィルタ処理のイメージを掴むために、図を用いて考え方を視覚化して直感的な流れを説明します。そして、背景にある厳密な原理について数式による詳細な定式化を行います。
３.１ 基本概念
下の図は、X線マイクロカロリメータで取得できるパルスの例で、左図（Ideal waveform）が理想的な波形、右図（Noise-laden waveform）がノイズの乗った波形となっています。パルスがいつも理想的な波形であれば、波高値は単純に$(\texttt{パルス電圧の最大値}) - (\texttt{パルスの立ち上がり前のオフセット電圧値})$で計算でき、同種のX線エネルギーなら必ず同じ波高値が得られますが、実際は本質的にノイズの影響を避けられないので、波高値が本来よりも誤差を持って算出されます。
そのため、試しに上の右図（Noise-laden waveform）のようなノイズが乗ったパルスを10000個作り、波高値を計算してヒストグラムにしてみると、下図（Histogram of PH）のようになります。2 この波高値こそがX線エネルギーの情報にあたりますから、ノイズの影響を受けた生の波高値をそのまま使ってしまうと、エネルギー分解能は悪く、X線マイクロカロリメータの本来の性能を引き出せません。
そこで、下の左図（Pulse and Noise FFT waveform）は上図（Pulse and Noise FFT waveforms）を作るにあたって作成したパルスとそれに乗せたノイズについて、周波数空間で表したものですが、このSNRを最大化するようなフィルタを取得される全パルスにかけます。これが最適フィルタ処理です。結果、波高値が下の右図（Histogram of PH and PHA）のPHAに変換され（後述：3.2章）、ヒストグラムを見れば分かるように、エネルギー分解能が良くなります。
３.２ 数学的定式化
それでは、具体的にどのようにして最適フィルタ処理を施すのか、数式を用いて説明していきましょう。
測定によって得られたパルスを$P(t)$として、これをFourier変換して周波数空間で表したものを$P(f)$とします。ここで、すべてのパルスが相似形であることを仮定すると、理想的な規格化パルス$M(f)$と振幅$A$によりパルスは$A\times M(f)$で書けます。さらに、ここにノイズ成分$N(f)$が加わると、$P(f)$は次のように表されます：
P(f) = A \times M(f) + N(f)\ . \label{eq1}\tag{1}
このとき、実際に得られたパルスと理想的なパルスの残差$\chi$が最小になるように、振幅$A$を最小二乗法で決定することを考えます。つまり、$\chi^2$を
\chi^2 \equiv \int \frac{|P(f) - A \times M(f)|^2}{|N(f)|^2}\ {\rm d}f \label{eq2}\tag{2}
と定義すると、これを最小にする$A$は
A = \frac{\displaystyle\int_{-\infty}^{\infty}\cfrac{P(f)M^\ast(f)}{|N(f)|^2}\ {\rm d}f}{\displaystyle\int_{-\infty}^{\infty}\cfrac{|M(f)|^2}{|N(f)|^2}\ {\rm d}f} = \frac{\displaystyle\int_{-\infty}^{\infty}\cfrac{P(f)}{M(f)}\left|\cfrac{M(f)}{N(f)}\right|^2\ {\rm d}f}{\displaystyle\int_{-\infty}^{\infty}\left|\cfrac{M(f)}{N(f)}\right|^2\ {\rm d}f} \label{eq3}\tag{3}
で与えられます。この式は、$|M(f)/N(f)|^2$を重みとしたときの、$P(f)/N(f)$の平均値になっていることがわかります。
途中式（折りたたみ）
まず、式$(\ref{eq2})$を$A$で偏微分します：
\begin{align}
\frac{\partial\chi^2}{\partial A} &= \frac{\partial}{\partial A}\int_{-\infty}^{\infty} \frac{|P(f)|^2 - A(P(f)M^\ast(f) + P^\ast(f)M(f)) + A^2|M(f)|^2}{|N(f)|^2}\ {\rm d}f \\
&= -\int_{-\infty}^{\infty}\frac{P(f)M^\ast(f) + P^\ast(f)M(f)}{|N(f)|^2}\ \mathrm{d}f + 2A\int_{-\infty}^{\infty}\frac{|M(f)|^2}{|N(f)|^2}\ \mathrm{d}f\ . \label{eq4}\tag{4}
\end{align}
ここで、上式が$0$のときに$\chi^2$が最小になるとすれば、$A$は容易に解くことができて
\begin{align}
A &= \frac{\displaystyle\int_{-\infty}^{\infty}\cfrac{P(f)M^\ast(f) + P^\ast(f)M(f)}{2|N(f)|^2}\ \mathrm{d}f}{\displaystyle\int_{-\infty}^{\infty}\cfrac{|M(f)|^2}{|N(f)|^2}\ \mathrm{d}f}\\
&= \frac{\displaystyle\int_{-\infty}^{\infty}\cfrac{P(f)M^\ast(f)}{2N(f)N^\ast(f)}\ \mathrm{d}f}{\displaystyle\int_{-\infty}^{\infty}\cfrac{|M(f)|^2}{|N(f)|^2}\ \mathrm{d}f} - \frac{\displaystyle\int_{\infty}^{-\infty}\cfrac{P^\ast(-f)M(-f)}{2N(-f)N^\ast(-f)}\ \mathrm{d}f}{\displaystyle\int_{-\infty}^{\infty}\cfrac{|M(f)|^2}{|N(f)|^2}\ \mathrm{d}f}\\
&= \frac{\displaystyle\int_{-\infty}^{\infty}\cfrac{P(f)M^\ast(f)}{2N(f)N^\ast(f)}\ \mathrm{d}f}{\displaystyle\int_{-\infty}^{\infty}\cfrac{|M(f)|^2}{|N(f)|^2}\ \mathrm{d}f} + \frac{\displaystyle\int_{-\infty}^{\infty}\cfrac{P(f)M^\ast(f)}{2N^\ast(f)N(f)}\ \mathrm{d}f}{\displaystyle\int_{-\infty}^{\infty}\cfrac{|M(f)|^2}{|N(f)|^2}\ \mathrm{d}f}\\
&= \frac{\displaystyle\int_{-\infty}^{\infty}\cfrac{P(f)M^\ast(f)}{|N(f)|^2}\ \mathrm{d}f}{\displaystyle\int_{-\infty}^{\infty}\cfrac{|M(f)|^2}{|N(f)|^2}\ \mathrm{d}f} \label{eq5}\tag{5}
\end{align}
となります。補足ですが、式$(\ref{eq5})$の変形では、実関数$g$のFourier変換$G(f)$について、$G^\ast(f)=G(-f)$が成り立つことを使いました。
最小二乗法
測定データの誤差が既知であるとき、標準的な$\chi^2$は次のように定義されました：
\chi^2 \equiv \sum_i \frac{(y_i - y(x_i))^2}{\sigma_i^2} \label{eq6}\tag{6}
- $y_i$：測定データ
- $y(x_i)$：理論モデル
- $\sigma_i$：測定データの誤差
さて、Parsevalの定理から式$(\ref{eq6})$はさらに
A = \frac{\displaystyle\int_{-\infty}^{\infty}\cfrac{P(f)M^\ast(f)}{|N(f)|^2}\ {\rm d}f}{\displaystyle\int_{-\infty}^{\infty}\cfrac{|M(f)|^2}{|N(f)|^2}\ {\rm d}f} = \frac{2\pi\displaystyle\int_{-\infty}^{\infty}P(t){\mathcal F^{-1}}\left[\cfrac{M(f)}{|N(f)|^2}\right]\ {\rm d}t}{\displaystyle\int_{-\infty}^{\infty}\cfrac{|M(f)|^2}{|N(f)|^2}\ {\rm d}f} \label{eq7}\tag{7}
と変形できます。ただし、$\mathcal F^{-1}$はFourier逆変換を表します。また、このときの
T(t) \equiv {\mathcal F^{-1}}\left[\cfrac{M(f)}{|N(f)|^2}\right] \label{eq8}\tag{8}
を最適フィルタの テンプレート（Template） と呼ぶことにします。
Parsevalの定理
Fourier変換された関数の積について、一般化されたParsevalの定理は以下のように記述されます：
\int_{-\infty}^{\infty} x(f)y^\ast(f)\ {\rm d}f = 2\pi\int_{-\infty}^{\infty} x(t)y^\ast(t)\ {\rm d}t = 2\pi\int_{-\infty}^{\infty} x(t){\mathcal F^{-1}}[y(f)]\ {\rm d}t \label{eq9}\tag{9}
以上から、式$(\ref{eq7})$の適当な規格化定数$Z$をとって 波高値（Pulse Height Analysis value; PHA） を書けば
{\rm PHA} = Z\int_{-\infty}^{\infty} P(t)T(t)\ {\rm d}t \label{eq10}\tag{10}
となり、離散的なデータに対しては
{\rm PHA} = Z\sum_{i} P(t_i)T(t_i) \label{eq11}\tag{11}
を得ます。このPHAがまさに、実際に得られたパルスと理想的なパルスの残差が最小になるように最適化された波高値になります。
３.３ エネルギー分解能の限界
最適フィルタ処理を施した場合のエネルギー分解能の限界（$1\sigma$エラー）を考えましょう。これは、式$(\ref{eq2})$の$\chi^2$が$1$だけ増える$A$の変化分で計算できます。よって、$\chi^2$を最小値$A_0$のまわりに二次展開します：
\chi^2(A_0+\delta A) \simeq \chi^2(A_0) + \frac{\partial\chi^2(A_0)}{\partial A}\delta A + \frac{1}{2}\frac{\partial^2\chi^2(A_0)}{\partial A^2}(\delta A)^2\ . \label{eq12}\tag{12}
ここで、$A_0$は$\chi^2(A)$を最小にする点、つまり一階微分が$0$となるように取るので、一次の項はゼロです。また、式$(\ref{eq4})$をさらに$A$で微分すれば、上式から
\chi^2(A_0+\delta A) - \chi^2(A_0) = \int_{-\infty}^{\infty} \frac{|M(f)|^2}{|N(f)|^2}\ {\rm d}f\ (\delta A)^2 \label{eq13}\tag{13}
を得ます。さて、上式の左辺が$1$のときに$1\sigma$エラーとなりますから
\delta A = \left(2\int_{0}^{\infty} \frac{|M(f)|^2}{|N(f)|^2}\ {\rm d}f\right)^{-\frac{1}{2}} \label{eq14}\tag{14}
です。ただし、$M(t)$および$N(t)$は実信号由来のため、$M(-f)=M^\ast(f),\ N(-f)=N^\ast(f)$が成り立ち、積分が対称になることを使いました。
$1\sigma$エラー
カイ二乗分布の自由度$1$における信頼区間の範囲$68\%$に対応します。その意味は以下の記事を参照ください。
式$(\ref{eq14})$はまだ振幅誤差の表式なので、エネルギーに較正し、さらに半値全幅（Full Width of Half Maximum; FWHM）にして書きます：
\Delta E_{\rm FWHM} = 2\sqrt{2\ln2}E\frac{\delta A}{A}
= 2\sqrt{2\ln2}E\left(2\int_{0}^{\infty} \frac{|A\times M(f)|^2}{|N(f)|^2}\ {\rm d}f\right)^{-\frac{1}{2}} \label{eq15}\tag{15}
離散的なデータに対しては
\Delta E_{\rm FWHM} = 2\sqrt{2\ln2}E\left(2\sum_{i} \frac{|A\times M(f_i)|^2}{|N(f_i)|^2}\right)^{-\frac{1}{2}} \tag{16}\label{eq16}
です。これにより期待されるエネルギー分解能を算出できます。
３.４ 要点整理
ここまでのおさらいです。最適フィルタ処理を施すことによって、最適化された波高値$\rm PHA$と期待されるエネルギー分解能$\Delta E_{\rm FWHM}$のそれぞれが以下で算出できます：
\begin{align}
\rm{PHA} &= Z\sum_{i} P(t_i){\mathcal F^{-1}}\left[\cfrac{\overline{P}(f)}{|N(f)|^2}\right]_i \label{eq17}\tag{17}\\
\Delta E_{\rm FWHM} &= 2\sqrt{2\ln2}E\left(2\sum_{i} \frac{|\overline{P}(f_i)|^2}{|N(f_i)|^2}\right)^{-\frac{1}{2}}\ . \label{eq18}\tag{18}
\end{align}
いずれも式$(\ref{eq11})$および式$(\ref{eq16})$を実用的なものに変形しました。これは、式$(\ref{eq1})$において測定信号のノイズが完璧にランダムであれば、足しこむほどノイズが下がるので、$M(f)$を全規格化パルスの平均$\overline{M}(f)$としたときに$\overline{P}(f)=A\times \overline{M}(f)$となり、式$(\ref{eq16})$は式$(\ref{eq18})$で表せます。しかし、式$(\ref{eq17})$と式$(\ref{eq11})$ではテンプレート波形部分の分子が異なり、不正確な変形に思えますが、これは規格化定数$Z$の中に$A$を吸収しているだけで、本質的な情報は変わりません。
以上から、測定信号に対して平均パルスとノイズのデータから、最適フィルタ処理を実行できることになります。
４ 解析方法
続いて、実際に最適フィルタ処理を行って波高値を補正する方法について説明していきます。最初に基本的な手順を説明し、それから実データに基づく具体的な手順を紹介します。実データには、 超伝導転移端センサ（Transition Edge Sensor; TES） を用いたX線マイクロカロリメータ（TESカロリメータ）に、$^{55}\rm Fe$線源を照射して得られた測定データを使用します。
なお、始めに紹介したXRISM搭載のResolveは、温度センサとしてシリコンサーミスタを用いた半導体カロリメータで、TESカロリメータとは異なります。違いの詳細については以下の記事を参照ください。
４.１ 解析の手順
基本的な手順
- 1. $N$個のパルスデータの平均をとったものをモデルパルスとして、それをFourier変換する。
- 2. $N$個のノイズデータをそれぞれFourier変換して、その平均をとる。
- 3. 式$(\ref{eq17})$に従って、テンプレートを作り、各パルスデータとテンプレートを数値積分してPHAに変換する。
実データの手順
手順は上の「基本的な手順」の通りです。具体的な進め方は以下のフローに従います。3 また、フローの最後にあるCalibrationは、$^{55}\rm Fe$線源から放射されるMnKα(5.9keV)とMnKβ(6.5keV)、および原点により、PHAからエネルギーへの較正曲線を作る工程のことです。
エネルギー較正曲線を作るにあたって、PHAの規格化定数は自由に決めて良いので、ここでは簡単に$Z=1$としました。
４.２ 解析コードの実装
Pythonによる実装例を示します。式$(\ref{eq17})$と式$(\ref{eq18})$をそのままコードとして書き起こしているだけなので、特に補足する点はありません。
最適フィルタ処理はavepulse, snr, template, pha = CalPHA(datp, datn)で実行します。2つのパラメータに、それぞれ持っている全パルスデータおよびノイズデータを入れてください。ただし、このときに代入するdatpとdatnは、いずれも行がデータ点、列が各イベントとするよう注意してください。
また、エネルギー分解能の限界はeres = cal_energy_resolution(energy, snr)で計算します。2つのパラメータに、それぞれ検出器に照射されたX線のエネルギー値（単位は$\rm eV$）と最適フィルタ処理によって得られたSNRを入れてください。
より詳細な解析については、2章のサンプルコードを参照ください。
import numpy as np
def cal_noise_pow(noises):
"""
取得した全ノイズデータを高速Fouier変換（Fast Fourier Transform; FFT）したのち、
ノイズスペクトルの平均を計算する関数。
Parameters:
noises (np.ndarray): ノイズデータ配列。データ点数 M 、イベント数 N（shape = [M, N])。
Returns:
npow (np.ndarray): ノイズスペクトルの平均値。イベント数 N（shape = [N]）。
"""
nfft = np.fft.fft(noises, axis=1)
npow = np.mean(np.abs(nfft)**2, axis=0)
return npow
def cal_pulse_pow(avepulse):
"""
取得した全パルスデータの平均を取り、それをFFTする関数。
Parameters:
pulses (np.ndarray): パルスデータ配列。データ点数 M 、イベント数 N（shape = [M, N])。
Returns:
ppow (np.ndarray): パルススペクトル。イベント数 N（shape = [N]）。
"""
pfft = np.fft.fft(avepulse)
ppow = np.abs(pfft)**2
return ppow
def cal_pha(pulses, noises):
"""
最適フィルタ処理を施し、PHをPHAに変換する関数。
Parameters:
pulses (np.ndarray): パルスデータ配列。データ点数 M 、イベント数 N（shape = [M, N])。
noises (np.ndarray): ノイズデータ配列。データ点数 M 、イベント数 N（shape = [M, N])。
Returns:
avepulse (np.ndarray): モデルパルス（平均パルス）。データ点数 M（shape = [M]）。
snr (np.ndarray) : SNR。データ点数 M（shape = [M]）。
template (np.ndarray): テンプレート波形。データ点数 M（shape = [M]）。
pha (np.ndarray) : PHA。イベント数 N（shape = [N]）。
"""
avepulse = np.mean(pulses, axis=0)
snr = np.sqrt(cal_pulse_pow(avepulse)/cal_noise_pow(noises))
template = np.fft.ifft(np.fft.fft(avepulse)/cal_noise_pow(noises)).real
pha = np.sum(pulses*template, axis=1)
return avepulse, snr, template, pha
def cal_energy_resolution(energy, snr):
"""
入射X線に対して、エネルギー分解能の限界（1σエラー）を計算する関数。
Parameters:
energy (float) : 入射X線のエネルギー値。
snr (np.ndarray): 最適フィルタ処理により得られたSNR。データ点数 M（shape = [M]）。
Returns:
eres (float): エネルギー分解能の限界値。
"""
eres = 2 * np.sqrt(2*np.log(2)) * energy / np.sqrt((snr**2).sum(axis=-1)*2)
return eres
# 最適フィルタ処理を実行
avepulse, snr, template, pha = CalPHA(datp, datn)
# 期待されるエネルギー分解能を計算
eres = cal_energy_resolution(energy, snr)
４.３ 実データにおける結果
下の左図（Pulse Height Histogram）が、実データの全イベントに対する波高値スペクトルになります。最適フィルタ処理を行って作成したテンプレート波形が中央図（Template waveform）で、これにより右図（PHA Histogram）の通りPHAに変換されました。
波高値とPHAのヒストグラムをエネルギーに較正して、MnKαにあたるエネルギー分解能を算出したところ、以下の結果が得られました。最適フィルタ処理によるエネルギー分解能の向上、つまりノイズによるばらつきの低減が確認されます。
| 評価項目 | エネルギー分解能（FWHM） |
|---|---|
| 単純な波高値（PH） | $32.63\pm0.69\ {\rm eV}$ |
| 最適フィルタ（PHA） | $10.13\pm0.34\ {\rm eV}$ |
| 限界値（$1\sigma$エラー） | $4.92\ {\rm eV}$ |
エネルギー分解能の算出
実データは、エネルギー分解能$\lesssim10\ {\rm eV}$を達成するTESカロリメータから得られたものです。このとき、Heisenbergの不確定性原理によって特性X線のエネルギー値は揺らぐため、その自然幅（$\lesssim4.5\ {\rm eV}$）を考慮する必要が出てきます。そこで、輝線スペクトルのフィッティングにはVoigt関数を用い、FWHMからエネルギー分解能を求めます。
詳細については、以下の記事を参照ください。
４.４ 補足：モデルパルスを作る際の注意点
最適フィルタ処理は、すべてのパルスが相似形である仮定の下での理論のため、モデルパルスをパルスデータの平均で作る場合、下図のダブルパルスのような無効イベントは除去するようにしてください。
また、理論モデルを使ってモデルパルスとする場合は、コード中のavepulseの変数にオリジナルの関数を入れてください。
５ まとめ
本記事では、X線マイクロカロリメータにおける最適フィルタ処理（Matched filter） の原理と、パルス信号の信号の波高値を求めつつノイズの影響を最小限に抑え、エネルギー分解能を向上させる方法について解説しました。具体的には、以下の内容を網羅的に説明しました。
- 最適フィルタ処理の原理：最小二乗法の考え方に基づいて、どのように信号を強調し、ノイズを抑制するのかを数式で示しました。
- 解析の具体的な手順：実データを用いた解析のフローチャートを用いて、具体的なステップを明確にしました。
- Pythonによる実装例：解析コードの主要な関数を提示し、実際のデータ処理への応用を可能にしました。
- 実データにおける結果：最適フィルタ処理がエネルギー分解能の向上に大きく寄与することを、実際の測定データと数値で示しました。
最適フィルタ処理による波高値の補正は、SNRが最大になるように定め、検出器の性能を最大限に引き出すことで有用な手法です。今回紹介したX線マイクロカロリメータの例だけでなく、微弱な信号を扱う分野で広く応用されています。ぜひ参考にしてみてください。
６ 補遺 ─ 高精度化のための更なる補正
実際には、最適フィルタ処理に加えて以下2つの補正も行います。
- Drift補正：時間経過とともに変化する信号のベースラインの補正をします。ベースライン、つまりオフセット電圧値は、装置由来の揺らぎなどで変動します。これにより、波高値が本来よりも誤差を持って算出されるため、必要な補正となります。
- Arrival Time補正：X線の正確な入射時刻を決定するための補正をします。取得される信号や最適フィルタ処理で作られるテンプレート波形は、実際には離散的です。これにより、パルスとテンプレートの内積を取る際、それぞれのピーク点が連続的に見たときのピークとずれ、波高値が本来よりも誤差を持って算出されるため、これも必須な補正となります。
この他にも解析方法に改善の余地はあり、ニューラルネットワークの使用やリアルタイム処理の実装など、様々なことが研究段階にあります。
６.１ Drift補正
先に紹介したXRISMのResolveでは、ACフィルタを通しておりDC成分が遮断され、Drift補正はできません。しかし、4章の実データの取得で用いたTESは、DC SQUID4による無信号時の値でDC成分を観測できるので、このDrift補正が可能になります。
以下のGIFは、信号のベースラインの変動を大袈裟にですが可視化したものです。グラフはTESの温度抵抗特性（$R-T$曲線）を表しており、超伝導と常伝導の状態間での急峻な抵抗変化が確認できます。また、グラフ左上はTESカロリメータの概念図で、X線が入射したときの応答について、本来のTESの動作点と、その動作点が揺らいだ場合の信号応答の違いをグラフ右下に比較として示しています。
このベースラインの変動は、本来同じエネルギーのX線に由来するはずの波高値にズレを生じさせます。これに対処するため、最適フィルタ処理により得られたPHAとオフセット電圧の値の相関を評価し、その相関がゼロとなるようにPHAを補正する処理をDrift補正と呼びます。
実際の解析では、最適フィルタ処理で得たPHAとオフセット電圧（パルス直前のノイズ成分の平均値）の散布図を作成し、その回帰分析から両者の相関を除去します。この際、通常は一つの輝線に注目して回帰直線を導き、これを全体のPHAに適用します。
解析コードの例は次の通りです。cpha = drift_correction(datp, pha, pha_min, pha_max)で実行します。4つのパラメータに、それぞれ取得した全パルスデータとPHA、さらにdrift補正を決定するために使用するPHA値の選択範囲の上下限値を入れてください。
import numpy as np
def cal_offset(pulses, window):
"""
取得した全パルスデータから、そのオフセット電圧値を計算する関数。
Parameters:
pulses (np.ndarray): パルスデータ配列。データ点数 M 、イベント数 N（shape = [M, N])。
window (float) : 信号全体に含まれるデータ点の総数。
Returns:
offset (np.ndarray): オフセット電圧値。イベント数 N（shape = [N]）。
"""
baselines = []
for pulse in pulses:
baseline = pulse[:int(window*5.e-2)] # 取得信号の形状によってベースラインの取り方は異なる
baselines.append(baseline)
offset = np.mean(baselines, axis=1)
return offset
def drift_correction(pulses, pha, pha_min, pha_max, window):
"""
取得した全パルスデータと、最適フィルタ処理で得たPHAから、drift補正を実行する関数。
Parameters:
pulses (np.ndarray): パルスデータ配列。データ点数 M 、イベント数 N（shape = [M, N])。
pha (np.ndarray) : 最適フィルタ処理から取得したPHA。イベント数 N（shape = [N]）。
pha_min (float) : drift補正を決定するために使用するPHA値の選択範囲の下限。
pha_max (float) : drift補正を決定するために使用するPHA値の選択範囲の上限。
window (float) : 信号全体に含まれるデータ点の総数。
Returns:
cpha (np.ndarray): drift補正をしたPHA。イベント数 N（shape = [N]）。
"""
offset = cal_offset(pulses, window)
mask = (pha >= pha_min) & (pha <= pha_max)
poly = np.polyfit(offset[mask], pha[mask], 1)
cpha = pha / np.polyval(poly, offset)
return cpha
# drift補正を実行
cpha = drift_correction(datp, pha, pha_min, pha_max)
4章の解析に対して、改めてDrift補正を行った結果は以下の通りで、エネルギー分解能は$9.58\pm0.34\ \rm{eV}$と算出されました。
６.２ Arrival Time補正
XRISMのResolveでは、TICK_SHIFTというパラメータで時間方向をシフトして最適化しています。より詳細な説明として、埼玉大学の武田さんの修士論文p.27から、該当部分の引用です：
実際のソフトウェア上でのクロスコリレーションは、テンプレートに対して生データを時間軸方向の前後に${\rm i}$ ADC sampleシフトさせることで計算される。${\rm i}$は$−8\leq{\rm i}\leq7$の範囲をとる。${\rm i}$番目のシフトでは波高値${\rm PHA}({\rm i})$が計算され、その最大値を探される。${\rm i}_{max}$番目のシフトで${\rm PHA}({\rm i}_{max})$が最大になると仮定する。この${\rm i}_{max}$は${\rm TICK\_SHIFT}$という$4\ {\rm bitの}$パラメータに保存される。ここから、本当のピークである${\rm PHA}_{peak}$が二次関数で補間して求められる。${\rm PHA}_{peak}$は${\rm PHA}(16\ {\rm bit})$、波高値が${\rm PHA}$となるときの時刻は${\rm TIME\_VERNIER}(4\ {\rm bit})$として保存される。
専門向けには The Practice of Pulse Processing, J. W. Fowler, B. K. Alpert, W. B. Doriese et al., (2016), Journal of Low Temperature Physics, 184, 374–381（arXivリンク）を参照ください。
この補正が必要な理由は、かなり大袈裟にしていますが下の図が示すように、X線の入射時刻の影響を受け、テンプレート波形と実際のパルス波形のピーク（ここでは最小値）部分のインデックスがズレることにあります。このズレにより、最適フィルタ処理におけるテンプレート波形を用いた信号の強調が不完全になります。そのため、テンプレートのデータ点を時間軸方向にシフトさせ、信号の強調を最大化させるArrival Time補正という処理を施します。
実際の解析では、上で引用した通り、テンプレート波形に対して取得したパルスを時間方向に$j$点ずつシフトさせながらクロスコリレーションを計算します。すなわち、まず式$(\ref{eq11})$に従ってシフト後のパルスとテンプレートの積分値を求めて、その結果得られる${\rm PHA}(j)$が最大となるシフト位置$j$を探索します。そして、この最大値を与える点${\rm PHA}(j)$と、その前後の${\rm PHA}(j-1)$と${\rm PHA}(j+1)$を用いて二次補間を行い、サンプリング点間のサブサンプル精度でより正確なピーク値を求めます。
二次補間とは、隣接する三点を通る二次関数を仮定し、その極値（ここでは最大値）を解析的に求めることで、離散データに対して連続的なピーク位置とピーク値を推定する手法です。
理論式（折りたたみ）
$j-1,\ j,\ j+1$の三点で二次関数$y(x) = ax^2 + bx + c$を仮定します。また、基準を$x=0$に置いて、$x=-1,\ 0,\ +1$に対応する値をそれぞれ$y_{-1}={\rm PHA}(j-1),\ y_0={\rm PHA}(j),\ y_{+1}={\rm PHA}(j+1)$とします。これらを代入すれば\begin{cases}
{\rm PHA}(j-1) &= a - b + c \\
{\rm PHA}(j) &= c \\
{\rm PHA}(j+1) &= a + b + c
\end{cases}
\label{eq19}\tag{19}
となるので、各定数は容易に解け
\begin{align}
a &= \frac{1}{2}\{{\rm PHA}(j-1)-2{\rm PHA}(j)+{\rm PHA}(j+1)\} \\
b &= \frac{1}{2}\{{\rm PHA}(j+1)-{\rm PHA}(j-1)\} \\
c &= {\rm PHA}(j)
\end{align}
\label{20}\tag{20}
を得ます。さて、二次関数$y=ax^2+bx+c$の頂点（極値）の位置は
x_{\rm peak} = -\frac{b}{2a} \tag{21}
なので、極大値そのものは
y_{\rm peak} = c - \frac{b^2}{4a} \tag{22}
です。したがって、この式に式$(\ref{20})$を代入すれば
{\rm PHA}_{\rm peak} = {\rm PHA}(j) - \frac{1}{8}\frac{({\rm PHA}(j+1)-{\rm PHA}(j-1))^2}{{\rm PHA}(j-1)-2{\rm PHA}(j)+{\rm PHA}(j+1)} \label{23}\tag{23}
となります。
解析コードの例は次の通りです。atpha, atpulse = arrival_time_correction(datp, template)で実行します。2つのパラメータに、それぞれ取得した全パルスデータとテンプレートを入れてください。ここでは、時間方向のシフトは$-20\leq j\leq 20$としています。
def arrival_time_correction(pulses, template):
"""
取得した全パルスデータから、Arrival Time補正を実行する関数。
Parameters:
pulses (np.ndarray) : パルスデータ配列。データ点数 M 、イベント数 N（shape = [M, N])。
template (np.ndarray): テンプレート波形。データ点数 M（shape = [M]）。
Returns:
pha_peak (np.ndarray) : Arrival Time補正をしたPHA。イベント数 N（shape = [N]）。
shift_pulse (np.ndarray): Arrival Timeお補正によりシフトしたパルスデータ配列。データ点数 M 、イベント数 N（shape = [M, N])。
"""
eventnums = pulses.shape[0]
SHIFT_RANGE = np.arange(-20, 20)
pha_shifted = np.zeros((eventnums, len(SHIFT_RANGE)))
for k, i in enumerate(SHIFT_RANGE):
pulse_shifted = np.roll(pulses, i, axis=1)
pha_shifted[:, k] = np.sum(template * pulse_shifted, axis=1)
idx_max = np.argmax(pha_shifted, axis=1)
pha_j = pha_shifted[np.arange(eventnums), idx_max]
tick_shift = SHIFT_RANGE[idx_max]
shift_pulse = np.array([np.roll(pulses[i], shift) for i, shift in enumerate(tick_shift)])
valid_mask = (idx_max > 0) & (idx_max < len(SHIFT_RANGE) - 1)
idx_valid = np.where(valid_mask)[0]
if len(idx_valid) > 0:
pre_pha_j = pha_shifted[idx_valid, idx_max[idx_valid] - 1]
next_pha_j = pha_shifted[idx_valid, idx_max[idx_valid] + 1]
denominator = 2 * pha_j[idx_valid] - next_pha_j - pre_pha_j
non_zero_mask = np.abs(denominator) > 1e-9
non_zero_idx = idx_valid[non_zero_mask]
final_pha_j = pha_j[non_zero_idx]
final_pre_pha_j = pre_pha_j[non_zero_mask]
final_next_pha_j = next_pha_j[non_zero_mask]
final_denominator = denominator[non_zero_mask]
final_numerator = final_next_pha_j - final_pre_pha_j
pha_peak = pha_j.copy()
pha_peak[non_zero_idx] = final_pha_j + (1. / 8.) * (final_numerator**2 / final_denominator)
else:
pha_peak = pha_j
return pha_peak, shift_pulse
atpha, atpulse = arrival_time_correction(datp, template)
実際の解析では、順序としては先に説明したDrift補正よりも先に行う補正です。また、Arrival Time補正はサンプリングが少ない場合に効果があるので，この補正の有効性はケースバイケースです。
関連記事
-
XRISMに搭載されている波形処理システム（Pulse Shape Processor; PSP）は、運用終了X線天文衛星「ひとみ」（ASTRO-H）と同等のものです。PSPの軌道上初期性能の結果については、In-flight performance of pulse-processing system of the ASTRO-H/Hitomi soft x-ray spectrometer, Y. Ishisaki, S. Yamada et al., (2018), JATIS, 4, 1, 011217 を参照ください。 ↩
-
ノイズの波形は、NumPyのnp.random.normalで作ったGaussianノイズと$50\ {\rm kHz}$および$100\ {\rm kHz}$のノイズを足し合わせたものです。また、ここでは波高値を$(\texttt{パルス電圧の最大値}) - (\texttt{パルスの立ち上がり前のオフセット電圧データ点の平均値})$で計算しました。 ↩
-
このフロー中のパルス信号は、３.１章で見たX線マイクロカロリメータから得られるパルスの形と異なり反転していますが、これはTESカロリメータで最終的に得られる電圧変化はマイナスがかかるためです。 ↩
-
TESカロリメータの微小な抵抗変化を読み出すための高感度磁気センサ（Superconducting QUantum Interference Device; SQUID）です。 ↩