相互相関関数の計算

西田究のホームページ (u-tokyo.ac.jp)
解説 (u-tokyo.ac.jp)

地震動に関する様々なシミュレーション結果をお手軽に見ることができる、素晴らしいサイトです。特に、地震波干渉法で相互相関関数が出来上がっていく過程を見たことがなかったので、非常に興味深く見ていました。

　　初期

数十分後

相互相関関数を求める際は、時間領域（直接法）か周波数領域（FFT）を利用します。配列の要素数が多ければ周波数領域での計算が早いでしょう。Python の scipy.signal.correlate では手法を指定できますし、指定しなくとも早い方を自動で選択してくれます（Ver.1.13.0）。 numpy.correlate は直接法のみです（Ver.1.26）。

以下、計算例です。設定次第なのですが、以下だと答えは①≠②=③
①直接法

sig1 = np.std(data1)
sig2 = np.std(data2)
data1 = data1-data1.mean()
data2 = data2-data2.mean()
cross_correlation = np.correlate(data1, data2, mode='same'
                                ) /sig1/sig2/N

②FFT1

f1 = np.fft.fft(data1)
f2 = np.fft.fft(data2)
cross_spectrum = f1 * np.conj(f2)
cross_correlation = np.real(np.fft.ifft(cross_spectrum)
                            ) /sig1/sig2/N

③FFT2

f1 = np.fft.rfft(data1)
f2 = np.fft.rfft(data2)
cross_spectrum2 = f1 * np.conj(f2)
cross_correlation = np.fft.irfft(cross_spectrum2
                                ) /sig1/sig2/N

0をパディング + シフトで①＝②となります。

①直接法

padded_data1 = np.pad(data1, (0, N), 'constant')
padded_data2 = np.pad(data2, (0, N), 'constant')
cross_correlation = np.correlate(padded_data1,
　　　　　　　　　　　　　　　　 padded_data2,
　　　　　　　　　　　　　　　　 mode='same'
                                ) /sig1/sig2/N

②FFT1

f1 = np.fft.fft(padded_data1)
f2 = np.fft.fft(padded_data2)
cross_spectrum = f1 * np.conj(f2)
cross_correlation = np.real(np.fft.ifft(cross_spectrum))
cross_correlation = np.fft.fftshift(cross_correlation
                                   ) /sig1/sig2/N