フィルタと前処理の基礎

FIR filterは畳み込み

\[ y[n]=\sum_{k=0}^{M}b_kx[n-k] \]

一般形のIIR filterは

\[ a_0y[n]=\sum_{k=0}^{M}b_kx[n-k]-\sum_{k=1}^{P}a_ky[n-k] \]

で表せます。\(a_0=1\)へ正規化した式がよく使われます。high-passはDC・ゆっくりしたdrift、low-passは高周波雑音、notchは電源周波数などを抑えます。ただし、filterは情報を追加せず、急峻な遮断はringingを生じ得ます。オフラインのforward-backward filterは位相ずれを相殺しますが非因果的で、振幅応答は片方向filterの \(|H(f)|^2\)、実効次数は概ね2倍になります。padding、初期条件、端点除外、pre-ringingも報告します。遮断周波数は「通過域端」「-3 dB点」など定義が異なるため、filter種、次数、通過域・阻止域、位相特性を併記します。

区間化、窓関数、時間周波数の基礎

矩形窓で有限区間を切ると端点の不連続によってspectral leakageが生じます。Hann、Hamming、Blackman、DPSSなどは漏れと周波数分解能の交換条件を変えます。STFTは

\[ X[k,m]=\sum_{n=0}^{L-1}x[n+mR]w[n]e^{-i2\pi kn/N_{FFT}} \]

で、\(L\) は窓長、\(R\) はhop長、通常 \(N_{FFT}\ge L\) です。長い窓は周波数分解能を上げる一方で急な時間変化をぼかし、短い窓はその逆です。重複率やzero-paddingを増やすと表示・bin samplingは細かくなりますが、intrinsicな時間・周波数分解能は上がりません。同じデータを重ねて使うため、統計上の独立な観測数とスペクトル推定の実効自由度は窓数より少なくなります。

連続wavelet変換（continuous wavelet transform: CWT）は、1つの固定窓ではなく、拡大・縮小したmother wavelet \(\psi\) と信号の一致度を計算します。scale \(a>0\)、時刻 \(b\) におけるCWTを

\[ W_x(a,b)=\frac{1}{\sqrt a}\int_{-\infty}^{\infty} x(t)\psi^*\!\left(\frac{t-b}{a}\right)dt \]

と定義します。連続式では \(a\) は \(t\) と同じ時間単位です。\(a\) が小さいほど短く高周波のwavelet、\(a\) が大きいほど長く低周波のwaveletになります。mother waveletの無次元中心角周波数を \(\omega_\psi\) とすると、連続scale \(a\) と、sample単位の離散scale \(a_d\) に対応する代表周波数は概ね

\[ f(a)=\frac{\omega_\psi}{2\pi a}, \qquad f(a_d)=\frac{\omega_\psi f_s}{2\pi a_d} \quad[\mathrm{Hz}] \]

です。ただし、scale–frequency変換はwaveletの定義とsoftwareの規約に依存します。scalogramは \(|W_x(a,b)|^2\) です。値が大きい点は、その時刻でそのwaveletに似た局所振動のenergyが大きいことを示します。これはPSDそのものではありません。上の連続energy normalizationで元信号の単位をU、時間をsとすると、\(W_x\) はU\(\sqrt{\mathrm{s}}\)、\(|W_x|^2\) はU²·sです。離散libraryではscaleをdimensionlessにして係数を元信号と同じ単位にするなど規約が異なるため、絶対値を比較するときはnormalizationを明記します。admissible waveletでは

\[ \lVert x\rVert_2^2= \frac{1}{C_\psi}\int_{-\infty}^{\infty}\int_0^\infty |W_x(a,b)|^2\frac{da\,db}{a^2} \]

というenergy保存関係を持ちます。離散実装ではscale・時間の積分重みを含めない単純な \(|W|^2\) の和は、元信号のenergyと一致しません。横方向に続くridgeは持続する振動、傾いたridgeは周波数の時間変化、短時間に多数scaleへ広がる縦方向の構造は急峻なtransientを表し得ます。ただし、瞬目、体動、spike、padding端も同じ形を作れます。

時間分解能と周波数分解能の交換条件はwaveletのcycle数・帯域幅で決まります。mother wavelet、中心周波数またはcycle数、scale/frequency grid、voices per octave、normalization、padding、平滑化を報告します。有限記録の端ではwaveletがデータ外へはみ出すため、cone of influence（COI）で端点影響が大きいと判定された領域はpaddingへの依存が強くなります。COIを表示し、主要解析から除く範囲を事前に決めます。COIはCWT固有の端点指標であり、STFTなどへ同じ名称を一般化しません。CWTの原型は Goupillaud, Grossmann & Morlet (1984) に示されています。

品質評価は解析と同じ重要度で扱います。少なくともvalid fraction、欠測率、飽和・clipping率、flatline率、artifact burden、除外チャネル・試行率を残します。SNRは定義が複数あるため、例えば

\[ SNR_{dB}=20\log_{10}\frac{RMS_{signal}}{RMS_{noise}} \]

のようにsignal区間とnoise区間を明示します。

2.5.1 EMD、EEMD、CEEMDANとHilbert–Huang変換

empirical mode decomposition（EMD、経験的モード分解）は、局所的に速い振動から順にintrinsic mode function（IMF）へ分けます。IMFは原則として、極値数とzero crossing数の差が1以下で、上下包絡の局所平均がほぼ0の成分です。残差を \(r_0[n]=x[n]\) とし、第 \(k\) 成分を求めるsiftingは次のように行います。

1. \(h_0[n]=r_{k-1}[n]\) とする。

2. \(h_j\) の極大値と極小値をそれぞれ補間し、上包絡 \(e_{\max,j}[n]\) と下包絡 \(e_{\min,j}[n]\) を作る。

