"""Finite-window spectral experiment. Python 3; no external packages."""
import cmath
import json
from math import cos, isfinite, pi


def spectrum(values, detrend="linear", window="hann"):
    x = [float(v) for v in values]
    n = len(x)
    if n < 4 or not all(isfinite(v) for v in x):
        raise ValueError("Need at least four finite real observations")
    if detrend not in ("none", "mean", "linear"):
        raise ValueError("Unknown detrending method")
    if window not in ("boxcar", "hann"):
        raise ValueError("Unknown window")
    center = (n - 1) / 2
    mean = sum(x) / n
    slope = (sum((i - center) * (v - mean) for i, v in enumerate(x))
             / sum((i - center)**2 for i in range(n)))
    y = [v if detrend == "none" else
         v - mean - (slope * (i - center) if detrend == "linear" else 0)
         for i, v in enumerate(x)]
    w = [1.0 if window == "boxcar" else
         0.5 - 0.5 * cos(2 * pi * i / n) for i in range(n)]
    amplitude = []
    for k in range(n // 2 + 1):
        z = sum(w[i] * y[i] * cmath.exp(-2j * pi * k * i / n)
                for i in range(n))
        multiplier = 1 if k == 0 or (n % 2 == 0 and k == n // 2) else 2
        amplitude.append(multiplier * abs(z) / sum(w))
    return amplitude


def peak(amplitude):
    # No non-DC peak is defined for a numerically flat residual.
    if max(amplitude[1:]) < 1e-10:
        return None
    return max(range(1, len(amplitude)), key=amplitude.__getitem__)


def experiment():
    n = 128
    x = [10 + 0.05*t + 2*cos(2*pi*10.5*t/n) for t in range(n)]
    finite = []
    for method, window in (("none", "boxcar"), ("mean", "boxcar"),
                           ("linear", "boxcar"), ("linear", "hann")):
        a = spectrum(x, method, window)
        k = peak(a)
        finite.append({"detrend": method, "window": window,
                       "DC": round(a[0], 6), "peak_bin": k,
                       "period_bars": round(n/k, 6),
                       "peak_amplitude": round(a[k], 6),
                       "bin_20_amplitude": round(a[20], 6)})
    changing = [20 + 0.02*t + cos(2*pi*t/(16 if t < 192 else 8))
                for t in range(384)]
    rolling = []
    for end in (127, 191, 255, 319, 383):
        a = spectrum(changing[end-n+1:end+1])
        k = peak(a)
        rolling.append({"start": end-n+1, "end": end, "peak_bin": k,
                        "period_bars": round(n/k, 6),
                        "amplitude_bin_8": round(a[8], 6),
                        "amplitude_bin_16": round(a[16], 6)})
    return {"finite_window": finite, "rolling": rolling}


if __name__ == "__main__":
    print(json.dumps(experiment(), indent=2))