code: def calc_hdev_phase(phase, rate, mj, stride):
      
  
      tau0 = 1.0 / float(rate)
      mj = int(mj)
      stride = int(stride)
      d3 = phase[3 * mj::stride]
      d2 = phase[2 * mj::stride]
      d1 = phase[1 * mj::stride]
      d0 = phase[::stride]
  
      n = min(len(d0), len(d1), len(d2), len(d3))
  
      v_arr = d3[:n] - 3 * d2[:n] + 3 * d1[:n] - d0[:n]
  
      s = np.sum(v_arr * v_arr)
  
      if n == 0:
          n = 1
  
      h = np.sqrt(s / 6.0 / float(n)) / float(tau0 * mj)
      e = h / np.sqrt(n)
      return h, e, n
docs: main calculation fungtion for HDEV and OHDEV
  
      Parameters
      ----------
      phase: np.array
          Phase data in seconds.
      rate: float
          The sampling rate for phase or frequency, in Hz
      mj: int
          M index value for stride
      stride: int
          Size of stride
  
      Returns
      -------
      (dev, deverr, n): tuple
          Array of computed values.
  
      Notes
      -----
      http://www.leapsecond.com/tools/adev_lib.c
                           1        N-3
           s2y(t) = --------------- sum [x(i+3) - 3x(i+2) + 3x(i+1) - x(i) ]^2
                    6*tau^2 (N-3m)  i=1
  
          N=M+1 phase measurements
          m is averaging factor
  
      NIST SP 1065 eqn (18) and (20) pages 20 and 21