code: def crosscov(x, y, axis=-1, all_lags=False, debias=True, normalize=True):
      
      if x.shape[axis] != y.shape[axis]:
          raise ValueError(
              'crosscov() only works on same-length sequences for now'
              )
      if debias:
          x = remove_bias(x, axis)
          y = remove_bias(y, axis)
      slicing = [slice(d) for d in x.shape]
      slicing[axis] = slice(None, None, -1)
      cxy = fftconvolve(x, y[tuple(slicing)].conj(), axis=axis, mode='full')
      N = x.shape[axis]
      if normalize:
          cxy /= N
      if all_lags:
          return cxy
      slicing[axis] = slice(N - 1, 2 * N - 1)
      return cxy[tuple(slicing)]
docs: Returns the crosscovariance sequence between two ndarrays.
      This is performed by calling fftconvolve on x, y[::-1]
  
      Parameters
      ----------
  
      x : ndarray
      y : ndarray
      axis : time axis
      all_lags : {True/False}
         whether to return all nonzero lags, or to clip the length of s_xy
         to be the length of x and y. If False, then the zero lag covariance
         is at index 0. Otherwise, it is found at (len(x) + len(y) - 1)/2
      debias : {True/False}
         Always removes an estimate of the mean along the axis, unless
         told not to (eg X and Y are known zero-mean)
  
      Returns
      -------
  
      cxy : ndarray
         The crosscovariance function
  
      Notes
      -----
  
      cross covariance of processes x and y is defined as
  
      .. math::
  
      C_{xy}[k]=E\{(X(n+k)-E\{X\})(Y(n)-E\{Y\})^{*}\}
  
      where X and Y are discrete, stationary (or ergodic) random processes
  
      Also note that this routine is the workhorse for all auto/cross/cov/corr
      functions.