function crossCorrelation(wave1, wave2, n) result(corr)
real :: wave1(n), wave2(n)
integer :: n, i
real :: sum1, sum2, ave1, ave2, dev1, dev2, var1, var2, cov
sum1 = 0.0
sum2 = 0.0
do i = 1, n
sum1 = sum1 + wave1(i)
sum2 = sum2 + wave2(i)
end do
ave1 = sum1 / real(n)
ave2 = sum2 / real(n)
var1 = 0.0
var2 = 0.0
cov = 0.0
do i = 1, n
dev1 = wave1(i) - ave1
dev2 = wave2(i) - ave2
var1 = var1 + dev1 * dev1
var2 = var2 + dev2 * dev2
cov = cov + dev1 * dev2
end do
var1 = var1 / real(n)
var2 = var2 / real(n)
cov = cov / real(n)
if (var1 == 0.0 .or. var2 == 0.0) then
corr = 0.0
else
corr = cov / (sqrt(var1) * sqrt(var2))
end if
end function crossCorrelation