Expected B1 ratio for given time-series length N and exponent mu
    
        FIXME: add reference (paper & link)
        
        The exponents are defined as
        S_y(f) = h_a f^alpha    (power spectrum of y)
        S_x(f) = g_b f^b        (power spectrum of x)
        bias = const * tau^mu
        
        and (b, alpha, mu) relate to eachother by:
        b    alpha   mu
        0    +2      -2
       -1    +1      -2   resolve between -2 cases with R(n)
       -2     0      -1
       -3    -1       0
       -4    -2      +1
       -5    -3      +2
       -6    -4      +3 for HDEV, by applying B1 to frequency data, and add +2 to resulting mu

def b1_theory(N, mu):
    """ Expected B1 ratio for given time-series length N and exponent mu
    
        FIXME: add reference (paper & link)
        
        The exponents are defined as
        S_y(f) = h_a f^alpha    (power spectrum of y)
        S_x(f) = g_b f^b        (power spectrum of x)
        bias = const * tau^mu
        
        and (b, alpha, mu) relate to eachother by:
        b    alpha   mu
        0    +2      -2
       -1    +1      -2   resolve between -2 cases with R(n)
       -2     0      -1
       -3    -1       0
       -4    -2      +1
       -5    -3      +2
       -6    -4      +3 for HDEV, by applying B1 to frequency data, and add +2 to resulting mu
    """
    
    # see Table 3 of Howe 2000
    if mu == 2:
        return float(N)*(float(N)+1.0)/6.0
        #up = N*(1.0-pow(N, mu))
        #down = 2*(N-1.0)*(1-pow(2.0, mu))
        #return up/down        
    elif mu == 1:
        return float(N)/2.0
    elif mu == 0:
        return N*np.log(N)/(2.0*(N-1.0)*np.log(2))
    elif mu == -1:
        return 1
    elif mu == -2:
        return (pow(N,2)-1.0)/(1.5*N*(N-1.0))
    else:
        up = N*(1.0-pow(N, mu))
        down = 2*(N-1.0)*(1-pow(2.0, mu))
        return up/down
        
    assert False