3. 局所平均

   \[ m_j[n]=\frac{e_{\max,j}[n]+e_{\min,j}[n]}{2} \]

   を引き、\(h_{j+1}[n]=h_j[n]-m_j[n]\) とする。

4. 指定したIMF停止条件を満たすまで2–3を反復し、\(c_k[n]=h_j[n]\) とする。

5. \(r_k[n]=r_{k-1}[n]-c_k[n]\) とし、残差の極値が2個未満になるまで次のIMFを求める。

したがって

\[ x[n]=\sum_{k=1}^{K}c_k[n]+r_K[n] \]

と再構成します。\(c_k\) と \(r_K\) の単位は元信号と同じです。小さい \(k\) は一般に速い成分、大きい \(k\) は遅い成分になりますが、IMF番号と周波数帯域は記録ごとに変わります。IMF 3を常にalpha、呼吸、LFなどとみなすことはできません。Huang et al. (1998)

EMDで異なるscaleが同じIMFへ入る、または同じscaleが複数IMFへ分かれる現象をmode mixingと呼びます。ensemble EMD（EEMD）は、標準偏差1の独立white noise \(w_i[n]\) を加えた

\[ x_i[n]=x[n]+\varepsilon\sigma_x w_i[n], \qquad i=1,\ldots,I \]

を各回EMDし、対応するIMFを

\[ \bar c_k[n]=\frac{1}{I}\sum_{i=1}^{I} \mathcal E_k\{x_i[n]\} \]

と平均します。\(\mathcal E_k\{\cdot\}\) はEMDで得た第 \(k\) IMFです。noiseがscale空間の参照を作り、mode mixingを軽減する狙いがあります。一方、有限ensembleでは加えたnoiseが概ね \(\varepsilon\sigma_x/\sqrt I\) のorderで残り、realizationごとにIMF数が異なる問題があります。Wu & Huang (2009)

complete ensemble EMD with adaptive noise（CEEMDAN）は、各段階の残差へ対応scaleのnoise modeを加え、1つの残差から次成分を作ります。以下では \(\varepsilon_k\) に各段階のnoise amplitudeと単位scaleを含めます。標準的な表記の一例は

\[ \widetilde c_1[n]=\frac{1}{I}\sum_{i=1}^{I} \mathcal E_1\{x[n]+\varepsilon_0w_i[n]\}, \qquad r_1[n]=x[n]-\widetilde c_1[n] \]

で、\(k\ge1\) に対して

\[ \widetilde c_{k+1}[n]=\frac{1}{I}\sum_{i=1}^{I} \mathcal E_1\!\left\{ r_k[n]+\varepsilon_k\mathcal E_k\{w_i[n]\} \right\}, \qquad r_{k+1}[n]=r_k[n]-\widetilde c_{k+1}[n] \]

と進めます。この定義では

\[ x[n]=\sum_{k=1}^{K}\widetilde c_k[n]+r_K[n] \]

が数値精度内で成立します。Torres et al. (2011) 改良版CEEMDANはnoise modeの正規化や局所平均operatorが異なるため、単に「CEEMDAN」とせず、algorithm、software/version、noise scalingを示します。Colominas, Schlotthauer & Torres (2014)

Hilbert–Huang変換（HHT）は、各IMFへHilbert変換を適用して解析信号を作ります。

\[ z_k(t)=c_k(t)+i\mathcal H\{c_k(t)\} =a_k(t)e^{i\phi_k(t)} \]

\[ a_k(t)=|z_k(t)|,\qquad f_k(t)=\frac{1}{2\pi}\frac{d}{dt} \operatorname{unwrap}\{\phi_k(t)\}\quad[\mathrm{Hz}] \]

\(a_k\) は元信号と同じ単位、\(f_k\) はHzです。energy表示のHilbert spectrumを

\[ S_H(t,f)=\sum_k a_k^2(t)\, \delta\{f-f_k(t)\} \]

とすると、連続周波数密度としての単位は信号単位²/Hzで、

\[ \int S_H(t,f)df=\sum_k a_k^2(t) \]

です。離散frequency binへ \(a_k^2\) を割り当てただけなら単位は信号単位²/binであり、Hz当たりにするにはbin幅で割ります。\(a_k^2\) は解析信号のenvelope squaredで、実信号の物理的な瞬時powerと常に同じではありません。例えば \(c_k(t)=A\cos(2\pi ft)\) では \(a_k^2=A^2\) ですが、実信号の1周期平均powerは \(A^2/2\) です。文献によって \(a_k\) を表示するamplitude spectrumと \(a_k^2\) を表示するenergy spectrumがあるため、係数とnormalizationを明記します。時間積分したmarginal Hilbert spectrum

\[ h(f)=\int S_H(t,f)dt \]

は、観測時間にわたり各瞬時周波数へ割り当てられたsquared-envelope量を表し、上の連続定義では単位は信号単位²·s/Hzです。

\(a_k^2(t)\) が大きい点は、そのIMFの局所振幅が大きいことを示します。横に安定した \(f_k(t)\) はほぼ一定周波数、上昇・下降する軌跡はfrequency modulationの候補です。軌跡の飛び、複数帯域への分裂、端での急変は真の非定常性だけでなく、mode mixing、spline端点、artifactでも生じます。\(a_k(t)\approx0\) では位相が定まらず、数値微分した瞬時周波数が発散しやすいため、振幅threshold、位相平滑化、除外ruleを定義します。Hilbert瞬時周波数は、1つのIMFが局所的に単一のAM–FM成分とみなせるときに最も解釈しやすくなります。

最低限、包絡補間法、端点延長、sifting停止条件と最大回数、EMD variant、\(\varepsilon_k\)、ensemble数 \(I\)、乱数seed、採用IMF、低振幅mask、位相微分法を報告します。再構成誤差

