Python で自己相関関数とかやる（途中・完成予定無し）
随分前に書こうと思って途中まで書いたけど面倒くさくなってやめたやつが下書きにあったからとりあえず公開。続きを書こうかと思ったけどびっくりするくらい覚えてないのでやめた。
どうも有末です。ご無沙汰しております。
研究関係でシグナルの自己相関関数やフーリエ変換について調べたので、未来の自分のためにメモを残そうと思います。 数学がそんなに得意ではない人間の解釈なのでご注意ください。訂正等のコメントは歓迎いたします。
前提知識というか数式
以下、特に断りがない限り入力として離散的な有限個数のデータを想定しています。
複素共役 (Complex conjugate)
今後式中に $\overline{f}$ の様な記載が出てきますが、これは次式で定義される複素共役 (Complex conjugate) を表しています。
$$
\overline{z} \equiv a - bi \qquad when ~~ z \equiv a + bi
$$
ただ、定義からもわかる通り虚数部がひっくり返っただけなので、実数しか扱わない実数値関数 (Real function) では $\overline{f} = f$ となります。 今回は実験データなどの実数値の解析を行うことを想定しているため定義上複素共役が出てきても、特に何かを行う必要はありません。
numpy では numpy.conjugate(x) として定義されています。
合成積 (Convolution)
今後式中に $f * g$ のような記載が出てきますが、これは次式で定義される合成積 (Convolution) を表しています。
$$
(f * g)(t) \equiv \int_{-\infty}^{\infty} f(\tau) g(t - \tau) d\tau
$$
定義上は連続関数である $f$ および $g$ を無限区間にて積分しますが、実験データなどは無限連続関数ではなく有限な離散値の集合になるため、集合 $a$ および $v$ の定義領域外（集合の領域外）を 0 と再定義し線形畳み込み (linear convolution)
$$
(a * v)(n) \equiv \sum_{m=-\infty}^{\infty} a_{m} v_{n - m}
$$
の計算を行います。なお、上式にて $m$ が集合 $a$ もしくは $v$ の定義域領域外に入った段階で積が 0 になり結果に影響を及ぼさないため事実上無限の計算は不要になります。
numpy では線形畳み込みを行う numpy.convolute(a, v) が定義されています。
この関数では渡された一次配列 a および v の値域がオーバラップする範囲 (N + M - 1 when N := len(a), M := len(v)）の線形畳み込み結果を返します。
期待値 (Expected value)
今後式中に $E[x]$ のような記載が出てきますが、これは次式で定義される期待値 (Expected value) を表しています。
$$
\begin{align}
E[X] &\equiv \int_{-\infty}^{\infty} xf(x)dx \\
&\sim \sum_{i=1}^{\infty} x_{i} p_{i},
\end{align}
$$
上式において $X$ は $f(x)$ を確立密度関数として持つ確率分布、$x_i$ は離散的な各値で $p_i$ はその確率です。
numpy では numpy.mean(x) が定義されています。
共分散 (Covariance)
2種のランダムデータがどの程度シンクロしているかの指標として共分散（Covariance）が以下のように定義されます。
$$ cov(X, Y) \equiv E[(X - E[X])(Y - E[Y])] $$
ここで $X$ および $Y$ は確率変数を表します。
numpy では numpy.cov(m, y) が定義されています。
相関という名前の数式
統計学と信号処理
In probability and statistics, the term cross-correlations is used for referring to the correlations between the entries of two random vectors X and Y, while the autocorrelations of a random vector X are considered to be the correlations between the entries of X itself, those forming the correlation matrix (matrix of correlations) of X. This is analogous to the distinction between autocovariance of a random vector and cross-covariance of two random vectors. One more distinction to point out is that in probability and statistics the definition of correlation always includes a standardising factor in such a way that correlations have values between −1 and +1.
上記は英語版 Wikipedia の Cross-correlation から引用したものですが、調べものをする際は基本的に読み飛ばす癖があるので、この部分を見つけるまでにかなりの時間を要してしまいました。
引用によりますと、自己相関（autocorrelation）や相互相関（cross-correlation）と言った場合に統計学（statistics）および信号処理（signal processing）にて定義に若干の違いが有るようです。 統計学においては「相関（correlation）」と言った場合は常に正規化されており、値は -1 ~ 1 を取るのに対し、信号処理では「相関」と言った場合は処理対象の値をそのまま返し、正規化されたものは特別に「相関係数（correlation coefficient）」と呼ばれるようです。
要約すると正規化の有無により以下のような呼称の違いが有るようです。
| 正規化 | 統計学 (statistics) | 信号処理 (signal processing) |
|---|---|---|
| × | 相互共分散 (cross-covariance) | 相互相関 (cross-correlation) |
| ○ | 相互相関 (cross-correlation) | 相互相関係数 (cross-correlation coefficient) |
| × | 自己共分散 (autocovariance) | 自己相関 (autocorrelation) |
| ○ | 自己相関 (autocorrelation) | 自己相関係数 (autocorrelation coefficient) |
また統計学の定義では「共分散（covariance）」および「相互共分散（cross-covariance）」に数式的な違いが見いだせませんでしたが、これは概念的な違いを明確にするための言葉のようでした（英語版 Wikipedia Cross-covariance より引用）。
In the case of two random vectors X=(X_1, X_2, ... , X_n) and Y=(Y_1, Y_2, ... , Y_n), the cross-covariance would be a square n by n matrix C{XY} with entries C{XY}(j,k) = cov(X_j, Y_k).\, Thus the term cross-covariance is used in order to distinguish this concept from the "covariance" of a random vector X, which is understood to be the matrix of covariances between the scalar components of X itself.
実際の解析時に使用する numpy では信号処理的な定義が用いられているため（関数の返り値が正規化されていない）今回は信号処理の定義および呼称を用いることとします。
相互相関 (Cross-correlation)
相互相関は $f \star g$ と表され、連続関数に対する定義は次式で表されます。
$$
\begin{align}
(f \star g)(t) &\equiv \overline{f}(-t) * g(t) \\
&\equiv \int_{-\infty}^{\infty} \overline{f}(\tau) g(t + \tau) d\tau
\end{align}
$$
合成積の場合と同様に、離散的な集合に対して相互相関関数を考える場合は以下のように無限積を取ります。
$$
(a \star v)(n) \equiv \sum_{m=-\infty}^{\infty} \overline{a}_{m} v_{n + m}
$$
numpy では合成積と同様に離散値に対応したものが numpy.correlate(a, v) として定義されています。
この関数も numpy.convolute(a, v) と同様に、渡された一次配列 a および v の値域がオーバラップする範囲 (N + M - 1 when N := len(a), M := len(v)）の相互相関関数の結果を返します。
自己相関 (Autocorrelation)
定義より自己相関は相互相関を用いて以下のように表すことができます。
$$
R(t) \equiv (f \star f)(t) \equiv \int_{-\infty}^{\infty} \overline{f}(\tau) f(t + \tau) d\tau
$$
また、離散的な集合に対しては
$$
R(n) \equiv (a \star a)(n) \equiv \sum_{m=-\infty}^{\infty} \overline{a}_{m} a_{n + m}
$$
とすることができます。
numpy では対応する関数は存在しませんが、上記の相互相関関数に同じ集合を渡すことで代用できます。
実装する
解析対象の定義
解析対象として 400, 800, 1280 Hz のサイン波を合成したシグナルを利用します。
import numpy as np import matplotlib.pyplot as pl # サンプルデータの定義 ======================================================== N = 250 # Signal length fs = 1000 # Sampling rate fn = fs * 0.5 # Nyquist frequency t = np.arange(N) / fs # 400 Hz sine_400 = np.sin(2. * np.pi * 400 * t) # 800 Hz sine_800 = 2 * np.sin(2. * np.pi * 800 * t) # 1280 Hz sine_1280 = 2 * np.sin(2. * np.pi * 1280 * t) y = sine_400 + sine_800 + sine_1280 pl.plot(t, y)
上記により以下のようなグラフが表示されます。
numpy.correlate を利用した場合
定義的には相互相関関数に対して同じ集合を割り当ててやれば自己相関関数になるはずです。 そこで、上記のサンプルデータを素直に割り当ててみます。
Rr = np.correlate(y, y, mode='full')
pl.plot(Rr)
先にも述べたとおり信号処理的な定義では正規化は行われていません。
そのため numpy.correlate 関数は
- 定義域: (-∞, ∞)
- 値域: (-∞, ∞)
となります。一方 Water survival probability など論文等で個人的によく見る自己相関関数は
- 定義域: (0, ∞)
- 値域: (-1;, 1)
となり（上記 Water survival probability に関して言えば）時間経過による相関の減衰を評価することができます。
解析ではこのようなグラフが欲しいことが多々有りますが、その場合は numpy.correlate() に対し
- データポイントの集合を平均値からの差分集合（偏差集合）にする
- データの散らばり具合（分散）で正規化する（※）
Rr(0)が定義上最大値なので、その値で正規化する
統計学的な定義では期待値（Expected value）が使われているため分散で正規化するだけで値域が -1 ~ 1 となりますが、信号処理的な定義では分散で （ここで日記は途切れている……
