Python 時系列分析 1,000本ノック– ノック81: Welch法 –

問題
答え
解説

次の Python コードの出力から何がわかりますか？

Python コード:

import numpy as np
from scipy.signal import welch

np.random.seed(42)
t = np.arange(500)
data = (
    100*np.sin(2*np.pi*t/50)
    +100*np.sin(2*np.pi*t/12)
    +np.random.randn(500))

freqs, psd = welch(data,fs=1,nperseg=256)
top = np.argsort(psd)[-5:][::-1]
print("Freq\tPSD\t Period")
for i in top:
    print(f"{freqs[i]:.3f}\t"
          f"{psd[i]:.1f}\t {1/freqs[i]:.2f}")

回答の選択肢:
(A) データに含まれるトレンドの傾きと切片がわかる
(B) データに含まれる周期成分とその強さがわかる
(C) データの定常性の有無がわかる
(D) データの自己相関が有意なラグ数がわかる

出力例:

Freq   PSD       Period
0.020  840112.0  51.20
0.082  738051.0  12.19
0.086  472590.2  11.64
0.023  297740.2  42.67
0.016  144390.3  64.00

正解: (B)
このコードは、Welch法によるパワースペクトル密度（PSD）推定を行い、時系列データに含まれる周期成分とその強度を周波数領域で分析しています。出力結果から、周期約51（周波数0.020）と周期約12（周波数0.082）に大きなピークが確認でき、これはデータ生成時に埋め込んだ2つのサイン波（周期50と周期12）に対応しています。

コードの解説

このコードは、複数の周期成分が混在する時系列データに対して、スペクトル解析で各周期成分の強さを特定しています。

import numpy as np
from scipy.signal import welch

np.random.seed(42)
t = np.arange(500)
data = (
    100*np.sin(2*np.pi*t/50)
    +100*np.sin(2*np.pi*t/12)
    +np.random.randn(500))

freqs, psd = welch(data,fs=1,nperseg=256)
top = np.argsort(psd)[-5:][::-1]
print("Freq\tPSD\t Period")
for i in top:
    print(f"{freqs[i]:.3f}\t"
          f"{psd[i]:.1f}\t {1/freqs[i]:.2f}")

簡単に説明します。

ライブラリのインポート

numpy：数値計算を行うためのライブラリ。ここでは乱数生成やサイン波の生成に使用しています。
scipy.signal.welch：Welch法によるパワースペクトル密度推定を行う関数です。

データの生成

100*np.sin(2*np.pi*t/50)：振幅100、周期50の正弦波です。
100*np.sin(2*np.pi*t/12)：振幅100、周期12の正弦波です。
np.random.randn(500)*0.1：標準偏差0.1のホワイトノイズです。
これら3つを足し合わせることで、2つの周期成分がノイズの中に埋め込まれたデータを生成しています。

Welch法の実行

welch(data, fs=1.0, nperseg=256)：Welch法でPSDを推定します。
- fs=1.0：サンプリング周波数を1.0に設定しています。これにより、返される周波数軸は0〜0.5の範囲になります。
- nperseg=256：各セグメントの長さを256に設定しています。セグメントごとにFFTを計算し、その平均を取ることでノイズの影響を抑えた推定を行います。
戻り値freqsは周波数の配列、psdは各周波数に対応するパワースペクトル密度の値です。

ピークの抽出

np.argsort(psd)[-5:][::-1]：PSDの値が大きい上位5つのインデックスを降順で取得しています。
1/freqs[i]：周波数を周期に変換しています（周期 = 1 / 周波数）。

選択肢の解説

(A) データに含まれるトレンドの傾きと切片がわかる → ×

スペクトル解析はデータを周波数成分に分解する手法であり、トレンドの傾きや切片を直接求めるものではありません。トレンドの推定には回帰分析やデトレンドなどの手法を使います。なお、トレンドがある場合は周波数0付近のPSDが大きくなることがありますが、傾きや切片を読み取ることはできません。

(B) データに含まれる周期成分とその強さがわかる → ○

正解です。PSDのピークの位置（周波数）から周期成分の周期がわかり、ピークの高さからその成分の強さ（パワー）がわかります。出力で最も大きいピーク（PSD=791.73、周期51.2）は振幅3の成分に、3番目のピーク（PSD=161.98、周期12.2）は振幅1.5の成分に対応しています。

(C) データの定常性の有無がわかる → ×

定常性の検定にはADF検定やKPSS検定を使います。スペクトル解析は定常性を仮定した上で周波数成分を分析する手法であり、定常性そのものを判定するものではありません。

(D) データの自己相関が有意なラグ数がわかる → ×

自己相関の分析にはACF（自己相関関数）やPACF（偏自己相関関数）を使います。ただし、スペクトル密度と自己相関関数はフーリエ変換を通じて理論的に等価な情報を含んでおり、この関係はウィーナー＝ヒンチンの定理として知られています。

Welch法とは

Welch法は、パワースペクトル密度を推定するための手法です。データ全体に対して一度にFFT（高速フーリエ変換）を計算するペリオドグラム法に対し、Welch法は以下の手順でノイズの影響を抑えます。

1. データをオーバーラップさせながら複数のセグメントに分割する（デフォルトで50%オーバーラップ）。
2. 各セグメントに窓関数（デフォルトでハニング窓）を適用する。
3. 各セグメントのペリオドグラムを計算する。
4. 全セグメントのペリオドグラムを平均する。

この平均化により、推定のばらつき（分散）が減少し、より信頼性の高いスペクトル推定が得られます。ただし、セグメントを短くするほど周波数分解能が低下するというトレードオフがあります。npersegを大きくすれば周波数分解能が上がりますが、平均に使えるセグメント数が減るため分散が増えます。

周波数と周期の関係

スペクトル解析の結果を読む際に重要な基本関係は以下の通りです。

$$
\text{周期} = \frac{1}{\text{周波数}} \quad , \quad \text{周波数} = \frac{1}{\text{周期}}
$$

今回のデータではfs=1.0としているため、周波数の単位は「サイクル/サンプル」です。例えば、周波数0.020は「約51サンプルで1回の振動」を意味し、データ生成時の周期50とほぼ一致します。

実際の時系列データでは、サンプリング周波数fsをデータの実際の収集頻度に合わせて設定します。例えば、月次データであればfs=12（年12回）と設定すると、周波数の単位が「サイクル/年」になり、ピークの位置から直接「何年周期か」を読み取ることができます。

また、解析可能な最大周波数はナイキスト周波数 $f_N = f_s / 2$ で決まります。これは「サンプリング間隔の2倍より短い周期は検出できない」ことを意味しています。