はじめに

2つの時系列データを見比べていると、

片方の変化が、もう片方より少し遅れて起きている気がする

という場面があります。

例えば、一般的なデータ解析では次のような状況です。

入力信号と出力信号の時間差を見たい
2つのセンサデータのズレを補正したい
ある波形が、別の波形に何秒遅れて似た形になるか知りたい
目で見るとズレていそうだが、何秒ズレているか数値で知りたい

このようなときに使える基本的な方法が、相互相関です。

この記事では、実データではなく、Pythonで作成したダミーデータを使います。
あらかじめ「真の時間遅れ」を入れておき、相互相関によってその遅れを推定できるかを確認します。

この記事のゴール

この記事では、次のことを目標にします。

時系列データとは何かを確認する
時間遅れ、ラグ、相関、相互相関の意味を理解する
片方の時系列を少しずつずらして、似ている位置を探す
Pythonでダミーデータを作り、時間遅れを推定する
推定した時間遅れを使って、2つの波形をそろえる

結論から言うと、この記事のコードでは、真の時間遅れを 0.80 秒 として作ったデータに対して、相互相関により 0.80 秒 と推定できます。

用語の整理

時系列データとは

時系列データとは、時間の順番に並んだデータです。

例えば、次のようなものです。

時刻:   0.00, 0.01, 0.02, 0.03, ...
値:     1.20, 1.25, 1.28, 1.22, ...

センサ信号、音声、株価、気温、心拍、装置のログなど、多くのデータは時系列として扱えます。

時間遅れとは

2つの時系列 $x(t)$ と $y(t)$ があるとします。

もし、$y(t)$ が $x(t)$ より少し遅れて同じように変化するなら、例えば次のように書けます。

$$
y(t) \approx x(t - \tau)
$$

ここで、$\tau$ が時間遅れです。

例えば $\tau = 0.8$ 秒なら、

$y$ は $x$ より 0.8 秒遅れている

という意味になります。

ラグとは

ラグとは、片方の時系列をどれだけずらすかを表す量です。

例えば、サンプリング間隔が 0.01 秒 のデータで、80点ずらすと、

$$
80 \times 0.01 = 0.8 \ \mathrm{s}
$$

なので、ラグは 0.8 秒 です。

相関とは

相関は、2つのデータがどれくらい一緒に変化するかを見る指標です。

ざっくり言うと、

一緒に増えたり減ったりする：相関が高い
片方が増えるともう片方が減る：負の相関
あまり関係なく動く：相関が低い

というイメージです。

相関係数は、よく $-1$ から $1$ の範囲で表されます。

 1 に近い  : よく似ている
 0 に近い  : あまり似ていない
-1 に近い  : 反対向きに似ている

相互相関とは

相互相関は、片方の時系列を少しずつずらしながら、もう片方とどれくらい似ているかを調べる方法です。

イメージとしては、次のような操作です。

yを -1.0 秒ずらして x と比べる
yを -0.9 秒ずらして x と比べる
yを -0.8 秒ずらして x と比べる
...
yを  0.0 秒ずらして x と比べる
...
yを +0.8 秒ずらして x と比べる
yを +0.9 秒ずらして x と比べる
yを +1.0 秒ずらして x と比べる

そして、最もよく似るラグを探します。

数式で書くと、連続時間では概念的に次のように表せます。

$$
C(\tau) = \int x(t)y(t+\tau)dt
$$

離散データでは、和で書けます。

$$
C(k) = \sum_t x(t)y(t+k)
$$

ただし、実際には平均値を引いたり、スケールの違いを補正したりして、正規化した相互相関を使うことが多いです。

この記事では、ラグごとに重なる部分を取り出し、その部分の相関係数を計算します。

今回の符号の約束

相互相関では、実装によってラグの符号の解釈が変わることがあります。
そのため、この記事では次のように約束します。

正のラグなら、$y$ が $x$ より遅れている

つまり、推定されたラグが +0.80 秒 なら、

$y$ は $x$ より 0.80 秒遅れている

と解釈します。

Google Colabで実行するコード

以下のコードは、Google Colabでそのまま実行できます。

このコードでは、次のことを行います。

基準となる時系列 $x$ を作る
$x$ を 0.80 秒 遅らせた時系列 $y$ を作る
$y$ にノイズを加える
相互相関で時間遅れを推定する
推定した遅れを使って $y$ を補正する
図を保存する

import os
import numpy as np
import matplotlib.pyplot as plt

# -----------------------------
# 1. 設定
# -----------------------------

SEED = 42
rng = np.random.default_rng(SEED)

