相互相関関数の計算方法の違いによる処理速度の差

8076 ワード
Python scipy numpy Python テキストリンク

はじめに

相互相関関数を計算するときに，時間領域で直接計算（直接法）する方法とフーリエ変換を利用して周波数領域で計算する方法（FFT法）があります．計算時間はFFT法のほうが速いと聞きましたが，Numpyのcorrelate関数は直接法でしか計算できません．そこで，勉強を兼ねてFFT法を実装して試してみようと思います．
作り終わってから気が付きましたが，Scipyのcorrelate関数はFFT法が使えます．

参考

こちらの資料を参考にさせていただきました．
コンボリューション関数（conv）と相関関数(xcorr)と離散フーリエ変換関数（fft）の関係

プログラム

ここからプログラム本体です．Juputer notebookで作成しました．

1. import

#coding:utf-8
%matplotlib notebook
import matplotlib.pyplot as plt
import numpy as np
import scipy.signal as sg

2. 試験信号

ランダムな信号を作成します．

N = 512
x, y = np.random.randn(2, N)
#正規化
x -= np.average(x)
x /= np.std(x)
y -= np.average(y)
y /= np.std(y)

3. 直接法

corr = np.correlate(x, y, mode="full") / N

4. FFT法 巡回問題対策済み

def xcorr(x, y):
    N = len(x) if len(x) > len(y) else len(y)
    #xとyの長さが同じになるように0を加える
    xz = np.append(x, np.zeros(len(y)-1))
    yz = np.append(np.zeros(len(x)-1), y)

    XY = np.fft.fft(xz) * np.conj(np.fft.fft(yz))
    xy = np.real(np.fft.ifft(XY)) / N
    return xy

5. 直接法とFFT法の計算結果の比較

FFT法では巡回問題を解決するために0を加えたため，インデックスが1つずれます．
1つずらして差分を取ると，両者の値は一致します．

diff = np.real(xy)[1:N*2-1]-np.append(0.0,corr)[1:N*2-1]

処理速度の比較

直接法とFFT法の処理速度を比較します．処理対象のデータ点数が少ないときはNumpyの関数が速いですが，点数が増えてくるとFFT法の利点が顕著に見れます．

まとめ

FFT法は直接法に比べて処理速度が速かったです．ただ，遅れ時間がそれほど大きくないことがわかっているときは，直接法で計算範囲を絞ったほうが速いと思います．