\[ e_{\mathrm{rec}}= \frac{\left\|x-\sum_kc_k-r_K\right\|_2}{\|x\|_2} \]

に加え、seedとparameterを変えたIMF spectraの重なり、端を除いた結果、既知のAM–FM模擬信号での回収率を確認します。

2.5.2 変分モード分解

variational mode decomposition（VMD、変分モード分解）は、信号をあらかじめ指定した \(K\) 個の狭帯域mode \(u_k(t)\) へ同時に分けます。各modeをHilbert変換でanalytic signalにし、中心角周波数 \(\omega_k\) を掛けて0 Hz付近へ移した後の帯域幅を最小化します。

\[ \min_{\{u_k\},\{\omega_k\}} \sum_{k=1}^{K} \left\| \partial_t \left[ \left(\delta(t)+\frac{i}{\pi t}\right)*u_k(t) \right]e^{-i\omega_kt} \right\|_2^2 \quad \text{subject to}\quad \sum_{k=1}^{K}u_k(t)=x(t) \]

実装では再構成constraintの二次penaltyとLagrange multiplierを加え、alternating direction method of multipliers（ADMM）で \(u_k\) と \(\omega_k\) を交互更新します。正周波数側のmode spectrum \(\widehat u_k(\omega)\) から中心周波数を

\[ \omega_k= \frac{\int_0^\infty\omega|\widehat u_k(\omega)|^2d\omega} {\int_0^\infty|\widehat u_k(\omega)|^2d\omega}, \qquad f_k=\frac{\omega_k}{2\pi}\quad[\mathrm{Hz}] \]

と更新します。\(u_k\) は元信号と同じ単位です。区間 \(n=0,\ldots,N-1\) におけるmodeのmean squareと構成比は

\[ P_k=\frac{1}{N}\sum_{n=0}^{N-1}u_k^2[n], \qquad \rho_k=\frac{P_k}{\sum_jP_j} \]

で、単位はそれぞれ信号単位²と無次元です。\(\Delta t\sum_nu_k^2[n]\) を用いる場合は信号単位²·sの区間energyになります。VMD modeは必ずしも互いに直交しないため、\(\sum_kP_k\) は元信号のmean squareと一致せず、\(\rho_k\) はmode間の二乗振幅構成比です。\(P_k\) や \(\rho_k\) が大きいmodeは、VMDの設定下でそのmodeの二乗振幅が相対的に大きいことを示します。生理的な活動量を直接表すわけではありません。modeへさらにHilbert変換を適用すれば局所振幅と瞬時周波数を求められますが、VMDが直接推定する \(\omega_k\) はmode全体の中心周波数です。

VMDは、信号が少数のcompact-band成分の和で近似できる場合に使いやすい方法です。大きいbandwidth penalty \(\alpha\) は各modeを狭くし、小さい \(\alpha\) は広い周波数変化を許します。\(K\) が多すぎると1つの成分が複数modeへ分裂し、少なすぎると異なる成分が混ざります。異なる初期値で中心周波数が入れ替わる、同じ成分へ複数modeが収束する、弱い成分が強い成分に吸収される場合があります。Fourier domainで解くため暗黙の周期性があり、記録端と急峻なonsetは特に不安定です。原著も、端点影響とmode数の手動指定を主要な限界として挙げています。Dragomiretskiy & Zosso (2014)

\(K\)、\(\alpha\)、noise toleranceまたはdual ascent parameter、DC modeの有無、初期中心周波数、収束許容値、最大反復、padding/mirror extensionを報告します。\(K\) と \(\alpha\) を変えたmodeの再現性、\(\|\sum_ku_k-x\|_2/\|x\|_2\)、mode spectraの重なり、端点除外後の結果を確認します。

2.5.3 Synchrosqueezing変換

synchrosqueezing transform（SST）は、CWTの各係数を、その局所位相から推定した瞬時周波数へ周波数方向に再配置します。2.1節の規約で、\(|W_x(a,b)|>\gamma\) の係数について

\[ \widehat f_x(a,b)= \frac{1}{2\pi} \operatorname{Im} \left\{ \frac{\partial_bW_x(a,b)}{W_x(a,b)} \right\} \quad[\mathrm{Hz}] \]

と推定します。閾値 \(\gamma\) は、ほぼ0の係数で比が発散することを避けます。frequency binを連続表記すると、synchrosqueezed CWTは

\[ T_x(b,f)= \int_{\{|W_x(a,b)|>\gamma\}} W_x(a,b)a^{-3/2} \delta\{f-\widehat f_x(a,b)\}\,da \]

です。離散実装では \(\delta\) を最寄りbinまたは有限幅kernelへの加算で置き換えます。CWTのぼやけたenergyを瞬時周波数付近へ集めるため、ridgeが細く見えます。ridge \(r_k(b)\) の近傍を積分すると

\[ \widehat x_k(b)\approx \frac{1}{C_\psi} \operatorname{Re} \int_{|f-r_k(b)|<\Delta f_k} T_x(b,f)\,df \]

として成分を再構成できます。\(C_\psi\) はwaveletとFourier規約に依存する定数です。\(\widehat x_k\) は元信号と同じ単位、\(r_k\) はHzです。\(T_x\) 自体の単位と絶対値はscale weighting、frequency-bin幅、normalizationに依存するため、software間で値を直接比較しません。

\(|T_x|\) または \(|T_x|^2\) が大きいridgeは、その時刻で対応する局所振動が強いことを表します。水平ridgeはほぼ一定周波数、傾いたridgeはchirp、途切れは振幅低下、成分消失、artifact rejection、ridge追跡失敗のいずれでも起こります。SSTの理論的な分離・再構成は、各成分が緩やかに変わるAM–FM信号で、瞬時周波数が十分離れていることを仮定します。近接または交差するridge、急な周波数jump、低SNR、\(|W|\approx0\)、COIで端点影響が大きい領域では成分の取り違えが起こり得ます。Daubechies, Lu & Wu (2011)