# サンプリング周波数 [Hz]
# 100 Hz なので、0.01秒ごとに1点のデータを取る
fs = 100

# データの長さ [s]
duration = 20

# 真の時間遅れ [s]
# 今回は y が x より 0.80 秒遅れるように作る
delay_true = 0.80

# 探索する最大ラグ [s]
# 今回は -3秒 から +3秒 の範囲で探す
max_lag_seconds = 3.0

# 時刻配列
t = np.arange(0, duration, 1/fs)

# 図の保存先
os.makedirs("figures", exist_ok=True)


# -----------------------------
# 2. ダミー信号を作る関数
# -----------------------------

def make_clean_signal(tt):
    """
    複数の山を持つダミー信号を作る関数。
    実データではなく、相互相関の説明用の人工データです。
    """
    centers = np.array([2.5, 6.0, 9.5, 14.0, 17.0])
    widths  = np.array([0.25, 0.45, 0.30, 0.60, 0.35])
    amps    = np.array([1.00, 0.70, 1.20, 0.80, 1.10])

    signal = np.zeros_like(tt, dtype=float)

    for c, w, a in zip(centers, widths, amps):
        signal += a * np.exp(-0.5 * ((tt - c) / w)**2)

    # ゆっくりした揺らぎも少し加える
    signal += 0.15 * np.sin(2*np.pi*0.35*tt)

    return signal


# -----------------------------
# 3. x と y を作る
# -----------------------------

# 基準信号 x
x_clean = make_clean_signal(t)

# x に少しノイズを加える
x = x_clean + 0.05 * rng.normal(size=t.size)

# y は x より delay_true 秒遅れた信号にする
# t - delay_true とすることで、山が右にずれる
y_clean = make_clean_signal(t - delay_true)

# y には少し強めのノイズと、ゆっくりした揺らぎを加える
y = y_clean + 0.12 * rng.normal(size=t.size)
y += 0.05 * np.sin(2*np.pi*0.05*t)


# -----------------------------
# 4. ラグごとの相関を計算する関数
# -----------------------------

def pearson_corr(a, b):
    """
    2つの配列 a, b の相関係数を計算する。
    """
    a = a - np.mean(a)
    b = b - np.mean(b)

    denom = np.sqrt(np.sum(a*a) * np.sum(b*b))

    if denom == 0:
        return np.nan

    return np.sum(a*b) / denom


def cross_correlation_by_lag(x, y, fs, max_lag_seconds):
    """
    x と y の相互相関を、ラグごとに計算する。

    この関数では、
    正のラグ = y が x より遅れている
    という符号の約束にしています。
    """
    n = len(x)

    max_lag_samples = int(max_lag_seconds * fs)
    lags_samples = np.arange(-max_lag_samples, max_lag_samples + 1)

    correlations = []

    for lag in lags_samples:
        if lag >= 0:
            # lag > 0 のとき：
            # y のほうが遅れていると仮定し、
            # x の前半と y の後半を比較する
            x_segment = x[:n-lag]
            y_segment = y[lag:]
        else:
            # lag < 0 のとき：
            # y のほうが進んでいると仮定し、
            # x の後半と y の前半を比較する
            x_segment = x[-lag:]
            y_segment = y[:n+lag]

        corr = pearson_corr(x_segment, y_segment)
        correlations.append(corr)

    correlations = np.array(correlations)
    lags_seconds = lags_samples / fs

    best_index = np.nanargmax(correlations)

    estimated_delay = lags_seconds[best_index]
    peak_corr = correlations[best_index]

    return lags_seconds, correlations, estimated_delay, peak_corr


# -----------------------------
# 5. 相互相関で時間遅れを推定する
# -----------------------------

lags_seconds, correlations, delay_est, peak_corr = cross_correlation_by_lag(
    x=x,
    y=y,
    fs=fs,
    max_lag_seconds=max_lag_seconds
)

print("===== Result =====")
print(f"True delay      : {delay_true:.3f} s")
print(f"Estimated delay : {delay_est:.3f} s")
print(f"Error           : {(delay_est - delay_true)*1000:.1f} ms")
print(f"Peak corr       : {peak_corr:.3f}")


# -----------------------------
# 6. 推定した遅れを使って y を補正する
# -----------------------------

# y は x より delay_est 秒遅れていると推定された。
# そこで、y を delay_est 秒ぶん前に戻す。
#
# y_aligned(t) = y(t + delay_est)
#
# と考えると、y の山が左に戻り、x と重なりやすくなる。
y_aligned = np.interp(
    t + delay_est,
    t,
    y,
    left=np.nan,
    right=np.nan
)


