相互相関関数と畳み込み

# はじめに

相互相関関数と畳み込みの関係についての備忘録.
数式も英語も碌に読めない素人の独学なので内容の正確性は保証しない.
例として挙げるコードは C# で書く.
コード中に出てくる Fourier.DFT、Fourier.IDFT は私の前の記事 ([短時間フーリエ変換とその逆変換](https://qiita.com/takayoshi1968/items/636372d464d02699ba4f)) で書いたものである.

相互相関関数や畳み込みは元来連続な関数に対して定義されるものであるが、この記事で扱うのはふたつの配列を受け取り、ひとつの配列を返すものとする(=> 定義域が離散的).
特に断りがない限り、相互相関関数や畳み込みに対する引数は double[] f, g であり、その結果は double[] h である.
f, g, h の長さはそれぞれ N, M, L で表すものとする.

# 相互相関関数 (Cross-Correlation)

名前の通り、ふたつの関数や数列についてその相関(≒類似度)を求めるためのものである.

一般に、ふたつの連続な関数 f(t), g(t) に対し、相互相関 f★g は次のように定義される

```math
(f \star g)(\tau) = \int_{-\infty}^\infty f(t)g(t-\tau) dt
```

この記事では配列 f[n], g[n] を扱いたいので離散的な場合の定義を書くと、

```math
(f \star g)[n] = \sum_{m=-\infty}^{\infty} f[m]g[m-n]
```

このままでは g に対して f が遅れている場合 n が負数になってしまう.
上の定義で $h = (f \star g)[n]$ について、値が 0 でない可能性のある区間は配列 g の長さ M を用いて (-M .. N) であるから、
これを [0 .. N+M-1) にするため以下のように改める.

```math
(f \star g)[n] = \sum_{m=-\infty}^{\infty} f[m]g[m-n+M-1]
```

プログラム上ではこの定義を採用することになるだろう.

## 手計算

数式を見てもピンとこないため例を挙げる.
f = { 1, 7, 5, 3 }, g = { 10, 6, 8, 4, 2 } の場合で手計算を示す.
このとき、(f の長さ =) N = 4, (g の長さ =) M = 5 である.

```text:手計算
h[0] = f[0] * g[4]
     = 1 * 2
     = 2

f     :              1  7  5  3
g << 4: 10  6  8  4  2
                     ~ 重なっている部分だけを使う

h[1] = f[0] * g[3] + f[1] * g[4]
     = 1 * 4 + 7 * 2
     = 18

f     :           1  7  5  3
g << 3: 10  6  8  4  2
                  ~~~~ 重なっている部分だけを使う

h[2] = f[0] * g[2] + f[1] * g[3] + f[2] * g[4]
     = 1 * 8 + 7 * 4 + 5 * 2
     = 46

f     :        1  7  5  3
g << 2: 10  6  8  4  2
               ~~~~~~~ 重なっている部分だけを使う

h[3] = f[0] * g[1] + f[1] * g[2] + f[2] * g[3] + f[3] * g[4]
     = 1 * 6 + 7 * 8 + 5 * 4 + 3 * 2
     = 88

f     :     1  7  5  3
g << 1: 10  6  8  4  2
            ~~~~~~~~~~ 重なっている部分だけを使う

h[4] = f[0] * g[0] + f[1] * g[1] + f[2] * g[2] + f[3] * g[3]
     = 1 * 10 + 7 * 6 + 5 * 8 + 3 * 4
     = 104

f     :  1  7  5  3
g << 0: 10  6  8  4  2
        ~~~~~~~~~~~ 重なっている部分だけを使う

h[5] = f[1] * g[0] + f[2] * g[1] + f[3] * g[2]
     = 7 * 10 + 5 * 6 + 3 * 8
     = 124

f     :  1  7  5  3
g >> 1:    10  6  8  4  2
           ~~~~~~~~ 重なっている部分だけを使う

h[6] = f[2] * g[0] + f[3] * g[1]
     = 5 * 10 + 3 * 6
     = 68

f     :  1  7  5  3
g >> 2:       10  6  8  4  2
              ~~~~~ 重なっている部分だけを使う

h[7] = f[3] * g[0]
     = 3 * 10
     = 30

f     :  1  7  5  3
g >> 3:          10  6  8  4  2
                 ~~ 重なっている部分だけを使う

以上から、h = { 2, 18, 46, 88, 104, 124, 68, 30 }, L = 8 (= N + M - 1)
```

## コード例

```C#:コード例
// h[0]  : f[0] と g[0] の一点のみが重なっている
// h[M-1]: f[0] と g[M-1] が重なっている
// h[N-1]: f[N-1] と g[0] が重なっている
// h[L-1]: f[N-1] と g[M-1] の一点のみが重なっている
public static double[] CrossCorrelation_Naive(double[] f, double[] g) {
  int N = f.Length;
  int M = g.Length;

  int L = N + M - 1;
  double[] h = new double[L];

  // n, m はそれぞれ i, j になっているが定義通り.
  for (int i = 0; i < L; ++i) {
    double val = 0;
    for (int j = 0; j < N; ++j) {
      int k = j - i + M - 1;
      if (0 <= k && k < M) {
        val += f[j] * g[k];
      }
    }// loop for j in [0 .. N)
    h[i] = val;
  }// loop for i in [0 .. L)

  return h;
}// CrossCorrelation_Naive(double[] f, double[] g)
```

## 可換性

相互相関関数は非可換である. 以下に図を示す.
上の段がオリジナルの配列、2、3段目が計算中の配列、最下段が計算結果である.
2、3段目の四角形は配列が重なった部分であり、最下段の縦線はその瞬間におけるカーソルである.
横線は y = 0 を表す.

$f \star g$

![cross-correlation_fg.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2329888/8d7521a9-a7b7-71fb-2dcf-5aa449925714.gif)

$g \star f$

![cross-correlation_gf.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2329888/d4a43fa2-816e-43f8-81ad-0cd15004aaa4.gif)

これらの結果はお互いの左右反転らしい.

# 線型畳み込み (Linear Convolution)

単に畳み込みと言った場合はおそらくこれを指すのだろう.
この記事では後述する循環畳み込みとの区別のため、頭に線型 (Linear) とつけている.

畳み込み演算について調べようとすると機械学習関連のページが上に出てきて2次元で説明していたり、相互相関関数と混同していたり、循環畳み込みの説明であったりと、極めて歯がゆい思いをさせられる.

一般に、ふたつの連続な関数 f(t), g(t) に対し、線型畳み込み f*g は次のように定義される

```math
(f \ast g)(\tau) = \int_{-\infty}^\infty f(t)g(\tau-t) dt
```

相互相関関数とよく似ているが、g の引数が異なる.

この記事では配列 f[n], g[n] を扱いたいので離散的な場合の定義を書くと、

```math
(f \ast g)[n] = \sum_{m=-\infty}^{\infty} f[m]g[n-m]
```

こちらは相互相関関数と異なり、このまま使っても定義域の問題がない.

## 手計算

相互相関関数の場合と同様に例を挙げる.
f = { 1, 7, 5, 3 }, g = { 10, 6, 8, 4, 2 } の場合で手計算を示す.
このとき、(f の長さ =) N = 4, (g の長さ =) M = 5 である.

自己相関関数で書いたものと同様に書くとややこしく見えるかもしれないが、畳み込みは f, g が正面衝突してすれ違うまでを表したものであると考えれば、こちらの方がむしろ直感的かもしれない.

reverse(xs) は配列を左右反転させるものであるとする.

```text:手計算
h[0] = f[0] * g[0]
     = 1 * 10
     = 10

        f      :             1  7  5  3
reverse(g) << 4: 2  4  8  6 10
                            ~~ 重なっている部分だけを使う

h[1] = f[0] * g[1] + f[1] * g[0]
     = 1 * 6 + 7 * 10
     = 76

        f      :          1  7  5  3
reverse(g) << 3: 2  4  8  6 10
                          ~~~~ 重なっている部分だけを使う

h[2] = f[0] * g[2] + f[1] * g[1] + f[2] * g[0]
     = 1 * 8 + 7 * 6 + 5 * 10
     = 100

        f      :       1  7  5  3
reverse(g) << 2: 2  4  8  6 10
                       ~~~~~~~ 重なっている部分だけを使う

h[3] = f[0] * g[3] + f[1] * g[2] + f[2] * g[1] + f[3] * g[0]
     = 1 * 4 + 7 * 8 + 5 * 6 + 3 * 10
     = 120

        f      :    1  7  5  3
reverse(g) << 1: 2  4  8  6 10
                    ~~~~~~~~~~ 重なっている部分だけを使う

h[4] = f[0] * g[4] + f[1] * g[3] + f[2] * g[2] + f[3] * g[1]
     = 1 * 2 + 7 * 4 + 5 * 8 + 3 * 6
     = 88

        f      : 1  7  5  3
reverse(g) << 0: 2  4  8  6 10
                 ~~~~~~~~~~ 重なっている部分だけを使う

h[5] = f[1] * g[4] + f[2] * g[3] + f[3] * g[2]
     = 7 * 2 + 5 * 4 + 3 * 8
     = 58

        f      : 1  7  5  3
reverse(g) >> 1:    2  4  8  6 10
                    ~~~~~~~ 重なっている部分だけを使う

h[6] = f[2] * g[4] + f[3] * g[3]
     = 5 * 2 + 3 * 4
     = 22

        f      : 1  7  5  3
reverse(g) >> 2:       2  4  8  6 10
                       ~~~~ 重なっている部分だけを使う

h[7] = f[3] * g[4]
     = 3 * 2
     = 6

        f      : 1  7  5  3
reverse(g) >> 2:          2  4  8  6 10
                         ~~ 重なっている部分だけを使う

以上から、h = { 10, 76, 100, 120, 88, 58, 22, 6 }, L = 8 (= N + M - 1)
```

## コード例

```C#:コード例
// h[0]  : f[0] と g[0] の一点のみが重なっている
// h[M-1]: f[0] と g[M-1] が重なっている
// h[N-1]: f[N-1] と g[0] が重なっている
// h[L-1]: f[N-1] と g[M-1] の一点のみが重なっている
public static double[] CrossCorrelation_Naive(double[] f, double[] g) {
  int N = f.Length;
  int M = g.Length;

  int L = N + M - 1;
  double[] h = new double[L];

  // n, m はそれぞれ i, j になっているが定義通り.
  for (int i = 0; i < L; ++i) {
    double val = 0;
    for (int j = 0; j < N; ++j) {
      int k = j - i + M - 1;
      if (0 <= k && k < M) {
        val += f[j] * g[k];
      }
    }// loop for j in [0 .. N)
    h[i] = val;
  }// loop for i in [0 .. L)

  return h;
}// CrossCorrelation_Naive(double[] f, double[] g)
```

## 可換性

畳み込みは可換である. 以下に図を示す.
上の段がオリジナルの配列、2、3段目が計算中の配列、最下段が計算結果である.
2、3段目の四角形は配列が重なった部分であり、最下段の縦線はその瞬間におけるカーソルである.
横線は y = 0 を表す.

$f \ast g$

![linear-convolution_fg.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2329888/15654083-ff9d-e3b7-d4c7-1727c49011f5.gif)

$g \ast f$

![linear-convolution_gf.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2329888/6cac1ae9-a811-934d-d469-21fa3b89fcc9.gif)

# 線型畳み込みと相互相関関数の関係

手計算の例も挙げたのであえて書かずとも良さそうな気はするが、一応明示的に書いておく.
配列の左右を反転させる関数 Reverse が存在するとき、以下のコードは正しい.

```C#:コード例
public static CrossCorrelation_WithLinearConvolution(double[] f, double[] g) {
  double[] rev_g = Reverse(g);
  return LinearConvolution(f, rev_g);

  // または以下:
  // double[] rev_f = Reverse(f);
  // return LinearConvolution(rev_f, g);
}// CrossCorrelation_WithLinearConvolution(double[] f, double[] g)

public static LinearConvolution_WithCrossCorrelation(double[] f, double[] g) {
  double[] rev_g = Reverse(g);
  return CrossCorrelation(f, rev_g);
}// LinearConvolution_WithCrossCorrelation(double[] f, double[] g)

private double[] Reverse(double[] xs) {
  int N = xs.Length;
  double[] ys = new double[N];
  for (int i = N / 2; i < N; ++i) {
    int k = N - 1 - i;
    ys[k] = xs[i];
    ys[i] = xs[k];
  }
  return ys;
}
```

# 循環畳み込み (Circular Convolution)

畳み込みの派生であり、同じ周期 T を持つ関数 f, g の間に定義される.

```math
(f \circledast g)(\tau) = \int_{t_0}^{t_0+T} f(t)g(\tau-t) dt
```

ここで、t_0 は任意の点である.

この記事では配列 f[n], g[n] を扱いたいので離散的な場合の定義を書くと、

```math
(f \circledast g)[n] = \sum_{m=-\infty}^{\infty} f[m]g[n-m]
```

配列の長さは有限であるから、f の長さ = g の長さ = T として以下のように改める.

```math
(f \circledast g)[n] = \sum_{m=0}^{T-1} f[m]g[(n-m)\ \ mod\ \ T]
```

相互相関関数の場合と同様に例を挙げる.
f = { 1, 7, 5, 3, 9 }, g = { 10, 6, 8, 4, 2 } の場合で手計算を示す.
このとき、f の長さ = g の長さ = T = 5 である.

rev_rep(xs) は配列を左右反転させた上でそれを繰り返すものであるとする.
| は繰り返しの区切りである.

```text:手計算
h[0] = f[0] * g[0] + f[1] * g[4] + f[2] * g[3] + f[3] * g[2] + f[4] * g[1]
     = 1 * 10 + 7 * 2 + 5 * 4 + 3 * 8 + 9 * 6
     = 122

        f      :             1   7  5  3  9
rev_rep(g) << 4: 2  4  8  6 10 | 2  4  8  6 10 
                            ~~~~~~~~~~~~~~~ 重なっている部分だけを使う

h[1] = f[0] * g[1] + f[1] * g[0] + f[2] * g[4] + f[3] * g[3] + f[4] * g[2]
     = 1 * 6 + 7 * 10 + 5 * 2 + 3 * 4 + 9 * 8
     = 170

        f      :          1  7   5  3  9
rev_rep(g) << 3: 2  4  8  6 10 | 2  4  8  6 10
                          ~~~~~~~~~~~~~~ 重なっている部分だけを使う

h[2] = f[0] * g[2] + f[1] * g[1] + f[2] * g[0] + f[3] * g[4] + f[4] * g[3]
     = 1 * 8 + 7 * 6 + 5 * 10 + 3 * 2 + 9 * 4
     = 142

        f      :       1  7  5   3  9
rev_rep(g) << 2: 2  4  8  6 10 | 2  4  8  6 10
                       ~~~~~~~~~~~~~~ 重なっている部分だけを使う

h[3] = f[0] * g[3] + f[1] * g[2] + f[2] * g[1] + f[3] * g[0] + f[4] * g[4]
     = 1 * 4 + 7 * 8 + 5 * 6 + 3 * 10 + 9 * 2
     = 138

        f      :    1  7  5  3   9
rev_rep(g) << 1: 2  4  8  6 10 | 2  4  8  6 10
                    ~~~~~~~~~~~~~~ 重なっている部分だけを使う

h[4] = f[0] * g[4] + f[1] * g[3] + f[2] * g[2] + f[3] * g[1] + f[4] * g[0]
     = 1 * 2 + 7 * 4 + 5 * 8 + 3 * 6 + 9 * 10
     = 178

        f      : 1  7  5  3  9
rev_rep(g) << 0: 2  4  8  6 10
                 ~~~~~~~~~~~~~ 重なっている部分だけを使う

以上から、h = { 122, 170, 142, 138, 178 }, L = 5 (= T)
```

## コード例

```C#:コード例
public static double[] CircularConvolution_Naive(double[] f, double[] g) {
  int T = f.Length;
  if (T != g.Length) throw new ArgumentException("f and g must have the same length");

  double[] h = new double[T];

  // n, m はそれぞれ i, j になっているが定義通り.
  for (int i = 0; i < T; ++i) {
    double val = 0;
    for (int j = 0; j < T; ++j) {
      int k = (i - j + T) % T;
      val += f[j] * g[k];
    }// loop for j in [0 .. T)
    h[i] = val;
  }// loop for i in [0 .. T)

  return h;
}// CircularConvolution_Naive(double[] f, double[] g)
```

## 可換性

循環畳み込みは可換である. 以下に図を示す.
上の段がオリジナルの配列、2、3段目が計算中の配列、最下段が計算結果である.
2、3段目の黒い四角形は配列が重なった部分であり、3段目の赤い四角形は配列ひとつ分であり、
また、最下段の縦線はその瞬間におけるカーソルである.
横線は y = 0 を表す.

$f \circledast g$

![circular-convolution_fg.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2329888/6f7b7d26-b269-3093-2b53-f1d209631cd5.gif)

$g \circledast f$

![circular-convolution_gf.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2329888/8f11de44-3bb4-7246-5159-f21a8ed7bb2b.gif)

# 循環畳み込みの高速化

これまでに見てきた素朴な実装では、相互相関関数、線型畳み込み、循環畳み込みいずれも2重ループになっていた.
データ数が多い場合、O(n<sup>2</sup>) では遅く現実的ではないこともあるだろう.

ところで、連続な関数 f, g に対する線型畳み込みについて、畳み込みの定理というものが知られている.

```math
\mathcal{F}[f \ast g] = \mathcal{F}[f] \cdot \mathcal{F}[g] \\
\text{ここで、} \mathcal{F}[f] \text{は} f \text{のフーリエ変換を表す.}
```

これの離散フーリエ変換版は以下のようになる.

```math
\mathcal{F}[f \circledast g] = DFT[f] \cdot DFT[g] \\
\text{ここで、} DFT[f] \text{は} f \text{の離散フーリエ変換を表す.}
```

このことから、循環畳み込みは2重ループを使わずに実装でき、そのオーダーは DFT のオーダーとなる.
DFT は Cooley-Tukey の高速フーリエ変換によって O(N logN) で計算可能であるから、これを用いれば十分高速となる.

## コード例

```C#:コード例
public static double[] CircularConvolution_FFT(double[] f, double[] g) {
  int T = f.Length;
  if (T != g.Length) throw new ArgumentException("f and g must have the same length");

  // 1. f, g を Complex 型に変換する
  Complex[] f_complex = new Complex[T];
  Complex[] g_complex = new Complex[T];
  for (int i = 0; i < T; ++i) f_complex[i] = f[i];
  for (int i = 0; i < T; ++i) g_complex[i] = g[i];

  // 2. FFT して、
  Complex[] f_spectrum = Fourier.FFT(f_complex);
  Complex[] g_spectrum = Fourier.FFT(g_complex);

  // 3. 得られたスペクトラムを掛け合わせて、
  Complex[] h_spectrum = new Complex[T];
  for (int i = 0; i < T; ++i) {
    h_spectrum[i] = f_spectrum[i] * g_spectrum[i];
  }

  // 4. iFFT する
  Complex[] h_complex = Fourier.IFFT(h_spectrum);

  // 5. 実数部分だけが必要なので取り出す
  double[] h = new double[T];
  for (int i = 0; i < T; ++i) {
    h[i] = h_complex[i].Real;
  }

  return h;
}// CircularConvolution_FFT(double[] f, double[] g)
```

# 線型畳み込みの高速化

循環畳み込みの計算量を O(N logN) にすることができたので、相互相関関数、線型畳み込みも循環畳み込みに帰着させることができれば高速化できるといえる.
同じ畳み込み仲間ということで、まずは線型畳み込みを循環畳み込みに帰着させる場合を示す.

```C#:コード例
public static double[] LinearConvolution_WithCircularConvolution(double[] f, double[] g) {
  int N = f.Length;
  int M = g.Length;

  int L = N + M - 1;

  // 1. f, g を長さ L になるようゼロパディングする(左パディング)
  double[] f_padded = new double[L];
  double[] g_padded = new double[L];
  for (int i = 0; i < N; ++i) f_padded[i] = f[i];
  for (int i = 0; i < M; ++i) g_padded[i] = g[i];

  // 2. 循環畳み込みを行う
  double[] h = CircularConvolution(f_padded, g_padded);

  return h;
}// LinearConvolution_WithCircularConvolution(double[] f, double[] g)
```

これをアニメーションにすると以下のようになる.

![linear-convolution_fg_as_circular.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2329888/6e4d8108-d27e-e2d1-2935-9fcef5c27c85.gif)

# 相互相関関数の高速化

相互相関関数を線型畳み込みに変換する方法を上の方(線型畳み込みと相互相関関数の関係)で示し、先ほど線型畳み込みを循環畳み込みに変換する方法を示した.
この2段階の変換によって相互相関関数を高速化できるだろう.

一方、次に示すコードも O(N logN) で相互相関関数を計算可能である.

```C#:コード例
public static double[] CrossCorrelation_DFT(double[] f, double[] g) {
  // 相互相関関数なので、線型畳み込みとはパディングの仕方と、一部計算に共役を用いる点で異なる
  // cf. http://frchick.blog129.fc2.com/blog-entry-717.html

  int N = f.Length;
  int M = g.Length;

  int L = N + M - 1;

  // 1. f, g を長さ L になるようゼロパディングする(f は左パディング、g は右パディング)
  Complex[] f_padded = new Complex[L];
  Complex[] g_padded = new Complex[L];
  for (int i = 0; i < N; ++i) f_padded[i + M - 1] = f[i];// 畳み込みとの違いそのいち
  for (int i = 0; i < M; ++i) g_padded[i] = g[i];

  // 2. DFT して、
  Complex[] f_spectrum = Fourier.DFT(f_padded);
  Complex[] g_spectrum = Fourier.DFT(g_padded);

  // 3. 得られたスペクトラムを、g 側は共役にした状態で掛け合わせて、
  Complex[] h_spectrum = new Complex[L];
  for (int i = 0; i < L; ++i) {
    h_spectrum[i] = f_spectrum[i] * Complex.Conjugate(g_spectrum[i]);// 畳み込みとの違いそのに
  }

  // 4. iDFT する
  Complex[] h_complex = Fourier.IDFT(h_spectrum);

  // 5. 実数部分だけが必要なので取り出す
  double[] h = new double[L];
  for (int i = 0; i < L; ++i) {
    h[i] = h_complex[i].Real;
  }

  return h;
}// CrossCorrelation_DFT(double[] f, double[] g)
```

数学的背景はよくわからないが、とにかくこのコードは機能しているように見える.
詳細は参考文献の[フラチキさんブログ 相互相関関数のFFTによる高速化](http://frchick.blog129.fc2.com/blog-entry-717.html)、およびそのページ内にリンクがある[法政大学、小林先生の資料](https://www.ikko.k.hosei.ac.jp/~matlab/xcorr.pdf)を参照.

# 参考文献

- Wikipedia
  日本語の相互相関関数のページは情報量が少なすぎるので英語版
  - [相互相関関数(英語)](https://en.wikipedia.org/wiki/Cross-correlation)
  - [畳み込み(英語)](https://en.wikipedia.org/wiki/Convolution)
  - [畳み込みの定理(英語)](https://en.wikipedia.org/wiki/Convolution_theorem)
    比較図やアニメーションがあってよい

- フラチキさんブログ
  [相互相関関数のFFTによる高速化](http://frchick.blog129.fc2.com/blog-entry-717.html)

- Qiita
  - [【改善版】C#でgifアニメを作る - ルパン三世のタイトルコール風](https://qiita.com/kob58im/items/d76ef98de48a5a17165e)
    記事の内容とは関係ないが、GIF画像を作るに際してコピペさせていただいた