mother wavelet、中心周波数/cycle数、scaleまたはfrequency grid、voices per octave、\(\gamma\)、数値微分法、ridge抽出法、ridge間隔・smoothness penalty、再構成幅 \(\Delta f_k\)、padding、COI除外を報告します。模擬chirpでridge誤差をHz単位で確認し、CWTとSSTを同じcolor normalizationで並べると、再配置前後を区別できます。

2.5.4 再配置spectrogramとWigner–Ville分布の注意

reassigned spectrogramは、STFT energyを計算した格子点のまま表示せず、各局所成分の群遅延と瞬時周波数へ移します。符号はSTFTの定義で変わるため、ここでは

\[ V_x^h(t,\omega)= \int_{-\infty}^{\infty} x(s)h^*(s-t)e^{-i\omega(s-t)}ds, \qquad \phi(t,\omega)=\arg V_x^h(t,\omega) \]

と定義します。\(|V_x^h(t,\omega)|>\gamma\) で

\[ \widehat t(t,\omega)= t-\frac{\partial\phi(t,\omega)}{\partial\omega}, \qquad \widehat\omega(t,\omega)= \frac{\partial\phi(t,\omega)}{\partial t} \]

です。この規約では、純音 \(e^{i\omega_0t}\) は \(\widehat\omega=\omega_0\)、impulse \(\delta(t-t_0)\) は \(\widehat t=t_0\) へ移ります。実装ではunwrapしたphaseの直接差分よりも、time-weighted windowとderivative windowによる追加STFTの比を使う方が安定します。その式の符号はlibraryのSTFT規約で確認します。

再配置後のenergy mapは

\[ S_R(t',\omega')= \iint |V_x^h(t,\omega)|^2 \delta\{t'-\widehat t(t,\omega)\} \delta\{\omega'-\widehat\omega(t,\omega)\} \,dt\,d\omega \]

です。値が大きい点は、元のspectrogramで周囲に分散していたenergyがその点へ多く集められたことを示します。単位は元のspectrogram normalizationに従います。瞬時周波数は \(\widehat f=\widehat\omega/(2\pi)\) Hzです。標準的なreassigned spectrogramは時間と周波数の両方へenergyを動かし、transientとchirpを見やすくしますが、energy mapだけから元の位相を一般に復元しません。SSTは通常、周波数方向だけに複素CWT係数を再配置して成分再構成を保つ点が異なります。Auger & Flandrin (1995)

\(|V|\) が小さい点ではphase derivativeが不安定なため、\(\gamma\) を設けます。同じ窓内で複数成分が重なると、局所重心が実在しない中間位置へ移ることがあります。鋭い表示はGabor–Heisenbergの制約を破って新しい情報を得たことを意味せず、既存energyの配置を変えた結果です。window、窓長、hop、FFT長、\(\gamma\)、補間・binning、phase convention、端で除外する半窓長を報告します。STFT端点にはCOIという名称を用いません。

Wigner–Ville分布（WVD）は

\[ W_x(t,f)=\int_{-\infty}^{\infty} x\!\left(t+\frac{\tau}{2}\right) x^*\!\left(t-\frac{\tau}{2}\right) e^{-i2\pi f\tau}d\tau \]

で、高い時間–周波数集中と

\[ \int W_x(t,f)df=|x(t)|^2, \qquad \int W_x(t,f)dt=|X(f)|^2 \]

というmarginalを持ちます。しかし \(x=x_1+x_2\) では \(W_x=W_{x_1}+W_{x_2}+2\operatorname{Re}W_{x_1x_2}\) となり、2成分の間に振動するcross-termと負値を作ります。raw WVDの中間ridgeや負値を生理成分として読みません。smoothed pseudo-WVDなどの平滑化はcross-termを減らしますが、時間・周波数集中をぼかします。WVDは発展的な比較対象に留め、多成分生体信号の第一選択にはしません。解析信号を用いた時間–周波数分布の原型は Ville (1948) に示されています。

2.5.5 Matching pursuit

matching pursuit（MP）は、候補波形の集合dictionary \(\mathcal D=\{g_\gamma\}\) から、現在の残差と最もよく一致するatomを1つずつ選ぶsparse分解です。各atomを \(\|g_\gamma\|_2=1\) とし、

\[ R^0x=x \]

から反復 \(m=0,\ldots,M-1\) で

\[ \gamma_m= \arg\max_{\gamma\in\mathcal D} \left|\langle R^mx,g_\gamma\rangle\right|, \qquad a_m=\langle R^mx,g_{\gamma_m}\rangle \]

\[ R^{m+1}x=R^mx-a_mg_{\gamma_m} \]

と更新します。\(M\) 回後は

\[ x=\sum_{m=0}^{M-1}a_mg_{\gamma_m}+R^Mx \]

で、各stepで選択atomと新残差が直交するため

\[ \|x\|_2^2= \sum_{m=0}^{M-1}|a_m|^2+\|R^Mx\|_2^2 \]

が成り立ちます。Gabor dictionaryのatom例は

\[ g_\gamma[n]= \frac{1}{\sqrt s} g\!\left(\frac{n-u}{s}\right) e^{i2\pi\xi n/f_s}, \qquad \gamma=(u,s,\xi) \]

です。\(u\) は中心sample、\(s\) はscale、\(\xi\) はHzです。real signalをcomplex atomで表す場合は、実部または共役atomの扱いを実装に合わせます。