# -----------------------------
# 7. 図1：元の2つの時系列
# -----------------------------

plt.figure(figsize=(10, 4))
plt.plot(t, x, label="x: reference signal")
plt.plot(t, y, label="y: delayed + noisy signal", alpha=0.8)
plt.title("Dummy time series: y is delayed relative to x")
plt.xlabel("time [s]")
plt.ylabel("amplitude [a.u.]")
plt.legend()
plt.tight_layout()
plt.savefig("figures/fig1_dummy_timeseries.png", dpi=160)
plt.show()


# -----------------------------
# 8. 図2：相互相関
# -----------------------------

plt.figure(figsize=(10, 4))
plt.plot(lags_seconds, correlations)
plt.axvline(delay_true, linestyle="--", label=f"true delay = {delay_true:.2f} s")
plt.axvline(delay_est, linestyle=":", label=f"estimated delay = {delay_est:.2f} s")
plt.title("Cross-correlation as a function of lag")
plt.xlabel("lag [s]")
plt.ylabel("normalized cross-correlation")
plt.legend()
plt.tight_layout()
plt.savefig("figures/fig2_cross_correlation.png", dpi=160)
plt.show()


# -----------------------------
# 9. 図3：補正後の2つの時系列
# -----------------------------

plt.figure(figsize=(10, 4))
plt.plot(t, x, label="x: reference signal")
plt.plot(t, y_aligned, label="y shifted earlier by estimated delay", alpha=0.8)
plt.title("After delay correction: signals are better aligned")
plt.xlabel("time [s]")
plt.ylabel("amplitude [a.u.]")
plt.legend()
plt.tight_layout()
plt.savefig("figures/fig3_aligned_timeseries.png", dpi=160)
plt.show()


# -----------------------------
# 10. 図4：一部を拡大して比較
# -----------------------------

plt.figure(figsize=(10, 4))

zoom = (t >= 5.0) & (t <= 8.0)

plt.plot(t[zoom], x[zoom], label="x: reference signal")
plt.plot(t[zoom], y[zoom], label="before correction", alpha=0.6)
plt.plot(t[zoom], y_aligned[zoom], label="after correction", alpha=0.9)

plt.title("Zoomed view before and after correction")
plt.xlabel("time [s]")
plt.ylabel("amplitude [a.u.]")
plt.legend()
plt.tight_layout()
plt.savefig("figures/fig4_zoom_before_after.png", dpi=160)
plt.show()

実行結果

上のコードを実行すると、次のような結果が表示されます。

===== Result =====
True delay      : 0.800 s
Estimated delay : 0.800 s
Error           : 0.0 ms
Peak corr       : 0.926

今回は、もともと $y$ を $x$ より 0.80 秒 遅らせて作りました。
相互相関による推定結果も 0.80 秒 になっているので、うまく時間遅れを推定できています。

図1：元の2つの時系列

図1では、青い波形が基準信号 $x$、オレンジ色の波形が遅れている信号 $y$ です。

$y$ は $x$ よりも全体的に右側にずれています。
つまり、$y$ のほうが遅れて変化しています。

この段階では、目で見ても「少し遅れていそう」と分かります。
しかし、目視だけでは正確に何秒遅れているかは分かりません。

そこで、相互相関を使います。

図2：ラグごとの相互相関

図2は、横軸がラグ、縦軸が相互相関です。

この図では、相互相関が最大になるラグが 0.80 秒 付近にあります。

つまり、

$y$ を 0.80 秒ぶん前に戻すと、$x$ と最もよく似る

ということです。

今回の符号の約束では、正のラグは「$y$ が $x$ より遅れている」という意味です。
したがって、推定結果は次のように解釈できます。

$y$ は $x$ より約 0.80 秒遅れている

図3：推定した遅れで補正した後

図3では、推定した時間遅れを使って、$y$ を前に戻しています。

具体的には、

$$
y_{\mathrm{aligned}}(t) = y(t + \hat{\tau})
$$

としています。

ここで、$\hat{\tau}$ は推定された時間遅れです。
今回は $\hat{\tau} = 0.80$ 秒です。

補正後の波形を見ると、$x$ と $y$ の山の位置がかなりよく重なっていることが分かります。

図4：補正前と補正後を拡大して見る

図4では、一部の時間範囲だけを拡大しています。

補正前の $y$ は、$x$ より右側にずれています。
一方で、補正後の $y$ は、$x$ とかなり近いタイミングで変化しています。

このように、相互相関を使うと、

2つの時系列がどれくらい時間的にずれているか

を定量的に推定できます。

なぜ相互相関で時間遅れが分かるのか

相互相関の考え方は、非常に素朴です。

