# Smoothed cepstral peak prominence (CPPS), the Chapter 8 primary outcome, from scratch on
# synthetic vowels. Cepstrum = spectrum of the dB power spectrum; a periodic voice gives a
# peak at the quefrency of its period. CPP = peak height (dB) above a regression line;
# CPPS first smooths over time and quefrency (after Hillenbrand et al. 1994; Hillenbrand &
# Houde 1996). The settings here are simplified ones.
import cmath, math, random

FS, N = 16000, 1024

def fft(x):
    if len(x) == 1:
        return x
    e, o, n = fft(x[0::2]), fft(x[1::2]), len(x)
    t = [cmath.exp(-2j * math.pi * k / n) * o[k] for k in range(n // 2)]
    return [e[k] + t[k] for k in range(n // 2)] + [e[k] - t[k] for k in range(n // 2)]

def vowel(breath, f0=200.0, jitter=0.002, dur=0.6, seed=3):
    rnd, x, ph = random.Random(seed), [], 0.0
    for _ in range(int(dur * FS)):
        ph += 2 * math.pi * f0 * (1 + rnd.gauss(0, jitter)) / FS
        harmonics = sum(math.sin(h * ph) / h ** 1.5 for h in range(1, 20))
        x.append(harmonics + rnd.gauss(0, 3 * breath))
    return x

def smooth(v, w):                                   # centred moving average, window w
    h = w // 2
    return [sum(v[max(0, i - h):i + h + 1]) / len(v[max(0, i - h):i + h + 1])
            for i in range(len(v))]

def cpp(x, smoothed=True, hop=320):
    win = [0.5 - 0.5 * math.cos(2 * math.pi * i / (N - 1)) for i in range(N)]
    ceps = []
    for s in range(0, len(x) - N, hop):
        spec = fft([x[s + i] * win[i] for i in range(N)])
        logp = [10 * math.log10(abs(v) ** 2 + 1e-12) for v in spec]
        ceps.append([20 * math.log10(abs(v) / N + 1e-9) for v in fft(logp)[:N // 2]])
    # 3 frames in time, 3 bins in quefrency
    if smoothed:
        cols = [smooth([c[q] for c in ceps], 3) for q in range(N // 2)]
        ceps = [smooth([cols[q][i] for q in range(N // 2)], 3) for i in range(len(ceps))]
    qs, vals = list(range(16, N // 2)), []          # regression line from 1 ms upward
    mq = sum(qs) / len(qs)
    for c in ceps:
        mc = sum(c[q] for q in qs) / len(qs)
        b = sum((q - mq) * (c[q] - mc) for q in qs) / sum((q - mq) ** 2 for q in qs)
        qp = max(range(FS // 330, FS // 60), key=lambda q: c[q])      # search 60-330 Hz
        vals.append(c[qp] - (mc + b * (qp - mq)))
    return sum(vals) / len(vals)

print("synthetic 200 Hz vowels | CPPS (smoothed) | CPP (unsmoothed)")
for name, b in (("clear", 0.02), ("mildly breathy", 0.08), ("breathy", 0.25)):
    x = vowel(b)
    print(f"  {name:21s} | {cpp(x):11.2f} dB  | {cpp(x, smoothed=False):12.2f} dB")
print("values depend on the settings: a trial must fix software, version and settings")
# synthetic 200 Hz vowels | CPPS (smoothed) | CPP (unsmoothed)
#   clear                 |       21.52 dB  |        23.03 dB
#   mildly breathy        |       15.54 dB  |        17.56 dB
#   breathy               |        9.06 dB  |        14.92 dB
# values depend on the settings: a trial must fix software, version and settings

# Sample size for the Chapter 8 trial: two parallel arms, 1:1, ANCOVA. Primary outcome:
# read-passage CPPS (dB) right after the session, adjusted for baseline (variance x
# 1-rho^2):
#   n per arm = 2 (z_a + z_b)^2 sd^2 (1 - rho^2) / delta^2 + z_a^2 / 4   (normal approx.)
# Assumptions, to be checked in a pilot used for nothing else: between-person SD 1.2 dB,
# pre-post correlation 0.7 (cautious), smallest effect worth acting on 0.4 dB.
import math
from statistics import NormalDist

z = NormalDist().inv_cdf

def per_arm(delta, sd, rho, alpha=0.05, power=0.90):
    za, zb = z(1 - alpha / 2), z(power)
    n = 2 * (za + zb) ** 2 * sd ** 2 * (1 - rho ** 2) / delta ** 2 + za ** 2 / 4
    return math.ceil(n)

SD, RHO, SESOI, LOSS = 1.2, 0.7, 0.4, 0.10
print(f"sd {SD} dB, alpha 0.05 two-sided, 90 % power; analysable participants per arm")
print("  SESOI (dB) | rho = 0.6 | rho = 0.7 | rho = 0.8")
for d in (0.3, 0.4, 0.5):
    cells = "   |   ".join(f"{per_arm(d, SD, r):5d}" for r in (0.6, 0.7, 0.8))
    print(f"     {d:.1f}     |   {cells}")
n = per_arm(SESOI, SD, RHO)
rec = math.ceil(n / (1 - LOSS))
print(f"chosen: SESOI {SESOI} dB, rho {RHO} -> {n} analysable per arm, {2 * n} in all")
print(f"allowing {LOSS:.0%} loss (unusable recordings, withdrawal): "
      f"recruit {rec} per arm, {2 * rec} in all")
print(f"if rho were only 0.6: {per_arm(SESOI, SD, 0.6)} per arm would be needed")
# sd 1.2 dB, alpha 0.05 two-sided, 90 % power; analysable participants per arm
#   SESOI (dB) | rho = 0.6 | rho = 0.7 | rho = 0.8
#      0.3     |     217   |     173   |     123
#      0.4     |     123   |      98   |      70
#      0.5     |      79   |      63   |      45
# chosen: SESOI 0.4 dB, rho 0.7 -> 98 analysable per arm, 196 in all
# allowing 10% loss (unusable recordings, withdrawal): recruit 109 per arm, 218 in all
# if rho were only 0.6: 123 per arm would be needed

# Decision rule for the Chapter 8 trial, fixed before any data are collected. diff =
# (straw + monitoring) minus (usual reading warm-up) in read-passage CPPS (dB),
# ANCOVA-adjusted, with 95 % CI; positive = more periodic, less noisy voice. SESOI 0.4 dB.
# Checked first: protocol fidelity and blinding of the acoustic analyst.
SESOI = 0.4

def verdict(diff, lo, hi, fidelity_ok=True, analyst_blind=True, harms=0):
    if not analyst_blind:
        return "NOT INTERPRETABLE: the acoustic analyst learned the allocation"
    if not fidelity_ok:
        return "NOT INTERPRETABLE: < 80 % of sessions delivered as scripted"
    if lo >= SESOI:
        v = "improves CPPS by at least the SESOI (0.4 dB)"
    elif lo > 0:
        v = "improves CPPS; size not shown to reach 0.4 dB"
    elif hi < 0:
        v = "LOWERS CPPS"
    elif lo >= -SESOI and hi <= SESOI:
        v = "no difference worth acting on (within +/- 0.4 dB)"
    else:
        v = "unclear: the interval is too wide to decide"
    if harms:
        v += f"; {harms} discomfort reports shown prominently"
    return v

for case in ((0.70, 0.45, 0.95), (0.30, 0.05, 0.55), (0.05, -0.20, 0.30),
             (-0.35, -0.60, -0.10), (0.20, -0.35, 0.75)):
    print(f"{case[0]:+.2f} dB ({case[1]:+.2f} to {case[2]:+.2f}): {verdict(*case)}")
print("+0.60 dB (+0.35 to +0.85), blinding broken:")
print("   " + verdict(0.6, 0.35, 0.85, analyst_blind=False))
print("+0.50 dB (+0.42 to +0.58), 3 reports of strain:")
print("   " + verdict(0.5, 0.42, 0.58, harms=3))
print("secondary outcomes (Borg CR10 effort, sustained-vowel CPPS, HNR): estimates with")
print("intervals only; they never change the verdict on the primary outcome")
# +0.70 dB (+0.45 to +0.95): improves CPPS by at least the SESOI (0.4 dB)
# +0.30 dB (+0.05 to +0.55): improves CPPS; size not shown to reach 0.4 dB
# +0.05 dB (-0.20 to +0.30): no difference worth acting on (within +/- 0.4 dB)
# -0.35 dB (-0.60 to -0.10): LOWERS CPPS
# +0.20 dB (-0.35 to +0.75): unclear: the interval is too wide to decide
# +0.60 dB (+0.35 to +0.85), blinding broken:
#    NOT INTERPRETABLE: the acoustic analyst learned the allocation
# +0.50 dB (+0.42 to +0.58), 3 reports of strain:
#    improves CPPS by at least the SESOI (0.4 dB); 3 discomfort reports shown prominently
# secondary outcomes (Borg CR10 effort, sustained-vowel CPPS, HNR): estimates with
# intervals only; they never change the verdict on the primary outcome