code: def mtotdev(data, rate=1.0, data_type="phase", taus=None):
      
      phase = input_to_phase(data, rate, data_type)
      (phase, ms, taus_used) = tau_generator(phase, rate, taus,
                                             maximum_m=float(len(phase))/3.0)
      devs = np.zeros_like(taus_used)
      deverrs = np.zeros_like(taus_used)
      ns = np.zeros_like(taus_used)
  
      for idx, mj in enumerate(ms):
          devs[idx], deverrs[idx], ns[idx] = calc_mtotdev_phase(phase, rate, mj)
  
      return remove_small_ns(taus_used, devs, deverrs, ns)
docs: PRELIMINARY - REQUIRES FURTHER TESTING.
          Modified Total deviation.
          Better confidence at long averages for modified Allan
  
          FIXME: bias-correction http://www.wriley.com/CI2.pdf page 6
  
          The variance is scaled up (divided by this number) based on the
          noise-type identified.
          WPM 0.94
          FPM 0.83
          WFM 0.73
          FFM 0.70
          RWFM 0.69
  
      Parameters
      ----------
      data: np.array
          Input data. Provide either phase or frequency (fractional,
          adimensional).
      rate: float
          The sampling rate for data, in Hz. Defaults to 1.0
      data_type: {'phase', 'freq'}
          Data type, i.e. phase or frequency. Defaults to "phase".
      taus: np.array
          Array of tau values, in seconds, for which to compute statistic.
          Optionally set taus=["all"|"octave"|"decade"] for automatic
          tau-list generation.
  
      NIST SP 1065 eqn (27) page 25