片方の波形を少しずつずらす
ずらすたびに、もう片方の波形とどれくらい似ているかを計算する
最も似ていたズレ量を採用する

今回のように、$y$ が $x$ の遅れたバージョンであれば、正しい遅れ量で比較したときに、2つの波形が最もよく重なります。

そのため、相互相関が最大になるラグを見れば、時間遅れを推定できます。

今回のコードでやっていること

今回のコードでは、ラグごとに次のような処理をしています。

例えば、ラグが正の値のときは、

x_segment = x[:n-lag]
y_segment = y[lag:]

としています。

これは、$y$ が $x$ より遅れていると仮定して、$x$ の前半部分と $y$ の後半部分を比較する操作です。

そして、重なっている部分について相関係数を計算します。

corr = pearson_corr(x_segment, y_segment)

この計算を、負のラグから正のラグまで繰り返します。

最後に、相関係数が最大になるラグを探します。

best_index = np.nanargmax(correlations)
estimated_delay = lags_seconds[best_index]

この estimated_delay が、推定された時間遅れです。

実データで使うときの注意点

相互相関は便利ですが、実データで使うときには注意が必要です。

1. 波形が似ていないと推定は難しい

相互相関は「似ている波形のズレ」を探す方法です。
そもそも2つの波形があまり似ていない場合、相互相関のピークは明確になりません。

2. ノイズが多いとピークが不安定になる

ノイズが大きいと、相互相関のピークがずれたり、複数のピークが出たりします。
必要に応じて、前処理や平滑化を考える必要があります。

ただし、フィルタをかけると波形のタイミングや振幅が変わることもあるため、前処理の影響も確認したほうがよいです。

3. 周期的な信号ではピークが複数出る

周期的な信号では、1周期ずれても似て見えることがあります。

そのため、相互相関に複数のピークが出ることがあります。
この場合は、探索するラグの範囲を現実的な範囲に制限することが重要です。

今回のコードでも、

max_lag_seconds = 3.0

として、探索範囲を -3 秒 から +3 秒 に制限しています。

4. トレンドやベースラインのずれに注意する

2つの時系列に大きなトレンドやベースラインのずれがあると、相互相関の結果に影響することがあります。

そのような場合は、

平均を引く
トレンドを除去する
比較したい周波数帯だけを見る
解析区間を限定する

といった工夫が必要になることがあります。

5. 相関が高いことは因果関係を意味しない

相互相関で時間遅れが推定できても、それだけで因果関係が証明されたわけではありません。

例えば、

$x$ が先に変化し、$y$ が遅れて変化した

ように見えても、実際には第三の要因が両方に影響している可能性もあります。

相互相関は、あくまで

時系列の類似性と時間的なズレを調べる方法

として理解するのがよいです。

まとめ

この記事では、相互相関を使って2つの時系列の時間遅れを推定しました。

ポイントは次の通りです。

相互相関は、片方の時系列を少しずつずらしながら似ている度合いを調べる方法
相互相関が最大になるラグから、時間遅れを推定できる
今回のダミーデータでは、真の遅れ 0.80 秒 に対して、推定値も 0.80 秒 になった
推定した遅れを使うと、2つの波形を時間的にそろえられる
実データでは、ノイズ、周期性、トレンド、前処理、探索範囲に注意が必要

相互相関は、時系列解析や信号処理の基本的な道具です。
難しい名前ですが、考え方は

ずらして、比べて、一番似ているところを探す

という、とても直感的なものです。

おまけ：試してみると理解が深まる変更

以下の値を変えて実行すると、相互相関の挙動をより理解しやすくなります。

真の遅れを変える

delay_true = 1.50

真の遅れを変えると、相互相関のピーク位置も変わります。

ノイズを強くする

y = y_clean + 0.30 * rng.normal(size=t.size)

ノイズを強くすると、相互相関のピークが少し分かりにくくなります。

探索範囲を狭くする

max_lag_seconds = 0.5

真の遅れが探索範囲の外にあると、正しい遅れを推定できません。

周期的な信号にする

山が少ない信号ではなく、ほぼ正弦波だけの信号にすると、相互相関に複数のピークが出やすくなります。

周期信号では、

どの周期のズレを正しい遅れと見るか

が問題になることがあります。

おわりに

相互相関は、時系列データの時間的なズレを調べるための基本的な方法です。

今回はダミーデータを使いましたが、考え方自体は、センサデータ、実験データ、信号処理、画像解析後の時系列など、さまざまな場面で使えます。

まずは今回のような人工データで動きを確認してから、実データに適用すると、結果の解釈がしやすくなります。