code: def b1_theory(N, 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
docs: 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