離散フーリエ変換(4) - 2次元DFT

今まで１次元のベクトルのみをあつかってきましたが、

相互相関関数で類似度を計算できる => 相互相関関数は畳み込みの定理でDFTして掛け合わせて逆DFTしたのとおんなじ

っていう話が２次元にも拡張できます。

今は素直な実装のDFTを使用していますが、あとでやる予定のFFTを使うとDFT、逆DFTの計算がはやく行えるので相互相関関数の計算がはやく出来ます。

まず、２次元の相互相関関数は以下です。
２つの $S \times T$ の行列FとGの相互相関関数は、

$C_{p, q} = \sum_{s = 0}^{S - 1} \sum_{t = 0}^{T - 1} F_{s, t} G_{s + p, t + q}$

となります。

$S \times T$ の行列のDFTは

$A_{p, q} = \sum_{s=0}^{S-1} \sum_{t=0}^{T-1}a_{s,t} \exp\{ \frac{-2\pi i p s}{S} \} \exp \{\frac{-2\pi i q t}{T} \}$

となります。

下のように変形すると分かりやすいですが、これは各行をDFTしたあと、その結果の行列の各列をDFTしたことと同じになります。

$A_{p, q} = \sum_{s=0}^{S-1} \left( \sum_{t=0}^{T-1}a_{s,t} \exp \{\frac{-2\pi i q t}{T} \} \right) \exp\{ \frac{-2\pi i p s}{S} \}$

これを踏まえて両者の結果が一致するのか確かめてみたいと思います。

#include <iostream>
#include <vector>
#include <complex>
#include <cmath>

using namespace std;

static const double PI = 3.14159265358979;

vector< complex<double> > dft(const vector< complex<double> > &x) {
  vector< complex<double> > result;
  for(size_t j = 0; j < x.size(); ++j) {
    complex<double> F_j(0.0, 0.0);
    for(size_t k = 0; k < x.size(); ++k) {
      double angle = -2.0 * PI * j * k / x.size();
      F_j += x.at(k) * complex<double>(cos(angle), sin(angle));
    }
    result.push_back(F_j);
  }

  return result;
}

void dft2d(const complex<double> orig[30][30], complex<double> result[30][30]) {
  for(int i = 0; i < 30; ++i) {
    vector< complex<double> > temp;
    for(int j = 0; j < 30; ++j)
      temp.push_back(orig[i][j]);
    
    vector< complex<double> > line = dft(temp);
    for(int j = 0; j < 30; ++j){
      result[i][j] = line.at(j);
    }
  }

  for(int i = 0; i < 30; ++i) {
    vector< complex<double> > temp;
    for(int j = 0; j < 30; ++j)
      temp.push_back(result[j][i]);

    vector< complex<double> > line = dft(temp);
    for(int j = 0; j < 30; ++j) 
      result[j][i] = line.at(j);
  }
}

int main() {
  complex<double> A[30][30];
  complex<double> B[10][10];
  for(int i = 0; i < 30; ++i) {
    for(int j = 0; j < 30; ++j) {
      A[i][j] = complex<double>(2, 0);
      if( 10 <= i && i < 20 && j >= 10 && j < 20 ) {
	A[i][j] = complex<double>((i-9) * (j-9), 0);
	B[i - 10][j - 10] = complex<double>((i-9) * (j-9), 0);
      }
    }
  }

  // 素直に相互相関を計算
  complex<double> result[20][20];
  for(int i = 0; i < 20; ++i) {
    for(int j = 0; j < 20; ++j) {
      complex<double> temp(0, 0);
      for(int k = 0; k < 10; ++k) {
	for(int l = 0; l < 10; ++l) {
	  temp += A[i + k][j + l] * B[k][l];
	}
      }
      result[i][j] = temp;
      cout << result[i][j].real() / 1.0e+05<< " ";
    }
    cout << endl;
  }

  complex<double> C[30][30];
  for(int i = 0; i < 30; ++i) {
    for(int j = 0; j < 30; ++j) {
      C[i][j] = complex<double>(0, 0);
      if( 0 <= i && i < 10 && j >= 0 && j < 10 ) {
	C[i][j] = complex<double>((i+1) * (j+1), 0);
      }
    }
  }

  complex<double> AF[30][30];
  complex<double> CF[30][30];
  complex<double> AC[30][30];
  complex<double> dft_result[30][30];
  // DFTしてから
  dft2d(A, AF);
  dft2d(C, CF);
  
  for(int i = 0; i < 30; ++i) {
    for(int j = 0; j < 30; ++j) {
      AC[i][j] = conj(AF[i][j]) * CF[i][j];
    }
  }
    
  // 逆変換
  dft2d(AC, dft_result);
  cout << endl << endl;
  for(int i = 0; i < 20; ++i) {
    for(int j = 0; j < 20; ++j) {
      cout << dft_result[i][j].real() / (1.0e+05 * 30.0 * 30.0) << " ";
    }
    cout << endl;
  }

  return 0;
}

結果

