エンジンの振動データを録ってみたら、雑音まみれで一見ただのギザギザにしか見えない。でもその奥には、回転に同期した周期成分が必ず潜んでいます。あるいは心電図から心拍の周期を測りたい、レーダーの受信波から自分が送ったパルスの反射を見つけたい――こうした「信号の中に隠れた繰り返し構造を、雑音ごしに掘り当てたい」という場面は、理工学のあらゆる分野に現れます。
このとき絶大な威力を発揮するのが自己相関関数（autocorrelation function）です。自己相関関数は、信号を時間方向に少しずつずらして「自分自身とどれだけ似ているか」を測ります。周期信号は1周期ぶんずらすと元とそっくりに戻るので、自己相関はそのラグ（ずらし量）で大きなピークを作ります。雑音は少しでもずらすと自分自身と似なくなるので、自己相関はラグ0以外でほぼ消えます。この非対称な振る舞いのおかげで、時間波形では埋もれて見えない周期が、自己相関の上でくっきり浮かび上がるのです。
自己相関関数が活きる応用先は数多くありますが、代表的なものを挙げると次の2つです。ひとつは基本周波数（ピッチ）の推定で、音声認識や楽器のチューナー、心拍・呼吸のモニタリングで周期を測るのに使われます。もうひとつはパワースペクトル密度（PSD）の推定で、ウィーナー・ヒンチンの定理を通じて自己相関とスペクトルが表裏一体であることを利用し、ランダム信号の周波数構造を調べます。さらにレーダー・ソナーのマッチトフィルタ、時系列解析のARモデル同定（ユール・ウォーカー方程式）など、応用は分野を横断して広がっています。
本記事の内容
- 自己相関関数の直感（信号を自分自身とずらして掛けて平均する）
- 確定信号と定常確率過程それぞれでの定義
- 偶関数性 $R(\tau)=R(-\tau)$、原点最大 $|R(\tau)|\le R(0)$ の証明
- 白色雑音の自己相関がデルタ関数になること
- ウィーナー・ヒンチンの定理（自己相関のフーリエ変換 = PSD）の導出
- 標本自己相関の推定とバイアス
- AR(1)過程の自己相関の導出
- Pythonでの周期検出とPSD一致の検証
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
自己相関関数とは
自己相関関数の気持ちは、実はとても素朴です。手元に信号 $x(t)$ があるとします。これを時間方向に $\tau$ だけずらしたコピー $x(t+\tau)$ を作り、元の信号と重ね合わせて、両者がどれだけ「似ているか」を数値化する――それだけです。
「似ている度合い」をどう測るかというと、2つの信号を各時刻で掛け算して、その結果を全時間にわたって平均します。掛け算は、両方が同時に正、あるいは両方が同時に負のときに大きな正の値を返し、符号が食い違うと負になります。ですから、ずらしたコピーが元とよく揃っているほど積の平均は大きくなり、バラバラなら打ち消し合ってゼロに近づきます。この「ずらし量 $\tau$ の関数として測った自己相似度」が自己相関関数 $R(\tau)$ です。
上の図がその手続きです。一番上が元の信号 $x(t)$、真ん中が元（青）と $\tau$ だけずらしたコピー（赤）を重ねたもの、一番下が両者の積 $x(t)x(t+\tau)$ です。この積を平均した値が $R(\tau)$ の1点になります。$\tau$ を少しずつ動かしながらこの操作を繰り返すと、$R(\tau)$ という1本の曲線が描けます。周期信号なら、ちょうど1周期ぶんずらしたところで元と再び揃うので、そのラグで $R(\tau)$ が大きな山を作ります。
ここで大事なのは、自己相関は相互相関の特別な場合だという視点です。相互相関 $R_{xy}(\tau)$ が2つの異なる信号 $x$ と $y$ の類似度を測るのに対し、自己相関は $y=x$、つまり自分自身との相互相関にほかなりません。相互相関を先に理解していれば、自己相関はその対角成分だと思えばよいわけです。
では、この「ずらして掛けて平均する」を数式に落とし込みましょう。信号が確定的（毎回同じ波形が得られる）な場合と、確率的（毎回違う実現値が出るランダム信号）な場合とで、平均の取り方が変わります。まず確定信号から定義します。
確定信号の自己相関の定義
信号が明確な数式で書けるような確定信号（deterministic signal）の場合、自己相関は「積の時間平均」として定義します。ただし信号が二乗可積分なエネルギー信号か、周期信号のようなパワー信号かで正規化が少し変わります。
エネルギー信号（$\int |x(t)|^2 dt < \infty$）に対しては、積の全時間積分をそのまま使います。
$$ \begin{equation} R_{xx}(\tau) = \int_{-\infty}^{\infty} x(t)\, x(t+\tau)\, dt \end{equation} $$
一方、正弦波のように永遠に続いてエネルギーが無限大になるパワー信号では、区間 $[-T, T]$ で平均してから $T\to\infty$ の極限を取ります。
$$ \begin{equation} R_{xx}(\tau) = \lim_{T\to\infty} \frac{1}{2T}\int_{-T}^{T} x(t)\, x(t+\tau)\, dt \end{equation} $$
離散時間信号 $x[n]$ の場合は積分が和に置き換わります。長さ $N$ のデータに対しては
$$ \begin{equation} R_{xx}[k] = \frac{1}{N}\sum_{n} x[n]\, x[n+k] \end{equation} $$
のように書けます。ここで $k$ が離散のラグ（サンプル数）です。積分・和のどちらでも、本質は「元の信号とラグだけずらしたコピーの積を、全域で足し合わせる」ことに尽きます。
具体的に、正弦波 $x(t) = A\cos(\omega_0 t)$ の自己相関を計算してみましょう。パワー信号の定義（式2）に代入すると
$$ R_{xx}(\tau) = \lim_{T\to\infty}\frac{1}{2T}\int_{-T}^{T} A\cos(\omega_0 t)\, A\cos(\omega_0 (t+\tau))\, dt $$
となります。ここで積和公式 $\cos\alpha\cos\beta = \tfrac{1}{2}[\cos(\alpha-\beta)+\cos(\alpha+\beta)]$ を使うと、被積分関数は
$$ A^2\cos(\omega_0 t)\cos(\omega_0 t + \omega_0\tau) = \frac{A^2}{2}\big[\cos(\omega_0\tau) + \cos(2\omega_0 t + \omega_0\tau)\big] $$
と2項に分かれます。第1項の $\cos(\omega_0\tau)$ は $t$ に依存しない定数なので、時間平均してもそのまま残ります。第2項 $\cos(2\omega_0 t + \omega_0\tau)$ は $t$ について振動する成分で、長時間平均するとプラスとマイナスが打ち消し合ってゼロに収束します。したがって
$$ \begin{equation} R_{xx}(\tau) = \frac{A^2}{2}\cos(\omega_0\tau) \end{equation} $$
が得られます。この結果は示唆に富んでいます。第一に、自己相関は元の正弦波と同じ周期を持っています（$\cos(\omega_0\tau)$ の周期は $2\pi/\omega_0$）。つまり自己相関を見れば元信号の周期がそのまま読み取れます。第二に、位相の情報が消えています。元の信号がどんな位相 $\phi$ を持っていても $R_{xx}(\tau)=\tfrac{A^2}{2}\cos(\omega_0\tau)$ となり、位相には依存しません。自己相関は「周期構造」だけを抽出し、「いつ始まったか」という位相は捨てるのです。
確定信号の定義が済んだところで、次に扱いたいのは雑音のようなランダム信号です。雑音は毎回違う波形になるので、1本の実現だけを見て時間平均するのではなく、確率過程としての期待値で捉える必要があります。
定常確率過程の自己相関の定義
雑音を含む現実の信号は、測るたびに違う値が出る確率過程 $\{X(t)\}$ です。各時刻 $t$ における値 $X(t)$ が確率変数であり、無数の実現（サンプルパス）の集まりとして信号を捉えます。このとき自己相関は、時間平均ではなく集団平均（アンサンブル平均）、すなわち期待値で定義します。
$$ \begin{equation} R_{XX}(t_1, t_2) = \mathbb{E}\big[X(t_1)\, X(t_2)\big] \end{equation} $$
一般にはこれは2つの時刻 $t_1, t_2$ の両方に依存する2変数関数です。しかしここで、多くの実用信号が満たす弱定常性（wide-sense stationary, WSS）を仮定します。弱定常とは、(1) 平均 $\mathbb{E}[X(t)]=\mu$ が時刻によらず一定で、(2) 自己相関が2つの時刻の差 $\tau = t_2 – t_1$ にしか依存しない、という性質です。定常性を仮定すると自己相関は1変数関数に簡約され、
$$ \begin{equation} R_{XX}(\tau) = \mathbb{E}\big[X(t)\, X(t+\tau)\big] \end{equation} $$
と書けます。「いつ測ったか」ではなく「どれだけ時間差があるか」だけで相関が決まる、というのが定常過程の便利さです。以降は平均を差し引いた $\mu=0$ の中心化された過程を考えることが多く、その場合 $R_{XX}(\tau)$ は自己共分散とも一致します。
確率過程では、時間平均と集団平均が一致するかどうかは自明ではありません。1本のサンプルパスを長時間平均した結果が、無数の実現にわたる期待値と一致する性質をエルゴード性（ergodicity）と呼びます。実務では真の期待値は手に入らないので、1本の観測データから時間平均で自己相関を推定します。この推定が意味を持つのは、対象がエルゴード的であるという暗黙の前提のおかげです。
$$ R_{XX}(\tau) \approx \lim_{T\to\infty}\frac{1}{2T}\int_{-T}^{T} x(t)\, x(t+\tau)\, dt $$
つまり、確定信号のところで書いた時間平均の式（式2）は、エルゴード的な定常過程では期待値（式6）の良い近似になっている、というわけです。この橋渡しがあるからこそ、1回の測定データから自己相関を計算する日々の作業が正当化されます。
定義が2通り（確定信号の時間平均、確率過程の期待値）出そろいました。どちらの定義でも、自己相関関数は同じ美しい性質を共有しています。次にその性質を証明していきましょう。
自己相関関数の性質
自己相関関数 $R(\tau)$ には、どんな信号に対しても必ず成り立つ普遍的な性質があります。特に重要なのが「偶関数であること」と「原点で最大値を取ること」の2つで、これらは自己相関の形を大きく制約します。
上の図はAR過程の標本自己相関ですが、ラグ $\tau=0$ を中心に左右対称（偶関数）で、しかも $\tau=0$ で最大値を取っていることがひと目でわかります。すべての値が $\pm R(0)$ の帯の中に収まっている点も見て取れます。以下ではこれらを数式で証明します。
性質1：偶関数性 $R(\tau) = R(-\tau)$
定常過程の定義（式6）から出発します。証明の鍵は、期待値が時刻の絶対位置によらない（定常性）ことです。
$$ R_{XX}(-\tau) = \mathbb{E}[X(t)\, X(t-\tau)] $$
ここで時刻を $t \to t+\tau$ とずらしても、定常性により期待値は変わりません。つまり $s = t+\tau$ とおくと
$$ R_{XX}(-\tau) = \mathbb{E}[X(s)\, X(s-\tau)]\Big|_{s=t+\tau} = \mathbb{E}[X(t+\tau)\, X(t)] $$
積は可換なので $X(t+\tau)X(t) = X(t)X(t+\tau)$ です。したがって
$$ \begin{equation} R_{XX}(-\tau) = \mathbb{E}[X(t)\, X(t+\tau)] = R_{XX}(\tau) \end{equation} $$
となり、偶関数性が示せました。直感的にも、「$\tau$ だけ先にずらして比べる」のと「$\tau$ だけ後ろにずらして比べる」のは、同じペアを比べているので当然同じ値になります。この偶関数性のおかげで、自己相関は非負のラグ $\tau\ge 0$ だけ計算すれば十分です。
性質2：原点最大 $|R(\tau)| \le R(0)$
次に、自己相関がラグ0で最大になることを示します。まず $R(0) = \mathbb{E}[X(t)^2]$ は信号の平均パワー（中心化していれば分散）そのものです。証明にはコーシー・シュワルツの不等式を使います。
任意の2つの確率変数 $U, V$ に対して、コーシー・シュワルツの不等式は
$$ \big(\mathbb{E}[UV]\big)^2 \le \mathbb{E}[U^2]\,\mathbb{E}[V^2] $$
を主張します。ここで $U = X(t)$、$V = X(t+\tau)$ と選びます。すると左辺は $\big(\mathbb{E}[X(t)X(t+\tau)]\big)^2 = R(\tau)^2$ です。右辺は $\mathbb{E}[X(t)^2]\,\mathbb{E}[X(t+\tau)^2]$ ですが、定常性により両方の因子が等しく $\mathbb{E}[X(t)^2] = \mathbb{E}[X(t+\tau)^2] = R(0)$ です。したがって
$$ R(\tau)^2 \le R(0)\cdot R(0) = R(0)^2 $$
となり、両辺の平方根を取れば
$$ \begin{equation} |R(\tau)| \le R(0) \end{equation} $$
が得られます。等号は $\tau=0$ で成立します（自分自身との相関だから完全一致）。この不等式は「どんな時間差でも、信号が自分自身と完全一致することを超える相似はありえない」という当たり前の事実を厳密に述べたものです。図の点線帯 $\pm R(0)$ からはみ出す値が存在しない理由がこれです。
性質3：$R(0)$ はパワー、ラグ増大で減衰
最後に振る舞いの傾向をまとめておきます。$R(0)$ は平均パワー（分散）で常に非負です。多くの現実の信号では、ラグ $\tau$ が大きくなるほど「昔の値との関係」は薄れていくので、$R(\tau) \to 0$（$\tau\to\infty$）と減衰します。減衰の速さがその信号の「記憶の長さ」を表します。ゆっくり減衰する信号は長い相関（強い予測可能性）を持ち、急に減衰する信号は短い記憶しか持ちません。
この「記憶がない極限」を突き詰めると、白色雑音という理想的な信号にたどり着きます。次はその自己相関を調べましょう。
白色雑音の自己相関はデルタ関数
白色雑音（white noise）は、各時刻の値が互いに完全に無相関で、平均0・分散 $\sigma^2$ を持つ理想的なランダム信号です。「白色」という名は、後で見るように全周波数に均等にパワーが分布する（白色光のアナロジー）ことに由来します。
白色雑音の自己相関を考えます。定義は $R(\tau) = \mathbb{E}[X(t)X(t+\tau)]$ でした。$\tau \ne 0$ のとき、$X(t)$ と $X(t+\tau)$ は異なる時刻の値であり、白色性の定義により無相関（平均0なので $\mathbb{E}[X(t)X(t+\tau)] = 0$）です。$\tau=0$ のときは $\mathbb{E}[X(t)^2] = \sigma^2$（分散）です。まとめると、離散時間では
$$ \begin{equation} R_{XX}[k] = \sigma^2\,\delta[k] = \begin{cases}\sigma^2 & k = 0 \\ 0 & k \ne 0\end{cases} \end{equation} $$
連続時間では
$$ \begin{equation} R_{XX}(\tau) = \sigma^2\,\delta(\tau) \end{equation} $$
となります。ここで $\delta$ はクロネッカーのデルタ（離散）またはディラックのデルタ（連続）です。白色雑音の自己相関は原点にだけ立つ1本の針、というわけです。
上の図は分散1の白色雑音の波形（左）と、その標本自己相関（右）です。$k=0$ で $R[0]\approx 0.96$（真の分散1にほぼ一致）と鋭く立ち上がり、$k\ge 1$ ではほぼ消えている（有限データの揺らぎで最大 $0.04$ 程度）ことがわかります。まさにデルタ関数状です。この「原点だけ」という性質は、白色雑音が過去といっさい相関を持たない=完全に予測不能であることの表れです。
デルタ関数状の自己相関は、後で見るウィーナー・ヒンチンの定理を通じて「平坦なスペクトル」に対応します。デルタ関数のフーリエ変換は定数だからです。この事実は、自己相関とスペクトルが表裏一体であることの最も単純で美しい例になっています。では、その表裏一体の関係――ウィーナー・ヒンチンの定理を導きましょう。
ウィーナー・ヒンチンの定理
自己相関関数の理論的なハイライトがウィーナー・ヒンチンの定理（Wiener–Khinchin theorem）です。定理は次のように述べます。
定常確率過程の自己相関関数のフーリエ変換は、パワースペクトル密度（PSD）に等しい。
$$ \begin{equation} S_{XX}(\omega) = \int_{-\infty}^{\infty} R_{XX}(\tau)\, e^{-j\omega\tau}\, d\tau \end{equation} $$
逆に、PSDを逆フーリエ変換すれば自己相関に戻ります。
$$ R_{XX}(\tau) = \frac{1}{2\pi}\int_{-\infty}^{\infty} S_{XX}(\omega)\, e^{j\omega\tau}\, d\omega $$
これは驚くべき主張です。時間領域で「信号がどれだけ自分と似ているか」を測った $R(\tau)$ と、周波数領域で「各周波数にどれだけパワーがあるか」を表す $S(\omega)$ が、フーリエ変換で1対1に結ばれている――時間の自己相似性と周波数のパワー分布は同じ情報の別表現だ、というのです。
なぜ成り立つのか（導出）
ランダム信号のフーリエ変換は数学的に微妙（サンプルパスが二乗可積分でない）ですが、ここでは有限区間 $[-T,T]$ で切り出した信号 $x_T(t)$ を使う標準的な導出をたどります。PSD は「有限区間のフーリエ変換の絶対値2乗を区間長で割り、期待値を取って区間長を無限大にした極限」として定義されます。
$$ S_{XX}(\omega) = \lim_{T\to\infty}\frac{1}{2T}\,\mathbb{E}\big[\,|X_T(\omega)|^2\,\big], \qquad X_T(\omega) = \int_{-T}^{T} x(t)\, e^{-j\omega t}\, dt $$
まず $|X_T(\omega)|^2 = X_T(\omega)\,X_T^*(\omega)$ を二重積分に書き下します。$X_T^*(\omega)$ の積分変数を $t’$ とすると
$$ |X_T(\omega)|^2 = \int_{-T}^{T}\!\int_{-T}^{T} x(t)\, x(t’)\, e^{-j\omega t}\, e^{+j\omega t’}\, dt\, dt’ $$
となります。ここで期待値を取ると、期待値と積分の順序を交換して $\mathbb{E}[x(t)x(t’)] = R_{XX}(t – t’)$（定常性より時刻差だけの関数）が現れます。
$$ \mathbb{E}\big[|X_T(\omega)|^2\big] = \int_{-T}^{T}\!\int_{-T}^{T} R_{XX}(t-t’)\, e^{-j\omega(t – t’)}\, dt\, dt’ $$
指数関数も $e^{-j\omega t}e^{+j\omega t’} = e^{-j\omega(t-t’)}$ とまとめました。次に変数変換 $\tau = t – t’$ を導入します。$t’$ を固定して $t$ で積分すると考え、内側の積分が $\tau$ の関数になるように書き換えると、$2T$ で割ったときに区間長で重み付けされた形が出てきます。$T\to\infty$ の極限では、この二重積分は1次元積分に collapse し、区間の端の寄与（三角窓）が消えて
$$ \lim_{T\to\infty}\frac{1}{2T}\,\mathbb{E}\big[|X_T(\omega)|^2\big] = \int_{-\infty}^{\infty} R_{XX}(\tau)\, e^{-j\omega\tau}\, d\tau $$
が残ります。これがまさに $S_{XX}(\omega)$ の定義（式12の左辺）と一致します。以上でウィーナー・ヒンチンの定理が導けました。要点は、「フーリエ変換の絶対値2乗の期待値」を展開すると必然的に $R_{XX}(t-t’)$ が湧き出し、それを1次元に畳み込むと自己相関のフーリエ変換になる、という構造です。
白色雑音で確かめる
定理を白色雑音でチェックしてみましょう。白色雑音の自己相関は $R(\tau)=\sigma^2\delta(\tau)$ でした（式11）。これをフーリエ変換すると、デルタ関数の変換は定数なので
$$ S_{XX}(\omega) = \int_{-\infty}^{\infty} \sigma^2\delta(\tau)\, e^{-j\omega\tau}\, d\tau = \sigma^2 $$
となり、全周波数で一定のPSDが得られます。これが「白色」の名の由来です。デルタ関数状の自己相関 ⟷ 平坦なスペクトル、という対応が定理からきれいに出てきました。
理論が出そろいました。ここからは実データ（有限長の観測）で自己相関をどう推定するかという実務の話に移り、最後にPythonで周期検出とPSD一致を確かめます。
標本自己相関の推定
実際に手に入るのは無限長の信号ではなく、有限個のサンプル $x[0], x[1], \dots, x[N-1]$ です。この有限データから自己相関を推定する式（標本自己相関）には、正規化の仕方で2つの流儀があります。
$$ \hat{R}_{\text{biased}}[k] = \frac{1}{N}\sum_{n=0}^{N-1-k} x[n]\, x[n+k], \qquad \hat{R}_{\text{unbiased}}[k] = \frac{1}{N-k}\sum_{n=0}^{N-1-k} x[n]\, x[n+k] $$
両者の違いは分母だけです。ラグ $k$ の項は、実際には $N-k$ 個の積しか足せません（データの端で相手がなくなるため）。不偏推定は正しく $N-k$ で割るので期待値が真値に一致しますが、$k$ が大きいと少数の積を割り増しするため分散が暴れます。バイアス付き推定は常に $N$ で割るので、大きなラグでは実際より小さめに見積もる（$\tau\to\infty$ で必ず0へ引き寄せられる）代わりに、分散が小さく安定します。加えて $1/N$ 正規化は、標本自己相関の系列が数学的に「半正定値」になる（有効なPSDを生む）という利点があり、実務ではバイアス付きが好まれます。本記事のPython実装も $1/N$ 正規化を使います。
なお numpy の np.correlate(x, x, mode="full") は正規化なしの生の相関和を返すので、自分で $1/N$ で割る必要があります。この関数は長さ $2N-1$ の配列を返し、その中央がラグ0、右半分が非負ラグ $k=0,1,\dots,N-1$ に対応します。偶関数性から非負ラグだけ見れば十分でした。
推定式が固まったので、いよいよ具体的な確率過程――AR(1)過程――の自己相関を理論と実測の両面から見ていきます。
具体例：AR(1)過程の自己相関
1次自己回帰過程（AR(1)）は、時系列解析で最も基本的な有色雑音モデルで、次式で定義されます。
$$ X[n] = \phi\, X[n-1] + \varepsilon[n], \qquad |\phi| < 1 $$
ここで $\varepsilon[n]$ は分散 $\sigma_\varepsilon^2$ の白色雑音（イノベーション）、$\phi$ は前の値をどれだけ引き継ぐかを表す係数です。$|\phi|<1$ は過程が定常になる条件です。白色雑音がデルタ状の自己相関を持つのに対し、AR(1) は過去を引きずるので、自己相関がなだらかに減衰する「有色」の信号になります。
この過程の自己相関 $R[k]=\mathbb{E}[X[n]X[n+k]]$ を導きましょう。定義式の両辺に $X[n-k]$（$k\ge 1$）を掛けて期待値を取ります。
$$ \mathbb{E}[X[n]X[n-k]] = \phi\,\mathbb{E}[X[n-1]X[n-k]] + \mathbb{E}[\varepsilon[n]X[n-k]] $$
ここで左辺は $R[k]$、右辺第1項は $\phi R[k-1]$ です。右辺第2項は、$\varepsilon[n]$ が「時刻 $n$ より前の $X[n-k]$（$k\ge1$）」と無相関なので0になります（イノベーションは過去に影響されないため）。したがって
$$ R[k] = \phi\, R[k-1], \qquad k \ge 1 $$
という漸化式が得られます。これを繰り返し適用すると $R[k] = \phi^k R[0]$（$k\ge0$）、偶関数性も合わせると一般に
$$ \begin{equation} R[k] = R[0]\,\phi^{|k|}, \qquad R[0] = \frac{\sigma_\varepsilon^2}{1-\phi^2} \end{equation} $$
となります。$R[0]$ は $k=0$ での式（分散）を、定常性 $\mathrm{Var}(X)=\phi^2\mathrm{Var}(X)+\sigma_\varepsilon^2$ を解いて求めたものです。要するにAR(1)の自己相関は公比 $\phi$ の幾何級数として減衰します。$\phi$ が1に近いほどゆっくり減り（長い記憶）、0に近いほど急に減ります（白色雑音に近い）。
上の図は $\phi=0.8$ のAR(1)過程の波形（左）と、その正規化自己相関 $\rho[k]=R[k]/R[0]$（右）です。標本値（青の棒）が理論曲線 $\phi^k$（赤の破線）にぴたりと重なっています。実測では $\rho[1]=0.804$（理論 $0.800$）、$\rho[3]=0.529$（理論 $0.512$）と、幾何減衰の理論とよく一致しました。波形が「なめらかに上下する（隣の値と似ている）」ことと、自己相関がゆっくり減衰することが対応しているのが見て取れます。
理論と具体例が揃いました。最後に、冒頭で掲げた「雑音に埋もれた周期の検出」と「自己相関のフーリエ変換がPSDになる」ことを、Pythonで実際に確かめましょう。
Pythonでの実装
周期検出：雑音に埋もれた正弦波
まず、10 Hz の正弦波に大きな白色雑音（振幅で正弦波の1.5倍）を加えた信号を作り、時間波形と自己相関を比べます。時間波形では周期が見えなくても、自己相関にはピークが出ることを確認します。
import numpy as np
import matplotlib.pyplot as plt
def acf_biased(x):
"""バイアス付き標本自己相関 R[k] (k>=0) を返す。1/N 正規化。"""
x = x - x.mean()
N = len(x)
full = np.correlate(x, x, mode="full") # 長さ 2N-1
return full[N - 1:] / N # 非負ラグ k=0..N-1
np.random.seed(42)
fs, f0, T = 200.0, 10.0, 5.0 # サンプリング200Hz, 正弦波10Hz, 5秒
N = int(fs * T)
t = np.arange(N) / fs
x = np.sin(2*np.pi*f0*t) + 1.5*np.random.randn(N) # 正弦波 + 強い雑音
r = acf_biased(x)
lag_ms = np.arange(N) / fs * 1000
period_ms = 1.0 / f0 * 1000 # 周期 = 100 ms
fig, ax = plt.subplots(2, 1, figsize=(9, 6))
ax[0].plot(t[:400]*1000, x[:400], lw=0.8)
ax[0].set_xlabel("時間 [ms]"); ax[0].set_ylabel("振幅")
ax[0].set_title("雑音に埋もれた10Hz正弦波")
ax[1].plot(lag_ms[:500], r[:500], color="red")
ax[1].set_xlabel("ラグ τ [ms]"); ax[1].set_ylabel("自己相関 R(τ)")
ax[1].set_title("自己相関に100ms周期のピーク")
plt.tight_layout(); plt.show()
このコードから、時間波形（上）は雑音でギザギザして周期がほとんど見えないのに対し、自己相関（下）には $\tau=100\,\text{ms}, 200\,\text{ms}, \dots$ と周期100 ms（=10 Hz）ごとに規則的なピークが並ぶことがわかります。実測ではラグ0で $R[0]\approx 2.62$（信号パワー0.5＋雑音分散2.25の和に近い）となり、第1ピークは約90〜100 msに現れました（雑音の揺らぎで真値100 msから1サンプルぶんずれる程度）。周期のピーク位置から基本周波数 $f_0 = 1/\tau_{\text{peak}}$ を読み取れる、というのがピッチ検出の原理です。
この図が示す通り、自己相関は雑音（ラグ0にしか相関を持たない）と周期信号（周期ごとにピーク）を自然に分離します。雑音のパワーはラグ0に集中して周期ピークを邪魔せず、周期成分だけが遠くのラグまで生き残るのです。これが「雑音ごしに周期を掘り当てる」自己相関の威力です。
ウィーナー・ヒンチンの定理の数値検証
次に、自己相関のフーリエ変換がPSD（ここでは信号のフーリエ変換の絶対値2乗＝ペリオドグラム）に一致することを確かめます。20 Hz と 50 Hz の正弦波に雑音を加えた信号を使います。
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(7)
fs, N = 256.0, 4096
t = np.arange(N) / fs
sig = np.sin(2*np.pi*20*t) + 0.7*np.sin(2*np.pi*50*t) + 0.8*np.random.randn(N)
sig = sig - sig.mean()
# 方法A: ペリオドグラム(信号のFFTの絶対値2乗)
freqs = np.fft.rfftfreq(N, 1/fs)
psd_direct = np.abs(np.fft.rfft(sig))**2 / N
# 方法B: ウィーナー・ヒンチン(自己相関のフーリエ変換)
r_full = np.correlate(sig, sig, mode="full") / N # 両側ACF, 長さ2N-1
Rk = np.zeros(N)
Rk[0] = r_full[N-1]
Rk[1:] = r_full[N:] + r_full[N-2::-1] # 正負ラグを畳む(実対称化)
psd_wh = np.fft.rfft(Rk).real
rel = np.max(np.abs(psd_wh[1:] - psd_direct[1:])) / np.max(psd_direct)
print(f"最大相対差 = {rel:.2e}") # -> 2.75e-16
plt.figure(figsize=(9, 4.5))
plt.plot(freqs, psd_direct, label="ペリオドグラム")
plt.plot(freqs, psd_wh, "r--", label="自己相関のフーリエ変換")
plt.xlim(0, 100); plt.xlabel("周波数 [Hz]"); plt.ylabel("PSD")
plt.legend(); plt.tight_layout(); plt.show()
このコードの出力では、2通りの方法で計算したPSDの最大相対差が $2.75\times 10^{-16}$、すなわち浮動小数点の丸め誤差レベルで完全に一致しました。ウィーナー・ヒンチンの定理が数値的にも厳密に成り立っていることの確認です。
図では青の実線（ペリオドグラム）と赤の破線（自己相関のフーリエ変換）が完全に重なり、20 Hz と 50 Hz に明瞭なピークが立っています。時間領域で計算した自己相関と、周波数領域のパワー分布が同じ情報を持つことが視覚的にも確認できました。この一致こそがウィーナー・ヒンチンの定理の内容であり、PSD推定に自己相関を経由する道（相関法／BT法）が使える根拠になっています。
まとめ
本記事では、自己相関関数について、直感から定義、性質の証明、ウィーナー・ヒンチンの定理、そしてPythonでの検証までを解説しました。
- 直感：自己相関 $R(\tau)$ は、信号を $\tau$ ずらして自分自身と掛けて平均した「自己相似度」。相互相関の $y=x$ の特別な場合
- 定義：確定信号では時間平均、定常確率過程では期待値 $R(\tau)=\mathbb{E}[X(t)X(t+\tau)]$。エルゴード性が両者を橋渡しする
- 性質：偶関数 $R(\tau)=R(-\tau)$、原点最大 $|R(\tau)|\le R(0)$（コーシー・シュワルツ）、$R(0)$ は平均パワー
- 白色雑音：$R(\tau)=\sigma^2\delta(\tau)$ とデルタ状。過去と無相関＝完全に予測不能
- ウィーナー・ヒンチンの定理：自己相関のフーリエ変換がPSD。時間の自己相似性と周波数のパワー分布は表裏一体
- 応用：正弦波の周期はピーク位置に現れ、雑音ごしの周期検出（ピッチ推定）に使える。AR(1)過程の自己相関は $\phi^{|k|}$ と幾何減衰する
自己相関は、時系列のARモデル同定（ユール・ウォーカー方程式）、PSD推定、レーダー・ソナーのマッチトフィルタ、そして機械学習の時系列特徴量まで、幅広い技術の土台になっています。次のステップとして、以下の記事も参考にしてください。