オシロスコープでアンプの出力雑音を眺めていると、画面はいつまでもランダムに揺れ続けます。この波形を止めることはできませんし、「10秒後の値」を予言することもできません。それでも私たちは平気で「この雑音は1 kHz付近にパワーが集中している」「このアンプの雑音は白い」と口にします。予測できない信号なのに、なぜ周波数ごとのパワー分布だけは語れるのでしょうか。
ここには、はっきりした技術的な問題が隠れています。フーリエ変換 $\int_{-\infty}^{\infty} x(t)e^{-j2\pi ft}dt$ が存在するには、少なくとも $x(t)$ が無限遠でおとなしく減衰していなければなりません。ところが定常な雑音は、いつまで待っても同じ強さで揺れ続けます。絶対可積分ではないし、二乗可積分でもない。つまりランダムな定常信号にフーリエ変換は直接使えないのです。にもかかわらずスペクトルアナライザは、雑音のスペクトルを何事もなく表示します。
この矛盾を解消するのが、本記事の主役であるウィーナー・ヒンチンの定理（Wiener–Khinchin theorem）です。定理の主張は驚くほど短くまとまります。「定常確率過程のパワースペクトル密度 $S(f)$ は、自己相関関数 $R(\tau)$ のフーリエ変換に等しい」。信号そのものはフーリエ変換できないけれど、信号から作った統計量である $R(\tau)$ ならフーリエ変換できる。そして $R(\tau)$ を変換したものが、まさに私たちがスペクトルと呼んでいたものに一致する — これが定理の中身です。
この定理は理屈の上での美しさだけの話ではありません。応用先を2つ挙げるだけでも、その射程の広さが見えてきます。ひとつは受信機の雑音設計です。アンテナで受けた微弱な信号に熱雑音がどう乗るかを見積もるとき、雑音の自己相関がほぼデルタ関数であること（＝白色）から、スペクトルが平坦であることを結論します。帯域幅 $B$ の受信機に入る雑音電力が $N_0B$ になるという基本式は、この定理の直接の帰結です。もうひとつはFFTによるスペクトル推定そのものです。測定器やライブラリが計算しているピリオドグラムは「データのFFTの二乗を長さで割った量」ですが、その期待値がなぜ真のスペクトルに近づくのか、なぜデータを増やしてもギザギザが消えないのかは、ウィーナー・ヒンチンの定理を有限長で書き直して初めて説明できます。振動計測の共振同定、脳波のα波・β波の分離、レーダーのクラッタ解析なども、すべて同じ土台の上に乗っています。
定理の全体像を1枚の絵にしておきます。数式に入る前に、「何を諦めて何に乗り換えたのか」という道筋だけ頭に入れてください。
赤い上のルートが「サンプルパスをそのままフーリエ変換する」道で、エネルギーが無限大になるため通行止めです。緑の下のルートが定理の道筋で、いったん自己相関 $R(\tau)$ という統計量に落としてからフーリエ変換します。そして黄色の下段は、スペクトルアナライザが実際にやっている測定手順（有限区間で切って二乗して $T$ で割って平均する）で、この測定量と緑ルートの終点が一致するというのが定理の主張です。
本記事の内容
- ランダム信号にフーリエ変換が直接使えない理由（絶対可積分でない）
- 定常確率過程と自己相関関数 $R(\tau)$ の整理
- 有限区間フーリエ変換のパワーの期待値からのPSDの定義
- ウィーナー・ヒンチンの定理の導出（連続時間・省略なし）
- 離散時間版の導出とピリオドグラムの期待値（バートレット窓・フェイエール核）
- $S(f)\ge 0$、$R(\tau)$ の正定値性、$R(0)=\int S(f)df$ の対応
- ボホナー・ヘルグロッツの定理との関係
- 白色雑音・AR(1)過程・正弦波＋雑音の3例での解析解
- Pythonで「標本自己相関のFFT」「ピリオドグラムの多数回平均」「理論スペクトル」の3本が一致することの検証
- ピリオドグラムが一致推定量にならない理由
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- 自己相関関数とは — $R(\tau)$ の定義、偶関数性、原点最大
- パワースペクトル密度（PSD）の理論 — PSDの位置づけとエネルギースペクトル密度との違い
- 確率過程入門 — 確率過程、定常性、アンサンブル平均
- フーリエ変換 — 変換対、パーシバルの等式
なぜランダム信号にフーリエ変換が使えないのか
まず、問題がどれくらい深刻なのかをきちんと見ておきましょう。ここを曖昧にしたまま定理だけ覚えると、「自己相関をフーリエ変換すればスペクトルになる」という呪文だけが残って、何が証明されたのか分からなくなります。
フーリエ変換 $\hat x(f)=\int_{-\infty}^{\infty}x(t)e^{-j2\pi ft}dt$ が普通の意味で存在するための十分条件は、$x$ が絶対可積分であること、つまり
$$ \int_{-\infty}^{\infty}|x(t)|\,dt < \infty $$
が成り立つことです。もう少し広げてプランシュレルの意味で考えるなら、$x$ が二乗可積分であること、すなわち $\int|x(t)|^2dt<\infty$（有限エネルギー）でも構いません。
ここで、定常なランダム信号を考えます。「定常」とは大雑把に言えば「統計的な性質が時間によらない」ということです。熱雑音は昨日も今日も同じ強さで揺れていますし、100年後も同じ強さで揺れているはずです。この過程の1本の実現（サンプルパス）$x(t)$ を取ってきて、そのエネルギーを計算してみましょう。二乗の期待値が定数 $E[X(t)^2]=P>0$ なので、
$$ E\left[\int_{-T/2}^{T/2}X(t)^2dt\right]=\int_{-T/2}^{T/2}E[X(t)^2]dt = P\,T $$
となり、$T\to\infty$ で発散します。エネルギーは無限大です。絶対可積分でもありません。したがって、定常確率過程のサンプルパスは、通常の意味でフーリエ変換できないのです。
この事実は、単なる数学的な神経質さではありません。実際にやってみれば分かります。長さ $N$ の白色雑音のFFTを取ると、確かに何らかの値は出てきます。しかし $N$ を2倍、4倍と増やしていっても、その値は特定の数に収束しません。バラつきの大きさが一向に減らないのです（この現象は本記事の後半で数値実験として確認します）。「収束しないもの」を信号の性質として語ることはできません。
実際に数値で見ておきましょう。定常なランダム信号を1本生成し、そのエネルギーとパワーを観測時間 $T$ の関数としてプロットしたものが次の図です。
中央のグラフで、累積エネルギー $\int_0^T x^2dt$ が観測時間に比例してまっすぐ伸び続けています。これが $PT$ が発散するという計算そのもので、フーリエ変換の存在条件が破れる理由です。一方、右のグラフのように $T$ で割ってパワーにすると、最初こそ暴れるものの $T=30$ 秒あたりでは平均パワーの水平線にぴたりと張りついています。発散を止める鍵は「時間で割る」ことだと、この2枚を見比べれば分かります。
では、どうすればよいでしょうか。突破口は2つの方向にあります。
方向1: エネルギーではなくパワーを見る。 エネルギーが無限に発散するのは、無限時間ぶんを足しているからです。ならば時間で割って、単位時間あたりの量、つまりパワーで測ればよい。$\frac{1}{T}\int_{-T/2}^{T/2}x(t)^2dt$ は $T\to\infty$ でも有限の値に落ち着くはずです。
方向2: 実現ごとではなくアンサンブル平均で見る。 1本のサンプルパスは永遠にランダムに揺れます。しかし多数の実現にわたって平均を取れば、ランダムさは平均化されて消えます。期待値 $E[\cdot]$ を取るのです。
この2つを両方使うのが、これから見る道筋です。そして重要なのは、フーリエ変換する対象を信号そのものから統計量に変えるという発想の転換です。定常過程の自己相関関数 $R(\tau)=E[X(t)X(t+\tau)]$ は、多くの場合 $|\tau|\to\infty$ で0に減衰するので、絶対可積分になり、素直にフーリエ変換できます。信号は変換できないが、信号の「相関構造」は変換できる。この一段のずらしが定理の核心です。
まずはその自己相関関数の性質を、定理で使う分だけ確認しておきましょう。
定常確率過程と自己相関関数
自己相関関数を、ここでいったん自分の言葉で言い直しておきます。信号を時間軸に沿って $\tau$ だけずらして、元の信号と掛けて平均する。$\tau=0$ なら自分自身との積なので必ず大きな値になります。$\tau$ を増やしていくと、信号が「$\tau$ 秒後の自分にどれだけ似ているか」を測ることになります。ゆっくり変化する信号は少しずらしても似ているので相関が長く残り、激しく暴れる信号はすぐに似なくなって相関が急落する。この「相関がどれだけ長く尾を引くか」が、そのままスペクトルの広がりに翻訳される — これがウィーナー・ヒンチンの定理を直感的に言い換えたものです。
言葉だけでは掴みにくいので、ゆっくり変化する信号と激しく暴れる信号を並べて、波形・自己相関・スペクトルの3点セットで比べてみます。
上段（$\phi=0.95$）は波形がゆったりうねり、自己相関はラグ60でもまだ0.05ほど残り、スペクトルは低周波に3桁近い差をつけて集中しています。下段（$\phi=0.20$）は波形がギザギザで、自己相関はラグ2でほぼ0に落ち、スペクトルはほぼ平坦です。自己相関の裾の長さとスペクトルの尖り方が完全に連動していることが、この2段の対比から読み取れます。定理はこの対応を「フーリエ変換対」という厳密な形で述べたものです。
正確な定義に進みます。確率過程 $\{X(t)\}$ が広義定常（wide-sense stationary, WSS）であるとは、次の2条件を満たすことを言います。
$$ \begin{equation} E[X(t)] = \mu \quad (\text{$t$ によらない定数}) \end{equation} $$
$$ \begin{equation} E[X(t)X(t+\tau)] = R(\tau) \quad (\text{$t$ によらず、ずれ $\tau$ だけの関数}) \end{equation} $$
以降、記述を簡単にするため $\mu = 0$ とします。平均が0でない場合は $X(t)-\mu$ に置き換えればよく、その効果はスペクトルに $\mu^2\delta(f)$ という直流成分のインパルスを足すだけです（実際の測定でも、まず平均を引いてからFFTするのが定石です）。
このとき $R(\tau)$ は次の性質を持ちます。証明は前提記事に譲り、結果だけ並べます。
- 偶関数性: $R(-\tau)=R(\tau)$。定常性から $E[X(t)X(t-\tau)]=E[X(t+\tau)X(t)]$ となるためです。
- 原点最大: $|R(\tau)|\le R(0)$。コーシー・シュワルツの不等式から出ます。
- 原点の値がパワー: $R(0)=E[X(t)^2]$ は信号の平均パワーそのものです。
この3つ目が特に重要です。$R(0)$ は「全周波数のパワーの合計」に対応するはずで、後で $R(0)=\int_{-\infty}^{\infty}S(f)df$ という形で定理から自動的に出てきます。
もうひとつ、定理の証明では使わないけれど概念的に大切な点があります。定常性は「時間シフトに対する不変性」であり、フーリエ解析はまさに時間シフトを対角化する道具だということです。時間シフト作用素の固有関数は複素指数 $e^{j2\pi ft}$ ですから、時間シフトで不変な構造（定常過程）を扱うのに周波数表現が自然に現れるのは、ある意味で必然です。ウィーナー・ヒンチンの定理は、この「必然」を具体的な等式として書き下したものだと見ることができます。
さて、$R(\tau)$ の準備ができました。次は、そもそもパワースペクトル密度をどう定義すべきかを詰めます。定義が決まらないと、証明すべきことも決まりません。
パワースペクトル密度をどう定義するか
「$S(f)$ とは $R(\tau)$ のフーリエ変換のことである」と定義してしまえば、定理は自明になってしまいます。それでは中身がありません。ウィーナー・ヒンチンの定理が定理として意味を持つのは、$S(f)$ を物理的に測定できる量として独立に定義したうえで、それが $R(\tau)$ のフーリエ変換に一致することを示すからです。
では物理的な定義とは何でしょうか。スペクトルアナライザや、私たちがFFTでやっていることを思い出してください。実際にやるのは次の3ステップです。
- 信号を有限の観測窓 $[-T/2, T/2]$ で切り出す（無限時間は測れないので）
- 切り出した波形をフーリエ変換して、その絶対値の二乗を取る
- 観測時間 $T$ で割って、単位時間あたりに直す
これをそのまま数式にします。まず、有限区間で切り出した信号のフーリエ変換を
$$ \begin{equation} \hat X_T(f) = \int_{-T/2}^{T/2} X(t)\,e^{-j2\pi ft}\,dt \end{equation} $$
と書きます。積分区間が有限なので、この積分は（サンプルパスが局所可積分である限り）確実に存在します。無限区間の発散問題はここで回避されました。
次に、そのパワーを観測時間で割ります。$\frac{1}{T}|\hat X_T(f)|^2$ は「周波数 $f$ 付近に、単位時間あたりどれだけのパワーがあるか」を表す量です。これはピリオドグラムと呼ばれます。ただしこの量は確率変数であり、実現ごとに値が変わります。そこでアンサンブル平均を取り、最後に $T\to\infty$ の極限を取ります。
$$ \begin{equation} S(f) \;\equiv\; \lim_{T\to\infty}\frac{1}{T}\,E\!\left[\,\bigl|\hat X_T(f)\bigr|^2\,\right] \end{equation} $$
これがパワースペクトル密度（PSD）の定義です。$|\hat X_T(f)|^2$ ではなく $\frac{1}{T}|\hat X_T(f)|^2$ にしたところが要点で、これによって「エネルギー密度」ではなく「パワー密度」になり、$T\to\infty$ でも有限に留まります。単位も確認しておくと、$X$ が電圧 [V] なら $\hat X_T$ は [V·s]、その二乗は [V²s²]、$T$ で割って [V²s] = [V²/Hz] となり、確かに「周波数あたりのパワー」の次元です。
なお、期待値と極限の順序は重要です。期待値を取らずに $\lim_{T\to\infty}\frac1T|\hat X_T(f)|^2$ としても、この極限は一般に存在しません（後半の数値実験で、$T$ を増やしてもバラつきが減らないことを見ます）。先にアンサンブル平均でランダム性をならしてから、時間を伸ばす。この順序でだけ、まともな極限が得られます。
この定義の4ステップを、実際の数値で追ってみます。
左は観測窓で切り出す操作、中央はその1回ぶんのピリオドグラムです。中央の橙色の線は真のスペクトル（黒破線）の周りで1桁以上も上下に暴れており、これを1本だけ見てスペクトルを語ることはできません。ところが右のように400回のアンサンブル平均を取ると、橙色の線は黒破線とほぼ完全に重なります。期待値 $E[\cdot]$ を定義に入れているのは飾りではなく、この暴れを消すために不可欠なのです。
定義が決まりました。ここからが本題です。この「測定手順をそのまま数式にした量」が、なぜ自己相関のフーリエ変換という、まったく別の顔をした量に一致するのか。導出していきましょう。
ウィーナー・ヒンチンの定理の導出
示すべきこと: 平均0の広義定常過程 $X(t)$ の自己相関 $R(\tau)$ が絶対可積分（$\int_{-\infty}^{\infty}|R(\tau)|d\tau<\infty$）ならば、上で定義した $S(f)$ は
$$ \begin{equation} S(f) = \int_{-\infty}^{\infty} R(\tau)\,e^{-j2\pi f\tau}\,d\tau \end{equation} $$
を満たす。すなわち $R(\tau)$ と $S(f)$ はフーリエ変換対をなす。
ステップ1: 絶対値の二乗を二重積分に開く
$X(t)$ は実数値なので $\overline{\hat X_T(f)} = \int_{-T/2}^{T/2}X(s)e^{+j2\pi fs}ds$ です。したがって
$$ \begin{align} \bigl|\hat X_T(f)\bigr|^2 &= \hat X_T(f)\,\overline{\hat X_T(f)} \\ &= \left(\int_{-T/2}^{T/2}X(t)e^{-j2\pi ft}dt\right)\left(\int_{-T/2}^{T/2}X(s)e^{+j2\pi fs}ds\right) \end{align} $$
2つの積分をまとめて二重積分にします。積分変数が $t$ と $s$ で独立なので、そのまま中に入れられます。
$$ \begin{equation} \bigl|\hat X_T(f)\bigr|^2 = \int_{-T/2}^{T/2}\!\!\int_{-T/2}^{T/2} X(t)X(s)\,e^{-j2\pi f(t-s)}\,dt\,ds \end{equation} $$
指数部が $e^{-j2\pi f(t-s)}$ と、$t$ と $s$ の差だけに依存する形になったことに注目してください。これが後で効いてきます。
ステップ2: 期待値を積分の中に入れる
両辺の期待値を取ります。積分区間が有限で、$E[|X(t)X(s)|]\le R(0)$ が有界なのでフビニの定理が使え、期待値と積分の順序を交換できます。
$$ \begin{align} E\!\left[\bigl|\hat X_T(f)\bigr|^2\right] &= \int_{-T/2}^{T/2}\!\!\int_{-T/2}^{T/2} E[X(t)X(s)]\,e^{-j2\pi f(t-s)}\,dt\,ds \\ &= \int_{-T/2}^{T/2}\!\!\int_{-T/2}^{T/2} R(t-s)\,e^{-j2\pi f(t-s)}\,dt\,ds \end{align} $$
ここで定常性を使いました。$E[X(t)X(s)]$ は本来2変数の関数ですが、定常過程では差 $t-s$ だけの関数 $R(t-s)$ になります。被積分関数全体が $t-s$ だけの関数になった — これがステップ3を可能にします。
ステップ3: 変数変換で二重積分を一重積分に潰す
被積分関数を $g(t-s) \equiv R(t-s)e^{-j2\pi f(t-s)}$ と略記します。示したいのは次の等式です。
$$ \begin{equation} \int_{-T/2}^{T/2}\!\!\int_{-T/2}^{T/2} g(t-s)\,dt\,ds = \int_{-T}^{T}\bigl(T-|\tau|\bigr)\,g(\tau)\,d\tau \end{equation} $$
変数変換 $\tau = t-s$、$v = s$ を行います。ヤコビアンは
$$ \frac{\partial(t,s)}{\partial(\tau,v)} = \begin{vmatrix} 1 & 1 \\ 0 & 1\end{vmatrix} = 1 $$
なので、測度はそのままです。問題は積分領域です。元の領域は $(t,s)$ 平面の一辺 $T$ の正方形 $[-T/2,T/2]^2$ でした。新しい変数では $t = \tau+v$、$s = v$ なので、条件は
$$ -\frac{T}{2}\le v\le \frac{T}{2}, \qquad -\frac{T}{2}\le \tau+v\le \frac{T}{2} $$
の2つです。後者を $v$ について解くと $-\frac{T}{2}-\tau\le v\le \frac{T}{2}-\tau$ となります。したがって $\tau$ を固定したとき、$v$ が動ける範囲は2つの区間の共通部分
$$ \max\left(-\frac{T}{2},\,-\frac{T}{2}-\tau\right)\;\le\; v \;\le\; \min\left(\frac{T}{2},\,\frac{T}{2}-\tau\right) $$
です。この区間の長さを計算します。$\tau\ge 0$ のときは下限が $-T/2$、上限が $T/2-\tau$ なので長さは $T-\tau$。$\tau<0$ のときは下限が $-T/2-\tau$、上限が $T/2$ なので長さは $T+\tau$。まとめると長さは $T-|\tau|$ であり、$|\tau|>T$ では共通部分が空になって長さ0です。
図形的に言えば、正方形を対角線の方向（$t-s$ が一定な方向）にスライスしているわけです。対角線上のスライスが最も長く（長さ $T$）、角に近づくほど短くなり、$|\tau|=T$ でちょうど点に縮む。この「長さの分布」が三角形になるのが、次に現れる係数 $T-|\tau|$ の正体です。図で見ると一目瞭然です。
左が $(t,s)$ 平面の積分領域で、色つきの直線がそれぞれ $\tau=t-s$ を固定したスライスです。$\tau=0$ の対角線は長さ $T=2$ ですが、$\tau=1.5$ まで離れると長さは $0.5$ しか残りません。右はその長さを $\tau$ の関数として並べたもので、きれいな三角形になります。この三角形は仮定したものではなく、正方形を斜めに切るという幾何から自動的に出てくるものだという点が重要です。
以上より、$g$ は $\tau$ だけの関数なので $v$ 積分はただ区間の長さを掛けるだけになり、
$$ \begin{equation} E\!\left[\bigl|\hat X_T(f)\bigr|^2\right] = \int_{-T}^{T}\bigl(T-|\tau|\bigr)R(\tau)\,e^{-j2\pi f\tau}\,d\tau \end{equation} $$
を得ます。
ステップ4: $T$ で割って三角窓の形にする
定義に合わせて両辺を $T$ で割ります。$\frac{T-|\tau|}{T}=1-\frac{|\tau|}{T}$ と書き直すと
$$ \begin{equation} \frac{1}{T}E\!\left[\bigl|\hat X_T(f)\bigr|^2\right] = \int_{-T}^{T}\left(1-\frac{|\tau|}{T}\right)R(\tau)\,e^{-j2\pi f\tau}\,d\tau \end{equation} $$
この形は非常に示唆的です。右辺は「$R(\tau)$ に三角形の窓 $w_T(\tau)=1-|\tau|/T$ を掛けてからフーリエ変換したもの」です。この三角窓はバートレット窓と呼ばれます。有限時間しか観測していないという事実が、自己相関に三角窓が掛かるという形で正確に現れているのです。観測時間が短いほど窓は鋭く尖り、遠いラグの相関情報が失われます。
この「窓が掛かる」という事実が、スペクトル側で何を引き起こすのかを見ておきます。
左のグラフでは、真の自己相関（黒）に対して $T=6$ の見え方（赤）が大きく削られ、$T=40$（青）ではほぼ黒に重なっています。右はそれをフーリエ変換した結果で、$T=6$ ではピーク高さが真値12に対し4.4程度まで潰れ、代わりに裾が広がっています。観測時間が短いほどスペクトルのピークは低く太くなる — 有限データのバイアスは、この段階ですでに完全に説明されています。
ステップ5: 極限を取る
いよいよ $T\to\infty$ とします。被積分関数を $h_T(\tau)$ と書くと
$$ h_T(\tau) = \left(1-\frac{|\tau|}{T}\right)\mathbb{1}_{\{|\tau|\le T\}}\,R(\tau)e^{-j2\pi f\tau} $$
です。ここで2つのことを確認します。第一に、各 $\tau$ を固定すると $T$ を十分大きく取れば $1-|\tau|/T\to 1$ かつ $\mathbb{1}_{\{|\tau|\le T\}}=1$ なので、$h_T(\tau)\to R(\tau)e^{-j2\pi f\tau}$ が各点で成り立ちます。第二に、$0\le 1-|\tau|/T\le 1$ なので
$$ |h_T(\tau)| \le |R(\tau)| $$
が $T$ によらず成り立ち、仮定より $|R|$ は可積分です。つまり可積分な優関数が存在します。ルベーグの優収束定理により、極限と積分を交換できて
$$ \begin{equation} S(f) = \lim_{T\to\infty}\frac{1}{T}E\!\left[\bigl|\hat X_T(f)\bigr|^2\right] = \int_{-\infty}^{\infty}R(\tau)\,e^{-j2\pi f\tau}\,d\tau \end{equation} $$
が示されました。$\square$
そして $R$ が可積分かつ連続なら、フーリエの反転公式から逆向きの関係も成り立ちます。
$$ \begin{equation} R(\tau) = \int_{-\infty}^{\infty} S(f)\,e^{+j2\pi f\tau}\,df \end{equation} $$
導出を振り返ると、本質的な仕事をしたのはステップ2とステップ3でした。ステップ2では定常性によって被積分関数が差だけの関数になり、ステップ3ではその構造のおかげで二重積分が一重積分に潰れ、そこで三角形の重み $T-|\tau|$ が自動的に生まれました。ステップ4以降は極限操作の始末です。定常性が二重積分を一重積分に落とす、これがウィーナー・ヒンチンの定理のエンジンだと理解しておくと、細部を忘れても再構成できます。
ところで、実際に私たちが扱うのは離散サンプル列です。同じ論法を離散時間に翻訳すると、有限データ長のFFTスペクトル推定の性質がそのまま出てきます。次はそれを見ましょう。
離散時間版とピリオドグラムの期待値
離散時間の広義定常過程 $\{x[n]\}$（平均0、自己相関 $R[k]=E[x[n]x[n+k]]$）を考えます。サンプリング周期を1に正規化し、正規化周波数 $f\in[-1/2,1/2)$ を使います。長さ $N$ のデータ $x[0],\dots,x[N-1]$ から作るピリオドグラムは
$$ \begin{equation} I_N(f) = \frac{1}{N}\left|\sum_{n=0}^{N-1}x[n]e^{-j2\pi fn}\right|^2 \end{equation} $$
です。連続時間と同じ手順を踏みます。まず絶対値の二乗を二重和に開き、期待値を取ります。
$$ \begin{align} E[I_N(f)] &= \frac{1}{N}\sum_{n=0}^{N-1}\sum_{m=0}^{N-1}E[x[n]x[m]]\,e^{-j2\pi f(n-m)} \\ &= \frac{1}{N}\sum_{n=0}^{N-1}\sum_{m=0}^{N-1}R[n-m]\,e^{-j2\pi f(n-m)} \end{align} $$
ここでも被加数が差 $n-m$ だけの関数になりました。連続時間では対角線に沿って積分領域をスライスしましたが、離散では$n-m=k$ となる格子点の個数を数えるだけです。$n,m\in\{0,\dots,N-1\}$ という $N\times N$ の格子で、差が $k$ になる点の個数はいくつでしょうか。$k\ge0$ なら $n=m+k$ かつ $0\le m\le N-1-k$ なので $N-k$ 個。$k<0$ なら同様に $N+k$ 個。まとめて $N-|k|$ 個（$|k|\le N-1$）です。したがって
$$ \begin{equation} E[I_N(f)] = \sum_{k=-(N-1)}^{N-1}\frac{N-|k|}{N}R[k]\,e^{-j2\pi fk} = \sum_{k=-(N-1)}^{N-1}\left(1-\frac{|k|}{N}\right)R[k]\,e^{-j2\pi fk} \end{equation} $$
連続時間のステップ4とまったく同じ形が出ました。バートレット窓 $w_B[k]=1-|k|/N$ が自己相関に掛かっています。$N\to\infty$ で $\sum_k|R[k]|<\infty$ なら、これは
$$ \begin{equation} S(f) = \sum_{k=-\infty}^{\infty}R[k]\,e^{-j2\pi fk}, \qquad R[k]=\int_{-1/2}^{1/2}S(f)e^{+j2\pi fk}df \end{equation} $$
に収束します。これが離散時間版のウィーナー・ヒンチンの定理です。
バイアスの正体はフェイエール核による平滑化
有限 $N$ でのズレをもう少し詳しく見ておきましょう。ラグ領域で「$R[k]$ と $w_B[k]$ の積」になっているということは、周波数領域では「$S(f)$ と $w_B$ のDTFTの畳み込み」になります。バートレット窓のDTFTは
$$ \mathcal{F}_N(\lambda) = \sum_{k=-(N-1)}^{N-1}\left(1-\frac{|k|}{N}\right)e^{-j2\pi\lambda k} = \frac{1}{N}\left|\sum_{n=0}^{N-1}e^{-j2\pi\lambda n}\right|^2 = \frac{1}{N}\left(\frac{\sin(\pi N\lambda)}{\sin(\pi\lambda)}\right)^2 $$
です。中央の等号は、バートレット列が矩形列の自己相関を $N$ で割ったものであることから従います（相関のDTFTは絶対値の二乗）。この $\mathcal{F}_N$ はフェイエール核と呼ばれ、非負で、1周期の積分が $w_B[0]=1$ になります。結局
$$ \begin{equation} E[I_N(f)] = \int_{-1/2}^{1/2}S(\nu)\,\mathcal{F}_N(f-\nu)\,d\nu \end{equation} $$
つまりピリオドグラムの期待値は、真のスペクトルを幅 $\sim 1/N$ の非負のカーネルでぼかしたものです。$N$ が大きくなるほどカーネルは細くなり、デルタ関数に近づいて $E[I_N(f)]\to S(f)$ となります。これがピリオドグラムのバイアスの正体で、スペクトルの鋭いピークが有限データでは丸まって見える理由でもあります。
左がフェイエール核そのもので、$N=8$ では幅が $0.1$ 以上ある鈍い山ですが、$N=128$ では高さ128の鋭い針になります。値がどこでも0以上であること（負に潜らないこと）も見て取れます。右はこの核で鋭いピークを持つスペクトルをぼかした結果で、$N=16$ ではピーク高さ3.20が0.54まで潰れ、$N=64$ で1.43、$N=256$ でようやく2.67まで戻ります。データ長 $N$ が周波数分解能を決めているのは、このカーネル幅 $\sim1/N$ そのものです。
有限データでの「厳密な」ウィーナー・ヒンチン恒等式
もうひとつ、実務で知っておくと得をする事実があります。ここまでは期待値の話でしたが、実は期待値を取らなくても成り立つ代数的な恒等式が存在します。バイアスありの標本自己相関を
$$ \hat R_b[k] = \frac{1}{N}\sum_{n=0}^{N-1-|k|}x[n]\,x[n+|k|], \qquad |k|\le N-1 $$
と定義すると、そのDTFTはピリオドグラムに厳密に一致します。
$$ \begin{equation} \sum_{k=-(N-1)}^{N-1}\hat R_b[k]\,e^{-j2\pi fk} = \frac{1}{N}\left|\sum_{n=0}^{N-1}x[n]e^{-j2\pi fn}\right|^2 = I_N(f) \end{equation} $$
理由は先ほどと同じ数え上げです。左辺を展開すると $\frac1N\sum_n\sum_m x[n]x[m]e^{-j2\pi f(n-m)}$ になり、これは右辺の二重和そのものです。つまり「標本自己相関をFFTする」ことと「データを直接FFTして二乗する」ことは、有限データでも同じ計算なのです。この恒等式は後のPython実験でも数値的に確認します（差は典型的に $10^{-14}$、最大でも $10^{-13}$ に収まります）。
左は $N=256$ のAR(1)データ1本について、$\hat R_b[k]$ のDTFTとピリオドグラムの差を全周波数でプロットしたものです。縦軸は対数で、値は $10^{-15}$ から $10^{-13}$ の範囲、つまり倍精度の丸め誤差そのものです。期待値を取る前の、たった1本のデータでも両者は同じ計算だということが数値的に確認できます。右のヒストグラムは後半の伏線で、ピリオドグラムの値が指数分布に従う（＝ばらつきが $N$ を増やしても減らない）ことを示しています。
また、$\hat R_b[k]$ の期待値は
$$ E[\hat R_b[k]] = \frac{N-|k|}{N}R[k] = \left(1-\frac{|k|}{N}\right)R[k] $$
です。ここでもバートレット窓が現れます。$1/N$ ではなく $1/(N-|k|)$ で割る「バイアスなし推定量」を使うと期待値は $R[k]$ に一致しますが、その代わりFFTしたときに負の値が出うるという困った性質を持ちます。有限データでスペクトル推定をするとき、あえてバイアスのある $\hat R_b$ を使うのは、非負性という物理的に必須な性質を守るためです。
ここまでで定理の中身と、有限データでの姿が見えました。次は、この定理から自動的に転がり出てくる帰結を整理します。
定理から出てくる4つの帰結
定理を証明すると、単に「変換対である」という以上のことが分かります。導出の途中式そのものが情報を持っているからです。
帰結1: $S(f)\ge 0$（非負性）
定義に戻ります。$S(f)=\lim_{T\to\infty}\frac1T E[|\hat X_T(f)|^2]$ の中身は絶対値の二乗の期待値ですから、どの $T$ でも非負です。非負な量の極限は非負なので
$$ S(f)\ge 0 \quad \text{for all } f $$
が直ちに言えます。これは物理的にも当然で、「ある周波数帯のパワーが負」ということはありえません。ただし数学的には決して自明ではないことを強調しておきます。$R(\tau)$ という時間領域の関数をフーリエ変換したら必ず非負になる、というのは $R$ に強い制約が課されていることを意味します。その制約こそが、次に述べる正定値性です。
帰結2: $S(f)$ は実数の偶関数
$R(\tau)$ が実数の偶関数なので、そのフーリエ変換は
$$ S(f)=\int_{-\infty}^{\infty}R(\tau)\bigl(\cos 2\pi f\tau – j\sin2\pi f\tau\bigr)d\tau = \int_{-\infty}^{\infty}R(\tau)\cos(2\pi f\tau)\,d\tau $$
となります。虚部の積分は「偶関数×奇関数＝奇関数」を対称区間で積分するので0です。したがって $S(f)$ は実数値であり、$\cos$ が偶関数なので $S(-f)=S(f)$ も従います。実信号のスペクトルが左右対称になり、片側スペクトルだけ表示すれば十分なのはこのためです。
帰結3: $R(0)=\int S(f)df$（パワーの保存）
逆変換の式 $R(\tau)=\int S(f)e^{j2\pi f\tau}df$ に $\tau=0$ を代入するだけです。
$$ \begin{equation} R(0) = E[X(t)^2] = \int_{-\infty}^{\infty}S(f)\,df \end{equation} $$
左辺は信号の全平均パワー、右辺は各周波数のパワー密度を全周波数で足し上げたもの。両者が等しいという、パーシバルの等式のパワー版です。$S(f)$ を「密度」と呼ぶ根拠がまさにこれで、$S(f)df$ が微小帯域 $[f,f+df]$ に含まれるパワーだと解釈できます。
この式は実務でそのまま使われます。たとえば白色雑音のPSDが $N_0/2$ [W/Hz]（両側）で、帯域幅 $B$ の理想フィルタを通したなら、通過後のパワーは $\int_{-B}^{B}\frac{N_0}{2}df = N_0B$ [W] です。受信機の雑音電力の基本式はこうして出てきます。
帰結4: フィルタを通したときのスペクトル変換
証明はここでは省きますが、同じ枠組みから重要な系が出ます。インパルス応答 $h(t)$、周波数応答 $H(f)$ のLTIシステムに定常過程 $X$ を入れると、出力 $Y$ のPSDは
$$ \begin{equation} S_Y(f) = |H(f)|^2\,S_X(f) \end{equation} $$
になります。時間領域では出力の自己相関が $R_Y = h(\tau)*h(-\tau)*R_X(\tau)$ という三重の畳み込みになって扱いにくいのですが、周波数領域では単なる掛け算です。フィルタ設計や雑音伝搬の計算がスペクトルで行われるのは、この単純さゆえです。ウィーナー・ヒンチンの定理は、この便利な公式を使うための前提を提供しています。
帰結3と帰結4を、数値で確かめておきましょう。
左はAR(1)（$\phi=0.7$）のスペクトルで、塗りつぶした面積を数値積分すると1.9608、これは $R[0]=\sigma^2/(1-\phi^2)=1.9608$ と小数点以下4桁まで一致します。右は白色雑音を通過域 $0.12\le|f|\le0.22$ のFIRフィルタに通した結果で、出力のピリオドグラム平均（橙）が理論の $|H(f)|^2$（黒破線）に4桁にわたって重なっています。阻止域の底が $10^{-4}$ 付近で止まって黒破線から離れるのは実測の誤りではなく、有限 $N$ のバートレット窓つき期待値（赤）がまさにその高さを予言しているためで、ここでもフェイエール核の裾が効いていることが分かります。
ここまでは定理から下る方向でした。逆に「どんな関数が自己相関になりうるか」を問うと、より深い構造が見えてきます。
自己相関の正定値性とボホナーの定理
適当に関数を書いて「これを自己相関に持つ定常過程はあるか」と問うと、答えはほとんどの場合ノーです。たとえば $R(\tau)$ が三角形（$1-|\tau|$ が $|\tau|\le1$、外は0）ならOKですが、$R(\tau)$ が矩形（$|\tau|\le1$ で1、外は0）だとダメです。なぜダメかは、フーリエ変換してみればすぐ分かります。矩形のフーリエ変換は $\mathrm{sinc}$ で、負の値を取るからです。帰結1に反します。
上段が矩形の候補で、右のフーリエ変換は $f\approx\pm0.72$ で $-0.44$ まで潜り込んでいます。赤く塗った領域は「パワーが負の周波数帯」であり、物理的にありえません。下段の三角形の候補は、フーリエ変換が $\mathrm{sinc}^2$ の形になって全周波数で0以上です。自己相関になれるかどうかは、フーリエ変換が非負かどうかで完全に判定できる — これが次に述べる正定値性とボホナーの定理の内容です。
この「自己相関になれる条件」を正面から定式化したのが正定値性です。
定義: 関数 $R:\mathbb{R}\to\mathbb{C}$ が正定値（positive definite）であるとは、任意の自然数 $n$、任意の時刻 $t_1,\dots,t_n$、任意の複素数 $a_1,\dots,a_n$ に対して
$$ \begin{equation} \sum_{i=1}^{n}\sum_{j=1}^{n} a_i\,\overline{a_j}\,R(t_i-t_j) \;\ge\; 0 \end{equation} $$
が成り立つことを言う。
定常過程の自己相関は必ず正定値です。証明は一行で済みます。$Z=\sum_{i=1}^n a_iX(t_i)$ という複素数値の確率変数を作ると
$$ \begin{align} 0 \le E\bigl[|Z|^2\bigr] &= E\left[\left(\sum_i a_iX(t_i)\right)\overline{\left(\sum_j a_jX(t_j)\right)}\right] \\ &= \sum_i\sum_j a_i\overline{a_j}\,E[X(t_i)X(t_j)] \\ &= \sum_i\sum_j a_i\overline{a_j}\,R(t_i-t_j) \end{align} $$
となるからです。二乗の期待値は非負、それだけです。
面白いのはここからで、この正定値性と $S(f)\ge0$ は完全に同じことを言っています。逆変換の式 $R(\tau)=\int S(f)e^{j2\pi f\tau}df$ を正定値性の左辺に代入してみましょう。
$$ \sum_i\sum_j a_i\overline{a_j}R(t_i-t_j) = \sum_i\sum_j a_i\overline{a_j}\int_{-\infty}^{\infty}S(f)e^{j2\pi f(t_i-t_j)}df $$
積分と有限和を交換し、指数を $e^{j2\pi ft_i}\cdot\overline{e^{j2\pi ft_j}}$ と分けます。すると和が「あるものとその共役の積」の形になり、絶対値の二乗にまとまります。
$$ \begin{equation} = \int_{-\infty}^{\infty}S(f)\left(\sum_i a_ie^{j2\pi ft_i}\right)\overline{\left(\sum_j a_je^{j2\pi ft_j}\right)}df = \int_{-\infty}^{\infty}S(f)\left|\sum_i a_ie^{j2\pi ft_i}\right|^2 df \end{equation} $$
被積分関数のうち $|\cdot|^2$ の部分は常に非負です。したがって $S(f)\ge0$ ならこの積分は必ず非負となり、正定値性が言えます。逆に、もしある区間で $S(f)<0$ なら、$\left|\sum_ia_ie^{j2\pi ft_i}\right|^2$ をその区間に集中させるように $a_i$ を選ぶことで積分を負にできるので、正定値性が破れます。
この同値性を厳密な形で述べたのがボホナーの定理です。「連続関数 $R(\tau)$ が正定値であることと、$R$ がある有限非負測度のフーリエ変換であることは同値」。離散時間版はヘルグロッツの定理と呼ばれます。ウィーナー・ヒンチンの定理は、この一般的な対応を「定常確率過程の自己相関」という具体的な文脈に落としたもの、と位置づけられます。
この視点には実用的な副産物もあります。第一に、$R$ が可積分でない場合（たとえば周期成分を含む過程）でも、$S$ を測度として扱えば定理が生き延びます。次節の「正弦波＋雑音」がまさにその例で、$S$ が連続な密度部分と離散的な線（デルタ）の和になります。第二に、ガウス過程回帰のカーネル設計もこの定理の応用です。定常カーネル $k(\tau)$ が有効（半正定値）であることを保証するには、そのフーリエ変換が非負であればよい。RBFカーネルが使えるのはガウス関数のフーリエ変換がガウス関数（非負）だからで、スペクトル混合カーネルは「スペクトル側で非負な関数を設計してから逆変換する」という発想そのものです。
理屈が揃いました。ここからは具体例で、$R$ と $S$ の対応を手で確かめていきます。
3つの具体例
例1: 白色雑音 — 相関がない＝スペクトルが平坦
離散時間の白色雑音 $w[n]$ は、平均0、分散 $\sigma^2$、異なる時刻では無相関という過程です。自己相関は
$$ R_w[k] = \sigma^2\delta[k] = \begin{cases}\sigma^2 & k=0\\ 0 & k\ne 0\end{cases} $$
これをフーリエ変換すると、和のうち $k=0$ の項だけが生き残って
$$ \begin{equation} S_w(f) = \sum_{k}\sigma^2\delta[k]e^{-j2\pi fk} = \sigma^2 \quad (\text{全ての } f) \end{equation} $$
完全に平坦です。「白色」という呼び名は、すべての周波数成分を等しく含む白色光からの類推です。逆にこの結果は直感的にも読めます。「隣のサンプルとの相関がまったくない」＝「どんなにゆっくりした構造もどんなに速い構造も特別扱いしない」＝「全周波数が均等」。
検算もしておきましょう。帰結3より $R_w[0]=\int_{-1/2}^{1/2}S_w(f)df=\sigma^2\cdot 1=\sigma^2$。確かに $R_w[0]=\sigma^2$ と一致します。正規化周波数の帯域幅が1だったので、密度と全パワーが同じ数値になっています。
例2: AR(1)過程 — 相関の減衰速度がスペクトルの傾きを決める
1次自己回帰過程
$$ x[n] = \phi\,x[n-1] + w[n], \qquad |\phi|<1,\quad w[n]\sim \text{WN}(0,\sigma^2) $$
を考えます。この過程の自己相関は（詳しい導出は前提記事に譲ります）
$$ \begin{equation} R[k] = \frac{\sigma^2}{1-\phi^2}\,\phi^{|k|} \end{equation} $$
です。指数的に減衰する形で、$\phi$ が1に近いほど相関が長く尾を引きます。これをフーリエ変換しましょう。$A=\sigma^2/(1-\phi^2)$ と置き、$k\ge0$ と $k<0$ に分けます。
$$ S(f) = A\sum_{k=-\infty}^{\infty}\phi^{|k|}e^{-j2\pi fk} = A\left(\sum_{k=0}^{\infty}\phi^ke^{-j2\pi fk} + \sum_{k=1}^{\infty}\phi^ke^{+j2\pi fk}\right) $$
$z=e^{-j2\pi f}$ と書けば、それぞれ公比 $\phi z$ と $\phi \bar z$ の等比級数です。$|\phi|<1$ なので収束して
$$ S(f) = A\left(\frac{1}{1-\phi z} + \frac{\phi\bar z}{1-\phi\bar z}\right) $$
通分します。分母は $(1-\phi z)(1-\phi\bar z)=|1-\phi z|^2$ で、分子は $(1-\phi\bar z)+\phi\bar z(1-\phi z)=1-\phi^2 z\bar z = 1-\phi^2$（$z\bar z=1$ を使いました）。よって
$$ S(f) = \frac{\sigma^2}{1-\phi^2}\cdot\frac{1-\phi^2}{|1-\phi e^{-j2\pi f}|^2} = \frac{\sigma^2}{|1-\phi e^{-j2\pi f}|^2} $$
分母を展開すると $|1-\phi e^{-j2\pi f}|^2 = (1-\phi\cos2\pi f)^2+(\phi\sin2\pi f)^2 = 1-2\phi\cos(2\pi f)+\phi^2$ なので、最終形は
$$ \begin{equation} S(f) = \frac{\sigma^2}{1-2\phi\cos(2\pi f)+\phi^2} \end{equation} $$
きれいな結果です。この式は、AR過程の伝達関数 $H(z)=1/(1-\phi z^{-1})$ を使った $S(f)=|H(e^{j2\pi f})|^2\sigma^2$（帰結4）とも一致しています。白色雑音をフィルタに通した出力、という見方がそのまま数式になっています。
具体的な数値を入れてみます。$\phi=0.7$、$\sigma^2=1$ とすると
- $S(0) = 1/(1-1.4+0.49) = 1/0.09 = 11.111$
- $S(0.25) = 1/(1-0+0.49) = 1/1.49 = 0.671$
- $S(0.5) = 1/(1+1.4+0.49) = 1/3.89 = 0.346$
低周波が高周波の32倍という、はっきりした低域集中型のスペクトルです。$\phi>0$ は「前の値を引きずる」動きなので、ゆっくりした変動が強調される — 直感と合っています。全パワーは $R[0]=\sigma^2/(1-\phi^2)=1/0.51=1.961$ で、これは $\int_{-1/2}^{1/2}S(f)df$ に等しいはずです。数値積分で確認すると 1.9608 となり、確かに一致します。
例3: 正弦波＋白色雑音 — 密度と「線」の共存
最後に、$R$ が絶対総和可能でない例を見ます。
$$ x[n] = A\cos(2\pi f_0 n + \Theta) + w[n], \qquad \Theta\sim \text{Uniform}(0,2\pi) $$
位相 $\Theta$ を一様乱数にするのは、過程を定常にするためです（位相を固定すると平均が時刻に依存してしまい定常でなくなります）。自己相関を計算します。積和公式 $\cos\alpha\cos\beta=\frac12[\cos(\alpha-\beta)+\cos(\alpha+\beta)]$ を使うと
$$ E\bigl[A^2\cos(2\pi f_0n+\Theta)\cos(2\pi f_0(n+k)+\Theta)\bigr] = \frac{A^2}{2}\cos(2\pi f_0k) + \frac{A^2}{2}E[\cos(2\pi f_0(2n+k)+2\Theta)] $$
第2項は $\Theta$ について一様分布で平均すると0になります（$\cos$ を1周期にわたって平均すると0）。したがって信号成分と雑音成分が独立であることも使って
$$ \begin{equation} R[k] = \frac{A^2}{2}\cos(2\pi f_0 k) + \sigma^2\delta[k] \end{equation} $$
ここで問題が起きます。第1項は $k\to\infty$ でも減衰せず、いつまでも振動し続けます。$\sum_k|R[k]|=\infty$ なので、これまでの証明の仮定（絶対総和可能）が破れます。しかしボホナー・ヘルグロッツの立場に立てば、$S$ を測度として扱うことで対応できます。$\cos(2\pi f_0k)=\frac12(e^{j2\pi f_0k}+e^{-j2\pi f_0k})$ を「デルタ関数のフーリエ変換」と読めば
$$ \begin{equation} S(f) = \frac{A^2}{4}\bigl[\delta(f-f_0)+\delta(f+f_0)\bigr] + \sigma^2 \end{equation} $$
（$[-1/2,1/2)$ で周期的に繰り返す）となります。連続な密度 $\sigma^2$（雑音）の上に、$\pm f_0$ に離散的な線スペクトル（正弦波）が立っている構造です。
この「線」は密度ではないので、ピリオドグラムでの現れ方が雑音とは根本的に違います。$Nf_0$ が整数になるようにビンを取ると、$f_0$ でのピリオドグラムの期待値は
$$ E[I_N(f_0)] \approx \frac{N A^2}{4} + \sigma^2 $$
となり、$N$ に比例して伸び続けます。密度部分（雑音の $\sigma^2$）は $N$ を増やしても一定なのに、線の高さだけが伸びる。これが「密度と線の違い」の数値的な現れです。後のPython実験で $A=1$、$f_0=0.15625$、$\sigma=0.5$、$N=512$ として実測すると、ピーク値は約128.25となり、予測値 $512\times1/4+0.25=128.25$ とぴたりと合います。逆にこの性質を使って、雑音に埋もれた正弦波を「$N$ を増やせば必ず検出できる」と保証するのが、コヒーレント積分の原理です。
赤い線（線スペクトルのピーク）は両対数軸で傾き1の直線に乗り、実測値は $N=64$ で16.28（予測16.25）、$N=1024$ で256.05（予測256.25）と、きっちり $N$ に比例して伸びています。一方、青い線（雑音床）は $N$ を16倍にしても0.24〜0.26の範囲にとどまり、$\sigma^2=0.25$ の水平線から動きません。「線」と「密度」はピリオドグラムの上でまったく違う振る舞いをするため、$N$ を伸ばすほど正弦波は雑音から相対的に浮き上がっていきます。
3つの例で $R$ と $S$ の対応を手計算しました。次はこれをPythonで数値的に確かめます。手計算した理論曲線と、データから推定した曲線が本当に重なるのかを見ましょう。
Pythonでの検証
以下では3種類の信号について、次の3本の曲線を重ねて比較します。
- (a) 標本自己相関のFFT — 各実現から $\hat R_b[k]$ を計算し、多数実現で平均してからDTFTしたもの
- (b) ピリオドグラムの多数回平均 — 各実現から $I_N(f)$ を計算し、多数実現で平均したもの
- (c) 理論スペクトル — 手計算した $S(f)$
定理が正しければ、(a)(b)(c)は重なるはずです。しかも(a)と(b)は先ほど示した恒等式によって厳密に一致し、(c)とはフェイエール核による平滑化のぶんだけわずかにずれる、という細かい予測も立ちます。
共通の準備
まず日本語フォントの設定と、共通で使う関数を定義します。
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
# 日本語フォント設定
for cand in ["Hiragino Sans", "Yu Gothic", "Noto Sans CJK JP", "IPAexGothic", "Meiryo"]:
if any(cand == f.name for f in matplotlib.font_manager.fontManager.ttflist):
plt.rcParams["font.family"] = cand
break
plt.rcParams["axes.unicode_minus"] = False
def ar1_process(n, phi, sigma, rng, burn=300):
"""AR(1)過程 x[k] = phi*x[k-1] + w[k] を生成（初期過渡を burn 分捨てる）"""
w = sigma * rng.standard_normal(n + burn)
y = np.zeros(n + burn)
for i in range(1, n + burn):
y[i] = phi * y[i - 1] + w[i]
return y[burn:]
def periodogram(x):
"""ピリオドグラム I_N(f) = |FFT(x)|^2 / N（周波数は fftfreq の順）"""
N = len(x)
return np.abs(np.fft.fft(x)) ** 2 / N
def biased_acf(x):
"""バイアスありの標本自己相関 R_b[k]（ラグ -(N-1)..(N-1)）"""
N = len(x)
return np.correlate(x, x, mode="full") / N
biased_acf で np.correlate の結果を $N$ で割っているのがポイントです。np.correlate はラグ $k$ で $N-|k|$ 個の積を足すので、$N$ で割ると本文の $\hat R_b[k]$ と一致します（$N-|k|$ で割ればバイアスなし推定量になります）。ここで $N$ を選んだのは、非負性を保証してくれるバイアスあり版だからです。
例1: 白色雑音での検証
もっとも単純な白色雑音から始めます。理論は $R[k]=\sigma^2\delta[k]$、$S(f)=\sigma^2$（平坦）です。
import numpy as np
import matplotlib.pyplot as plt
sigma, N, M = 1.0, 256, 2000 # 標準偏差 / データ長 / 実現数
rng = np.random.default_rng(0)
acc_per = np.zeros(N) # ピリオドグラムの累積
acc_acf = np.zeros(2 * N - 1) # 標本自己相関の累積
for _ in range(M):
x = sigma * rng.standard_normal(N)
acc_per += periodogram(x)
acc_acf += biased_acf(x)
avg_per = acc_per / M
avg_acf = acc_acf / M
f = np.fft.fftfreq(N) # 正規化周波数
lags = np.arange(-(N - 1), N)
# (a) 平均標本自己相関のDTFT
acf_spec = np.array([np.sum(avg_acf * np.exp(-2j * np.pi * ff * lags)).real
for ff in f])
theory = np.full_like(f, sigma ** 2) # (c) 理論スペクトル
order = np.argsort(f)
plt.figure(figsize=(10, 4.5))
plt.plot(f[order], acf_spec[order], lw=3.0, alpha=0.45, label="(a) 標本自己相関のFFT")
plt.plot(f[order], avg_per[order], lw=1.2, label=f"(b) ピリオドグラム{M}回平均")
plt.plot(f[order], theory[order], "k--", lw=1.8, label="(c) 理論スペクトル $\\sigma^2$")
plt.xlabel("正規化周波数 $f$"); plt.ylabel("パワースペクトル密度 $S(f)$")
plt.title("白色雑音：自己相関のフーリエ変換とピリオドグラム平均は一致する")
plt.legend(); plt.grid(alpha=0.3); plt.tight_layout(); plt.show()
print("平均パワー R[0] =", avg_acf[N - 1], " 理論 =", sigma ** 2)
print("スペクトルの平均 =", avg_per.mean(), " 理論 =", sigma ** 2)
左は標本自己相関で、$k=0$ に1.0024が立ち、それ以外のラグはすべて0近傍に貼りついています。$R[k]=\sigma^2\delta[k]$ という理論形そのものです。右のスペクトルでは縦軸を0.85〜1.15に拡大していますが、それでも3本の曲線は互いに区別できないほど重なり、$S(f)=1$ の水平線の周りに±5%程度で収まっています。
このグラフでは3本の線が完全に重なり、$S(f)=1$ の水平線になります。(a)の太い半透明の線の上に(b)の細い線がぴったり乗っているのは偶然ではなく、本文で示した恒等式（標本自己相関のDTFT＝ピリオドグラム）が実現ごとに厳密に成り立っているためです。理論値との一致も良好で、$R[0]$ の推定値は1.00前後、スペクトルの周波数平均も1.00前後になります。白色雑音では $R[k]$ が $k=0$ の一点にしか値を持たないためバートレット窓の影響を受けず、有限 $N$ のバイアスがそもそも生じないのも特徴です。
例2: AR(1)過程での検証
次は相関のある過程です。ここではバイアスの効果もはっきり見えます。
import numpy as np
import matplotlib.pyplot as plt
phi, sigma, N, M = 0.7, 1.0, 256, 2000
rng = np.random.default_rng(11)
acc_per = np.zeros(N)
acc_acf = np.zeros(2 * N - 1)
for _ in range(M):
x = ar1_process(N, phi, sigma, rng)
acc_per += periodogram(x)
acc_acf += biased_acf(x)
avg_per, avg_acf = acc_per / M, acc_acf / M
f = np.fft.fftfreq(N)
lags = np.arange(-(N - 1), N)
acf_spec = np.array([np.sum(avg_acf * np.exp(-2j * np.pi * ff * lags)).real
for ff in f])
S_th = sigma ** 2 / (1 - 2 * phi * np.cos(2 * np.pi * f) + phi ** 2) # (c)
# 有限Nのバイアス（バートレット窓をかけた理論値）
R_th = sigma ** 2 * phi ** np.abs(lags) / (1 - phi ** 2)
w_bartlett = 1 - np.abs(lags) / N
S_biased = np.array([np.sum(w_bartlett * R_th * np.exp(-2j * np.pi * ff * lags)).real
for ff in f])
order = np.argsort(f)
plt.figure(figsize=(10, 4.8))
plt.plot(f[order], acf_spec[order], lw=3.2, alpha=0.40, label="(a) 標本自己相関のFFT")
plt.plot(f[order], avg_per[order], lw=1.2, label=f"(b) ピリオドグラム{M}回平均")
plt.plot(f[order], S_th[order], "k--", lw=1.8, label="(c) 理論スペクトル")
plt.plot(f[order], S_biased[order], ":", lw=1.8, label="(c') バートレット窓つき理論値")
plt.yscale("log")
plt.xlabel("正規化周波数 $f$"); plt.ylabel("パワースペクトル密度 $S(f)$（対数軸）")
plt.title(f"AR(1)過程（$\\phi$={phi}）：3本の曲線が重なる")
plt.legend(); plt.grid(alpha=0.3, which="both"); plt.tight_layout(); plt.show()
print("f=0 推定 %.4f / 理論 %.4f / 窓つき理論 %.4f"
% (avg_per[0], S_th[0], S_biased[0]))
i = np.argmin(np.abs(f - 0.25))
print("f=0.25 推定 %.4f / 理論 %.4f" % (avg_per[i], S_th[i]))
print("f=0.5 推定 %.4f / 理論 %.4f" % (avg_per[N // 2], S_th[N // 2]))
print("R[0] 推定 %.4f / 理論 %.4f" % (avg_acf[N - 1], sigma ** 2 / (1 - phi ** 2)))
print("平均相対誤差（推定 vs 理論） = %.3f" % np.mean(np.abs(avg_per - S_th) / S_th))
左の自己相関は $R[0]$ の実測1.943（理論1.961）から指数的に減衰し、ラグ25ではほぼ0です。この減衰こそが $\sum_k|R[k]|<\infty$ を保証し、フーリエ変換を可能にしています。右のスペクトルは対数軸で、4本の曲線がほぼ完全に重なっています。$f=0$ 付近だけをよく見ると、実測（橙）は真値（黒破線）よりわずかに低く、バートレット窓つき理論値（緑点線）の側に寄っています。
実行すると、$f=0$ で推定10.73に対し理論11.11、$f=0.25$ で推定0.700に対し理論0.671、$f=0.5$ で推定0.337に対し理論0.346となり、全周波数での平均相対誤差は約2%に収まります。低周波が高周波の30倍以上あるという、手計算で予想した低域集中の形がきちんと再現されています。$R[0]$ の推定値1.943も理論値1.961に近い値です。
さらに細かく見ると、残った誤差の正体もはっきりします。$f=0$ での推定値10.73は、真のスペクトル11.11よりも、バートレット窓をかけた理論値10.99のほうに近い。これは偶然ではなく、$E[I_N(f)]$ が真の $S(f)$ ではなくフェイエール核で平滑化されたものに等しいという本文の結果そのものです。鋭いピークは有限データでは必ず低く丸められる — このバイアスは実現数 $M$ をいくら増やしても消えず、データ長 $N$ を伸ばして初めて減ります。
例3: 正弦波＋白色雑音での検証
最後に、密度と線が共存する例を見ます。
import numpy as np
import matplotlib.pyplot as plt
A, f0, sigma = 1.0, 0.15625, 0.5 # f0*N=80 がちょうど整数になるよう選ぶ
N, M = 512, 2000
rng = np.random.default_rng(3)
n = np.arange(N)
acc_per = np.zeros(N)
acc_acf = np.zeros(2 * N - 1)
for _ in range(M):
theta = rng.uniform(0, 2 * np.pi) # ランダム位相（定常化）
x = A * np.cos(2 * np.pi * f0 * n + theta) + sigma * rng.standard_normal(N)
acc_per += periodogram(x)
acc_acf += biased_acf(x)
avg_per, avg_acf = acc_per / M, acc_acf / M
f = np.fft.fftfreq(N)
lags = np.arange(-(N - 1), N)
acf_spec = np.array([np.sum(avg_acf * np.exp(-2j * np.pi * ff * lags)).real
for ff in f])
order = np.argsort(f)
fig, ax = plt.subplots(1, 2, figsize=(13, 4.6))
ax[0].plot(lags, avg_acf, lw=1.0, label="標本自己相関（平均）")
R_th = (A ** 2 / 2) * np.cos(2 * np.pi * f0 * lags) + sigma ** 2 * (lags == 0)
ax[0].plot(lags, (1 - np.abs(lags) / N) * R_th, "k--", lw=1.0, label="理論×バートレット窓")
ax[0].set_xlim(-80, 80); ax[0].set_xlabel("ラグ $k$"); ax[0].set_ylabel("$R[k]$")
ax[0].set_title("自己相関：減衰しない振動＋原点の突起"); ax[0].legend(); ax[0].grid(alpha=0.3)
ax[1].semilogy(f[order], acf_spec[order], lw=3.0, alpha=0.45, label="(a) 標本自己相関のFFT")
ax[1].semilogy(f[order], avg_per[order], lw=1.0, label="(b) ピリオドグラム平均")
ax[1].axhline(sigma ** 2, color="k", ls="--", lw=1.5, label="(c) 雑音の密度 $\\sigma^2$")
ax[1].axhline(N * A ** 2 / 4 + sigma ** 2, color="r", ls=":", lw=1.5,
label="線のピーク予測 $NA^2/4+\\sigma^2$")
ax[1].set_xlabel("正規化周波数 $f$"); ax[1].set_ylabel("$S(f)$（対数軸）")
ax[1].set_title("スペクトル：平坦な密度の上に線が立つ"); ax[1].legend(fontsize=8)
ax[1].grid(alpha=0.3, which="both")
plt.tight_layout(); plt.show()
i0 = np.argmin(np.abs(f - f0))
print("ピーク値 %.3f / 予測 %.3f" % (avg_per[i0], N * A ** 2 / 4 + sigma ** 2))
print("雑音床 %.4f / 予測 %.4f" % (avg_per[np.argmin(np.abs(f - 0.35))], sigma ** 2))
print("全パワー R[0] %.4f / 理論 %.4f" % (avg_acf[N - 1], A ** 2 / 2 + sigma ** 2))
左では標本自己相関（青）が理論×バートレット窓（黒破線）に完全に重なり、ラグ80まで振幅0.4前後の振動が続いています。ラグ0だけが0.75に跳ね上がっているのが雑音の寄与です。右のスペクトルは対数軸で、$\pm f_0$ に幅1ビンの鋭い線が128.25まで立ち、それ以外は $\sigma^2=0.25$ の平坦な床になっています。線と床の高さの比が約500倍もあることが、対数軸でこそはっきり見て取れます。
左のグラフでは、自己相関がラグを大きくしても減衰せずに振動し続ける様子が見えます（三角窓のぶんだけ振幅がゆっくり細っていくのは、有限データの効果です）。ラグ0だけが $A^2/2+\sigma^2=0.75$ まで跳ね上がっているのは、そこにだけ雑音の $\sigma^2$ が乗るからです。$\sum_k|R[k]|$ が発散するというのは、この「いつまでも減衰しない振動」を指しています。
右のグラフでは、$\pm f_0$ に鋭い線が立ち、それ以外は $\sigma^2=0.25$ の平坦な床になります。実測のピーク値は128.248で、予測 $NA^2/4+\sigma^2 = 512/4+0.25 = 128.25$ とほぼ完全に一致しました。雑音床の実測も0.2533で予測0.25とよく合います。全パワー $R[0]$ の実測0.7499も理論0.75と一致します。ここでも(a)と(b)は重なったままで、$R$ が絶対総和可能でなくても有限データの恒等式は破れないことが分かります。
なぜピリオドグラムは滑らかにならないのか
ここまでは「多数回平均したピリオドグラム」を使ってきました。では、平均せずに1回だけ計算したらどうなるでしょうか。$N$ を増やせば滑らかになるでしょうか。
import numpy as np
import matplotlib.pyplot as plt
sigma = 1.0
Ns = [64, 256, 1024, 4096]
rng = np.random.default_rng(7)
# (1) 単一実現のピリオドグラム（N を変えて比較）
fig, ax = plt.subplots(1, 2, figsize=(13, 4.6))
for N in Ns:
x = sigma * rng.standard_normal(N)
P = periodogram(x)
f = np.fft.fftfreq(N)
m = f >= 0
ax[0].plot(f[m], P[m], lw=0.7, alpha=0.8, label=f"N={N}")
ax[0].axhline(sigma ** 2, color="k", ls="--", lw=2, label="真の $S(f)=\\sigma^2$")
ax[0].set_xlabel("正規化周波数 $f$"); ax[0].set_ylabel("$I_N(f)$")
ax[0].set_title("データ長を増やしてもギザギザは消えない")
ax[0].legend(fontsize=8); ax[0].grid(alpha=0.3)
# (2) 特定周波数でのピリオドグラム値の標準偏差
means, stds = [], []
for N in Ns:
rng2 = np.random.default_rng(7)
vals = [periodogram(sigma * rng2.standard_normal(N))[N // 4] for _ in range(5000)]
vals = np.array(vals)
means.append(vals.mean()); stds.append(vals.std())
print(f"N={N:5d} 平均={vals.mean():.3f} 標準偏差={vals.std():.3f}")
ax[1].plot(Ns, means, "o-", label="平均（$\\to S(f)$ に収束）")
ax[1].plot(Ns, stds, "s-", label="標準偏差（減らない）")
ax[1].axhline(sigma ** 2, color="k", ls="--", lw=1.2, label="$\\sigma^2$")
ax[1].set_xscale("log", base=2); ax[1].set_ylim(0, 1.5)
ax[1].set_xlabel("データ長 $N$"); ax[1].set_ylabel("$I_N(0.25)$ の統計量")
ax[1].set_title("ピリオドグラムは一致推定量ではない")
ax[1].legend(fontsize=9); ax[1].grid(alpha=0.3)
plt.tight_layout(); plt.show()
左の4枚は同じ縦軸で並べてあります。$N=64$ でも $N=4096$ でも、値が0から7あたりまで暴れる範囲はまったく変わらず、変わるのは点の密度だけです。右のグラフでは、平均（青）が $N$ によらず1.00前後に張りつく一方、標準偏差（赤）も1.00前後から下がりません。平均だけ収束して分散が残るという、一致推定量でない推定量の典型的な姿です。
$N=64$ でも $N=4096$ でもギザギザの激しさがまったく変わらないことが見て取れます。増えるのは点の密度だけで、真の値 $S(f)=1$ の周りでの暴れ方は同じままです。標準偏差の数値がそれを裏づけます。実測では $N=64$ で平均1.025・標準偏差1.016、$N=4096$ で平均1.000・標準偏差0.991となり、平均は真値に収束するのに標準偏差は1のまま減らないのです。
この振る舞いには理論的な説明がつきます。ガウス性の白色雑音の場合、$I_N(f)$ は $f\ne0,\pm1/2$ において漸近的に $\frac{S(f)}{2}\chi^2_2$ 分布（自由度2のカイ二乗の定数倍、すなわち指数分布）に従います。自由度は $N$ によらず常に2です。データを増やしても推定に使える「独立な情報の数」は増えず、周波数分解能が上がって推定点が増えるだけ。だから分散が減らない。標準偏差が $S(f)$ と同じ大きさになるのも、指数分布の平均と標準偏差が等しいことの帰結です。
このことは、ウィーナー・ヒンチンの定理が期待値についての主張であって、1回の測定についての主張ではないことを改めて示しています。実務でスペクトルを推定するときは、定理をそのまま素朴に使うのではなく、分散を減らす工夫が必須です。代表的には、データを区間に分けて各区間のピリオドグラムを平均するウェルチ法、周波数方向に隣接ビンを平均するダニエル法、標本自己相関にラグ窓をかけてから変換するブラックマン・チューキー法などがあります。いずれも「分解能を犠牲にして自由度を稼ぐ」という同じトレードオフの上に乗っています。
まとめ
本記事では、ウィーナー・ヒンチンの定理を出発点の問題設定から導出まで丁寧にたどりました。
- 問題の所在: 定常確率過程のサンプルパスは絶対可積分でも二乗可積分でもないので、フーリエ変換が直接定義できない。エネルギーは $T\to\infty$ で発散する
- 突破口: 信号そのものではなく、統計量である自己相関 $R(\tau)$ をフーリエ変換する。$R$ は多くの場合減衰するので可積分になる
- PSDの定義: 測定手順をそのまま数式化して $S(f)=\lim_{T\to\infty}\frac1T E[|\hat X_T(f)|^2]$ とする。期待値を先に、極限を後に取る順序が本質的
- 定理: $\int|R|<\infty$ のとき $S(f)=\int R(\tau)e^{-j2\pi f\tau}d\tau$。導出の要は、定常性により二重積分の被積分関数が差だけの関数になり、対角方向のスライスの長さ $T-|\tau|$（バートレット窓）が現れること
- 有限区間の姿: $\frac1T E[|\hat X_T|^2]=\int(1-|\tau|/T)R(\tau)e^{-j2\pi f\tau}d\tau$。離散版では $E[I_N(f)]=\sum(1-|k|/N)R[k]e^{-j2\pi fk}$ となり、フェイエール核による平滑化としても書ける
- 厳密な恒等式: バイアスありの標本自己相関のDTFTは、期待値を取らなくてもピリオドグラムに厳密に一致する（数値実験で誤差 $10^{-13}$ 以下を確認）
- 4つの帰結: $S(f)\ge0$、$S(-f)=S(f)$、$R(0)=\int S(f)df$（全パワー）、$S_Y(f)=|H(f)|^2S_X(f)$
- 正定値性との対応: $R$ の正定値性と $S\ge0$ は同値。この一般形がボホナーの定理（離散版はヘルグロッツの定理）で、ガウス過程のカーネル設計にも直結する
- 3つの例: 白色雑音は $R[k]=\sigma^2\delta[k]\leftrightarrow S(f)=\sigma^2$、AR(1)は $R[k]\propto\phi^{|k|}\leftrightarrow S(f)=\sigma^2/(1-2\phi\cos2\pi f+\phi^2)$、正弦波＋雑音は連続な密度の上に線スペクトルが立つ
- ピリオドグラムの限界: 期待値は真値に収束するが分散は $N$ を増やしても減らない（漸近的に自由度2のカイ二乗）。実用には平均化が必須
この定理を押さえると、スペクトル解析の風景が一段クリアになります。「なぜFFTの二乗がスペクトルなのか」「なぜ窓関数が必要なのか」「なぜ平均しないと使い物にならないのか」が、すべて同じ一本の式から説明できるようになるからです。
次のステップとして、以下の記事も参考にしてください。
- ピリオドグラムによるスペクトル解析 — 本記事の推定量を実践的に使う
- ウェルチ法によるPSD推定 — 分散を減らす標準的な手法
- スペクトル漏れと窓関数の選び方 — フェイエール核の裾が引き起こす問題への対処
- ARモデル（自己回帰モデル） — 例2のAR(1)を一般次数へ拡張する
- ウィーナーフィルタの理論 — PSDを使った最適フィルタ設計
- 短時間フーリエ変換（STFT） — 定常性が成り立たない信号への拡張