import math, random

FS, F0 = 44100, 220.0

def synth(jitter, shimmer, noise, dur=0.6, seed=1):
    """A sung vowel-like tone: each cycle has its own period and amplitude, plus breath noise."""
    rnd = random.Random(seed)
    x, t0 = [], 1 / F0
    while len(x) < dur * FS:
        T = t0 * (1 + rnd.gauss(0, jitter))
        A = 1 + rnd.gauss(0, shimmer)
        n = round(T * FS)
        x += [A * sum(math.sin(2 * math.pi * h * k / n) / h ** 1.2 for h in range(1, 9)) for k in range(n)]
    return [v + rnd.gauss(0, noise) for v in x]

def smooth(x, w=9):
    h = w // 2
    out, acc = [], sum(x[:w])
    out = [0.0] * h
    for k in range(h, len(x) - h - 1):
        out.append(acc / w)
        acc += x[k + h + 1] - x[k - h]
    return out + [0.0] * (len(x) - len(out))

def measure(x):
    n0 = round(FS / F0)
    xs = smooth(x)                       # light low-pass so noise does not move the peaks
    # cycle peaks: the largest sample in successive windows that follow the last peak
    peaks, i = [], max(range(10, n0 + 10), key=lambda k: xs[k])
    while i + int(1.3 * n0) < len(x) - 10:
        j = max(range(i + int(0.7 * n0), i + int(1.3 * n0)), key=lambda k: xs[k])
        a, b, c = xs[j - 1], xs[j], xs[j + 1]
        peaks.append((j + 0.5 * (a - c) / (a - 2 * b + c), b))   # parabolic interpolation
        i = j
    P = [(q[0] - p[0]) / FS for p, q in zip(peaks, peaks[1:])]
    A = [v for _, v in peaks]
    jit = sum(abs(a - b) for a, b in zip(P, P[1:])) / (len(P) - 1) / (sum(P) / len(P)) * 100
    shim = sum(abs(a - b) for a, b in zip(A, A[1:])) / (len(A) - 1) / (sum(A) / len(A)) * 100
    # harmonics-to-noise ratio from the normalised autocorrelation at the period lag
    m = sum(x) / len(x)
    y = [v - m for v in x]
    e = sum(v * v for v in y)
    r = max(sum(y[k] * y[k + L] for k in range(len(y) - L)) / e * len(y) / (len(y) - L)
            for L in range(int(0.8 * n0), int(1.2 * n0)))
    hnr = 10 * math.log10(r / (1 - r))
    return jit, shim, hnr

for name, j, s, nz in (("tense, rough", 0.012, 0.06, 0.25), ("relaxed, clear", 0.003, 0.02, 0.05)):
    jit, shim, hnr = measure(synth(j, s, nz))
    print(f"{name:15s} measured jitter {jit:4.2f} %, shimmer {shim:4.2f} %, HNR {hnr:5.1f} dB "
          f"(built with jitter ≈ {1.128 * j * 100:.2f} %, shimmer ≈ {1.128 * s * 100:.1f} %)")
# tense, rough    measured jitter 2.23 %, shimmer 9.50 %, HNR   9.7 dB (built with jitter ≈ 1.35 %, shimmer ≈ 6.8 %)
# relaxed, clear  measured jitter 0.65 %, shimmer 2.79 %, HNR  22.7 dB (built with jitter ≈ 0.34 %, shimmer ≈ 2.3 %)