マイクで拾った音声に、別のマイクで録った同じ音声が混ざっているとします。2つの録音を重ね合わせて比較すれば「似ている」ことは直感的にわかりますが、どれくらい似ているか、そしてどれだけ時間がずれているかを数値として取り出すにはどうすればよいでしょうか。
GPS受信機は、衛星から届く微弱な信号の中から既知のコード列を探し出し、到達時刻の差から自分の位置を算出します。レーダーは送信パルスの反射波が何マイクロ秒遅れて戻ってきたかを計測し、目標までの距離を求めます。これらに共通するのは「2つの信号のずれ（遅延）を正確に測りたい」という要求です。
この問いに答えるのが相互相関関数（cross-correlation function）です。相互相関関数を理解すると、次のような場面で強力な武器になります。
- 遅延推定: GPS・レーダー・超音波センサなど、信号の到達時間差から距離や位置を求める
- テンプレートマッチング: 長い信号の中から既知のパターンを検出する（音声認識、画像処理）
- システム同定: 入力と出力の相互相関からインパルス応答を推定する
- ノイズ低減: 相関のある成分だけを取り出し、無相関なノイズを除去する
本記事の内容
- 相関の直感的理解 — 2つの信号がどれだけ「似ているか」の尺度
- 自己相関関数の定義と性質
- 相互相関関数の定義と性質
- フーリエ変換との関係（相互スペクトル密度）
- 正規化相関とピアソン相関係数の関係
- 遅延推定への応用（GPS・レーダーの原理）
- テンプレートマッチングの仕組み
- 離散相関とFFTによる高速計算
- Python実装：遅延推定とテンプレートマッチング
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- 畳み込み定理の理論と証明 — 畳み込みとフーリエ変換の関係
- フーリエ変換の理論 — 連続フーリエ変換の定義と性質
相関の直感 — 「ずらして掛けて足す」
相関（correlation）の本質は、驚くほどシンプルです。2つの信号を少しずつずらしながら重ね合わせ、各点で値を掛け算してから全部足し合わせる — それだけです。
もう少し具体的にイメージしましょう。信号 $x(t)$ のグラフが手元にあるとして、もう1つの信号 $y(t)$ のグラフを透明シートに印刷します。透明シートを $x(t)$ のグラフの上に重ね、左右にスライドさせます。2つの波形がぴったり一致する位置では、山と山、谷と谷が重なるので、掛け算の結果は常に正の値になり、合計は大きくなります。逆に、山と谷がずれて重なると掛け算は負になり、合計は小さく（または負に）なります。
この「ずらし量」を $\tau$（タウ）と書き、合計値を $\tau$ の関数として描いたものが相関関数です。相関関数のピーク位置が、2つの信号が最もよく一致するずらし量、すなわち時間遅延を教えてくれるのです。
この操作は畳み込み（convolution）と非常に似ています。畳み込みでは一方の信号を時間反転してからスライドさせますが、相関では時間反転せずにそのままスライドさせます。この違いが数学的にどう表れるかは、後のセクションで詳しく見ていきます。
まずは、自分自身との相関、すなわち「自己相関関数」から定義を確認しましょう。
自己相関関数の定義と性質
自己相関関数とは
相互相関関数を理解する前に、まず自分自身との相関を考えると見通しがよくなります。自己相関関数（autocorrelation function）は、「信号が自分自身のずらしたコピーとどれだけ似ているか」を遅延 $\tau$ の関数として表したものです。
たとえば、正弦波は周期的に繰り返すので、1周期分ずらすと元の波形に完全に重なります。つまり、自己相関関数は周期 $T$ ごとにピークを持つはずです。一方、ホワイトノイズはどんなにずらしても自分自身と似た形にはなりません。自己相関は $\tau = 0$ でのみピークを持ち、それ以外ではほぼゼロになります。このように、自己相関関数は信号の「時間構造」を特徴づける道具です。
エネルギー有限信号（$\int_{-\infty}^{\infty} |x(t)|^2 dt < \infty$）に対して、自己相関関数 $R_{xx}(\tau)$ は次のように定義されます。
$$ \begin{equation} R_{xx}(\tau) = \int_{-\infty}^{\infty} x(t) \, x(t + \tau) \, dt \end{equation} $$
$\tau = 0$ を代入すると $R_{xx}(0) = \int |x(t)|^2 dt$ となり、信号の全エネルギーに一致します。
パワー有限信号（周期信号や確率過程など、エネルギーが無限だが平均パワーが有限な信号）に対しては、時間平均を取った形で定義します。
$$ \begin{equation} R_{xx}(\tau) = \lim_{T \to \infty} \frac{1}{2T} \int_{-T}^{T} x(t) \, x(t + \tau) \, dt \end{equation} $$
自己相関関数の主な性質
自己相関関数は、以下の重要な性質を持ちます。
性質1: 原点で最大値を取る
$$ |R_{xx}(\tau)| \leq R_{xx}(0) \quad \text{（すべての } \tau \text{ に対して）} $$
これは直感に合います。信号を全くずらさないとき（$\tau = 0$）が最も自分自身と似ているのは当然です。コーシー・シュワルツの不等式から厳密に証明できます。
$$ \left| \int x(t) x(t+\tau) dt \right|^2 \leq \int |x(t)|^2 dt \cdot \int |x(t+\tau)|^2 dt = R_{xx}(0)^2 $$
性質2: 偶関数
$$ R_{xx}(\tau) = R_{xx}(-\tau) $$
$\tau$ だけ右にずらして重ねることと、$\tau$ だけ左にずらして重ねることは同じ一致度を与えます。証明は、積分の変数を $u = t + \tau$ と置換すれば直ちに得られます。
性質3: フーリエ変換はパワースペクトル密度
自己相関関数のフーリエ変換は、信号のパワースペクトル密度 $S_{xx}(f)$ に等しくなります。これはウィーナー・ヒンチンの定理として知られ、時間領域と周波数領域を結ぶ基本的な関係式です。
$$ \begin{equation} S_{xx}(f) = \int_{-\infty}^{\infty} R_{xx}(\tau) \, e^{-j2\pi f \tau} \, d\tau \end{equation} $$
パワースペクトル密度は $S_{xx}(f) = |X(f)|^2 \geq 0$ であるため、常に非負です。この非負性は自己相関関数の原点最大性と対応しています。
自己相関関数が信号の「自分自身との類似度」を測るのに対し、次に定義する相互相関関数は「2つの異なる信号の類似度」を測る道具です。
相互相関関数の定義と性質
定義
いよいよ本題の相互相関関数です。2つのマイクでそれぞれ録音した信号 $x(t)$ と $y(t)$ があるとします。音源が片方のマイクに近い場合、一方の信号は他方より少し遅れて到着します。この遅延を知りたいとき、$y(t)$ を $\tau$ だけずらしながら $x(t)$ との掛け算の合計を計算すれば、最もよく一致する $\tau$ が遅延量を教えてくれます。
エネルギー有限信号に対する相互相関関数 $R_{xy}(\tau)$ は、次のように定義されます。
$$ \begin{equation} R_{xy}(\tau) = \int_{-\infty}^{\infty} x(t) \, y(t + \tau) \, dt \end{equation} $$
この式が言っていることを日本語に翻訳すると、「$y$ を $\tau$ だけ時間シフトし、$x$ と掛け合わせて全時刻にわたって積算する」です。
パワー有限信号の場合は、自己相関関数と同様に時間平均を取ります。
$$ \begin{equation} R_{xy}(\tau) = \lim_{T \to \infty} \frac{1}{2T} \int_{-T}^{T} x(t) \, y(t + \tau) \, dt \end{equation} $$
注意: 教科書によって $x(t) \, y(t + \tau)$ と $x(t + \tau) \, y(t)$ のどちらを採用するか、あるいは $y$ に複素共役を付けるかが異なります。本記事では実数信号を前提とし、$x(t) \, y(t + \tau)$ を採用します。複素信号の場合は $x^*(t) \, y(t + \tau)$ とするのが一般的です。
畳み込みとの関係
相互相関関数は、畳み込みと密接に関係しています。畳み込みの定義は
$$ (x * y)(\tau) = \int_{-\infty}^{\infty} x(t) \, y(\tau – t) \, dt $$
です。相互相関関数の定義式で $u = t + \tau$ と置換すると
$$ R_{xy}(\tau) = \int_{-\infty}^{\infty} x(u – \tau) \, y(u) \, du $$
と書き直せます。これを畳み込みの形と比較すると、相互相関は $x$ の時間反転 $x(-t)$ との畳み込みに他なりません。
$$ \begin{equation} R_{xy}(\tau) = [x(-t) * y(t)](\tau) \end{equation} $$
畳み込みでは一方の信号を「裏返してスライドさせる」のに対し、相関では「そのままスライドさせる」という直感的な違いが、数式上では $x(t)$ を $x(-t)$ に置き換えるという形で現れているのです。
相互相関関数の性質
相互相関関数は自己相関関数とは異なるいくつかの特徴を持ちます。
性質1: 偶関数ではない（一般に非対称）
$$ R_{xy}(\tau) = R_{yx}(-\tau) $$
$x$ と $y$ を入れ替えると、遅延の符号が反転します。自己相関は偶関数でしたが、相互相関はそうとは限りません。これは物理的に自然です。「$x$ に対して $y$ が $\tau_0$ だけ遅れている」ことと「$y$ に対して $x$ が $\tau_0$ だけ遅れている」ことは、遅延の符号が反転する関係にあります。
証明は定義に従って確認できます。
$$ R_{yx}(-\tau) = \int y(t) \, x(t – \tau) \, dt $$
ここで $u = t – \tau$ と置換すると $t = u + \tau$、$dt = du$ となるので
$$ R_{yx}(-\tau) = \int y(u + \tau) \, x(u) \, du = R_{xy}(\tau) $$
が得られます。
性質2: 上界の存在
$$ |R_{xy}(\tau)| \leq \sqrt{R_{xx}(0) \cdot R_{yy}(0)} $$
コーシー・シュワルツの不等式から導かれます。右辺は $x$ と $y$ のエネルギーの幾何平均です。$R_{xy}(\tau)$ の絶対値がこの値に近いほど、2つの信号の類似度が高いことを意味します。
性質3: 無相関な信号の相互相関はゼロ
$x(t)$ と $y(t)$ が統計的に無相関であれば、$R_{xy}(\tau) = 0$ です。ノイズと信号の相互相関を計算すると、信号成分だけが残りノイズが消える — これが相互相関を使ったノイズ低減の原理です。
ここまでで、相互相関関数の時間領域での定義と性質を把握しました。次に、周波数領域の視点からこの関数を眺めてみましょう。フーリエ変換を使うと、相互相関の計算が格段に効率的になるだけでなく、周波数ごとの相関構造という新しい知見も得られます。
フーリエ変換との関係 — 相互スペクトル密度
相互相関定理
畳み込み定理を知っていれば、相互相関のフーリエ変換は自然に導けます。前のセクションで示したように、相互相関は $x(-t)$ と $y(t)$ の畳み込みです。畳み込み定理によれば、畳み込みのフーリエ変換は各信号のフーリエ変換の積になります。
$x(t)$ のフーリエ変換を $X(f)$、$y(t)$ のフーリエ変換を $Y(f)$ とします。時間反転のフーリエ変換に関して、$x(-t) \xleftrightarrow{\mathcal{F}} X^*(f)$（実数信号の場合）という性質を使います。
畳み込み定理を適用すると
$$ \mathcal{F}\{R_{xy}(\tau)\} = \mathcal{F}\{x(-t) * y(t)\} $$
右辺にフーリエ変換の乗算則を適用して
$$ = X^*(f) \cdot Y(f) $$
が得られます。この量を相互スペクトル密度（cross-spectral density）あるいはクロスパワースペクトルと呼び、$S_{xy}(f)$ と書きます。
$$ \begin{equation} S_{xy}(f) = X^*(f) \cdot Y(f) \end{equation} $$
これが相互相関定理（cross-correlation theorem）です。
$$ \begin{equation} R_{xy}(\tau) \xleftrightarrow{\mathcal{F}} S_{xy}(f) = X^*(f) Y(f) \end{equation} $$
相互スペクトル密度の意味
相互スペクトル密度 $S_{xy}(f)$ は一般に複素数です。その大きさ $|S_{xy}(f)|$ は周波数 $f$ において2つの信号が共有するパワーを表し、偏角 $\angle S_{xy}(f)$ は周波数 $f$ における位相差を表します。
もし $y(t) = x(t – \tau_0)$（つまり $y$ が $x$ の単純な遅延コピー）であれば、$Y(f) = X(f) e^{-j2\pi f \tau_0}$ なので
$$ S_{xy}(f) = X^*(f) \cdot X(f) e^{-j2\pi f \tau_0} = |X(f)|^2 e^{-j2\pi f \tau_0} $$
位相が $-2\pi f \tau_0$ と周波数に比例するため、位相スペクトルの傾きから遅延 $\tau_0$ を読み取ることもできます。
自己相関の場合との比較
$y = x$ の場合、$S_{xy}(f) = X^*(f) X(f) = |X(f)|^2 = S_{xx}(f)$ となり、パワースペクトル密度に帰着します。パワースペクトル密度は実数かつ非負ですが、相互スペクトル密度は一般に複素数です。この違いは、自己相関が偶関数であるのに対し相互相関が非対称であることと対応しています（偶関数のフーリエ変換は実数になるため）。
フーリエ変換との関係がわかると、相互相関の計算を周波数領域で効率的に行えるようになります。しかしその前に、相互相関関数のもう1つの重要な変種である「正規化相関」を見ておきましょう。これにより、信号の振幅に依存しない、純粋な「形の類似度」を測ることができます。
正規化相関とピアソン相関係数
なぜ正規化が必要か
相互相関関数 $R_{xy}(\tau)$ の値は、信号の振幅に依存します。$x(t)$ の振幅を2倍にすれば $R_{xy}(\tau)$ も2倍になります。つまり、$R_{xy}$ の値が大きいからといって「形がよく似ている」とは限りません — 単に信号が大きいだけかもしれないのです。
これは日常的な場面でも身に覚えがあるはずです。2つの株価チャートが似た形をしているかを比べたいとき、一方が100円台、他方が10,000円台なら、単純に値を掛けた合計を比較しても意味がありません。形の類似度を比べるには、まずスケールを揃える必要があります。
正規化相互相関関数の定義
正規化相互相関関数（normalized cross-correlation, NCC）は、相互相関をエネルギーで割ることで $[-1, 1]$ の範囲に正規化した量です。
$$ \begin{equation} \rho_{xy}(\tau) = \frac{R_{xy}(\tau)}{\sqrt{R_{xx}(0) \cdot R_{yy}(0)}} \end{equation} $$
分母は $x$ と $y$ のエネルギーの幾何平均です。コーシー・シュワルツの不等式により $|\rho_{xy}(\tau)| \leq 1$ が保証されます。
$\rho_{xy}(\tau) = 1$ のとき、$x(t)$ と $y(t+\tau)$ は定数倍の違いを除いて完全に一致します。$\rho_{xy}(\tau) = -1$ なら完全に逆位相、$\rho_{xy}(\tau) = 0$ なら無相関です。
ピアソン相関係数との関係
統計学で頻出するピアソン相関係数は、2つの確率変数 $X, Y$ に対して
$$ r = \frac{\mathrm{Cov}(X, Y)}{\sigma_X \sigma_Y} = \frac{E[(X – \mu_X)(Y – \mu_Y)]}{\sqrt{E[(X-\mu_X)^2] \cdot E[(Y-\mu_Y)^2]}} $$
と定義されます。分子は共分散、分母は標準偏差の積です。
正規化相互相関関数は、これを「時系列の各時刻」を確率変数のサンプルと見なし、「遅延 $\tau$ を導入したもの」と解釈できます。平均値を差し引いたバージョン（ゼロ平均信号に対する正規化相互相関）はピアソン相関係数の自然な一般化です。
信号処理では、実用上、局所的な正規化を行うことも多いです。テンプレートマッチングでは、テンプレートの長さに対応する区間ごとにエネルギーを計算して正規化します。この局所正規化については、テンプレートマッチングのセクションで具体的に扱います。
ここまでで相互相関関数の理論的な土台が整いました。ここからは、この道具が実際にどのように使われるか、代表的な2つの応用を見ていきましょう。
遅延推定への応用 — GPS・レーダーの原理
基本アイデア
相互相関関数の最も重要な応用の1つが時間遅延推定（time delay estimation, TDE）です。既知の信号 $s(t)$ を送信し、ターゲットで反射して戻ってきた受信信号 $r(t)$ との相互相関を計算します。受信信号が
$$ r(t) = \alpha \, s(t – \tau_0) + n(t) $$
と表せるとき（$\alpha$ は減衰係数、$\tau_0$ は遅延、$n(t)$ はノイズ）、相互相関関数 $R_{sr}(\tau)$ は $\tau = \tau_0$ でピークを持ちます。
これを直感的に理解しましょう。$R_{sr}(\tau)$ の計算では、送信信号 $s(t)$ を $\tau$ だけずらしながら受信信号 $r(t)$ に重ねています。$\tau = \tau_0$ のとき、$s(t+\tau_0)$ と $r(t)$ 中の $s(t-\tau_0+\tau_0) = s(t)$ が完全に位相が合い、積の合計が最大になります。一方、ノイズ $n(t)$ は $s(t)$ と無相関なので、積分すると打ち消し合います。
数式での確認
具体的に計算してみましょう。$r(t) = \alpha s(t – \tau_0) + n(t)$ を相互相関の定義に代入します。
$$ R_{sr}(\tau) = \int s(t) \, r(t + \tau) \, dt $$
$r(t+\tau)$ を展開して
$$ = \int s(t) \left[ \alpha \, s(t + \tau – \tau_0) + n(t + \tau) \right] dt $$
線形性から2つの積分に分離すると
$$ = \alpha \int s(t) \, s(t + \tau – \tau_0) \, dt + \int s(t) \, n(t + \tau) \, dt $$
第1項は送信信号の自己相関 $\alpha R_{ss}(\tau – \tau_0)$ であり、第2項は信号とノイズの相互相関 $R_{sn}(\tau)$ です。
$$ \begin{equation} R_{sr}(\tau) = \alpha \, R_{ss}(\tau – \tau_0) + R_{sn}(\tau) \end{equation} $$
$s(t)$ と $n(t)$ が無相関であれば $R_{sn}(\tau) \approx 0$ となり
$$ R_{sr}(\tau) \approx \alpha \, R_{ss}(\tau – \tau_0) $$
自己相関関数は $\tau = 0$ で最大値を取る（性質1）ので、$R_{ss}(\tau – \tau_0)$ は $\tau = \tau_0$ で最大値を取ります。したがって、$R_{sr}(\tau)$ のピーク位置から $\tau_0$ を読み取ることができるのです。
GPSでの応用
GPS衛星はPRN符号（Pseudo-Random Noise code）と呼ばれる疑似乱数系列を送信しています。受信機は、自分が生成した同じPRN符号と受信信号の相互相関を計算し、ピーク位置から信号の到達時間を求めます。4つ以上の衛星からの到達時間差を連立すると、受信機の3次元位置と時計のずれが算出できます。
PRN符号が使われる理由は、自己相関特性が優れているからです。PRN符号の自己相関は $\tau = 0$ で鋭いピークを持ち、それ以外ではほぼゼロになります。このインパルス状の自己相関があるおかげで、相互相関のピークも鋭くなり、遅延推定の精度が高まります。
レーダーでの応用
レーダーもGPSと同じ原理を利用しています。送信パルスと受信信号の相互相関を計算し、ピーク位置から反射体までの往復時間 $\tau_0$ を求めます。距離 $d$ は光速 $c$ を用いて
$$ d = \frac{c \, \tau_0}{2} $$
で算出されます（往復なので2で割る）。パルス圧縮レーダーでは、チャープ信号（周波数が時間変化する信号）を使い、送信信号との相互相関で高い距離分解能を実現します。これはマッチドフィルタの理論と深く関わっています。
遅延推定では「信号全体」の時間ずれを求めましたが、次に見るテンプレートマッチングでは、「長い信号の中から短いパターンを探す」という問題に相互相関を適用します。
テンプレートマッチング
テンプレートマッチングとは
テンプレートマッチングとは、長い信号（あるいは画像）の中に、探したいパターン（テンプレート）がどこに存在するかを見つけ出す手法です。たとえば、数時間分の音声データの中から特定のキーワードが発話された時刻を検出する、工場のラインで製品画像の中から欠陥パターンを探す、といった場面で使われます。
考え方は遅延推定と本質的に同じです。テンプレート $h(t)$（長さ $T_h$ の短い信号）を長い信号 $x(t)$ に沿ってスライドさせ、各位置で相互相関を計算します。相互相関が大きい位置に、テンプレートと似た波形が存在しています。
局所正規化の重要性
テンプレートマッチングでは、信号の局所的な振幅変動の影響を除去するために局所正規化が不可欠です。たとえば、信号の中に振幅が大きい区間があると、その区間ではテンプレートとの形が違っていても単純な相互相関の値が大きくなってしまいます。
局所正規化相互相関（locally normalized cross-correlation）は、テンプレートの長さに対応する区間ごとに正規化を行います。離散信号の場合、位置 $m$ における値は
$$ \begin{equation} \rho(m) = \frac{\sum_{n=0}^{N_h – 1} h(n) \, x(m+n)}{\sqrt{\sum_{n=0}^{N_h-1} h(n)^2 \cdot \sum_{n=0}^{N_h-1} x(m+n)^2}} \end{equation} $$
と定義されます。ここで $N_h$ はテンプレートの長さです。分母が位置 $m$ ごとに変わるため、信号の振幅変動の影響が抑えられます。
$\rho(m)$ が閾値（たとえば 0.8）を超える位置を「マッチ」と判定します。この閾値はタスクに応じて調整します。偽陽性を減らしたければ閾値を上げ、見逃しを減らしたければ下げるというトレードオフがあります。
2次元への拡張
画像処理でのテンプレートマッチングは、1次元の相互相関を2次元に拡張したものです。テンプレート画像を対象画像上で上下左右にスライドさせ、各位置で正規化相互相関を計算します。OpenCVの matchTemplate 関数がまさにこの処理を行っています。
ここまでで理論と応用の概要を把握しました。次に、計算の実装に目を向けましょう。相互相関の定義式をそのまま計算すると、信号長 $N$ に対して $O(N^2)$ の計算量がかかります。FFTを使えば、これを $O(N \log N)$ に高速化できます。
離散相関とFFTによる高速計算
離散相互相関関数の定義
連続信号に対する相互相関関数の定義を、離散信号に適用します。長さ $N$ の離散信号 $x[n]$ と $y[n]$ に対して、離散相互相関関数は
$$ \begin{equation} R_{xy}[m] = \sum_{n=0}^{N-1} x[n] \, y[n + m] \end{equation} $$
と定義されます。遅延（ラグ） $m$ は $-(N-1) \leq m \leq N-1$ の範囲を取り、出力の長さは $2N – 1$ です。$n + m$ が $[0, N-1]$ の範囲外になる場合は $y[n+m] = 0$（ゼロパディング）とするか、あるいは和の範囲を適切に制限します。
直接計算の計算量
定義式をそのまま実装すると、各ラグ $m$ に対して $O(N)$ の掛け算と足し算が必要です。ラグの総数が $2N – 1$ なので、全体の計算量は $O(N^2)$ です。$N = 10^6$（サンプリングレート 1 MHz で 1 秒間の信号）の場合、$10^{12}$ 回の演算が必要になり、リアルタイム処理は困難です。
FFTによる高速化
相互相関定理を離散版に適用すれば、FFTで高速に計算できます。アルゴリズムは以下の手順です。
ステップ1: $x[n]$ と $y[n]$ をそれぞれ長さ $2N – 1$ 以上にゼロパディングします。これは円状相関と線形相関の違いを吸収するためです。FFTは周期的な信号を前提としているため、ゼロパディングなしでは巻き込みが発生します。
ステップ2: ゼロパディングした信号のFFTを計算します。
$$ X[k] = \text{FFT}(x_{\text{padded}}), \quad Y[k] = \text{FFT}(y_{\text{padded}}) $$
ステップ3: 相互スペクトル密度を計算します。
$$ S_{xy}[k] = X^*[k] \cdot Y[k] $$
ステップ4: 逆FFTで時間領域に戻します。
$$ R_{xy}[m] = \text{IFFT}(S_{xy}[k]) $$
各FFTの計算量は $O(N \log N)$ であり、乗算は $O(N)$ なので、全体の計算量は $O(N \log N)$ です。$N = 10^6$ の場合でも約 $2 \times 10^7$ 回の演算で済み、直接計算の5万分の1になります。
NumPy/SciPyでの実装
Pythonでは numpy.correlate や scipy.signal.correlate がこの計算を行います。scipy.signal.correlate はデフォルトで信号長に応じて直接法とFFT法を自動選択する method='auto' を使用します。method='fft' を指定すればFFT法を強制できます。
また、scipy.signal.fftconvolve を使って相互相関を計算することもできます。相関が「時間反転 + 畳み込み」であることを利用して
$$ R_{xy}[m] = (x[-n] * y[n])[m] $$
を計算します。scipy.signal.correlate の内部実装も、FFTモードではこの方法を使っています。
それでは、ここまでの理論をPythonコードで実際に動かして確認しましょう。まず遅延推定の実装から始めます。
Python実装1: 相互相関による遅延推定
問題設定
ここでは、以下のシナリオをシミュレーションします。
- 既知の送信信号（チャープ信号）を生成する
- 送信信号を一定時間遅延させ、ノイズを加えて受信信号とする
- 相互相関を計算して遅延量を推定し、真の値と比較する
チャープ信号を使う理由は、広帯域にわたってエネルギーを持つため自己相関のピークが鋭く、遅延推定の精度が高いためです。これはレーダーのパルス圧縮と同じ原理です。
まず、信号の生成と相互相関の計算を行います。
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
# パラメータ設定
fs = 10000 # サンプリング周波数 [Hz]
T = 0.1 # 送信信号の長さ [s]
true_delay = 0.035 # 真の遅延 [s]
snr_db = 5 # 信号対雑音比 [dB]
# 時間軸
t_sig = np.arange(0, T, 1/fs)
# 送信信号: チャープ（100Hzから2000Hzに線形掃引）
f0, f1 = 100, 2000
tx_signal = signal.chirp(t_sig, f0, T, f1)
# 受信信号: 遅延 + ノイズ
delay_samples = int(true_delay * fs)
total_len = len(tx_signal) + delay_samples + 500 # 余裕を持たせる
rx_signal = np.zeros(total_len)
rx_signal[delay_samples:delay_samples + len(tx_signal)] = tx_signal
# ノイズ付加
sig_power = np.mean(tx_signal**2)
noise_power = sig_power / (10**(snr_db / 10))
noise = np.sqrt(noise_power) * np.random.randn(total_len)
rx_signal += noise
# ゼロパディングした送信信号
tx_padded = np.zeros(total_len)
tx_padded[:len(tx_signal)] = tx_signal
上のコードでは、チャープ信号を生成し、350サンプル（0.035秒）遅延させた上でSNR 5 dBのノイズを加えています。SNR 5 dBは信号パワーがノイズパワーの約3.2倍という比較的ノイジーな環境です。
次に、相互相関を計算してピーク位置から遅延を推定します。
# 相互相関の計算（FFT法）
corr = signal.correlate(rx_signal, tx_signal, mode='full', method='fft')
lags = signal.correlation_lags(len(rx_signal), len(tx_signal), mode='full')
# ピーク検出
peak_idx = np.argmax(corr)
estimated_delay_samples = lags[peak_idx]
estimated_delay = estimated_delay_samples / fs
print(f"真の遅延: {true_delay:.4f} s ({delay_samples} samples)")
print(f"推定遅延: {estimated_delay:.4f} s ({estimated_delay_samples} samples)")
print(f"推定誤差: {abs(estimated_delay - true_delay)*1e6:.1f} μs")
上のコードでは、SciPyの signal.correlate を用いてFFT法で相互相関を計算し、ピーク位置のラグから遅延を推定しています。SNR 5 dBという比較的ノイジーな環境でも、相互相関のピーク位置は真の遅延とほぼ一致するはずです。これは、相関演算が信号成分をコヒーレントに加算する一方で、ノイズ成分は位相がランダムなため積分により打ち消し合うためです。
続いて、結果を可視化します。
fig, axes = plt.subplots(3, 1, figsize=(12, 10))
# 送信信号
t_tx = np.arange(len(tx_padded)) / fs
axes[0].plot(t_tx * 1000, tx_padded, color='#00bcd4', linewidth=0.8)
axes[0].set_title('Transmitted Signal (Chirp)', fontsize=13)
axes[0].set_xlabel('Time [ms]')
axes[0].set_ylabel('Amplitude')
axes[0].set_xlim([0, total_len / fs * 1000])
axes[0].grid(True, alpha=0.3)
# 受信信号
t_rx = np.arange(len(rx_signal)) / fs
axes[1].plot(t_rx * 1000, rx_signal, color='#ff9800', linewidth=0.5, alpha=0.8)
axes[1].axvline(true_delay * 1000, color='red', linestyle='--', label=f'True delay = {true_delay*1000:.1f} ms')
axes[1].set_title('Received Signal (Delayed + Noise, SNR=5dB)', fontsize=13)
axes[1].set_xlabel('Time [ms]')
axes[1].set_ylabel('Amplitude')
axes[1].legend()
axes[1].grid(True, alpha=0.3)
# 相互相関
lag_time = lags / fs * 1000
axes[2].plot(lag_time, corr, color='#4caf50', linewidth=0.8)
axes[2].axvline(estimated_delay * 1000, color='red', linestyle='--',
label=f'Estimated delay = {estimated_delay*1000:.1f} ms')
axes[2].axvline(true_delay * 1000, color='blue', linestyle=':',
label=f'True delay = {true_delay*1000:.1f} ms', alpha=0.7)
axes[2].set_title('Cross-Correlation', fontsize=13)
axes[2].set_xlabel('Lag [ms]')
axes[2].set_ylabel('Correlation')
axes[2].legend()
axes[2].grid(True, alpha=0.3)
axes[2].set_xlim([-20, 120])
plt.tight_layout()
plt.savefig('cross_correlation_delay_estimation.png', dpi=150, bbox_inches='tight')
plt.show()
上のグラフから、以下の重要なポイントが読み取れます。
-
受信信号（中段）ではノイズに埋もれて遅延を目視で特定するのは困難ですが、相互相関（下段）には明確なピークが現れます。これが相互相関のノイズ耐性の本質です。信号成分は $\tau = \tau_0$ で位相が揃って強め合い、ノイズ成分は位相がランダムなので積分で打ち消し合います。
-
ピーク位置は真の遅延とほぼ一致しています。サンプリング周波数の逆数（$1/f_s = 0.1$ ms）が遅延推定の分解能の限界です。サブサンプル精度が必要な場合は、ピーク付近を放物線でフィッティングするなどの補間手法が使われます。
-
ピークの幅が狭いことに注目してください。これはチャープ信号の広帯域特性によるものです。帯域幅 $B$ のチャープ信号の自己相関のメインローブ幅は約 $1/B$ です。$B = 2000 – 100 = 1900$ Hz なので、ピーク幅は約 $0.53$ ms となり、高い時間分解能が得られています。
Python実装2: テンプレートマッチング
問題設定
次に、テンプレートマッチングの実装を行います。長い信号の中にテンプレートが複数箇所に埋め込まれている状況を想定します。
- 5秒間のノイズ混じりの信号を生成する
- 3箇所にテンプレート（短いパルス列）を埋め込む
- 正規化相互相関を計算してテンプレートの位置を検出する
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
# パラメータ
fs = 8000 # サンプリング周波数 [Hz]
duration = 5 # 信号長 [s]
N = fs * duration
# テンプレート: 特徴的なパルス列（AM変調された正弦波バースト）
t_template = np.arange(0, 0.05, 1/fs) # 50ms
template = np.sin(2 * np.pi * 800 * t_template) * np.hanning(len(t_template))
N_h = len(template)
# 長い信号の生成（背景ノイズ）
np.random.seed(42)
long_signal = 0.3 * np.random.randn(N)
# テンプレートを3箇所に埋め込む（異なる振幅で）
true_positions = [0.8, 2.3, 3.9] # 秒
amplitudes = [1.0, 0.6, 1.2] # 各挿入箇所の振幅
for pos, amp in zip(true_positions, amplitudes):
idx = int(pos * fs)
long_signal[idx:idx + N_h] += amp * template
このコードでは、800 Hzの正弦波にハニング窓をかけた50 msのバースト信号をテンプレートとし、5秒間の背景ノイズの中の3箇所に異なる振幅で埋め込んでいます。テンプレートの振幅が場所によって異なるため、正規化が重要になります。
次に、正規化相互相関を計算します。
# 正規化相互相関の計算
def normalized_cross_correlation(x, template):
"""局所正規化相互相関を計算する"""
N_x = len(x)
N_h = len(template)
N_out = N_x - N_h + 1
# テンプレートのエネルギー（固定値）
template_energy = np.sqrt(np.sum(template**2))
ncc = np.zeros(N_out)
for m in range(N_out):
segment = x[m:m + N_h]
segment_energy = np.sqrt(np.sum(segment**2))
if segment_energy > 1e-10: # ゼロ除算防止
ncc[m] = np.sum(template * segment) / (template_energy * segment_energy)
return ncc
ncc = normalized_cross_correlation(long_signal, template)
上の関数は定義に忠実な直接実装です。各位置でテンプレート長の区間を切り出し、テンプレートとの内積を局所エネルギーの幾何平均で割っています。$O(N \cdot N_h)$ の計算量がかかりますが、原理の理解には十分です。
検出結果を可視化します。
fig, axes = plt.subplots(3, 1, figsize=(14, 10))
# 元の信号
t_long = np.arange(N) / fs
axes[0].plot(t_long, long_signal, color='#78909c', linewidth=0.4, alpha=0.8)
for pos in true_positions:
axes[0].axvline(pos, color='red', linestyle='--', alpha=0.7, linewidth=1)
axes[0].set_title('Signal with Embedded Templates (3 locations)', fontsize=13)
axes[0].set_xlabel('Time [s]')
axes[0].set_ylabel('Amplitude')
axes[0].grid(True, alpha=0.3)
# テンプレート
axes[1].plot(t_template * 1000, template, color='#e91e63', linewidth=1.5)
axes[1].set_title('Template (50ms burst)', fontsize=13)
axes[1].set_xlabel('Time [ms]')
axes[1].set_ylabel('Amplitude')
axes[1].grid(True, alpha=0.3)
# 正規化相互相関
t_ncc = np.arange(len(ncc)) / fs
axes[2].plot(t_ncc, ncc, color='#4caf50', linewidth=0.8)
# 閾値とピーク検出
threshold = 0.6
axes[2].axhline(threshold, color='orange', linestyle='--', label=f'Threshold = {threshold}')
# ピーク検出
peaks, properties = signal.find_peaks(ncc, height=threshold, distance=int(0.1*fs))
for p in peaks:
axes[2].axvline(t_ncc[p], color='red', linestyle=':', alpha=0.7)
axes[2].annotate(f'{t_ncc[p]:.2f}s\n(NCC={ncc[p]:.3f})',
xy=(t_ncc[p], ncc[p]),
xytext=(t_ncc[p]+0.15, ncc[p]),
fontsize=9, color='red',
arrowprops=dict(arrowstyle='->', color='red', alpha=0.7))
axes[2].set_title('Normalized Cross-Correlation', fontsize=13)
axes[2].set_xlabel('Time [s]')
axes[2].set_ylabel('NCC')
axes[2].legend(loc='upper right')
axes[2].grid(True, alpha=0.3)
axes[2].set_ylim([-0.5, 1.1])
plt.tight_layout()
plt.savefig('template_matching.png', dpi=150, bbox_inches='tight')
plt.show()
# 検出結果の表示
print("=== テンプレート検出結果 ===")
print(f"真の位置: {true_positions}")
print(f"検出位置: {[f'{t_ncc[p]:.3f}' for p in peaks]}")
print(f"NCC値: {[f'{ncc[p]:.3f}' for p in peaks]}")
このグラフから、以下のことが確認できます。
-
3箇所すべてのテンプレートが正しく検出されています。正規化相互相関のピーク位置は真の埋め込み位置とほぼ一致しています。
-
振幅が異なっても検出できています。2つ目のテンプレート（振幅0.6）はNCC値がやや低めになりますが、それでも閾値を十分に超えています。これが正規化の効果です。正規化しない相互相関を使った場合、振幅の大きい区間で偽検出が発生する可能性があります。
-
背景ノイズの区間ではNCC値が低く抑えられています。ノイズはテンプレートと形が似ていないため、正規化相互相関は小さな値にとどまります。閾値を適切に設定すれば、偽陽性なくテンプレートを検出できています。
Python実装3: FFTによる高速相互相関
直接法とFFT法の比較
最後に、FFTによる高速化の効果を実際に計測してみましょう。
import numpy as np
import time
from scipy import signal
def cross_corr_direct(x, y):
"""直接法による相互相関（NumPyのcorrelateを使用）"""
return np.correlate(x, y, mode='full')
def cross_corr_fft(x, y):
"""FFT法による相互相関"""
n = len(x) + len(y) - 1
fft_size = 2**int(np.ceil(np.log2(n))) # 2のべき乗に切り上げ
X = np.fft.fft(x, fft_size)
Y = np.fft.fft(y, fft_size)
corr = np.fft.ifft(np.conj(X) * Y).real
return corr[:n]
# 異なる信号長で計測
lengths = [256, 512, 1024, 2048, 4096, 8192, 16384, 32768]
times_direct = []
times_fft = []
times_scipy = []
for N in lengths:
x = np.random.randn(N)
y = np.random.randn(N)
# 直接法
start = time.perf_counter()
for _ in range(5):
cross_corr_direct(x, y)
times_direct.append((time.perf_counter() - start) / 5)
# FFT法（自作）
start = time.perf_counter()
for _ in range(5):
cross_corr_fft(x, y)
times_fft.append((time.perf_counter() - start) / 5)
# SciPy（自動選択）
start = time.perf_counter()
for _ in range(5):
signal.correlate(x, y, method='auto')
times_scipy.append((time.perf_counter() - start) / 5)
上のコードでは、信号長 $N$ を256から32768まで変化させながら、直接法（numpy.correlate）、自作FFT法、SciPyの自動選択の3手法で計算時間を計測しています。各条件で5回の計測を行い平均をとることで、計測のばらつきを抑えています。信号長が大きくなるにつれ、直接法の $O(N^2)$ とFFT法の $O(N \log N)$ の計算量の差が顕著に現れるはずです。
計測結果を可視化します。
fig, ax = plt.subplots(figsize=(10, 6))
ax.loglog(lengths, times_direct, 'o-', color='#f44336', linewidth=2, label='Direct (numpy.correlate)')
ax.loglog(lengths, times_fft, 's-', color='#2196f3', linewidth=2, label='FFT (manual)')
ax.loglog(lengths, times_scipy, '^-', color='#4caf50', linewidth=2, label='SciPy auto')
# 理論的なスケーリング（参考線）
N_ref = np.array(lengths, dtype=float)
ax.loglog(N_ref, times_direct[0] * (N_ref / N_ref[0])**2, '--', color='#f44336', alpha=0.3, label=r'$O(N^2)$ reference')
ax.loglog(N_ref, times_fft[0] * (N_ref * np.log2(N_ref)) / (N_ref[0] * np.log2(N_ref[0])),
'--', color='#2196f3', alpha=0.3, label=r'$O(N \log N)$ reference')
ax.set_xlabel('Signal Length N', fontsize=12)
ax.set_ylabel('Computation Time [s]', fontsize=12)
ax.set_title('Cross-Correlation: Direct vs FFT Computation Time', fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3, which='both')
plt.tight_layout()
plt.savefig('fft_speedup.png', dpi=150, bbox_inches='tight')
plt.show()
# 速度比の表示
print("=== 計算時間の比較 ===")
print(f"{'N':>8} | {'Direct':>10} | {'FFT':>10} | {'Speedup':>8}")
print("-" * 45)
for i, N in enumerate(lengths):
speedup = times_direct[i] / times_fft[i]
print(f"{N:>8} | {times_direct[i]:>10.6f} | {times_fft[i]:>10.6f} | {speedup:>7.1f}x")
このグラフと数値結果から、以下のことがわかります。
-
直接法の計算時間は $O(N^2)$ に従って増大しています。信号長を2倍にすると計算時間はおよそ4倍になっています。赤い破線（$O(N^2)$ の参考線）との一致がこれを裏付けています。
-
FFT法は $O(N \log N)$ で増大しており、信号長が大きくなるほど直接法との差が開きます。$N = 32768$ では数十倍以上の高速化が得られています。
-
SciPyの
method='auto'は信号長に応じて適切な方法を自動選択しています。短い信号では直接法のオーバーヘッドが小さいためFFTとの差は僅かですが、長い信号ではFFT法が選択され高速に動作しています。
計算結果の一致確認
FFT法で得られた結果が直接法と一致することも確認しておきます。
# 結果の一致を確認
N_test = 1024
x_test = np.random.randn(N_test)
y_test = np.random.randn(N_test)
corr_direct = cross_corr_direct(x_test, y_test)
corr_fft = cross_corr_fft(x_test, y_test)
max_error = np.max(np.abs(corr_direct - corr_fft))
print(f"直接法とFFT法の最大誤差: {max_error:.2e}")
最大誤差は浮動小数点演算の丸め誤差（$10^{-10}$ 程度）に収まるはずです。これにより、FFT法が直接法と数学的に等価であることが数値的にも確認できます。FFTの内部では浮動小数点の加減算が直接法とは異なる順序で行われるため、完全にゼロにはなりませんが、実用上は無視できる差です。
相互相関の実用上の注意点
ここまでの理論と実装を踏まえ、相互相関を実務で使う際に気をつけるべきポイントを整理します。
ゼロパディングの重要性
FFTで相互相関を計算する際、ゼロパディングを忘れると円状相関（circular correlation）が計算されてしまいます。円状相関は信号が周期的に繰り返すと仮定した場合の相関であり、信号の端が反対側に巻き込まれてしまいます。線形相関を得るには、2つの信号の長さの和から1を引いた長さ以上にゼロパディングする必要があります。
サブサンプル精度の遅延推定
離散相関のピーク位置はサンプル単位でしか得られないため、遅延推定の精度はサンプリング間隔 $\Delta t = 1/f_s$ に制限されます。GPS受信機やレーダーではこれでは不十分な場合があり、以下の手法でサブサンプル精度を実現します。
放物線補間: ピーク位置とその両隣の3点に放物線をフィッティングし、放物線の頂点をサブサンプル精度のピーク位置とします。ピークの位置を $m_0$、相関値を $R[m_0-1], R[m_0], R[m_0+1]$ とすると
$$ \hat{m} = m_0 + \frac{R[m_0-1] – R[m_0+1]}{2(R[m_0-1] – 2R[m_0] + R[m_0+1])} $$
でサブサンプル精度の推定が得られます。
ゼロパディングFFT法: 周波数領域でゼロパディング（FFTのサイズを信号長より大きくする）することで、逆FFT後の時間分解能を上げる方法もあります。
正規化方法の選択
前述のように、正規化の方法は用途によって使い分けます。
| 方法 | 定義 | 用途 |
|---|---|---|
| 非正規化 | $R_{xy}(\tau)$ そのまま | 遅延推定（SNRが既知の場合） |
| エネルギー正規化 | $R_{xy}/\sqrt{R_{xx}(0) R_{yy}(0)}$ | 全体的な類似度比較 |
| 局所正規化 | 区間ごとにエネルギーで割る | テンプレートマッチング |
| 係数正規化 | 平均除去 + 標準偏差で割る | 統計的相関分析（ピアソン相関） |
目的に応じて適切な正規化を選ぶことが、誤検出を防ぎ、ロバストな結果を得るための鍵です。
多次元への拡張
ここまでは1次元の信号を扱いましたが、相互相関は任意の次元に拡張できます。2次元画像のテンプレートマッチング、3次元ボリュームデータの位置合わせなど、基本原理は全く同じです。計算量の増大に対しても、多次元FFTによる高速化が有効です。
まとめ
本記事では、相互相関関数の理論を基礎から解説し、Pythonで実装しました。
- 相関の直感: 2つの信号をずらしながら重ね合わせ、一致度を遅延の関数として表したもの。畳み込みとは「一方を時間反転するかしないか」の違いがある
- 自己相関関数: 信号自身との相関。$\tau = 0$ で最大値を取り、偶関数。フーリエ変換するとパワースペクトル密度になる（ウィーナー・ヒンチンの定理）
- 相互相関関数: 2つの異なる信号の類似度を遅延の関数として表す。一般に非対称（$R_{xy}(\tau) = R_{yx}(-\tau)$）
- フーリエ変換との関係: 相互相関のフーリエ変換は $X^*(f)Y(f)$（相互スペクトル密度）。この関係によりFFTで $O(N \log N)$ の高速計算が可能
- 正規化相関: 振幅に依存しない類似度の尺度。ピアソン相関係数の時系列版
- 遅延推定: GPS・レーダーの根幹技術。相互相関のピーク位置が時間遅延を教える
- テンプレートマッチング: 局所正規化相互相関により、振幅の異なるパターンも検出可能
相互相関関数は、信号処理の中で最も基本的かつ実用的な道具の1つです。本記事で扱った遅延推定とテンプレートマッチングは、より高度な手法 — マッチドフィルタ、適応フィルタ、MUSIC法 — の出発点にもなっています。
次のステップとして、以下の記事も参考にしてください。
- パワースペクトル密度の理論 — 自己相関とスペクトルの関係をさらに深掘り
- マッチドフィルタの理論 — 相互相関を最適化の観点から見た「SNR最大化フィルタ」