编辑代码

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