def autocorr_noise_id(x, af, data_type="phase", dmin=0, dmax=2):
    """ Lag-1 autocorrelation based noise identification
        
    Parameters
    ----------
    x: numpy.array
        phase or fractional frequency time-series data
        minimum recommended length is len(x)>30 roughly.
    af: int
        averaging factor
    data_type: string {'phase', 'freq'}
        "phase" for phase data in seconds
        "freq" for fractional frequency data
    dmin: int 
        minimum required number of differentiations in the algorithm
    dmax: int
        maximum number of differentiations
        defaults to 2 for ADEV
        set to 3 for HDEV
        
    Returns
    -------
    alpha_int: int
        noise-slope as integer
    alpha: float
        noise-slope as float
    d: int
        number of differentiations of the time-series performed
    
    Notes
    -----
        http://www.stable32.com/Auto.pdf
        http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.503.9864&rep=rep1&type=pdf
        
            
    Power law noise identification using the lag 1 autocorrelation
    Riley,W.J. et al.
    18th European Frequency and Time Forum (EFTF 2004)
    https://ieeexplore.ieee.org/document/5075021
        
    """
    d = 0 # number of differentiations
    lag = 1
    if data_type is "phase":
        
        if af>1:
            #x = scipy.signal.decimate(x, af, n=1, ftype='fir')
            x = x[0:len(x):af] # decimate by averaging factor
        x = detrend(x, deg=2) # remove quadratic trend (frequency offset and drift)
    elif data_type is "freq":
        # average by averaging factor
        y_cut = np.array(x[:len(x)-(len(x)%af)]) # cut to length
        assert len(y_cut)%af == 0
        y_shaped = y_cut.reshape((int(len(y_cut)/af), af))
        x = np.average(y_shaped, axis=1) # average
        x = detrend(x, deg=1) # remove frequency drift
    
    # require minimum length for time-series
    if len(x)<30:
        print("autocorr_noise_id() Don't know how to do noise-ID for time-series length= %d"%len(x))
        raise NotImplementedError
    
    while True:
        r1 = lag1_acf(x)
        rho = r1/(1.0+r1)
        if d >= dmin and ( rho < 0.25 or d >= dmax ):
            p = -2*(rho+d)
            #print r1
            #assert r1 < 0
            #assert r1 > -1.0/2.0
            phase_add2 = 0
            if data_type is "phase":
                phase_add2 = 2
            alpha = p+phase_add2
            alpha_int = int( -1.0*np.round(2*rho) - 2.0*d )+phase_add2 
            #print "d=",d,"alpha=",p+2
            return alpha_int, alpha, d, rho
        else:
            x = np.diff(x)
            d = d + 1
    assert False