はじめに
相互相関関数を計算するときに,時間領域で直接計算(直接法)する方法とフーリエ変換を利用して周波数領域で計算する方法(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法は直接法に比べて処理速度が速かったです.ただ,遅れ時間がそれほど大きくないことがわかっているときは,直接法で計算範囲を絞ったほうが速いと思います.