VSB Power Line Blog - Detrended Fluctuation Analysis FD

dfa_fd.py

def logarithmic_n(min_n, max_n, factor):
  # stop condition: min * f^x = max
  # => f^x = max/min
  # => x = log(max/min) / log(f)
  max_i = int(np.floor(np.log(1.0 * max_n / min_n) / np.log(factor)))
  ns = [min_n]
  for i in range(max_i + 1):
    n = int(np.floor(min_n * (factor ** i)))
    if n > ns[-1]:
      ns.append(n)
  return ns

def dfa(data):
  data = np.asarray(data)
  total_N = len(data)
  # nvals contains the different number of subseries to be used
  nvals = logarithmic_n(4, 0.1 * total_N, 1.5)
  walk = np.cumsum(data - np.mean(data))
  fluctuations = []

  for n in nvals:
    # total_N - (total_N % n) will ensure we get whole number of subseries of data
    d = walk[:total_N - (total_N % n)]
    d = d.reshape((total_N // n, n))

    # Calculate trends in each subseries (window)
    x = np.arange(n)
    # tpoly contains the coefficients obtained after best-fitting a polynomial in each subseries (window)
    tpoly = [np.polyfit(x, d[i], 2) for i in range(len(d))]
    tpoly = np.array(tpoly)
    # Predicted values is obtained by fitting the coefficients to the subseries
    trend = np.array([np.polyval(tpoly[i], x) for i in range(len(d))])
    # Calculate the standard deviations ("fluctuations") of walks in d around trend
    flucs = np.sqrt(np.sum((d - trend)**2, axis=1) / n)
    # Calculate the mean fluctuation over all the subseries
    f_n = np.sum(flucs)/len(flucs)
    fluctuations.append(f_n)
  # filter zeros from fluctuations
  # nonzero = np.where(fluctuations != 0)
  # nvals = np.array(nvals)[nonzero]
  # fluctuations = fluctuations[nonzero]
  if len(fluctuations) == 0:
    # all fluctutations are zero => we cannot fit a line
    poly = [np.nan, np.nan]
  else:
    poly = np.polyfit(np.log10(nvals), np.log10(fluctuations), deg=1)
  return poly[0]