atomがdimensionlessかつunit normなら \(a_m\) は元信号と同じ単位、\(|a_m|^2\) は離散二乗normへの寄与で信号単位²です。物理時間を含むenergyは \(\Delta t|a_m|^2\) で信号単位²·sです。大きい \(|a_m|^2\) は、そのdictionary内で選ばれた時刻・scale・周波数のatomが残差をよく説明したことを意味します。時間–周波数上でatomが群を作ればburstやtransientの候補ですが、1つの生理源を表すとは限りません。

MPは少数の局在波形で記録を近似できる場合に有用です。一方、greedy選択は全atomの組合せに対するglobal最適解ではありません。似たatomが多いcoherent dictionaryでは選択が非一意になり、わずかなnoiseで時刻・scaleが入れ替わります。dictionaryにない波形は多数atomへ分裂し、停止を遅らせるとnoiseまで説明します。同じデータでatomを探索し、選ばれたatomだけに通常の未補正検定を行うとselection biasが生じます。Mallat & Zhang (1993)

dictionaryのwindow形、scale・frequency・time grid、atom normalization、MP/orthogonal MPの別、最大atom数 \(M\)、残差energyまたは係数threshold、境界atomの扱い、高速近似法を報告します。停止曲線 \(\|R^Mx\|_2^2/\|x\|_2^2\)、模擬atomの回収誤差、dictionaryを変えた安定性、hold-out dataでの再現性を確認します。

2.5.6 特異スペクトル解析

singular spectrum analysis（SSA、特異スペクトル解析）は、1次元時系列を遅延座標のHankel matrixへ埋め込み、singular value decomposition（SVD）で低rank成分へ分けます。\(x_0,\ldots,x_{N-1}\)、window長 \(L\)、\(2\le L\le N-1\)、\(K=N-L+1\) とし、

\[ \mathbf X_j= \begin{bmatrix} x_j&x_{j+1}&\cdots&x_{j+L-1} \end{bmatrix}^{\mathsf T}, \qquad j=0,\ldots,K-1 \]

\[ \mathbf X= \begin{bmatrix} \mathbf X_0&\mathbf X_1&\cdots&\mathbf X_{K-1} \end{bmatrix} =\mathbf U\boldsymbol\Sigma\mathbf V^{\mathsf T} =\sum_{i=1}^{d}\sigma_i\mathbf u_i\mathbf v_i^{\mathsf T} \]

とします。\(\mathbf X\) の各anti-diagonalは同じ時刻の値を持つためHankel matrixです。選んだindex集合 \(G\) の成分を

\[ \mathbf X^{(G)}= \sum_{i\in G}\sigma_i\mathbf u_i\mathbf v_i^{\mathsf T} \]

とgroup化し、anti-diagonal averagingで1次元へ戻します。\(n=0,\ldots,N-1\) に対して

