相互相関関数のお遊び (python)

ちょっと必要があって簡単な相互相関関数の計算をした。

一応、離散データ $x(i)$ , $y(i)$ から相互相関関数 $R_{xy}(i)$ を

$R_{xy}(i)=\frac{1}{N-i}\sum_{j=0}^{N-1-i}x(j) \cdot y(i+j)$

で計算させてみた*1。

#!/usr/bin/env python

import array    
import math
import ROOT
ROOT.gStyle.SetTitleBorderSize(0)

def cross_correlation(array1, array2):
    nx = len(array1)
    ny = len(array2)
    if nx != ny:
        print "Error: lengths of two lists does not math each other."

    cross_corr = array.array('f')
    for i in range(nx):
        sum = 0.0
        for j in range(nx - i):
            sum = sum + float(array1[j]) * float(array2[i+j])
        cross_corr.append(sum/float(nx-i))
    return cross_corr

resolution = 0.2
x  = array.array('f')
y1 = array.array('f')
y2 = array.array('f')
for i in range(80):
    x.append(i)
    y1.append(math.cos(resolution*i))
    y2.append(math.sin(resolution*i))

r = cross_correlation(y1, y2)

can = ROOT.TCanvas("can", "can", 0, 0, 600, 500)

gra1 = ROOT.TGraph(len(y1), x, y1)
gra1.SetTitle('Cross-Correlation;X;Y')
gra1.SetMarkerColor(1)
gra1.SetLineColor(1)
gra1.SetMarkerSize(0.8)
gra1.SetMarkerStyle(8)
gra1.Draw("APL")

gra2 = ROOT.TGraph(len(y2), x, y2)
gra2.SetMarkerColor(3)
gra2.SetLineColor(3)
gra2.SetMarkerSize(0.8)
gra2.SetMarkerStyle(8)
gra2.Draw("PL same")

gra3 = ROOT.TGraph(len(r), x, r)
gra3.SetMarkerColor(2)
gra3.SetLineColor(2)
gra3.SetMarkerSize(0.8)
gra3.SetMarkerStyle(8)
gra3.Draw("PL same")

*1:http://www4.plala.or.jp/tamo/dsp/07cross-corr.pdf