0.0605 0.0605 0.0605 0.0605 0.0605 0.0605 0.0605 0.0605 0.0605 0.0605 0.0605 0.0605 0.0605 0.0605 0.0605 0.0605 0.0605 0.0605 0.0605 0.0605
0.0605 0.0595 0.0596 0.0607 0.0627 0.0655 0.069 0.0731 0.0777 0.0827 0.088 0.0845 0.0809 0.0773 0.0738 0.0705 0.0675 0.0649 0.0628 0.0613
0.0605 0.0596 0.06169 0.06648 0.07368 0.083 0.09415 0.10684 0.12078 0.13568 0.15125 0.1391 0.12686 0.11482 0.10327 0.0925 0.0828 0.07446 0.06777 0.06302
0.0605 0.0607 0.06648 0.07728 0.09254 0.1117 0.1342 0.15948 0.18698 0.21614 0.2464 0.221 0.19562 0.17082 0.14716 0.1252 0.1055 0.08862 0.07512 0.06556
0.0605 0.0627 0.07368 0.09254 0.11838 0.1503 0.1874 0.22878 0.27354 0.32078 0.3696 0.3269 0.28442 0.24306 0.20372 0.1673 0.1347 0.10682 0.08456 0.06882
0.0605 0.0655 0.083 0.1117 0.1503 0.1975 0.252 0.3125 0.3777 0.4463 0.517 0.4535 0.3905 0.3293 0.2712 0.2175 0.1695 0.1285 0.0958 0.0727
0.0605 0.069 0.09415 0.1342 0.1874 0.252 0.32625 0.4084 0.4967 0.5894 0.68475 0.5975 0.5111 0.4273 0.34785 0.2745 0.209 0.1531 0.10855 0.0771
0.0605 0.0731 0.10684 0.15948 0.22878 0.3125 0.4084 0.51424 0.62778 0.74678 0.869 0.7556 0.64346 0.53482 0.43192 0.337 0.2523 0.18006 0.12252 0.08192
0.0605 0.0777 0.12078 0.18698 0.27354 0.3777 0.4967 0.62778 0.76818 0.91514 1.0659 0.9245 0.78482 0.64962 0.52166 0.4037 0.2985 0.20882 0.13742 0.08706
0.0605 0.0827 0.13568 0.21614 0.32078 0.4463 0.5894 0.74678 0.91514 1.09118 1.2716 1.1009 0.93242 0.76946 0.61532 0.4733 0.3467 0.23882 0.15296 0.09242
0.0605 0.088 0.15125 0.2464 0.3696 0.517 0.68475 0.869 1.0659 1.2716 1.48225 1.2815 1.0835 0.8921 0.71115 0.5445 0.396 0.2695 0.16885 0.0979
0.0605 0.0845 0.1391 0.221 0.3269 0.4535 0.5975 0.7556 0.9245 1.1009 1.2815 1.109 0.9389 0.7745 0.6191 0.476 0.3485 0.2399 0.1535 0.0926
0.0605 0.0809 0.12686 0.19562 0.28442 0.3905 0.5111 0.64346 0.78482 0.93242 1.0835 0.9389 0.79634 0.65858 0.52838 0.4085 0.3017 0.21074 0.13838 0.08738
0.0605 0.0773 0.11482 0.17082 0.24306 0.3293 0.4273 0.53482 0.64962 0.76946 0.8921 0.7745 0.65858 0.54658 0.44074 0.3433 0.2565 0.18258 0.12378 0.08234
0.0605 0.0738 0.10327 0.14716 0.20372 0.2712 0.34785 0.43192 0.52166 0.61532 0.71115 0.6191 0.52838 0.44074 0.35793 0.2817 0.2138 0.15598 0.10999 0.07758
0.0605 0.0705 0.0925 0.1252 0.1673 0.2175 0.2745 0.337 0.4037 0.4733 0.5445 0.476 0.4085 0.3433 0.2817 0.225 0.1745 0.1315 0.0973 0.0732
0.0605 0.0675 0.0828 0.1055 0.1347 0.1695 0.209 0.2523 0.2985 0.3467 0.396 0.3485 0.3017 0.2565 0.2138 0.1745 0.1395 0.1097 0.086 0.0693
0.0605 0.0649 0.07446 0.08862 0.10682 0.1285 0.1531 0.18006 0.20882 0.23882 0.2695 0.2399 0.21074 0.18258 0.15598 0.1315 0.1097 0.09114 0.07638 0.06598
0.0605 0.0628 0.06777 0.07512 0.08456 0.0958 0.10855 0.12252 0.13742 0.15296 0.16885 0.1535 0.13838 0.12378 0.10999 0.0973 0.086 0.07638 0.06873 0.06334
0.0605 0.0613 0.06302 0.06556 0.06882 0.0727 0.0771 0.08192 0.08706 0.09242 0.0979 0.0926 0.08738 0.08234 0.07758 0.0732 0.0693 0.06598 0.06334 0.06148


