ピアノのA4の音は440Hzの振動です。しかし実際にピアノの音を録音して波形を見ると、単純な正弦波ではなく複雑な形をしています。この複雑な波形の中にどのような周波数成分がどのくらいの強さで含まれているかを分析する — これがスペクトル解析の本質です。
時系列分析は大きく分けて時間領域（time domain）のアプローチと周波数領域（frequency domain）のアプローチがあります。これまで見てきたAR、MA、ARIMAモデルは時間領域のアプローチですが、スペクトル解析は周波数領域のアプローチに属します。
スペクトル解析を理解することで、以下のような場面で威力を発揮します。
- 信号処理 — 通信信号やセンサーデータに含まれる特定の周波数成分の検出とフィルタリング
- 機械振動の診断 — 回転機械の故障診断において、特定の周波数に現れるピークから異常を検出する
- 気象データの周期性分析 — 気温や降水量の長期データに含まれるエルニーニョ周期や太陽活動周期の同定
本記事の内容
- 時間領域と周波数領域の関係
- パワースペクトル密度（PSD）の定義
- ウィーナー・ヒンチンの定理（自己相関↔PSD）
- ピリオドグラムの導出と統計的性質
- 平均化法（Bartlett法、Welch法）
- ARスペクトル推定（パラメトリック法）
- Python（scipy.signal）での実装
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
時間領域と周波数領域
なぜ周波数領域で考えるのか
時間領域の分析では、自己相関関数（ACF）やARモデルの係数を通じて、「過去の値と現在の値の依存関係」を捉えていました。これは時間軸に沿った逐次的な関係です。
一方、周波数領域の分析では、時系列を「様々な周波数の正弦波の重ね合わせ」として捉えます。各周波数成分がどのくらいの強さ（パワー）を持っているかを調べることで、時系列の周期的な構造を明らかにします。
日常的なアナロジーで言えば、音楽を聴いているとき、耳は自然と「低音がずっしりしている」「高音がキラキラしている」と感じます。これは脳が無意識にスペクトル分析を行っているようなものです。イコライザーの表示で低音域から高音域までの音量がバー表示されるのは、まさにスペクトル解析の視覚化です。
フーリエ解析との関係
フーリエ解析は、任意の関数を正弦波と余弦波の和として分解する数学的手法です。連続関数に対するフーリエ変換は次のように定義されます。
$$ X(f) = \int_{-\infty}^{\infty} x(t) e^{-i2\pi ft} dt $$
離散時系列 $\{x_0, x_1, \dots, x_{N-1}\}$ に対しては、離散フーリエ変換（DFT）を使います。
$$ \begin{equation} X_k = \sum_{n=0}^{N-1} x_n e^{-i2\pi kn/N}, \quad k = 0, 1, \dots, N-1 \end{equation} $$
$X_k$ は周波数 $f_k = k/N$（サンプリング周波数で正規化）における複素振幅を表します。$|X_k|^2$ がその周波数におけるパワー（エネルギー）の指標になります。
フーリエ解析が確定的な信号の周波数分解であるのに対し、スペクトル解析は確率過程の周波数構造を扱います。確率過程は実現するたびに異なる時系列を生成するため、個々のフーリエ係数ではなく、その統計的な性質（期待値）に着目するのがスペクトル解析の特徴です。
この統計的な観点から定義されるのが、次に見るパワースペクトル密度です。
パワースペクトル密度（PSD）の定義
直感的な理解
パワースペクトル密度（Power Spectral Density, PSD）$S(f)$ は、定常過程の「エネルギーが周波数軸上にどのように分布しているか」を表す関数です。
イメージとしては、PSDは時系列のエネルギーの「密度関数」です。確率密度関数が「確率が値の軸上にどう分布しているか」を示すように、PSDは「パワーが周波数軸上にどう分布しているか」を示します。
定義
弱定常過程 $\{x_t\}$ の自己共分散関数を $\gamma(h) = \text{Cov}(x_t, x_{t+h})$ とするとき、パワースペクトル密度は次のように定義されます。
$$ \begin{equation} S(f) = \sum_{h=-\infty}^{\infty} \gamma(h) e^{-i2\pi fh} \end{equation} $$
ここで $f \in [-1/2, 1/2]$ は正規化周波数です（サンプリング周波数の半分を1/2とする）。
この式は、自己共分散関数 $\gamma(h)$ の離散時間フーリエ変換（DTFT）に他なりません。つまり、PSDは自己共分散関数を周波数領域に変換したものです。
PSDの性質
PSDは以下の重要な性質を持ちます。
- 非負性: $S(f) \geq 0$（全ての $f$ に対して）
- 対称性: 実数値過程の場合 $S(f) = S(-f)$
- 全パワー: $\int_{-1/2}^{1/2} S(f) df = \gamma(0) = \text{Var}(x_t)$
性質3は特に重要です。PSDを全周波数にわたって積分すると、時系列の分散（= 全パワー）が得られます。これはパーセバルの定理の時系列版です。
PSDの定義は自己共分散関数のフーリエ変換として与えられました。逆に、PSDから自己共分散関数を復元できるでしょうか。その答えがウィーナー・ヒンチンの定理です。
ウィーナー・ヒンチンの定理
定理の内容
ウィーナー・ヒンチンの定理は、パワースペクトル密度と自己共分散関数がフーリエ変換対であることを述べます。
$$ \begin{equation} S(f) = \sum_{h=-\infty}^{\infty} \gamma(h) e^{-i2\pi fh} \end{equation} $$
$$ \begin{equation} \gamma(h) = \int_{-1/2}^{1/2} S(f) e^{i2\pi fh} df \end{equation} $$
つまり、$\gamma(h)$ と $S(f)$ は一対一に対応しています。時間領域の情報（自己共分散）と周波数領域の情報（PSD）は、フーリエ変換を通じて等価なのです。
AR(1)過程のPSD
具体例としてAR(1)過程 $x_t = \phi x_{t-1} + \varepsilon_t$（$|\phi| < 1$, $\varepsilon_t \sim \text{WN}(0, \sigma^2)$）のPSDを導出してみましょう。
AR(1)の自己共分散関数は $\gamma(h) = \sigma^2 \phi^{|h|} / (1 – \phi^2)$ です。これをPSDの定義式に代入します。
$$ S(f) = \sum_{h=-\infty}^{\infty} \frac{\sigma^2 \phi^{|h|}}{1 – \phi^2} e^{-i2\pi fh} $$
$h \geq 0$ と $h < 0$ の和に分けます。
$$ S(f) = \frac{\sigma^2}{1-\phi^2}\left[\sum_{h=0}^{\infty} (\phi e^{-i2\pi f})^h + \sum_{h=1}^{\infty} (\phi e^{i2\pi f})^h\right] $$
各等比級数の公式を適用すると（$|\phi| < 1$ なので収束保証あり）、
$$ S(f) = \frac{\sigma^2}{1-\phi^2}\left[\frac{1}{1 – \phi e^{-i2\pi f}} + \frac{\phi e^{i2\pi f}}{1 – \phi e^{i2\pi f}}\right] $$
2つの項を通分して整理すると、
$$ \begin{equation} S(f) = \frac{\sigma^2}{|1 – \phi e^{-i2\pi f}|^2} = \frac{\sigma^2}{1 – 2\phi\cos(2\pi f) + \phi^2} \end{equation} $$
$\phi > 0$ の場合、$f = 0$（低周波数）でPSDが最大となり、高周波数に向かって減少します。これはAR(1)過程が「ゆっくりした変動（低周波成分）が支配的」であることを意味し、正の自己相関を持つことと整合しています。
逆に $\phi < 0$ の場合、$f = 1/2$（最高周波数）でPSDが最大となり、高周波成分（急激な振動）が支配的です。
MA(q)過程のPSD
一般のMA(q)過程 $x_t = \sum_{j=0}^{q} \theta_j \varepsilon_{t-j}$（$\theta_0 = 1$）のPSDは、
$$ \begin{equation} S(f) = \sigma^2 \left|\sum_{j=0}^{q} \theta_j e^{-i2\pi fj}\right|^2 = \sigma^2 |\Theta(e^{-i2\pi f})|^2 \end{equation} $$
MA多項式 $\Theta(z) = 1 + \theta_1 z + \cdots + \theta_q z^q$ を単位円上（$z = e^{-i2\pi f}$）で評価した絶対値の2乗に $\sigma^2$ をかけたものです。
ウィーナー・ヒンチンの定理により、PSDの推定は自己共分散の推定と等価であることがわかりました。しかし、実際にPSDを推定するにはどうすればよいでしょうか。最も直接的な方法がピリオドグラムです。
ピリオドグラムの導出と統計的性質
ピリオドグラムの定義
長さ $N$ の時系列データ $\{x_0, x_1, \dots, x_{N-1}\}$ に対して、ピリオドグラムは次のように定義されます。
$$ \begin{equation} I(f_k) = \frac{1}{N} |X_k|^2 = \frac{1}{N} \left|\sum_{n=0}^{N-1} x_n e^{-i2\pi kn/N}\right|^2 \end{equation} $$
ここで $f_k = k/N$（$k = 0, 1, \dots, N-1$）はフーリエ周波数です。
ピリオドグラムは「データのDFTの絶対値の2乗を $N$ で割ったもの」であり、各周波数におけるパワーの推定量です。FFT（高速フーリエ変換）を用いると $O(N \log N)$ で効率的に計算できます。
ピリオドグラムとPSDの関係
ピリオドグラム $I(f_k)$ はPSD $S(f_k)$ の推定量として期待されますが、重要な統計的性質があります。
$N \to \infty$ のとき、ピリオドグラムの期待値は、
$$ \mathbb{E}[I(f_k)] \to S(f_k) $$
つまり、ピリオドグラムは漸近的に不偏です。しかし、ピリオドグラムの分散は $N$ が大きくなっても減少しません。
$$ \text{Var}[I(f_k)] \approx S(f_k)^2 \quad (N \to \infty) $$
分散がPSDの2乗に比例し、$N$ に依存しないのです。これは、ピリオドグラムが不一致推定量（inconsistent estimator）であることを意味します。サンプルサイズを増やしても推定の精度が向上しないのです。
直感的には、サンプルサイズ $N$ を増やすとフーリエ周波数の解像度が上がり、推定すべき点の数も $N$ に比例して増えるため、各点あたりの情報量は変わらないのです。
ピリオドグラムの問題点の可視化
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(42)
# AR(1)過程の生成
phi = 0.8
sigma2 = 1.0
N = 512
eps = np.random.normal(0, np.sqrt(sigma2), N)
x = np.zeros(N)
x[0] = eps[0]
for t in range(1, N):
x[t] = phi * x[t-1] + eps[t]
# ピリオドグラムの計算
freqs = np.fft.rfftfreq(N)
X = np.fft.rfft(x)
periodogram = (1.0 / N) * np.abs(X) ** 2
# 理論PSD
S_theory = sigma2 / (1 - 2*phi*np.cos(2*np.pi*freqs) + phi**2)
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
# 線形スケール
axes[0].plot(freqs, periodogram, linewidth=0.5, alpha=0.7, label='Periodogram')
axes[0].plot(freqs, S_theory, linewidth=2, color='red', label='True PSD')
axes[0].set_xlabel('Frequency')
axes[0].set_ylabel('Power')
axes[0].set_title('Periodogram vs True PSD (Linear)')
axes[0].legend()
# 対数スケール
axes[1].semilogy(freqs, periodogram, linewidth=0.5, alpha=0.7, label='Periodogram')
axes[1].semilogy(freqs, S_theory, linewidth=2, color='red', label='True PSD')
axes[1].set_xlabel('Frequency')
axes[1].set_ylabel('Power (log)')
axes[1].set_title('Periodogram vs True PSD (Log)')
axes[1].legend()
plt.tight_layout()
plt.show()
ピリオドグラムの結果を見ると、理論PSD（赤線）のまわりにピリオドグラム（青線）が大きくばらついていることがわかります。
- 全体的な傾向は捉えている — 低周波にパワーが集中する（$\phi > 0$ のAR(1)の特徴）傾向は見えています。
- ばらつきが非常に大きい — ピリオドグラムは真のPSDのまわりに激しく振動しており、個々の周波数でのPSD推定としては信頼性が低いです。
- $N$ を増やしても改善しない — これがピリオドグラムが不一致推定量であることの帰結です。
この問題を解決するために、次のセクションではピリオドグラムの平滑化手法を見ていきます。
平均化法（Bartlett法・Welch法）
Bartlett法
ダニエル・バートレットは、ピリオドグラムの分散を減らすための単純かつ効果的な方法を提案しました。アイデアは「データを複数のセグメントに分割し、各セグメントのピリオドグラムを平均する」というものです。
長さ $N$ のデータを $L$ 個の重なりのないセグメント（各長さ $M = N/L$）に分割し、各セグメントのピリオドグラム $I_l(f)$（$l = 1, \dots, L$）を計算して平均します。
$$ \begin{equation} \hat{S}_{\text{Bartlett}}(f) = \frac{1}{L} \sum_{l=1}^{L} I_l(f) \end{equation} $$
各 $I_l(f)$ は独立な推定量なので、平均により分散が $1/L$ に減少します。
$$ \text{Var}[\hat{S}_{\text{Bartlett}}(f)] \approx \frac{S(f)^2}{L} $$
ただし、セグメント長 $M$ が短くなると周波数解像度が $1/M$ に低下するため、分散の低減と周波数解像度のトレードオフが生じます。
Welch法
ピーター・ウェルチはバートレット法をさらに改良し、以下の2つの工夫を加えました。
- セグメントのオーバーラップ — 隣接するセグメントを50%程度重ねることで、同じデータ長からより多くのセグメントを取り出す（独立性はやや低下するが、データの利用効率が向上）
- 窓関数の適用 — 各セグメントにハニング窓やハミング窓を適用し、スペクトルの漏れ（spectral leakage）を抑制
$$ \begin{equation} \hat{S}_{\text{Welch}}(f) = \frac{1}{L} \sum_{l=1}^{L} \frac{1}{M \cdot U} \left|\sum_{n=0}^{M-1} w_n x_{l,n} e^{-i2\pi fn/M}\right|^2 \end{equation} $$
ここで $w_n$ は窓関数、$U = (1/M) \sum_{n=0}^{M-1} w_n^2$ は窓関数のパワー正規化定数です。
Welch法の実装と確認
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
np.random.seed(42)
# AR(1)過程の生成
phi = 0.8
sigma2 = 1.0
N = 2048
eps = np.random.normal(0, np.sqrt(sigma2), N)
x = np.zeros(N)
x[0] = eps[0]
for t in range(1, N):
x[t] = phi * x[t-1] + eps[t]
# 理論PSD
freqs_theory = np.linspace(0, 0.5, 500)
S_theory = sigma2 / (1 - 2*phi*np.cos(2*np.pi*freqs_theory) + phi**2)
# 1. ピリオドグラム（素のFFT）
freqs_raw = np.fft.rfftfreq(N)
X = np.fft.rfft(x)
periodogram = (1.0 / N) * np.abs(X) ** 2
# 2. Welch法（異なるセグメント長）
nperseg_values = [64, 256, 512]
fig, axes = plt.subplots(2, 2, figsize=(12, 8))
# 素のピリオドグラム
axes[0, 0].semilogy(freqs_raw, periodogram, linewidth=0.3, alpha=0.7)
axes[0, 0].semilogy(freqs_theory, S_theory, linewidth=2, color='red')
axes[0, 0].set_title(f'Raw Periodogram (N={N})')
axes[0, 0].set_xlabel('Frequency')
axes[0, 0].set_ylabel('PSD')
# Welch法 各セグメント長
for idx, nperseg in enumerate(nperseg_values):
freqs_w, psd_w = signal.welch(x, nperseg=nperseg, noverlap=nperseg//2)
ax = axes[(idx + 1) // 2, (idx + 1) % 2]
ax.semilogy(freqs_w, psd_w, linewidth=1.5, color='steelblue')
ax.semilogy(freqs_theory, S_theory, linewidth=2, color='red', alpha=0.7)
n_segments = (N - nperseg) // (nperseg // 2) + 1
ax.set_title(f'Welch (nperseg={nperseg}, ~{n_segments} segments)')
ax.set_xlabel('Frequency')
ax.set_ylabel('PSD')
plt.tight_layout()
plt.show()
Welch法の結果から、分散と解像度のトレードオフが視覚的に理解できます。
- 素のピリオドグラム（左上）: 周波数解像度は最高ですが、激しいばらつきのため真のPSDの形が読み取りにくいです。
- nperseg=64（右上）: 多数のセグメントに分割するため分散が大幅に減少し、滑らかな推定が得られます。ただし、周波数解像度は低く、細かい周波数構造が見えなくなっています。
- nperseg=256（左下）: 中間的なバランスで、分散の低減と解像度の維持のバランスが良好です。
- nperseg=512（右下）: 解像度は高いですが、セグメント数が少ないため分散がやや大きくなっています。
実務では、nperseg=256程度（オーバーラップ50%）が良いバランスとして使われることが多いです。
Welch法は非パラメトリックなPSD推定法ですが、データがARモデルに従うことが事前にわかっている場合、パラメトリックなアプローチがより効率的です。次にARスペクトル推定を見ていきましょう。
ARスペクトル推定
アイデア
AR(p)過程のPSDは、先ほど導出したAR(1)の結果を一般化すると、
$$ \begin{equation} S(f) = \frac{\sigma^2}{|\Phi(e^{-i2\pi f})|^2} = \frac{\sigma^2}{|1 – \phi_1 e^{-i2\pi f} – \phi_2 e^{-i4\pi f} – \cdots – \phi_p e^{-i2\pi pf}|^2} \end{equation} $$
という閉じた形で表されます。
したがって、データからARモデルのパラメータ $\phi_1, \dots, \phi_p$ と $\sigma^2$ を推定すれば、PSDが直ちに得られます。これがARスペクトル推定の基本的なアイデアです。
利点と注意点
ARスペクトル推定の利点は以下の通りです。
- 少ないデータでも滑らかなPSD推定が得られる — パラメータ数が $p + 1$ と少ないため、分散が小さい
- 周波数解像度がデータ長に制約されない — 任意の細かさで $S(f)$ を評価できる
- 鋭いピークの検出に優れる — ピリオドグラムやWelch法ではぼやけてしまう鋭い周波数ピークを検出できる
注意点としては、データが実際にAR過程に従わない場合（例えばMA成分が含まれる場合）、推定にバイアスが生じる可能性があることです。
実装と比較
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from statsmodels.tsa.ar_model import AutoReg
np.random.seed(42)
# 2つの正弦波 + ノイズ（AR過程ではないデータ）
N = 512
t = np.arange(N)
f1, f2 = 0.1, 0.15 # 近接した2つの周波数
x = 2 * np.sin(2 * np.pi * f1 * t) + 1.5 * np.sin(2 * np.pi * f2 * t)
x += np.random.normal(0, 1, N)
# 1. ピリオドグラム
freqs_raw = np.fft.rfftfreq(N)
X = np.fft.rfft(x)
periodogram = (1.0 / N) * np.abs(X) ** 2
# 2. Welch法
freqs_w, psd_w = signal.welch(x, nperseg=128, noverlap=64)
# 3. ARスペクトル推定
ar_order = 30
ar_model = AutoReg(x, lags=ar_order).fit()
ar_coeffs = ar_model.params[1:] # AR係数
sigma2_ar = ar_model.sigma2
freqs_ar = np.linspace(0, 0.5, 1000)
z = np.exp(-1j * 2 * np.pi * freqs_ar[:, np.newaxis] * np.arange(1, ar_order + 1))
Phi = 1 - np.sum(ar_coeffs * z, axis=1)
psd_ar = sigma2_ar / np.abs(Phi) ** 2
fig, axes = plt.subplots(3, 1, figsize=(10, 10))
axes[0].plot(freqs_raw, 10 * np.log10(periodogram), linewidth=0.5, alpha=0.7)
axes[0].axvline(x=f1, color='red', linestyle='--', alpha=0.5, label=f'f={f1}')
axes[0].axvline(x=f2, color='green', linestyle='--', alpha=0.5, label=f'f={f2}')
axes[0].set_title('Periodogram')
axes[0].set_ylabel('PSD (dB)')
axes[0].legend()
axes[1].plot(freqs_w, 10 * np.log10(psd_w), linewidth=1.5, color='steelblue')
axes[1].axvline(x=f1, color='red', linestyle='--', alpha=0.5, label=f'f={f1}')
axes[1].axvline(x=f2, color='green', linestyle='--', alpha=0.5, label=f'f={f2}')
axes[1].set_title('Welch Method (nperseg=128)')
axes[1].set_ylabel('PSD (dB)')
axes[1].legend()
axes[2].plot(freqs_ar, 10 * np.log10(psd_ar), linewidth=1.5, color='darkorange')
axes[2].axvline(x=f1, color='red', linestyle='--', alpha=0.5, label=f'f={f1}')
axes[2].axvline(x=f2, color='green', linestyle='--', alpha=0.5, label=f'f={f2}')
axes[2].set_title(f'AR Spectral Estimation (order={ar_order})')
axes[2].set_ylabel('PSD (dB)')
axes[2].set_xlabel('Frequency')
axes[2].legend()
plt.tight_layout()
plt.show()
3つの手法の比較から、各推定法の特徴が明確にわかります。
- ピリオドグラム（上段）: 2つの正弦波の周波数 $f_1 = 0.1$ と $f_2 = 0.15$ に鋭いピークが見られますが、ノイズによるスプリアスなピークも多数存在しています。
- Welch法（中段）: セグメント長を128とすると、ピークのばらつきは減りますが、周波数解像度が低下して2つのピークが分離しにくくなる場合があります。
- ARスペクトル推定（下段）: 2つの周波数ピークが滑らかに分離されています。ARモデルのパラメトリックな仮定を活用することで、高い周波数解像度と滑らかさを両立しています。
実データでのスペクトル解析
最後に、実際の時系列データのスペクトル特性を分析してみましょう。
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import signal
from statsmodels.datasets import get_rdataset
# AirPassengersデータ
data = get_rdataset("AirPassengers").data
data.columns = ['time', 'passengers']
passengers = data['passengers'].values
# 対数変換 + トレンド除去（1次差分）
log_pass = np.log(passengers)
detrended = np.diff(log_pass)
N = len(detrended)
# Welch法
freqs_w, psd_w = signal.welch(detrended, nperseg=min(64, N//2),
noverlap=min(32, N//4))
fig, axes = plt.subplots(2, 2, figsize=(12, 8))
# 元データ
axes[0, 0].plot(passengers, linewidth=1.5)
axes[0, 0].set_title('AirPassengers (Original)')
axes[0, 0].set_xlabel('Month')
# トレンド除去後
axes[0, 1].plot(detrended, linewidth=1)
axes[0, 1].set_title('Δ log(Passengers)')
axes[0, 1].set_xlabel('Month')
# ピリオドグラム
freqs_raw = np.fft.rfftfreq(N)
X = np.fft.rfft(detrended)
periodogram = (1.0 / N) * np.abs(X) ** 2
axes[1, 0].plot(freqs_raw, periodogram, linewidth=0.8)
axes[1, 0].axvline(x=1/12, color='red', linestyle='--', alpha=0.7,
label='1/12 (yearly)')
axes[1, 0].axvline(x=1/6, color='green', linestyle='--', alpha=0.7,
label='1/6 (semi-annual)')
axes[1, 0].set_title('Periodogram')
axes[1, 0].set_xlabel('Frequency (cycles/month)')
axes[1, 0].legend(fontsize=8)
# Welch法
axes[1, 1].plot(freqs_w, psd_w, linewidth=1.5, color='steelblue')
axes[1, 1].axvline(x=1/12, color='red', linestyle='--', alpha=0.7,
label='1/12 (yearly)')
axes[1, 1].axvline(x=1/6, color='green', linestyle='--', alpha=0.7,
label='1/6 (semi-annual)')
axes[1, 1].set_title('Welch PSD Estimate')
axes[1, 1].set_xlabel('Frequency (cycles/month)')
axes[1, 1].legend(fontsize=8)
plt.tight_layout()
plt.show()
AirPassengersデータのスペクトル解析結果から、重要な周期的構造が読み取れます。
- 周波数 $1/12$ に明確なピーク — これは12か月周期、すなわち年次の季節パターンに対応します。AirPassengersデータが毎年同じ月に同様のパターンを示すことの周波数領域での反映です。
- 周波数 $1/6$ にもピークが見られる — これは6か月の半年周期の成分で、年次パターンの第2高調波（倍音）に対応します。
- 低周波成分は差分によりほぼ除去されている — 1次差分をとったことでトレンドが除去され、低周波のパワーが抑制されています。
まとめ
本記事では、時系列のスペクトル解析の理論と実装を詳しく解説しました。
- パワースペクトル密度（PSD）: 時系列のエネルギーが周波数軸上にどのように分布しているかを示す関数。自己共分散関数のフーリエ変換として定義される
- ウィーナー・ヒンチンの定理: PSDと自己共分散関数がフーリエ変換対であるという基本定理。時間領域と周波数領域の解析が等価であることを保証する
- ピリオドグラム: PSDの最も直接的な推定量だが、不一致推定量であり分散が大きい
- Welch法: セグメント分割・オーバーラップ・窓関数により分散を低減する標準的な非パラメトリック手法。分散と周波数解像度のトレードオフがある
- ARスペクトル推定: ARモデルの仮定を活用したパラメトリック手法。少ないデータでも滑らかな推定が可能で、鋭いピーク検出に優れる
スペクトル解析は時系列の「見えない周期構造」を可視化する強力なツールです。時間領域の分析と組み合わせることで、データのより深い理解が得られます。
次のステップとして、以下の記事も参考にしてください。