ラジオを聴いているとき、チューニングを回すと異なる放送局の音声が聞こえてきます。これは、電波のエネルギーが周波数ごとに異なる「場所」に分布しているからです。では、ある信号のエネルギーが具体的にどの周波数にどれだけ集中しているのか — これを定量的に記述するのがパワースペクトル密度（Power Spectral Density, PSD）です。
PSD はあらゆる信号処理の基盤となる概念であり、その応用範囲は驚くほど広いです。
- 通信工学: 受信信号のスペクトルを調べて帯域外雑音を除去するフィルタを設計する
- 振動解析: 構造物の振動データから共振周波数を特定し、疲労破壊を予防する
- 天文学: 宇宙背景放射のスペクトルを解析して宇宙の温度分布を推定する
- 脳科学: 脳波（EEG）のアルファ波、ベータ波をPSD で分離して認知状態を判定する
PSD の理論の中核を成すのがウィーナー・ヒンチンの定理です。この定理は「PSD は自己相関関数のフーリエ変換に等しい」という驚くべきシンプルな関係を述べており、時間領域と周波数領域を橋渡しする最も重要な結果の一つです。
本記事の内容
- エネルギースペクトル密度からPSDへの導入
- 確定信号とランダム信号に対するPSDの定義
- 自己相関関数の定義と性質
- ウィーナー・ヒンチンの定理の厳密な導出
- 白色雑音のPSD
- ピリオドグラムとウェルチ法によるPSD推定
- Pythonによる各種信号のPSD計算と比較
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- フーリエ変換の定義と性質を完全解説 — フーリエ変換の定義、パーシバルの等式
- 確率過程入門 — 確率過程、定常性、期待値
エネルギースペクトル密度 — PSDへの足がかり
信号のエネルギーとパーシバルの等式
PSD を理解するための第一歩は、信号の「エネルギー」を周波数の世界で表現することです。オーディオ信号であれば、エネルギーが大きいほど音が大きいというイメージが自然でしょう。
連続時間信号 $x(t)$ のエネルギー $E$ は次のように定義されます。
$$ E = \int_{-\infty}^{\infty} |x(t)|^2 \, dt $$
この量が有限であるとき、$x(t)$ をエネルギー信号（energy signal）と呼びます。パルス信号や減衰振動のように「いつかは消える」信号が典型例です。
一方、フーリエ変換の世界ではパーシバルの等式（Parseval’s theorem）が成り立ちます。信号 $x(t)$ のフーリエ変換を $X(f)$ とすると
$$ \int_{-\infty}^{\infty} |x(t)|^2 \, dt = \int_{-\infty}^{\infty} |X(f)|^2 \, df $$
ここで $f$ は周波数（Hz）です。左辺は時間領域でのエネルギー、右辺は周波数領域でのエネルギーを表しています。パーシバルの等式は「エネルギーの総量は時間領域で計算しても周波数領域で計算しても同じ」ということを意味します。物理的に言えば、信号をフーリエ変換しても、エネルギーが増えたり減ったりすることはないのです。
エネルギースペクトル密度の定義
パーシバルの等式の右辺に着目しましょう。被積分関数 $|X(f)|^2$ は「周波数 $f$ あたりのエネルギーの密度」と解釈できます。これをエネルギースペクトル密度（Energy Spectral Density, ESD）と呼びます。
$$ S_E(f) = |X(f)|^2 $$
ESD を $f_1$ から $f_2$ まで積分すれば、その周波数帯域に含まれるエネルギーが得られます。
$$ E_{[f_1, f_2]} = \int_{f_1}^{f_2} S_E(f) \, df $$
たとえば、短いパルス信号のフーリエ変換は広い周波数帯域にわたって広がります。これは直感にも合います。短いパルスはあらゆる周波数の波の重ね合わせで構成されるため、エネルギーが周波数方向に分散するのです。
ESD は「エネルギー信号」に対しては完全に定義が明確です。しかし、現実の信号処理で扱いたい信号の多くは、ESD だけでは不十分です。永遠に続く正弦波や、終わりのないランダム雑音のことを考えましょう。こうした信号はエネルギーが無限大になるため、フーリエ変換が通常の意味では存在しません。このギャップを埋めるのが「パワースペクトル密度」です。
パワースペクトル密度（PSD）の定義
エネルギー信号からパワー信号へ
永遠に鳴り続けるサイレンの音を考えてみましょう。この信号のエネルギーは明らかに無限大です。しかし、1秒あたりのエネルギー — つまり平均パワーならば有限値になります。エネルギー信号とパワー信号の違いは、バケツの水の総量と蛇口から出る水の流量の違いに似ています。バケツの水は有限ですが、蛇口を開けっぱなしにすれば水量は無限に増えます。しかし、毎秒の流量は一定です。
信号 $x(t)$ の平均パワー $P$ は次のように定義されます。
$$ P = \lim_{T \to \infty} \frac{1}{2T} \int_{-T}^{T} |x(t)|^2 \, dt $$
この量が $0 < P < \infty$ であるとき、$x(t)$ をパワー信号（power signal）と呼びます。定常的な正弦波やランダム雑音がこれに該当します。
確定信号のPSD
パワー信号に対しても、ESD のように「周波数ごとのパワーの分布」を知りたいわけです。そこで、信号を時間区間 $[-T, T]$ で切り出した窓掛け信号 $x_T(t)$ を考えます。
$$ x_T(t) = \begin{cases} x(t) & |t| \leq T \\ 0 & |t| > T \end{cases} $$
$x_T(t)$ は有限区間の信号なのでエネルギーが有限であり、フーリエ変換 $X_T(f)$ が存在します。このとき、$x_T(t)$ のESD は $|X_T(f)|^2$ です。時間区間 $2T$ で割って正規化し、$T \to \infty$ の極限をとれば、「単位時間あたりのスペクトル密度」が得られます。これがPSD です。
$$ \begin{equation} S_x(f) = \lim_{T \to \infty} \frac{|X_T(f)|^2}{2T} \end{equation} $$
この定義が意味するところを確認しましょう。$|X_T(f)|^2$ は時間幅 $2T$ 内でのエネルギースペクトル密度であり、$2T$ で割ることで1秒あたりのパワーの周波数分布に変換しています。$T \to \infty$ の極限は、十分に長い時間にわたって観測したときの「定常的な」スペクトル分布を取り出す操作です。
PSD を全周波数にわたって積分すると、平均パワーが得られます。
$$ P = \int_{-\infty}^{\infty} S_x(f) \, df $$
これはパーシバルの等式の「パワー版」であり、PSD が確かにパワーの周波数分布を表していることを裏付けています。
ランダム信号のPSD
通信やセンサ計測で扱う雑音信号は確定的ではなく確率的です。毎回異なる実現値をとるランダム過程 $X(t)$ に対して、先ほどの確定信号の定義をそのまま適用すると、実現値ごとにPSD が変わってしまいます。
ランダム過程に対するPSD は、アンサンブル平均（期待値）を導入して定義します。
$$ \begin{equation} S_X(f) = \lim_{T \to \infty} \frac{E\left[|X_T(f)|^2\right]}{2T} \end{equation} $$
ここで $E[\cdot]$ は確率的な期待値であり、$X_T(f)$ は確率過程 $X(t)$ を $[-T, T]$ で窓掛けした信号のフーリエ変換です。期待値をとることで、個々の実現値に依存しないスペクトル特性が得られます。
この定義は、広義定常過程（Wide-Sense Stationary, WSS）の場合に特に扱いやすくなります。WSS 過程とは、平均値が時間に依存せず、自己相関関数が時間差 $\tau$ のみの関数であるような過程です。WSS 過程ではウィーナー・ヒンチンの定理が成り立ち、PSD を自己相関関数のフーリエ変換として直接計算できます。
PSD の定義を押さえたところで、次に自己相関関数を導入しましょう。PSD と自己相関関数の間の深い関係こそが、ウィーナー・ヒンチンの定理の核心です。
自己相関関数の定義と性質
自己相関関数の直感的な意味
ある信号を少しだけ時間方向にずらして、元の信号と重ね合わせることを想像してみてください。ずらした量が小さければ、信号はほとんど同じ形なので「よく似ている」でしょう。ずらす量が大きくなると、信号の形が異なってきて「似ていない」と感じるはずです。この「時間ずれに対する類似度」を定量化するのが自己相関関数です。
ゆっくりと変化する低周波信号は、少し時間をずらしてもあまり形が変わらないため、自己相関が広い範囲で大きい値をとります。逆に、高周波成分が多い信号は、わずかな時間ずれでも急速に自己相関が減衰します。このように、自己相関の「幅」と周波数成分の「帯域幅」は反比例の関係にあります — これは後に導出するウィーナー・ヒンチンの定理を通じて厳密に理解できます。
確定信号の自己相関関数
確定的なパワー信号 $x(t)$ に対する自己相関関数は次のように定義されます。
$$ R_x(\tau) = \lim_{T \to \infty} \frac{1}{2T} \int_{-T}^{T} x(t) \, x^*(t – \tau) \, dt $$
ここで $x^*$ は複素共役を表します。実信号の場合は $x^* = x$ です。$\tau$ は時間のずれ（ラグ）を表すパラメータです。
$\tau = 0$ のとき、自己相関関数は信号自身との内積になるため
$$ R_x(0) = \lim_{T \to \infty} \frac{1}{2T} \int_{-T}^{T} |x(t)|^2 \, dt = P $$
となり、信号の平均パワーに等しくなります。
ランダム過程の自己相関関数
WSS 確率過程 $X(t)$ に対する自己相関関数は
$$ R_X(\tau) = E\left[X(t)\, X^*(t – \tau)\right] $$
と定義されます。WSS 過程の仮定により、この値は $t$ に依存せず $\tau$ のみの関数です。
自己相関関数の基本性質
自己相関関数にはいくつかの重要な性質があります。これらは後の導出やPSD の解釈で繰り返し使われます。
性質1: 対称性
$$ R_x(-\tau) = R_x^*(\tau) $$
実信号の場合は $R_x(-\tau) = R_x(\tau)$ であり、自己相関関数は偶関数です。
この性質を示しましょう。定義から
$$ R_x(-\tau) = E[X(t) X^*(t + \tau)] $$
$s = t + \tau$ と変数変換すると
$$ R_x(-\tau) = E[X(s – \tau) X^*(s)] = \left(E[X(s) X^*(s – \tau)]\right)^* = R_x^*(\tau) $$
性質2: 原点で最大値
$$ |R_x(\tau)| \leq R_x(0) \quad \text{（全ての } \tau \text{ に対して）} $$
これはコーシー・シュワルツの不等式から直ちに従います。信号は時間ずれ $\tau = 0$ のとき自分自身と最も相関が強く、$\tau$ が大きくなるにつれて相関が弱まるのが一般的です。
性質3: 原点の値は平均パワー
$$ R_x(0) = E[|X(t)|^2] = P $$
これは平均パワーの定義そのものです。
自己相関関数の性質を踏まえると、PSD との関係を導く準備が整いました。次節では、PSD と自己相関関数を結びつけるウィーナー・ヒンチンの定理を導出します。
ウィーナー・ヒンチンの定理
定理の主張
ウィーナー・ヒンチンの定理は、WSS 過程 $X(t)$ のパワースペクトル密度 $S_X(f)$ と自己相関関数 $R_X(\tau)$ がフーリエ変換対であることを述べます。
$$ \begin{equation} S_X(f) = \int_{-\infty}^{\infty} R_X(\tau) \, e^{-j2\pi f \tau} \, d\tau = \mathcal{F}\{R_X(\tau)\} \end{equation} $$
$$ \begin{equation} R_X(\tau) = \int_{-\infty}^{\infty} S_X(f) \, e^{j2\pi f \tau} \, df = \mathcal{F}^{-1}\{S_X(f)\} \end{equation} $$
この定理の威力は絶大です。PSD を直接計算しようとすると「無限時間の信号のフーリエ変換」という扱いにくい極限操作が必要ですが、自己相関関数さえ分かれば、そのフーリエ変換をとるだけで PSD が得られるのです。
導出
ここから、ウィーナー・ヒンチンの定理を導出します。出発点はランダム信号のPSD の定義（式(2)）です。
WSS 過程 $X(t)$ を時間区間 $[-T, T]$ で窓掛けした信号のフーリエ変換は
$$ X_T(f) = \int_{-T}^{T} X(t) \, e^{-j2\pi ft} \, dt $$
PSD の定義から
$$ S_X(f) = \lim_{T \to \infty} \frac{E[|X_T(f)|^2]}{2T} $$
$|X_T(f)|^2 = X_T(f) \cdot X_T^*(f)$ を展開します。
$$ |X_T(f)|^2 = \left(\int_{-T}^{T} X(t) \, e^{-j2\pi ft} \, dt\right) \left(\int_{-T}^{T} X^*(s) \, e^{j2\pi fs} \, ds\right) $$
2つの積分を1つにまとめると
$$ |X_T(f)|^2 = \int_{-T}^{T} \int_{-T}^{T} X(t) X^*(s) \, e^{-j2\pi f(t-s)} \, dt \, ds $$
この二重積分に期待値をとります。$X(t)$ は確率変数ですが、積分と期待値の交換が可能であるとすると（可積分条件を仮定）
$$ E[|X_T(f)|^2] = \int_{-T}^{T} \int_{-T}^{T} E[X(t) X^*(s)] \, e^{-j2\pi f(t-s)} \, dt \, ds $$
WSS の仮定により、$E[X(t) X^*(s)] = R_X(t – s)$ です。ここで変数変換 $\tau = t – s$ を導入します。$s$ を固定して $t$ を $\tau$ に置き換えると $t = s + \tau$, $dt = d\tau$ です。
$$ E[|X_T(f)|^2] = \int_{-T}^{T} \left(\int_{-T-s}^{T-s} R_X(\tau) \, e^{-j2\pi f\tau} \, d\tau \right) ds $$
$\tau$ についての積分範囲は $s$ に依存しますが、$\tau$ の積分を先に行う形に順序を交換します。$\tau$ の取りうる範囲は $-2T \leq \tau \leq 2T$ です。$\tau$ を固定したとき、$s$ の範囲は $\max(-T, -T+\tau) \leq s \leq \min(T, T+\tau)$ を満たす必要があります。
$|\tau| \leq 2T$ のとき、$s$ の積分範囲の長さは $2T – |\tau|$ です。したがって
$$ E[|X_T(f)|^2] = \int_{-2T}^{2T} (2T – |\tau|) \, R_X(\tau) \, e^{-j2\pi f\tau} \, d\tau $$
$2T$ で割ると
$$ \frac{E[|X_T(f)|^2]}{2T} = \int_{-2T}^{2T} \left(1 – \frac{|\tau|}{2T}\right) R_X(\tau) \, e^{-j2\pi f\tau} \, d\tau $$
$T \to \infty$ の極限をとります。被積分関数中の $\left(1 – \frac{|\tau|}{2T}\right)$ は、$\tau$ が有限のとき $1$ に収束します。$R_X(\tau)$ が絶対可積分（$\int_{-\infty}^{\infty} |R_X(\tau)| d\tau < \infty$）であれば、ルベーグの優収束定理を適用できます。
$$ S_X(f) = \lim_{T \to \infty} \frac{E[|X_T(f)|^2]}{2T} = \int_{-\infty}^{\infty} R_X(\tau) \, e^{-j2\pi f\tau} \, d\tau $$
これがウィーナー・ヒンチンの定理です。PSD は自己相関関数のフーリエ変換に等しいことが示されました。逆変換の存在はフーリエ変換の反転公式から保証されます。
導出の要点の整理
導出のポイントを振り返りましょう。
- 出発点: PSD の定義 — 窓掛け信号のフーリエ変換の二乗の期待値を時間幅で正規化して極限をとる
- WSS の活用: $E[X(t)X^*(s)] = R_X(t-s)$ と書けることで、二重積分を $\tau = t – s$ の1変数に帰着させた
- 三角窓の出現: $s$ 積分を実行すると、重み $(1 – |\tau|/2T)$ が自然に現れる。これはバートレット窓（三角窓）と呼ばれ、有限時間窓による切り出しの効果を反映している
- 極限操作: $T \to \infty$ で三角窓が消え、純粋な自己相関関数のフーリエ変換が残る
PSD の非負性
ウィーナー・ヒンチンの定理から直ちに導かれる重要な性質として、PSD は常に非負です。
$$ S_X(f) \geq 0 \quad \text{（全ての } f \text{ に対して）} $$
これは定義から明らかです。$|X_T(f)|^2 \geq 0$ であり、期待値も非負なので、その極限も非負です。PSD が「パワーの密度」であることを考えれば、負の値をとらないのは物理的にも自然な要請です。
PSD と自己相関の対応関係の直感
ウィーナー・ヒンチンの定理が意味する対応関係を直感的に整理しておきましょう。
自己相関関数 $R_X(\tau)$ が広く裾を引く（ゆっくり減衰する）とき、信号は「過去の値を長く記憶している」ことを意味します。フーリエ変換の性質として、時間領域で幅が広い関数は周波数領域では幅が狭くなります。したがって、長い相関をもつ信号のPSD は狭い周波数帯域に集中します。低周波成分が支配的な信号を思い浮かべれば納得できるでしょう。
逆に、自己相関が急激に減衰する信号は「過去の値をすぐに忘れる」ランダム性の高い信号であり、PSD は広い帯域にわたって分布します。この極端な場合が次に議論する白色雑音です。
白色雑音のPSD
白色雑音とは
白色雑音（white noise）とは、あらゆる周波数に均一にパワーが分布している理想的なランダム信号です。名前の由来は光の白色光にあります。白色光がすべての可視光の周波数を均等に含むように、白色雑音はすべての周波数成分を均等に含みます。
白色雑音 $W(t)$ の自己相関関数はデルタ関数で表されます。
$$ R_W(\tau) = \frac{N_0}{2} \delta(\tau) $$
ここで $N_0/2$ はパワースペクトル密度の値（定数）であり、$\delta(\tau)$ はディラックのデルタ関数です。
この自己相関関数が意味するのは、白色雑音は「時間ずれがゼロのときだけ相関をもち、それ以外の瞬間とは全く無相関」ということです。過去の値を一切記憶しない、究極的にランダムな信号です。
白色雑音のPSDの導出
ウィーナー・ヒンチンの定理を適用すると
$$ S_W(f) = \int_{-\infty}^{\infty} R_W(\tau) \, e^{-j2\pi f\tau} \, d\tau = \int_{-\infty}^{\infty} \frac{N_0}{2} \delta(\tau) \, e^{-j2\pi f\tau} \, d\tau $$
デルタ関数の性質 $\int_{-\infty}^{\infty} \delta(\tau) g(\tau) d\tau = g(0)$ を使うと
$$ S_W(f) = \frac{N_0}{2} \cdot e^{0} = \frac{N_0}{2} $$
PSD が周波数 $f$ に依存しない定数 $N_0/2$ となりました。これこそが「全周波数に均一にパワーが分布する」ことの数学的表現です。
白色雑音の平均パワー
白色雑音の平均パワーは
$$ P = \int_{-\infty}^{\infty} S_W(f) \, df = \int_{-\infty}^{\infty} \frac{N_0}{2} \, df = \infty $$
パワーが無限大になります。これは白色雑音が数学的な理想化であり、現実には厳密に存在しないことを示しています。実際の雑音は、ある周波数帯域を超えると減衰するため「帯域制限白色雑音」としてモデル化されます。それでも、白色雑音は解析上非常に扱いやすく、かつ実用的にも良い近似となるため、通信工学や制御工学の基本モデルとして広く使われています。
白色雑音の議論により、ウィーナー・ヒンチンの定理の威力を実感できたのではないでしょうか。自己相関がデルタ関数ならばPSD は定数、自己相関が指数減衰ならばPSD はローレンツ型 — このように自己相関の形とPSD の形がフーリエ変換を通じて1対1に対応するのです。
次に、実際の有限長データからPSD をどのように推定するかという実践的な問題に進みます。
ピリオドグラムによるPSD推定
理論から推定へ
ウィーナー・ヒンチンの定理は、自己相関関数が既知であれば PSD を厳密に計算できることを保証します。しかし現実には、有限個の離散的なサンプルデータしか手に入りません。無限時間の信号を観測することは不可能なので、有限長データから PSD を「推定」する必要があります。
PSD 推定法には大きく分けてノンパラメトリック法（モデルを仮定しない）とパラメトリック法（信号モデルを仮定する）がありますが、ここではノンパラメトリック法の代表であるピリオドグラムとウェルチ法を取り上げます。
ピリオドグラムの定義
ピリオドグラム（periodogram）は、PSD の最も直接的な推定量です。$N$ 点のデータ $x[0], x[1], \dots, x[N-1]$ に対して、離散フーリエ変換（DFT）を用いて次のように定義されます。
$$ \begin{equation} \hat{S}_{\text{per}}(f) = \frac{1}{N f_s} \left| \sum_{n=0}^{N-1} x[n] \, e^{-j2\pi fn/f_s} \right|^2 = \frac{|\text{DFT}\{x[n]\}|^2}{N f_s} $$
ここで $f_s$ はサンプリング周波数です。この定義は、確定信号のPSD の定義（式(1)）の離散版と見なすことができます。
ピリオドグラムの統計的性質
ピリオドグラムには残念ながら大きな弱点があります。
漸近不偏性: $N \to \infty$ で $E[\hat{S}_{\text{per}}(f)] \to S_X(f)$ が成り立ちます。つまり、データ数を無限にすれば期待値は真のPSD に収束します。
非一致性: しかし、ピリオドグラムの分散は $N$ を増やしても減少しません。具体的には、白色ガウス雑音の場合
$$ \text{Var}[\hat{S}_{\text{per}}(f)] \approx S_X^2(f) $$
分散が推定量の二乗のオーダーであり、$N$ に依存しないのです。これはデータ数をいくら増やしてもピリオドグラムのばらつきが小さくならないことを意味します。統計推定量として「一致推定量ではない」（inconsistent）というのは致命的な欠点です。
直感的には次のように理解できます。データ長 $N$ を増やすと、DFT の周波数分解能は $f_s / N$ まで細かくなります。分解能が上がるほど各周波数ビンに入るデータ点が減るため、個々のビンの推定精度は改善されません。分解能の向上と推定精度の間にトレードオフが存在するのです。
ピリオドグラムのこのばらつきの大きさは、グラフを見ると一目瞭然です。PSD の推定値がギザギザに激しく振動し、真の PSD の滑らかな形状を把握しにくくなります。このギザギザを抑える方法が必要です。
ウェルチ法によるPSD推定
ウェルチ法の着想
ピリオドグラムの分散を減らすためのアイデアは非常にシンプルです。長いデータをいくつかの短い区間（セグメント）に分割し、各セグメントのピリオドグラムを計算して平均をとれば、分散は平均化によって小さくなるはずです。
P.D. ウェルチが1967年に提案した方法（ウェルチ法, Welch’s method）は、この着想に2つの工夫を加えています。
- 窓関数の適用: 各セグメントに窓関数（ハニング窓、ハミング窓など）を掛けて、スペクトル漏れ（spectral leakage）を抑制する
- セグメントのオーバーラップ: 隣接セグメントを部分的に重ねることで、データを有効活用する
ウェルチ法のアルゴリズム
ウェルチ法は以下の手順で実行されます。
ステップ1: セグメント分割
データ $x[0], \dots, x[N-1]$ をセグメント長 $L$、オーバーラップ量 $D$ で分割します。$k$ 番目のセグメントは
$$ x_k[n] = x[n + kD], \quad n = 0, 1, \dots, L-1 $$
セグメント数は $K = \lfloor (N – L) / D \rfloor + 1$ です。一般的にはオーバーラップ率 50%（$D = L/2$）がよく使われます。
ステップ2: 窓関数の適用
各セグメントに窓関数 $w[n]$ を掛けます。
$$ y_k[n] = w[n] \cdot x_k[n] $$
窓関数を適用する理由は、有限長の信号を切り出す操作自体が矩形窓を掛けることに相当し、周波数領域で sinc 関数による裾の広がり（スペクトル漏れ）を引き起こすからです。ハニング窓やハミング窓は裾を滑らかにすることで、この漏れを軽減します。
ステップ3: 修正ピリオドグラムの計算
各セグメントの修正ピリオドグラムを計算します。
$$ \hat{S}_k(f) = \frac{1}{L f_s U} \left| \sum_{n=0}^{L-1} y_k[n] \, e^{-j2\pi fn/f_s} \right|^2 $$
ここで $U$ は窓関数のパワーの正規化定数です。
$$ U = \frac{1}{L} \sum_{n=0}^{L-1} |w[n]|^2 $$
ステップ4: 平均化
$K$ 個の修正ピリオドグラムの算術平均をとります。
$$ \hat{S}_{\text{Welch}}(f) = \frac{1}{K} \sum_{k=0}^{K-1} \hat{S}_k(f) $$
ウェルチ法の統計的性質
$K$ 個のセグメントが独立であれば、分散は $1/K$ に減少します。50% オーバーラップの場合、隣接セグメントは完全には独立ではありませんが、窓関数の効果により実効的な自由度は増加し、分散は大幅に低減されます。
$$ \text{Var}[\hat{S}_{\text{Welch}}(f)] \approx \frac{S_X^2(f)}{K_{\text{eff}}} $$
ここで $K_{\text{eff}}$ はオーバーラップと窓関数に依存する実効セグメント数です。50% オーバーラップのハニング窓では $K_{\text{eff}} \approx 0.75 K$ 程度になります。
ただし、セグメント長 $L$ を短くするほど周波数分解能が $f_s / L$ に粗くなります。つまりウェルチ法にも分散の低減と周波数分解能の間のトレードオフが存在します。短いセグメントで多く平均すればギザギザは減りますが、近接した2つの周波数成分を区別する能力が下がるのです。実用上はデータ長やアプリケーションに応じてセグメント長とオーバーラップ率を調整します。
ここまでの理論を踏まえて、Python で実際にPSD を計算してみましょう。理論的な予測が実際のデータで確認できることを実感できるはずです。
Pythonでの実装
各種信号の生成とPSDの理論値
まず、PSD の理論値が既知のいくつかの信号を生成し、理論と推定の対応を確認します。対象とする信号は以下の3つです。
- 正弦波: 周波数 $f_0$ の正弦波のPSD は $f = \pm f_0$ にデルタ関数的なピークをもつ
- 白色雑音: PSD は全周波数で一定
- 色付き雑音（AR(1)過程）: 低周波側にパワーが偏る
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
# --- パラメータ設定 ---
fs = 1000 # サンプリング周波数 [Hz]
N = 8192 # データ点数
t = np.arange(N) / fs # 時間軸
# --- 信号生成 ---
np.random.seed(42)
# (1) 正弦波 (50 Hz, 振幅1) + 弱い白色雑音
f0 = 50
x_sine = np.sin(2 * np.pi * f0 * t) + 0.5 * np.random.randn(N)
# (2) 白色ガウス雑音
sigma2 = 1.0 # 分散
x_white = np.sqrt(sigma2) * np.random.randn(N)
# (3) AR(1) 色付き雑音: x[n] = a*x[n-1] + w[n]
a = 0.95
x_ar1 = np.zeros(N)
w = np.random.randn(N)
for n in range(1, N):
x_ar1[n] = a * x_ar1[n - 1] + w[n]
print(f"データ点数: {N}, サンプリング周波数: {fs} Hz")
print(f"周波数分解能: {fs/N:.4f} Hz")
print(f"正弦波信号の分散: {np.var(x_sine):.4f}")
print(f"白色雑音の分散: {np.var(x_white):.4f}")
print(f"AR(1)雑音の分散: {np.var(x_ar1):.4f}")
上のコードでは3種類の信号を生成しています。正弦波信号には弱い白色雑音を加えて、現実の計測状況に近づけています。AR(1) 過程は $x[n] = a \cdot x[n-1] + w[n]$ という差分方程式で生成されるランダム過程で、係数 $a = 0.95$ により強い正の相関をもちます。$a$ が 1 に近いほど低周波成分が強まり、0 に近づくと白色雑音に近づきます。
AR(1) 過程のPSD の理論値は、$z$ 変換を用いて導出でき
$$ S_{\text{AR1}}(f) = \frac{\sigma_w^2 / f_s}{|1 – a e^{-j2\pi f/f_s}|^2} = \frac{\sigma_w^2 / f_s}{1 – 2a\cos(2\pi f/f_s) + a^2} $$
となります。ここで $\sigma_w^2$ は駆動白色雑音の分散です。
ピリオドグラムの計算
次に、3つの信号に対してピリオドグラムを計算します。
def periodogram(x, fs):
"""ピリオドグラムを計算する"""
N = len(x)
X = np.fft.rfft(x)
freqs = np.fft.rfftfreq(N, d=1/fs)
# 片側PSD（正の周波数のみ、両側分のパワーを反映）
Pxx = (np.abs(X)**2) / (N * fs)
Pxx[1:-1] *= 2 # DC成分とナイキスト以外は2倍
return freqs, Pxx
# 各信号のピリオドグラム
freq_s, psd_sine_per = periodogram(x_sine, fs)
freq_w, psd_white_per = periodogram(x_white, fs)
freq_a, psd_ar1_per = periodogram(x_ar1, fs)
# 可視化
fig, axes = plt.subplots(3, 1, figsize=(10, 10))
axes[0].semilogy(freq_s, psd_sine_per, alpha=0.7, linewidth=0.5)
axes[0].set_title("Periodogram: Sine (50 Hz) + Noise")
axes[0].set_ylabel("PSD [V²/Hz]")
axes[0].set_xlim([0, fs/2])
axes[1].semilogy(freq_w, psd_white_per, alpha=0.7, linewidth=0.5)
axes[1].axhline(y=sigma2/fs, color='r', linestyle='--', label=f"Theory: {sigma2/fs:.4f}")
axes[1].set_title("Periodogram: White Noise")
axes[1].set_ylabel("PSD [V²/Hz]")
axes[1].set_xlim([0, fs/2])
axes[1].legend()
# AR(1) 理論PSD
f_theory = np.linspace(0, fs/2, 1000)
psd_ar1_theory = (1.0 / fs) / (1 - 2*a*np.cos(2*np.pi*f_theory/fs) + a**2)
axes[2].semilogy(freq_a, psd_ar1_per, alpha=0.7, linewidth=0.5, label="Periodogram")
axes[2].semilogy(f_theory, psd_ar1_theory, 'r-', linewidth=2, label="Theory")
axes[2].set_title("Periodogram: AR(1) Process (a=0.95)")
axes[2].set_ylabel("PSD [V²/Hz]")
axes[2].set_xlabel("Frequency [Hz]")
axes[2].set_xlim([0, fs/2])
axes[2].legend()
plt.tight_layout()
plt.savefig("periodogram_three_signals.png", dpi=150, bbox_inches="tight")
plt.show()
ピリオドグラムの結果から、3つの重要な観察ができます。
- 正弦波信号: 50 Hz 付近に鋭いピークが現れ、それ以外の周波数帯域ではノイズフロアが広がっています。これは理論どおり、正弦波のPSD がデルタ関数的なピークをもつことと整合します。
- 白色雑音: 理論値 $\sigma^2 / f_s$（赤の破線）の周りに推定値が激しく振動しています。ピリオドグラムの「非一致性」がはっきり見て取れます。平均的には理論値に近いのですが、個々の周波数ビンの値は大きくばらつきます。
- AR(1) 過程: 低周波領域でパワーが大きく、高周波に向かって減衰する傾向は理論曲線（赤線）と一致しますが、やはりギザギザが目立ちます。
いずれの信号でもピリオドグラムのばらつきが大きく、真のPSD の形状を滑らかに把握するのは困難です。これがウェルチ法を使う動機となります。
ウェルチ法の実装と比較
SciPy の signal.welch を使ってウェルチ法を適用し、ピリオドグラムと比較します。さらに、セグメント長を変えたときの分散と分解能のトレードオフも確認します。
# --- ウェルチ法によるPSD推定 ---
# セグメント長を変えて比較
nperseg_list = [256, 1024, 4096]
fig, axes = plt.subplots(3, 1, figsize=(10, 12))
# --- 正弦波 + 雑音 ---
axes[0].semilogy(freq_s, psd_sine_per, alpha=0.3, linewidth=0.5, label="Periodogram")
for nperseg in nperseg_list:
f_welch, psd_welch = signal.welch(x_sine, fs=fs, nperseg=nperseg,
noverlap=nperseg//2, window='hann')
axes[0].semilogy(f_welch, psd_welch, linewidth=1.5,
label=f"Welch (L={nperseg})")
axes[0].set_title("Sine (50 Hz) + Noise")
axes[0].set_ylabel("PSD [V²/Hz]")
axes[0].set_xlim([0, 200])
axes[0].legend(fontsize=8)
# --- 白色雑音 ---
axes[1].semilogy(freq_w, psd_white_per, alpha=0.3, linewidth=0.5, label="Periodogram")
for nperseg in nperseg_list:
f_welch, psd_welch = signal.welch(x_white, fs=fs, nperseg=nperseg,
noverlap=nperseg//2, window='hann')
axes[1].semilogy(f_welch, psd_welch, linewidth=1.5,
label=f"Welch (L={nperseg})")
axes[1].axhline(y=sigma2/fs, color='k', linestyle='--', linewidth=1, label="Theory")
axes[1].set_title("White Noise")
axes[1].set_ylabel("PSD [V²/Hz]")
axes[1].set_xlim([0, fs/2])
axes[1].legend(fontsize=8)
# --- AR(1) 過程 ---
axes[2].semilogy(freq_a, psd_ar1_per, alpha=0.3, linewidth=0.5, label="Periodogram")
for nperseg in nperseg_list:
f_welch, psd_welch = signal.welch(x_ar1, fs=fs, nperseg=nperseg,
noverlap=nperseg//2, window='hann')
axes[2].semilogy(f_welch, psd_welch, linewidth=1.5,
label=f"Welch (L={nperseg})")
axes[2].semilogy(f_theory, psd_ar1_theory, 'k--', linewidth=1.5, label="Theory")
axes[2].set_title("AR(1) Process (a=0.95)")
axes[2].set_ylabel("PSD [V²/Hz]")
axes[2].set_xlabel("Frequency [Hz]")
axes[2].set_xlim([0, fs/2])
axes[2].legend(fontsize=8)
plt.tight_layout()
plt.savefig("welch_comparison.png", dpi=150, bbox_inches="tight")
plt.show()
ウェルチ法の比較結果から、以下のことが読み取れます。
- 分散の低減: セグメント長 $L = 256$ のウェルチ法（セグメント数 $K \approx 63$）は、ピリオドグラムと比べて圧倒的に滑らかです。白色雑音のPSD がほぼ水平な直線に見えるようになり、理論値との一致が明確です。
- 分解能とのトレードオフ: 正弦波信号を見ると、$L = 256$ ではピークの幅が広がっている一方、$L = 4096$ ではピークが鋭くなっています。これはセグメント長が短いほど周波数分解能（$f_s / L$）が粗くなるためです。$L = 256$ では分解能が約 3.9 Hz、$L = 4096$ では約 0.24 Hz です。
- AR(1) の理論曲線との一致: $L = 1024$ のウェルチ法は分散の低減と分解能のバランスが良く、理論曲線にかなり近い滑らかな推定結果が得られています。実用上は $L = 1024$ 程度が多くの場面で良い出発点となります。
自己相関関数からのPSD計算
ウィーナー・ヒンチンの定理を直接使い、自己相関関数のフーリエ変換としてPSD を計算する方法も確認しましょう。
def psd_from_autocorr(x, fs):
"""自己相関関数のFFTによるPSD推定"""
N = len(x)
# バイアス付き自己相関推定
Rxx = np.correlate(x, x, mode='full') / N
# 自己相関のフーリエ変換
S = np.fft.rfft(Rxx[N-1:]) # 正のラグのみ使用
freqs = np.fft.rfftfreq(N, d=1/fs)
# 実部をとる（理論上は実数）
Pxx = np.real(S) / fs
Pxx[1:-1] *= 2
return freqs, Pxx
# AR(1) 過程で検証
freq_ac, psd_ar1_ac = psd_from_autocorr(x_ar1, fs)
fig, ax = plt.subplots(figsize=(10, 5))
ax.semilogy(freq_a, psd_ar1_per, alpha=0.3, linewidth=0.5, label="Periodogram")
ax.semilogy(freq_ac, psd_ar1_ac, 'g-', alpha=0.7, linewidth=1.0,
label="Autocorrelation method")
f_welch, psd_welch_1024 = signal.welch(x_ar1, fs=fs, nperseg=1024,
noverlap=512, window='hann')
ax.semilogy(f_welch, psd_welch_1024, 'b-', linewidth=1.5, label="Welch (L=1024)")
ax.semilogy(f_theory, psd_ar1_theory, 'r--', linewidth=2, label="Theory")
ax.set_xlabel("Frequency [Hz]")
ax.set_ylabel("PSD [V²/Hz]")
ax.set_title("PSD Estimation Comparison for AR(1) Process")
ax.set_xlim([0, fs/2])
ax.legend()
plt.tight_layout()
plt.savefig("psd_autocorr_comparison.png", dpi=150, bbox_inches="tight")
plt.show()
このグラフでは、4つの手法を重ね合わせています。自己相関法（緑線）は、ウィーナー・ヒンチンの定理を直接適用した結果であり、全データの自己相関を使っているためピリオドグラムと本質的に同じスペクトル形状を示します。一方、ウェルチ法（青線）は平均化の効果により滑らかな推定結果が得られ、理論曲線（赤破線）に最も近い推定となっています。自己相関法とピリオドグラムの形状がほぼ一致することは、ウィーナー・ヒンチンの定理が「PSD は自己相関のフーリエ変換」と述べていることの数値的な裏付けです。
窓関数の影響を可視化する
ウェルチ法で使う窓関数の選択がPSD推定にどのように影響するかを確認します。
# --- 窓関数によるスペクトル漏れの違い ---
# 近接した2つの正弦波（50 Hz と 55 Hz）
x_two_sine = (np.sin(2 * np.pi * 50 * t)
+ 0.5 * np.sin(2 * np.pi * 55 * t)
+ 0.3 * np.random.randn(N))
windows = ['boxcar', 'hann', 'blackman']
window_labels = ['Rectangular', 'Hann', 'Blackman']
fig, ax = plt.subplots(figsize=(10, 5))
for win, label in zip(windows, window_labels):
f_w, psd_w = signal.welch(x_two_sine, fs=fs, nperseg=1024,
noverlap=512, window=win)
ax.semilogy(f_w, psd_w, linewidth=1.5, label=label)
ax.set_xlabel("Frequency [Hz]")
ax.set_ylabel("PSD [V²/Hz]")
ax.set_title("Window Function Effect on Two Close Sinusoids (50 Hz + 55 Hz)")
ax.set_xlim([30, 80])
ax.legend()
plt.tight_layout()
plt.savefig("window_comparison.png", dpi=150, bbox_inches="tight")
plt.show()
窓関数の比較からは、2つの重要な知見が得られます。
- 矩形窓（Rectangular）: メインローブが最も狭く周波数分解能は高いですが、サイドローブが大きいため、50 Hz と 55 Hz の間の谷が浅くなり、弱い成分（55 Hz）が強い成分（50 Hz）のサイドローブに埋もれやすくなります。
- ブラックマン窓（Blackman）: サイドローブが最も小さく、近接周波数の分離が良好ですが、メインローブが広がるため、ピークの幅がやや広くなります。ハニング窓はこの両者の中間的な特性をもち、多くの実用場面で良い選択となります。
窓関数の選択は「分解能」と「漏れ抑制」のトレードオフであり、アプリケーションの要求に応じて使い分ける必要があります。
実践的なPSD解析の例
最後に、複数の周波数成分をもつ実用的な信号に対して、ウェルチ法でPSD 解析を行う一連の流れを示します。
# --- 複合信号のPSD解析 ---
# 3つの正弦波 + 帯域制限雑音
np.random.seed(123)
N_long = 32768
t_long = np.arange(N_long) / fs
# 信号成分
s1 = 1.0 * np.sin(2 * np.pi * 30 * t_long) # 30 Hz, 振幅 1.0
s2 = 0.5 * np.sin(2 * np.pi * 100 * t_long) # 100 Hz, 振幅 0.5
s3 = 0.2 * np.sin(2 * np.pi * 250 * t_long) # 250 Hz, 振幅 0.2
# 帯域制限雑音（50-150 Hz のバンドパスフィルタ）
noise = np.random.randn(N_long)
b_bp, a_bp = signal.butter(4, [50, 150], btype='band', fs=fs)
bandlimited_noise = 0.3 * signal.filtfilt(b_bp, a_bp, noise)
x_complex = s1 + s2 + s3 + bandlimited_noise
# ウェルチ法でPSD推定
f_psd, psd_complex = signal.welch(x_complex, fs=fs, nperseg=2048,
noverlap=1024, window='hann')
fig, axes = plt.subplots(2, 1, figsize=(10, 8))
# 時間波形
axes[0].plot(t_long[:2000], x_complex[:2000], linewidth=0.5)
axes[0].set_xlabel("Time [s]")
axes[0].set_ylabel("Amplitude")
axes[0].set_title("Composite Signal: Time Domain")
axes[0].set_xlim([0, 2])
# PSD
axes[1].semilogy(f_psd, psd_complex, linewidth=1.5)
axes[1].set_xlabel("Frequency [Hz]")
axes[1].set_ylabel("PSD [V²/Hz]")
axes[1].set_title("Composite Signal: PSD (Welch, L=2048)")
axes[1].set_xlim([0, fs/2])
# ピーク位置にアノテーション
for freq_peak in [30, 100, 250]:
idx = np.argmin(np.abs(f_psd - freq_peak))
axes[1].annotate(f"{freq_peak} Hz",
xy=(f_psd[idx], psd_complex[idx]),
xytext=(f_psd[idx]+20, psd_complex[idx]*5),
arrowprops=dict(arrowstyle='->', color='red'),
fontsize=10, color='red')
plt.tight_layout()
plt.savefig("composite_signal_psd.png", dpi=150, bbox_inches="tight")
plt.show()
この解析結果から、ウェルチ法によるPSD推定の実用的な能力がよくわかります。
- 3つの正弦波成分の検出: 30 Hz、100 Hz、250 Hz に明確なピークが現れています。ピークの高さは振幅の二乗に比例しており、30 Hz（振幅 1.0）のピークが最も高く、250 Hz（振幅 0.2）のピークが最も低くなっています。これは PSD がパワー（振幅の二乗）の密度であることと一致します。
- 帯域制限雑音の検出: 50 Hz から 150 Hz の帯域で PSD のフロアが持ち上がっており、バンドパスフィルタを通した雑音成分が検出されています。100 Hz のピークはこの雑音帯域内にありますが、正弦波成分が十分に強いため、雑音フロアの上に明確に突出しています。
- 時間波形との対比: 時間領域の波形（上段）を見ただけでは、3つの周波数成分を視覚的に分離するのは困難です。PSD 解析（下段）を通じて初めて、信号に含まれる周波数構造が明確に可視化されます。これがスペクトル解析の価値です。
PSD の単位と片側・両側スペクトル
単位について
PSD の単位を整理しておきます。信号 $x(t)$ の単位が $[\text{V}]$（電圧）であるとき
- 自己相関関数 $R_X(\tau)$ の単位: $[\text{V}^2]$
- PSD $S_X(f)$ の単位: $[\text{V}^2 / \text{Hz}]$
PSD を帯域幅 $\Delta f$ で積分するとパワー $[\text{V}^2]$ が得られることからも、この単位は整合しています。dB 表記ではよく $10 \log_{10}(S_X(f))$ が使われ、単位は $[\text{dB/Hz}]$ となります。
片側スペクトルと両側スペクトル
理論的な PSD の定義では、$-\infty < f < \infty$ の両側スペクトルを考えています。実信号の場合、$S_X(-f) = S_X(f)$ が成り立つため、正の周波数のみで表示する片側スペクトルがよく使われます。
$$ G_X(f) = 2 S_X(f), \quad f \geq 0 $$
片側スペクトルを使えば $P = \int_0^{\infty} G_X(f) df$ でパワーが得られます。SciPy の signal.welch はデフォルトで片側スペクトルを返します（return_onesided=True）。文献によって両側・片側のどちらを使っているかが異なるため、PSD の値を比較する際は注意が必要です。
この区別は些細に見えますが、通信工学では $N_0/2$（両側）と $N_0$（片側）の混同がしばしば設計ミスにつながるため、常に意識しておくべきです。
具体例 — 指数減衰する自己相関のPSD
問題設定
ウィーナー・ヒンチンの定理の計算例として、自己相関関数が指数減衰する場合のPSD を解析的に求めてみましょう。
$$ R_X(\tau) = \sigma^2 e^{-\alpha |\tau|} $$
ここで $\sigma^2 = R_X(0)$ は平均パワー、$\alpha > 0$ は減衰率です。$\alpha$ が大きいほど自己相関の減衰が速く、信号の「記憶」が短いことを意味します。この形の自己相関は、RC回路を通した白色雑音や、オルンシュタイン・ウーレンベック過程など多くの物理系で現れます。
PSD の導出
ウィーナー・ヒンチンの定理を適用します。
$$ S_X(f) = \int_{-\infty}^{\infty} \sigma^2 e^{-\alpha|\tau|} e^{-j2\pi f\tau} \, d\tau $$
絶対値を場合分けして、$\tau < 0$ と $\tau \geq 0$ に分けます。
$$ S_X(f) = \sigma^2 \left[ \int_{-\infty}^{0} e^{\alpha\tau} e^{-j2\pi f\tau} d\tau + \int_{0}^{\infty} e^{-\alpha\tau} e^{-j2\pi f\tau} d\tau \right] $$
第1項の積分では、指数部が $(\alpha – j2\pi f)\tau$ となります。$\tau \to -\infty$ で $e^{(\alpha – j2\pi f)\tau} \to 0$（$\alpha > 0$ なので収束）より
$$ \int_{-\infty}^{0} e^{(\alpha – j2\pi f)\tau} d\tau = \left[\frac{e^{(\alpha – j2\pi f)\tau}}{\alpha – j2\pi f}\right]_{-\infty}^{0} = \frac{1}{\alpha – j2\pi f} $$
第2項も同様に計算します。
$$ \int_{0}^{\infty} e^{-(\alpha + j2\pi f)\tau} d\tau = \frac{1}{\alpha + j2\pi f} $$
2つの項を足し合わせ、通分すると
$$ S_X(f) = \sigma^2 \left(\frac{1}{\alpha – j2\pi f} + \frac{1}{\alpha + j2\pi f}\right) = \sigma^2 \cdot \frac{2\alpha}{\alpha^2 + (2\pi f)^2} $$
得られたPSD はローレンツ型（ローレンツ分布の形）です。
$$ \begin{equation} S_X(f) = \frac{2\alpha\sigma^2}{\alpha^2 + (2\pi f)^2} \end{equation} $$
結果の検証
この結果を検証しましょう。
パワーの確認: PSD を全周波数にわたって積分すると
$$ \int_{-\infty}^{\infty} S_X(f) df = \sigma^2 \cdot 2\alpha \int_{-\infty}^{\infty} \frac{df}{\alpha^2 + (2\pi f)^2} $$
$u = 2\pi f/\alpha$ と変数変換すると
$$ = \sigma^2 \cdot 2\alpha \cdot \frac{1}{2\pi\alpha} \int_{-\infty}^{\infty} \frac{du}{1 + u^2} = \sigma^2 \cdot \frac{1}{\pi} \cdot \pi = \sigma^2 $$
確かに $R_X(0) = \sigma^2$ に一致しており、パーシバルの関係が成り立っています。
極限の確認: $\alpha \to \infty$（自己相関が瞬時に減衰）の極限では、$R_X(\tau) \to \sigma^2 \delta(\tau)$ に近づき、PSD は全帯域にわたってフラットになるはずです。実際、$\alpha \to \infty$ で $S_X(f) \to 2\sigma^2/\alpha \cdot (\alpha^2/(\alpha^2 + (2\pi f)^2))$ の広がりが無限大になる様子が確認でき、白色雑音の極限に整合します。
逆に $\alpha \to 0$ の極限では相関が無限に長くなり、PSD は $f = 0$ に集中する鋭いピークとなります。これは信号がほぼ定数（DC成分のみ）になることに対応します。
このように、自己相関の時間スケール $1/\alpha$ と PSD の帯域幅が反比例する関係は、フーリエ変換の「時間幅と帯域幅の不確定性関係」の一つの表れです。
次のPython コードで、この解析解と数値推定の対応を可視化しましょう。
# --- 指数減衰自己相関のPSD: 解析解と数値推定の比較 ---
alpha_values = [5, 20, 100]
sigma2_exp = 1.0
N_ou = 32768
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
for alpha in alpha_values:
# オルンシュタイン・ウーレンベック過程の離散近似
dt = 1 / fs
x_ou = np.zeros(N_ou)
for n in range(1, N_ou):
x_ou[n] = np.exp(-alpha * dt) * x_ou[n-1] + np.sqrt(
sigma2_exp * (1 - np.exp(-2 * alpha * dt))) * np.random.randn()
# 数値的な自己相関
max_lag = 500
lags = np.arange(max_lag)
Rxx_num = np.array([np.mean(x_ou[:N_ou-k] * x_ou[k:]) for k in lags])
# 理論自己相関
tau_theory = lags * dt
Rxx_theory = sigma2_exp * np.exp(-alpha * tau_theory)
axes[0].plot(tau_theory, Rxx_num, alpha=0.7, linewidth=1,
label=f"Numerical (α={alpha})")
axes[0].plot(tau_theory, Rxx_theory, '--', linewidth=1.5,
label=f"Theory (α={alpha})")
# PSD: ウェルチ法
f_w, psd_w = signal.welch(x_ou, fs=fs, nperseg=2048,
noverlap=1024, window='hann')
# PSD: 理論
psd_theory = 2 * alpha * sigma2_exp / (alpha**2 + (2*np.pi*f_w)**2)
axes[1].semilogy(f_w, psd_w, alpha=0.7, linewidth=1,
label=f"Welch (α={alpha})")
axes[1].semilogy(f_w, psd_theory, '--', linewidth=1.5,
label=f"Theory (α={alpha})")
axes[0].set_xlabel("Lag τ [s]")
axes[0].set_ylabel("R_X(τ)")
axes[0].set_title("Autocorrelation Function")
axes[0].set_xlim([0, 0.3])
axes[0].legend(fontsize=7)
axes[1].set_xlabel("Frequency [Hz]")
axes[1].set_ylabel("PSD [V²/Hz]")
axes[1].set_title("Power Spectral Density (Lorentzian)")
axes[1].set_xlim([0, fs/2])
axes[1].legend(fontsize=7)
plt.tight_layout()
plt.savefig("exponential_autocorr_psd.png", dpi=150, bbox_inches="tight")
plt.show()
このグラフから、ウィーナー・ヒンチンの定理が実データでも成立していることが確認できます。
- 左図（自己相関）: 数値計算による自己相関推定値（実線）が理論値（破線）とよく一致しています。$\alpha$ が大きいほど自己相関の減衰が速いことが視覚的に明確です。
- 右図（PSD）: ウェルチ法による推定（実線）がローレンツ型の理論 PSD（破線）とよく一致しています。$\alpha = 5$ のとき PSD は低周波に集中し急峻に減衰しますが、$\alpha = 100$ では広い帯域にわたってほぼフラットになり、白色雑音に近づいていることがわかります。
- 自己相関とPSD の対応: $\alpha$ が小さいほど自己相関が広く（左図）、PSD が狭い帯域に集中する（右図）。$\alpha$ が大きいほど自己相関が狭く、PSD が広帯域に広がる。時間幅と帯域幅の反比例関係がはっきり見えています。
応用上の注意点
ゼロパディングの効果
FFT の計算でゼロパディング（データの末尾にゼロを追加してDFT の点数を増やす）を行うと、PSD のプロットが滑らかになります。しかし、これは周波数分解能が向上したわけではなく、既存の情報を補間して表示しているだけです。ゼロパディングは「見た目の滑らかさ」を改善しますが、新たな周波数情報を生み出すわけではない点に注意が必要です。
非定常信号への対応
ウィーナー・ヒンチンの定理は WSS 過程を前提としています。非定常信号（時間とともに統計的性質が変化する信号）に対しては、短時間フーリエ変換（STFT）によるスペクトログラムやウェーブレット変換など、時間−周波数解析の手法が必要になります。ウェルチ法は短いセグメントごとに PSD を計算するため、ある程度の非定常性には頑健ですが、信号の統計的性質が急激に変化する場合には限界があります。
サンプリング定理との関係
離散信号のPSD はナイキスト周波数 $f_s / 2$ までしか定義されません。$f_s / 2$ を超える周波数成分が存在する場合、エイリアシングによって PSD が歪みます。PSD 推定を正しく行うためには、対象信号がサンプリング定理を満たしていることが前提条件です。
まとめ
本記事では、パワースペクトル密度（PSD）の理論をエネルギースペクトル密度から出発して体系的に解説し、ウィーナー・ヒンチンの定理を厳密に導出しました。
- PSD の定義: エネルギー信号のフーリエ変換の二乗（ESD）を出発点とし、パワー信号に対しては窓掛け信号のフーリエ変換の二乗を時間幅で正規化して極限をとることで PSD を定義しました
- ウィーナー・ヒンチンの定理: WSS 過程の PSD が自己相関関数のフーリエ変換に等しいことを、PSD の定義から変数変換と極限操作を通じて導出しました
- 白色雑音: 自己相関がデルタ関数であるとき PSD が定数になることを示し、白色雑音が理想化であることを確認しました
- PSD 推定: ピリオドグラムが非一致推定量であることを確認し、ウェルチ法による分散低減を理論と実装の両面で検証しました
- 解析例: 指数減衰自己相関に対するローレンツ型 PSD を解析的に導出し、数値推定との一致を確認しました
PSD は信号処理のあらゆる分野の基礎となる概念です。今後は、PSD の知識を活かして以下のトピックに進むことができます。
次のステップとして、以下の記事も参考にしてください。
- 熱雑音の理論 — 白色雑音モデルの物理的起源と、ナイキストの定理による熱雑音 PSD の導出
- ディジタルフィルタの基礎 — PSD を用いたフィルタの周波数応答設計と、信号の帯域制限