\[ \widetilde x_n^{(G)}= \frac{1}{q_n} \sum_{\substack{1\le i\le L,\;1\le j\le K\\i+j=n+2}} X_{ij}^{(G)}, \qquad q_n= \#\{(i,j):i+j=n+2\} \]

です。\(\widetilde x_n^{(G)}\) は元信号と同じ単位です。\(\sigma_i\) も元信号と同じ単位ですが、scaleはwindow長と標本数に依存します。trajectory matrix内のenergy構成比

\[ \eta_i=\frac{\sigma_i^2}{\sum_j\sigma_j^2} \]

は無次元です。\(\sum_i\sigma_i^2=\|\mathbf X\|_F^2\) では中央のsampleほど複数回現れるため、\(\eta_i\) は埋め込みmatrixのenergy割合であり、元時系列の単純な分散割合と完全には同じではありません。

大きい \(\eta_i\) は、選んだ \(L\) の遅延空間でその成分が強いことを表します。低周波の滑らかなreconstructed componentはtrend候補、近い大きさの2つの特異値と位相が約90°ずれたeigenvector pairは振動候補です。しかしcolored noiseや短い自己相関系列も、滑らかな成分やpairを作ります。特異値の大きさだけでsignalとnoiseを分けず、surrogateまたはMonte Carlo null、screeの安定性、成分間weighted correlationを確認します。

SSAはtrend、準周期成分、noiseが低rank構造として分離できる系列に向きます。\(L\) が短すぎると遅い周期を表せず、長すぎると推定varianceと端点影響が増えます。近い周期、急に変わる周波数、短い記録、欠測は成分を混ぜます。端では \(q_n\) が小さく、diagonal averagingに使う要素が少ないため不確かさが大きくなります。window長とgroupingに普遍的な正解はありません。Vautard, Yiou & Ghil (1992) SSAにつながる遅延座標SVDの基礎は Broomhead & King (1986) に示されています。

\(L\)、demean/detrendとscaling、欠測処理、SVD solver、rankとgroupの選択法、surrogate nullと反復数、端点除外を報告します。複数の \(L\) とgroupingでreconstructed component、周期、\(\eta_i\) が保たれるかを感度分析します。

2.5.7 共通の検証と報告

1. 解析前に対象成分の周波数、持続時間、振幅変化、許容する時間変化を定義する。
2. sampling frequency、解析区間、detrend/filter、欠測・artifact maskを全手法で揃える。
3. 既知成分、colored noise、transient、端点を含む模擬信号で、周波数誤差、時間誤差、再構成誤差、偽成分を測る。
4. mode数、noise amplitude、wavelet、ridge幅、dictionary、SSA windowなど主要parameterを変え、結論が保たれる範囲を示す。
5. paddingされた範囲を結果から除き、CWT/SSTではCOI、STFT系では半窓長、EMDでは包絡端点、VMDでは周期境界、SSAでは小さい \(q_n\) を別々に扱う。
6. energy、振幅、PSD、割合を混同せず、normalization、片側/両側、bin幅、物理単位を示す。
7. 原信号、標準的なPSD/STFT/CWT、分解成分、残差または再構成を同じ時間軸で確認する。
8. 生理学的な名称は、独立した測定、既知の周波数・時刻、介入、再現性で裏付けられた場合だけ付ける。

2.5.8 LMS・NLMS・RLSによる適応雑音除去

適応雑音除去は、主sensor \(d[n]\) に混ざったartifact \(a[n]\) を、加速度、EOG、ECG、近傍noise電極等の参照signal \(u[n]\) から逐次予測して引く方法です。[一般的実務]。tap数 \(M\) の参照vectorを

\[ \mathbf u[n]=[u[n],u[n-1],\ldots,u[n-M+1]]^T \]

とすると、artifact推定と残差信号は

\[ \widehat a[n]=\mathbf w^T[n]\mathbf u[n],\qquad e[n]=d[n]-\widehat a[n] \]

です。least mean squares（LMS）は、目的関数を \(E[e^2]/2\) とする規約で

\[ \mathbf w[n+1]=\mathbf w[n]+\mu e[n]\mathbf u[n] \]

と更新します。理想的な定常入力で平均収束する目安は \(0<\mu<2/\lambda_{max}(R_u)\) ですが、実データで信号保存を保証する推奨値ではありません。normalized LMS（NLMS）は

\[ \mathbf w[n+1]=\mathbf w[n]+ \frac{\mu_N}{\epsilon+|\mathbf u[n]|_2^2} e[n]\mathbf u[n] \]

で、参照振幅の変化へ対応しやすくします。理想条件の数学的範囲は通常 \(0<\mu_N<2\) です。\(\mu_N\) が大きいほど急な伝達変化へ速く追従しますが、係数揺れと目的信号の除去も増え得ます。

recursive least squares（RLS）は、forgetting factor \(0<\lambda\le1\) を用いた

\[ J_n=\sum_{i=0}^{n}\lambda^{n-i}e^2[i] \]

を逐次最小化します。

\[ \xi[n]=d[n]-\mathbf w^T[n-1]\mathbf u[n] \]

\[ \mathbf k[n]=\frac{\mathbf P[n-1]\mathbf u[n]} {\lambda+\mathbf u^T[n]\mathbf P[n-1]\mathbf u[n]} \]

\[ \mathbf w[n]=\mathbf w[n-1]+\mathbf k[n]\xi[n],\qquad \mathbf P[n]=\lambda^{-1} \left[\mathbf P[n-1]-\mathbf k[n]\mathbf u^T[n]\mathbf P[n-1]\right] \]

\(\lambda\) が1に近いほど長い過去を保ち、近似的記憶長は \(1/(1-\lambda)\) sample、時間定数は \(-\Delta t/\ln\lambda\) です。小さい \(\lambda\) は速い変化へ追従しますがvarianceと数値的不安定性を増やします。

\(d,e,\widehat a\) は主signalと同じ単位、\(\mathbf w\) は主signal単位／参照signal単位、MSEは主signal単位²です。係数 \(\mathbf w\) は参照sensorから主sensorへ混入する統計的伝達関係の推定であり、生理結合係数ではありません。参照へ目的の脳反応、心拍反応、眼球運動等が入ると、それも除去されます。したがって、除去前後のevent波形・PSD、residual–reference coherence、係数軌跡、模擬混入での回収率を確認します。\(f_s\)、tap数、参照と単位、遅延、\(\mu/\mu_N,\epsilon,\lambda,P_0\)、更新停止規則、初期過渡除外を報告します（Widrow et al., 1975）。

2.5.9 状態空間model・Kalman filter・時変AR

状態空間modelは、直接見えない状態 \(\mathbf x_k\) が時間発展し、sensor値 \(\mathbf y_k\) を生むと表します。線形Gaussian形は

\[ \mathbf x_k=F_k\mathbf x_{k-1}+B_k\mathbf u_k+\mathbf w_k, \qquad \mathbf w_k\sim\mathcal N(0,Q_k) \]

\[ \mathbf y_k=H_k\mathbf x_k+\mathbf v_k, \qquad \mathbf v_k\sim\mathcal N(0,R_k) \]

です。Kalman filterの予測は

\[ \widehat{\mathbf x}_{k|k-1}=F_k\widehat{\mathbf x}_{k-1|k-1}+B_k\mathbf u_k, \qquad P_{k|k-1}=F_kP_{k-1|k-1}F_k^T+Q_k \]

で、innovationと更新は

\[ \boldsymbol\nu_k=\mathbf y_k-H_k\widehat{\mathbf x}_{k|k-1}, \quad S_k=H_kP_{k|k-1}H_k^T+R_k, \quad K_k=P_{k|k-1}H_k^TS_k^{-1} \]

\[ \widehat{\mathbf x}_{k|k}=\widehat{\mathbf x}_{k|k-1}+K_k\boldsymbol\nu_k \]

です。covarianceは数値的に安定なJoseph形

\[ P_{k|k}=(I-K_kH_k)P_{k|k-1}(I-K_kH_k)^T+K_kR_kK_k^T \]

で更新できます。\(Q\) を大きく置くと状態が速く変化し得ると仮定して観測へ追従し、\(R\) を大きく置くと観測を不確かとみなして予測を重くします。filterは現在までのdataだけを使うcausal推定、smootherは未来の観測も使うoffline推定です。normalized innovation squared

\[ NIS_k=\boldsymbol\nu_k^TS_k^{-1}\boldsymbol\nu_k \]

が連続して大きい場合、artifactだけでなく \(F,H,Q,R\) の不整合も疑います。出力はmodel仮定下の推定状態と不確かさであり、覚醒、stress、自律状態等そのものではありません。非線形modelにはEKF、UKF、particle filter等があります（Kalman, 1960）。

時間で変わるスペクトルを追う一例として、AR係数を状態にした

\[ x[n]+\sum_{p=1}^{P}a_p[n]x[n-p]=\varepsilon[n],\qquad \mathbf a[n]=\mathbf a[n-1]+\boldsymbol\eta[n] \]

をKalman/RLS等で推定し、

\[ S^{(2)}(f,n)=\frac{\sigma_\varepsilon^2[n]} {f_s\left|1+\sum_{p=1}^{P}a_p[n]e^{-i2\pi fp/f_s}\right|^2} \]

とします。これはinnovations varianceを信号単位²、\(f\)をHzとした両側PSDで、単位は信号単位²/Hzです。実数信号の片側PSDは、DCとNyquist以外を2倍します。\(Q\) またはforgettingを大きくするとpeak移動へ速く追従する一方、PSDが揺れます。AR次数、安定性制約、初期値、process-noise共分散\(Q\)とobservation-noise共分散\(R\)をそれぞれ、更新・表示間隔、片側／両側、credible/confidence interval、残差白色性を報告します。\(Q\)と\(R\)は次元・形が異なり得るため、無条件に比\(Q/R\)とは書きません。固定window法より滑らかな時間変化を得ても、急変を生理eventと自動同定しません。

2.5.10 PELT変化点検出とHidden Markov model

変化点検出は、指定した特徴の平均、variance、spectrum等が変わる時刻を探します。PELT等のpenalized segmentationは

\[ 0=\tau_0<\tau_1<\cdots<\tau_m<\tau_{m+1}=N \]

として

\[ \min_{\tau_{1:m}} \left[\sum_{j=0}^{m}\mathcal C \left(y_{\tau_j+1:\tau_{j+1}}\right)+\beta m\right] \]

を解きます。平均変化用の単純なcostは

\[ \mathcal C(y_{a:b})=\sum_{t=a}^{b}(y_t-\bar y_{a:b})^2 \]

です。この場合 \(\beta\) を大きくすると少数の長いsegment、小さくすると多数の短いsegmentになります。検出時刻は「このcostで統計特性が変わった時刻」であり、生理eventの時刻や介入の因果効果そのものではありません。filter端、artifact、missing、自己相関も変化点を作ります。特徴・単位、cost、penalty、最小segment長、offline/online、端点・missing処理、penalty感度とfalse-alarm検証を報告します（Killick et al., 2012）。

Hidden Markov model（HMM）は、観測 \(\mathbf y_t\) が見えない離散状態 \(z_t\) の間を遷移するとします。

\[ P(z_1=k)=\pi_k,\qquad P(z_t=j\mid z_{t-1}=i)=A_{ij},\qquad p(\mathbf y_t\mid z_t=k)=b_k(\mathbf y_t) \]

\[ p(\mathbf y_{1:T},z_{1:T})=\pi_{z_1}b_{z_1}(\mathbf y_1) \prod_{t=2}^{T}A_{z_{t-1},z_t}b_{z_t}(\mathbf y_t) \]

forward–backwardでstate posterior、Viterbiで一つの最尤state列を得ます。指標はposterior probability、fractional occupancy、transition probability、switching rate、dwell timeです。通常HMMではstate \(k\) の滞在長が幾何分布となり

\[ E[D_k]=\frac{1}{1-A_{kk}} \]

sampleです。window stepを掛ければ時間になります。滞在時間分布を直接表したい場合はhidden semi-Markov modelを検討します。state番号は計算上のclusterで、覚醒、眠気、交感優位等のlabelではありません。state数、feature、window/step、emission、標準化、初期値・seed、model選択、decoding、group間state対応、外部annotationでの妥当性を報告します（Baum et al., 1970）。

2.5.11 Robust PCAによる低rank・sparse分解

多channel×timeまたはepoch×feature行列を

\[ D=L+S+E \]

とし、反復的な低rank成分 \(L\)、まれな大振幅成分 \(S\)、小残差 \(E\) へ分ける方法です。stable principal component pursuitの代表形は

\[ \min_{L,S}\|L\|_*+\lambda\|S\|_1 \quad\text{subject to}\quad \|D-L-S\|_F\le\delta \]

です。\(\|L\|_*\) は特異値和、\(\|S\|_1\) は要素絶対値和です。\(S\) の絶対値やsparse率が大きいほど、この行列配置と \(\lambda,\delta\) の下で局所的な外れ成分を多く割り当てたことを表します。しかし \(L=\)生体信号、\(S=\)artifactとは自動的に対応しません。広く続くartifactは低rank側へ、短い真の脳・心臓反応はsparse側へ入ることがあります。行・列の意味、中心化・scale、\(\lambda,\delta\)、solver、停止条件、rank、sparse率、残差、成分採用規則と模擬回収を報告します（Candès et al., 2011）。

2.6.1 MF-DFA

multifractal detrended fluctuation analysis（MF-DFA）は、DFAを「小さい揺らぎと大きい揺らぎでscale則が違うか」まで拡張します。profile \(Y(k)=\sum_{i=1}^{k}(x_i-\bar x)\) を長さ \(s\) の区間へ分け、各区間 \(v\) で\(m\)次多項式trend \(p_{v,m}\) を除いて

\[ F^2(v,s)=\frac1s\sum_{j=1}^{s}\{Y_v(j)-p_{v,m}(j)\}^2, \qquad F_q(s)=\left[\frac1{2N_s}\sum_{v=1}^{2N_s}\{F^2(v,s)\}^{q/2}\right]^{1/q} \]

を計算します（余りを捨てないよう両端から区切るため \(2N_s\)、\(q=0\) はlog平均の極限）。\(F_q(s)\propto s^{h(q)}\) の傾き \(h(q)\)、\(\tau(q)=qh(q)-1\)、\(\alpha=d\tau/dq\)、\(f(\alpha)=q\alpha-\tau(q)\) が出力です。\(F_q\) はprofileと同じ信号単位、\(h,\tau,\alpha,f\) と幅 \(\Delta\alpha\) は無次元です。正の\(q\)は大振幅、負の\(q\)は小振幅を強調し、\(h(q)\)が一定ならmonofractal、\(\Delta\alpha\)が広いほどscale特性の多様性が大きい候補です。定常増分（fGn型）を想定できる場合の \(h(2)>0.5\)／\(<0.5\) は、採用scaleでの持続／反持続と整合しますが、生理機能の良否ではありません。短系列、外れ値、非多項式trend、crossover、負の\(q\)でのほぼ0の分散は偽の幅を作ります。\(N\)、標本化、\(s\)範囲と各scaleの区間数、\(m\)、\(q\) grid、fit法・\(R^2\)、crossover処理、端点、surrogateと有限長biasを必ず報告します（Kantelhardt et al., 2002）。

2.6.2 Cross-recurrence / joint recurrence

cross-recurrence plot（CRP）は、埋込み状態 \(\mathbf x_i=(x_i,x_{i+\tau_x},\ldots)\) と \(\mathbf y_j\) の近さを

\[ CR_{ij}=\Theta\!\left(\varepsilon-\|\mathbf x_i-\mathbf y_j\|\right) \]

で表します。対してjoint recurrence plot（JRP）は各信号内の再帰 \(R^x_{ij},R^y_{ij}\) を別々に作り、同期した同一時刻対について \(JR_{ij}=R^x_{ij}R^y_{ij}\) とします。CRPは「xの状態とyの状態が似る時刻対応」、JRPは「両系が同時に各自の過去状態へ戻ること」を問うため、交換できません。recurrence rate（1の割合）とDET（最小長 \(l_{min}\) 以上の対角線上にある点の割合）は無次元、線長とlagはsample、\(\Delta t\) を掛ければ秒です。高いrateは近い状態が多く、高いDET・長い対角線は似た時間発展が続く方向を示しますが、結合や因果を証明しません。CRPには比較可能な座標・埋込み次元、JRPには同期時刻が必要です。scale差、閾値、自己相関、drift、共通刺激だけでも線が増えます。埋込み次元・delay、norm、標準化、\(\varepsilon\)または固定recurrence rate、Theiler window、\(l_{min}\)、時間合わせ、欠測、surrogateを報告します（Marwan & Kurths, 2002、Romano et al., 2004）。

2.6.3 DCCA / MF-DCCA

detrended cross-correlation analysis（DCCA）は、同時系列 \(x,y\) の累積profileをscale \(s\) ごとに区切り、各区間でtrendを除いた残差 \(\tilde X_v(j),\tilde Y_v(j)\) の共分散

\[ C_v(s)=\frac1s\sum_{j=1}^{s}\tilde X_v(j)\tilde Y_v(j), \qquad F_{xy}^2(s)=\frac1{2N_s}\sum_vC_v(s), \qquad \rho_{DCCA}(s)=\frac{F_{xy}^2(s)}{\sqrt{F_{xx}^2(s)F_{yy}^2(s)}} \]

を求めます。\(F_{xy}^2\) はx単位×y単位、\(\rho_{DCCA}\in[-1,1]\) は無次元で、正負はscale \(s\) の同方向／逆方向の共変動、絶対値が大きいほど強い方向です。MF-DCCAでは一例として \(G_q(s)=(2N_s)^{-1}\sum_v\operatorname{sign}(C_v)|C_v|^{q/2}\) を作り、符号を保存して \(|G_q(s)|^{1/q}\propto s^{h_{xy}(q)}\) をfitします。\(|G_q|^{1/q}\) の単位は \(\sqrt{\text{x単位}\times\text{y単位}}\)、\(h_{xy}\) は無次元です。その\(q\)依存は、小振幅（負\(q\)）と大振幅（正\(q\)）でcross-scale則が違う候補で、直接的な生理結合量ではありません。両profileに共通するtrend、時刻ずれ、狭いscale範囲、共分散の符号反転、\(C_v\approx0\) での負\(q\)は推定を壊します。\(N\)、同期・resampling、scaleと区間法、detrend次数、\(q\) grid、符号・\(q=0\)規約、fit範囲・不確かさ、surrogate、多重比較を報告します（Podobnik & Stanley, 2008、Zhou, 2008）。

2.6.4 Dynamic time warping

形は似るが速度やdurationが異なるeventを比較するdynamic time warping（DTW）では、両端を対応させる標準的な規約なら \(D(0,0)=0\)、\(D(i,0)=D(0,j)=+\infty\)（\(i,j>0\)）と初期化し、局所距離 \(d(x_i,y_j)\) から \(i=1,\ldots,N_x\)、\(j=1,\ldots,N_y\) について

\[ D(i,j)=d(x_i,y_j)+\min\{D(i-1,j),D(i,j-1),D(i-1,j-1)\} \]

を動的計画法で計算します。最終cost \(D(N_x,N_y)\) またはpath長で規格化したcostが小さいほど、許した時間伸縮の範囲で二波形が似ています。\(d=|x_i-y_j|\) ならcostの単位は信号単位、\(d=(x_i-y_j)^2\) なら信号単位²、z-score後なら無次元です。path長で割ってもこの単位は変わりません。warping pathは、どの時点同士が対応したかを示します。open-begin/open-endを許す場合は上の境界条件が変わるため、Sakoe–Chiba band等の最大ずれ、step pattern、振幅標準化、端点条件、距離、path normalizationを報告します。制約が緩すぎると無関係な波形も似せ、実時間の遅延や因果関係を保存しません。EOG blink/saccade、ECG morphology、呼吸cycle等のtemplate照合に使えますが、標本化とfilterを揃え、学習templateはvalidation dataから独立に作ります。