0.0605 0.0605 0.0605 0.0605 0.0605 0.0605 0.0605 0.0605 0.0605 0.0605 0.0605 0.0605 0.0605 0.0605 0.0605 0.0605 0.0605 0.0605 0.0605 0.0605
0.0605 0.0595 0.0596 0.0607 0.0627 0.0655 0.069 0.0731 0.0777 0.0827 0.088 0.0845 0.0809 0.0773 0.0738 0.0705 0.0675 0.0649 0.0628 0.0613
0.0605 0.0596 0.06169 0.06648 0.07368 0.083 0.09415 0.10684 0.12078 0.13568 0.15125 0.1391 0.12686 0.11482 0.10327 0.0925 0.0828 0.07446 0.06777 0.06302
0.0605 0.0607 0.06648 0.07728 0.09254 0.1117 0.1342 0.15948 0.18698 0.21614 0.2464 0.221 0.19562 0.17082 0.14716 0.1252 0.1055 0.08862 0.07512 0.06556
0.0605 0.0627 0.07368 0.09254 0.11838 0.1503 0.1874 0.22878 0.27354 0.32078 0.3696 0.3269 0.28442 0.24306 0.20372 0.1673 0.1347 0.10682 0.08456 0.06882
0.0605 0.0655 0.083 0.1117 0.1503 0.1975 0.252 0.3125 0.3777 0.4463 0.517 0.4535 0.3905 0.3293 0.2712 0.2175 0.1695 0.1285 0.0958 0.0727
0.0605 0.069 0.09415 0.1342 0.1874 0.252 0.32625 0.4084 0.4967 0.5894 0.68475 0.5975 0.5111 0.4273 0.34785 0.2745 0.209 0.1531 0.10855 0.0771
0.0605 0.0731 0.10684 0.15948 0.22878 0.3125 0.4084 0.51424 0.62778 0.74678 0.869 0.7556 0.64346 0.53482 0.43192 0.337 0.2523 0.18006 0.12252 0.08192
0.0605 0.0777 0.12078 0.18698 0.27354 0.3777 0.4967 0.62778 0.76818 0.91514 1.0659 0.9245 0.78482 0.64962 0.52166 0.4037 0.2985 0.20882 0.13742 0.08706
0.0605 0.0827 0.13568 0.21614 0.32078 0.4463 0.5894 0.74678 0.91514 1.09118 1.2716 1.1009 0.93242 0.76946 0.61532 0.4733 0.3467 0.23882 0.15296 0.09242
0.0605 0.088 0.15125 0.2464 0.3696 0.517 0.68475 0.869 1.0659 1.2716 1.48225 1.2815 1.0835 0.8921 0.71115 0.5445 0.396 0.2695 0.16885 0.0979
0.0605 0.0845 0.1391 0.221 0.3269 0.4535 0.5975 0.7556 0.9245 1.1009 1.2815 1.109 0.9389 0.7745 0.6191 0.476 0.3485 0.2399 0.1535 0.0926
0.0605 0.0809 0.12686 0.19562 0.28442 0.3905 0.5111 0.64346 0.78482 0.93242 1.0835 0.9389 0.79634 0.65858 0.52838 0.4085 0.3017 0.21074 0.13838 0.08738
0.0605 0.0773 0.11482 0.17082 0.24306 0.3293 0.4273 0.53482 0.64962 0.76946 0.8921 0.7745 0.65858 0.54658 0.44074 0.3433 0.2565 0.18258 0.12378 0.08234
0.0605 0.0738 0.10327 0.14716 0.20372 0.2712 0.34785 0.43192 0.52166 0.61532 0.71115 0.6191 0.52838 0.44074 0.35793 0.2817 0.2138 0.15598 0.10999 0.07758
0.0605 0.0705 0.0925 0.1252 0.1673 0.2175 0.2745 0.337 0.4037 0.4733 0.5445 0.476 0.4085 0.3433 0.2817 0.225 0.1745 0.1315 0.0973 0.0732
0.0605 0.0675 0.0828 0.1055 0.1347 0.1695 0.209 0.2523 0.2985 0.3467 0.396 0.3485 0.3017 0.2565 0.2138 0.1745 0.1395 0.1097 0.086 0.0693
0.0605 0.0649 0.07446 0.08862 0.10682 0.1285 0.1531 0.18006 0.20882 0.23882 0.2695 0.2399 0.21074 0.18258 0.15598 0.1315 0.1097 0.09114 0.07638 0.06598
0.0605 0.0628 0.06777 0.07512 0.08456 0.0958 0.10855 0.12252 0.13742 0.15296 0.16885 0.1535 0.13838 0.12378 0.10999 0.0973 0.086 0.07638 0.06873 0.06334
0.0605 0.0613 0.06302 0.06556 0.06882 0.0727 0.0771 0.08192 0.08706 0.09242 0.0979 0.0926 0.08738 0.08234 0.07758 0.0732 0.0693 0.06598 0.06334 0.06148

完全に一致！

次回はDFTして掛け合わせて逆DFTするとなんで相互相関と同じになるのかを簡単に証明したいと思います