// organ_dsp.el — the AFFERENT DSP organ, own-core, ported from peripheral/src/periph.swift.
//
// The mirror of speech.el: where speech.el RENDERS meaning out through a voice
// (efferent), this module HEARS — it takes a raw 16-bit PCM RIFF/WAVE capture
// and metabolizes it into GEOMETRY: a compact 8-number audio descriptor, a
// voice-signature (F0 + formants F1-F5 by LPC), and an LPC analysis-resynthesis
// that speaks the heard voice back from its own signature.
//
// Everything here is float physics (math_sin/math_cos/math_sqrt), not the
// fixed-point integer path speech.el uses for synthesis — the analysis side
// needs the dynamic range that autocorrelation and Levinson-Durbin demand.
//
// Binary I/O note: an El String truncates at the first NUL, so a WAV can never
// be read through fs_read(). We read it through fs_read_b64_chunk(), which
// hands back plain-ASCII base64 of a byte window, and decode that base64 HERE,
// in El, into a [Int] of byte values. Chunks are a multiple of 3 bytes so each
// base64 window decodes cleanly with no padding except the final one.
//
// Every function is prefixed dsp_ so nothing here can collide with speech.el
// (write_wav / wav_le16 / wav_le32 / sp_* / voice_* all live there and are NOT
// redefined). imitate returns [Int] samples — hand them to speech.el's write_wav.

// ---------------------------------------------------------------------------
// Base64 -> bytes (own-core; the only way binary reaches El intact)
// ---------------------------------------------------------------------------

// Standard RFC 4648 alphabet A-Za-z0-9+/ -> 0..63. Padding '=' and any other
// character -> -1 (a sentinel; there is no exception handling in El).
fn dsp_b64_val(c: Int) -> Int {
    if c >= 65 {
        if c <= 90 {
            return c - 65
        }
    }
    if c >= 97 {
        if c <= 122 {
            return c - 71
        }
    }
    if c >= 48 {
        if c <= 57 {
            return c + 4
        }
    }
    if c == 43 {
        return 62
    }
    if c == 47 {
        return 63
    }
    return 0 - 1
}

// Read a whole file as a list of byte values 0..255. Empty list on failure.
fn dsp_read_bytes(path: String) -> [Int] {
    let bytes: [Int] = native_list_empty()
    let size: Int = fs_size(path)
    if size <= 0 {
        return bytes
    }
    let chunk: Int = 60000                 // multiple of 3 -> no interior padding
    let off: Int = 0
    while off < size {
        let s: String = fs_read_b64_chunk(path, off, chunk)
        let sl: Int = str_len(s)
        if sl < 4 {
            return bytes
        }
        let i: Int = 0
        while i + 3 < sl {
            let c0: Int = dsp_b64_val(str_char_code(s, i))
            let c1: Int = dsp_b64_val(str_char_code(s, i + 1))
            let c2: Int = dsp_b64_val(str_char_code(s, i + 2))
            let c3: Int = dsp_b64_val(str_char_code(s, i + 3))
            if c0 < 0 {
                return bytes
            }
            if c1 < 0 {
                return bytes
            }
            let b0: Int = c0 * 4 + c1 / 16
            bytes = native_list_append(bytes, b0)
            if c2 >= 0 {
                let lo1: Int = c1 - (c1 / 16) * 16
                let b1: Int = lo1 * 16 + c2 / 4
                bytes = native_list_append(bytes, b1)
                if c3 >= 0 {
                    let lo2: Int = c2 - (c2 / 4) * 4
                    let b2: Int = lo2 * 64 + c3
                    bytes = native_list_append(bytes, b2)
                }
            }
            i = i + 4
        }
        off = off + chunk
    }
    return bytes
}

fn dsp_rd16(b: [Int], o: Int) -> Int {
    let b0: Int = native_list_get(b, o)
    let b1: Int = native_list_get(b, o + 1)
    return b0 + b1 * 256
}

fn dsp_rd32(b: [Int], o: Int) -> Int {
    let b0: Int = native_list_get(b, o)
    let b1: Int = native_list_get(b, o + 1)
    let b2: Int = native_list_get(b, o + 2)
    let b3: Int = native_list_get(b, o + 3)
    return b0 + b1 * 256 + b2 * 65536 + b3 * 16777216
}

// Four bytes at o compared against a 4-char ASCII chunk id.
fn dsp_chunk_is(b: [Int], o: Int, id: String) -> Bool {
    let k: Int = 0
    while k < 4 {
        let got: Int = native_list_get(b, o + k)
        let want: Int = str_char_code(id, k)
        if got != want {
            return false
        }
        k = k + 1
    }
    return true
}

// ---------------------------------------------------------------------------
// readWavSamples — 16-bit PCM RIFF/WAVE -> normalized samples.
// Walks chunks to find 'fmt ' and 'data', so JUNK/FLLR padding (which
// AVAudioRecorder emits) is stepped over rather than mistaken for audio.
// Channel 0 only if stereo.
//
// The returned list is PACKED: [ sr, ch, n, s0, s1, ... s(n-1) ] with the three
// header numbers carried as Floats (El has no tuples). Use dsp_wav_sr /
// dsp_wav_ch / dsp_wav_n / dsp_wav_pcm to open it. Empty list on failure.
// ---------------------------------------------------------------------------
fn dsp_read_wav(path: String) -> [Float] {
    let out: [Float] = native_list_empty()
    let d: [Int] = dsp_read_bytes(path)
    let dn: Int = native_list_len(d)
    if dn <= 44 {
        return out
    }
    let ch: Int = 0
    let sr: Int = 0
    let bits: Int = 0
    let o: Int = 12
    while o + 8 <= dn {
        let sz: Int = dsp_rd32(d, o + 4)
        if dsp_chunk_is(d, o, "fmt ") {
            if o + 24 <= dn {
                ch = dsp_rd16(d, o + 10)
                sr = dsp_rd32(d, o + 12)
                bits = dsp_rd16(d, o + 22)
            }
        }
        if dsp_chunk_is(d, o, "data") {
            if bits != 16 {
                return out
            }
            if ch <= 0 {
                return out
            }
            // Faithful to periph.swift, including its `d.count - 1` bound and
            // the `while i + 1 < end` test — the last sample of a file whose
            // data chunk runs to EOF is dropped there, so it is dropped here.
            let start: Int = o + 8
            let end: Int = start + sz
            if end > dn - 1 {
                end = dn - 1
            }
            let samples: [Float] = native_list_empty()
            let step: Int = 2 * ch
            let i: Int = start
            let count: Int = 0
            while i + 1 < end {
                let v: Int = dsp_rd16(d, i)
                if v >= 32768 {
                    v = v - 65536
                }
                let f: Float = int_to_float(v) / 32768.0
                samples = native_list_append(samples, f)
                count = count + 1
                i = i + step
            }
            let srf: Float = int_to_float(sr)
            let chf: Float = int_to_float(ch)
            let nf: Float = int_to_float(count)
            out = native_list_append(out, srf)
            out = native_list_append(out, chf)
            out = native_list_append(out, nf)
            let j: Int = 0
            while j < count {
                let sv: Float = native_list_get(samples, j)
                out = native_list_append(out, sv)
                j = j + 1
            }
            return out
        }
        let adv: Int = 8 + sz + (sz - (sz / 2) * 2)
        if adv <= 0 {
            return out
        }
        o = o + adv
    }
    return out
}

fn dsp_wav_sr(w: [Float]) -> Int {
    if native_list_len(w) < 3 {
        return 0
    }
    let v: Float = native_list_get(w, 0)
    return float_to_int(v)
}

fn dsp_wav_ch(w: [Float]) -> Int {
    if native_list_len(w) < 3 {
        return 0
    }
    let v: Float = native_list_get(w, 1)
    return float_to_int(v)
}

fn dsp_wav_n(w: [Float]) -> Int {
    if native_list_len(w) < 3 {
        return 0
    }
    let v: Float = native_list_get(w, 2)
    return float_to_int(v)
}

// The bare sample list, unpacked from the header-prefixed form.
fn dsp_wav_pcm(w: [Float]) -> [Float] {
    let out: [Float] = native_list_empty()
    let n: Int = dsp_wav_n(w)
    let i: Int = 0
    while i < n {
        let v: Float = native_list_get(w, i + 3)
        out = native_list_append(out, v)
        i = i + 1
    }
    return out
}

// ---------------------------------------------------------------------------
// small float helpers (El has no unary minus on Float in every position, and
// no min/max builtin — so they are written out)
// ---------------------------------------------------------------------------
fn dsp_fabs(v: Float) -> Float {
    if v < 0.0 {
        return 0.0 - v
    }
    return v
}

// ---------------------------------------------------------------------------
// computeAudio — the compact audio descriptor.
// Returns the 8-number vector [seconds, sr, ch, rms, peak, zcr, centroid, f0].
// Empty list if the WAV cannot be read (no exceptions in El — sentinels only).
// ---------------------------------------------------------------------------
fn dsp_compute_audio(path: String) -> [Float] {
    let vec: [Float] = native_list_empty()
    let w: [Float] = dsp_read_wav(path)
    let n: Int = dsp_wav_n(w)
    if n <= 0 {
        return vec
    }
    let sr: Int = dsp_wav_sr(w)
    let ch: Int = dsp_wav_ch(w)
    let s: [Float] = dsp_wav_pcm(w)
    let nf: Float = int_to_float(n)
    let srf: Float = int_to_float(sr)
    let seconds: Float = nf / srf

    // energy / peak / zero crossings
    let sumsq: Float = 0.0
    let peak: Float = 0.0
    let zc: Float = 0.0
    let i: Int = 0
    while i < n {
        let v: Float = native_list_get(s, i)
        sumsq = sumsq + v * v
        let av: Float = dsp_fabs(v)
        if av > peak {
            peak = av
        }
        if i > 0 {
            let pv: Float = native_list_get(s, i - 1)
            let a: Bool = pv < 0.0
            let b: Bool = v < 0.0
            if a != b {
                zc = zc + 1.0
            }
        }
        i = i + 1
    }
    let rms: Float = math_sqrt(sumsq / nf)
    let zcr: Float = zc / nf * srf              // ~2*dominant freq for tonal

    // Spectral centroid via a coarse 64-bin DFT on a centered 2048 window.
    let ww: Int = 2048
    if n < ww {
        ww = n
    }
    let off: Int = (n - ww) / 2
    if off < 0 {
        off = 0
    }
    let two: Float = 2.0
    let pi: Float = math_pi()
    let num: Float = 0.0
    let den: Float = 0.0
    let bins: Int = 64
    let binsf: Float = 128.0                    // Double(2*bins)
    let k: Int = 1
    while k < bins {
        let kf: Float = int_to_float(k)
        let f: Float = kf * srf / binsf
        let re: Float = 0.0
        let im: Float = 0.0
        let j: Int = 0
        while j < ww {
            let jf: Float = int_to_float(j)
            let ang: Float = (0.0 - two) * pi * kf * jf / binsf
            let xv: Float = native_list_get(s, off + j)
            let cv: Float = math_cos(ang)
            let sv: Float = math_sin(ang)
            re = re + xv * cv
            im = im + xv * sv
            j = j + 1
        }
        let mag: Float = math_sqrt(re * re + im * im)
        num = num + f * mag
        den = den + mag
        k = k + 1
    }
    let centroid: Float = 0.0
    if den > 0.0 {
        centroid = num / den
    }

    // F0 by autocorrelation over the plausible speech range 70-400 Hz.
    let lag_min: Int = sr / 400
    let lag_max: Int = sr / 70
    if lag_max > n - 1 {
        lag_max = n - 1
    }
    let bound: Int = off + ww
    if bound > n {
        bound = n
    }
    let best_lag: Int = 0
    let best_corr: Float = 0.0
    if lag_max > lag_min {
        let lag: Int = lag_min
        while lag <= lag_max {
            let c: Float = 0.0
            let p: Int = 0
            while p + lag < bound {
                let a1: Float = native_list_get(s, off + p)
                let a2: Float = native_list_get(s, off + p + lag)
                c = c + a1 * a2
                p = p + 1
            }
            if c > best_corr {
                best_corr = c
                best_lag = lag
            }
            lag = lag + 1
        }
    }
    let f0: Float = 0.0
    if best_lag > 0 {
        f0 = srf / int_to_float(best_lag)
    }

    vec = native_list_append(vec, seconds)
    vec = native_list_append(vec, srf)
    let chf: Float = int_to_float(ch)
    vec = native_list_append(vec, chf)
    vec = native_list_append(vec, rms)
    vec = native_list_append(vec, peak)
    vec = native_list_append(vec, zcr)
    vec = native_list_append(vec, centroid)
    vec = native_list_append(vec, f0)
    return vec
}

// ---------------------------------------------------------------------------
// LPC core: hamming, autocorr, Levinson-Durbin, formant peak-pick, pitch.
// ---------------------------------------------------------------------------

fn dsp_hamming(x: [Float]) -> [Float] {
    let n: Int = native_list_len(x)
    if n < 2 {
        return x
    }
    let out: [Float] = native_list_empty()
    let pi: Float = math_pi()
    let dn: Float = int_to_float(n - 1)
    let i: Int = 0
    while i < n {
        let v: Float = native_list_get(x, i)
        let ang: Float = 2.0 * pi * int_to_float(i) / dn
        let cv: Float = math_cos(ang)
        let wv: Float = 0.54 - 0.46 * cv
        out = native_list_append(out, v * wv)
        i = i + 1
    }
    return out
}

// r[lag] = sum_i x[i]*x[i-lag], lag = 0..p. Returns p+1 numbers.
fn dsp_autocorr(x: [Float], p: Int) -> [Float] {
    let n: Int = native_list_len(x)
    let r: [Float] = native_list_empty()
    let lag: Int = 0
    while lag <= p {
        let acc: Float = 0.0
        let i: Int = lag
        while i < n {
            let a: Float = native_list_get(x, i)
            let b: Float = native_list_get(x, i - lag)
            acc = acc + a * b
            i = i + 1
        }
        r = native_list_append(r, acc)
        lag = lag + 1
    }
    return r
}

// Levinson-Durbin -> LPC coeffs a[0..p] with A(z) = 1 + sum a[k] z^-k.
// Returns p+2 numbers: a[0..p] followed by the residual energy at index p+1.
fn dsp_levinson(r: [Float], p: Int) -> [Float] {
    let a: [Float] = native_list_empty()
    a = native_list_append(a, 1.0)
    let z: Int = 1
    while z <= p {
        a = native_list_append(a, 0.0)
        z = z + 1
    }
    let err: Float = native_list_get(r, 0)
    if err <= 0.0 {
        a = native_list_append(a, 0.0)
        return a
    }
    let i: Int = 1
    let stopped: Bool = false
    while i <= p {
        if stopped {
            i = p + 1
        } else {
            let acc: Float = native_list_get(r, i)
            if i > 1 {
                let j: Int = 1
                while j < i {
                    let aj: Float = native_list_get(a, j)
                    let rij: Float = native_list_get(r, i - j)
                    acc = acc + aj * rij
                    j = j + 1
                }
            }
            let k: Float = (0.0 - acc) / err
            // na = a, then na[i] = k, then na[j] = a[j] + k*a[i-j] for 1<=j<i
            let na: [Float] = native_list_empty()
            na = native_list_append(na, 1.0)
            let m: Int = 1
            while m <= p {
                let v: Float = native_list_get(a, m)
                if m == i {
                    v = k
                }
                if m < i {
                    let am: Float = native_list_get(a, m)
                    let aim: Float = native_list_get(a, i - m)
                    v = am + k * aim
                }
                na = native_list_append(na, v)
                m = m + 1
            }
            a = na
            err = err * (1.0 - k * k)
            if err <= 0.0 {
                stopped = true
            }
            i = i + 1
        }
    }
    a = native_list_append(a, err)
    return a
}

// Formant peaks from the LPC all-pole spectral envelope, 512 steps over 0..sr/2,
// kept in 150..5200 Hz, at most 5. Returns FLAT PAIRS [f1, bw1, f2, bw2, ...]
// with the crude -3 dB (peak/sqrt2) bandwidth estimate.
fn dsp_formants(a: [Float], sr: Int) -> [Float] {
    let peaks: [Float] = native_list_empty()
    let p: Int = native_list_len(a) - 1
    if p < 1 {
        return peaks
    }
    let steps: Int = 512
    let stepsf: Float = 512.0
    let srf: Float = int_to_float(sr)
    let pi: Float = math_pi()
    let mag: [Float] = native_list_empty()
    let s: Int = 0
    while s < steps {
        let wq: Float = pi * int_to_float(s) / stepsf     // 0..pi -> 0..sr/2
        let re: Float = 0.0
        let im: Float = 0.0
        let k: Int = 0
        while k <= p {
            let ak: Float = native_list_get(a, k)
            let ang: Float = wq * int_to_float(k)
            let cv: Float = math_cos(ang)
            let sv: Float = math_sin(ang)
            re = re + ak * cv
            im = im - ak * sv
            k = k + 1
        }
        let d: Float = math_sqrt(re * re + im * im)
        if d < 0.000000001 {
            d = 0.000000001
        }
        mag = native_list_append(mag, 1.0 / d)
        s = s + 1
    }
    let found: Int = 0
    let t: Int = 1
    while t < steps - 1 {
        if found < 5 {
            let m0: Float = native_list_get(mag, t - 1)
            let m1: Float = native_list_get(mag, t)
            let m2: Float = native_list_get(mag, t + 1)
            let rise: Bool = m1 > m0
            let fall: Bool = m1 >= m2
            if rise {
                if fall {
                    let f: Float = int_to_float(t) * srf / 2.0 / stepsf
                    if f > 150.0 {
                        if f < 5200.0 {
                            let thr: Float = m1 / 1.4142
                            let lo: Int = t
                            let scan_lo: Bool = true
                            while scan_lo {
                                if lo > 0 {
                                    let mv: Float = native_list_get(mag, lo)
                                    if mv > thr {
                                        lo = lo - 1
                                    } else {
                                        scan_lo = false
                                    }
                                } else {
                                    scan_lo = false
                                }
                            }
                            let hi: Int = t
                            let scan_hi: Bool = true
                            while scan_hi {
                                if hi < steps - 1 {
                                    let mv2: Float = native_list_get(mag, hi)
                                    if mv2 > thr {
                                        hi = hi + 1
                                    } else {
                                        scan_hi = false
                                    }
                                } else {
                                    scan_hi = false
                                }
                            }
                            let bw: Float = int_to_float(hi - lo) * srf / 2.0 / stepsf
                            peaks = native_list_append(peaks, f)
                            peaks = native_list_append(peaks, bw)
                            found = found + 1
                        }
                    }
                }
            }
        }
        t = t + 1
    }
    return peaks
}

// Autocorrelation pitch over 70-400 Hz with the 0.30 voicing threshold.
// Returns 0.0 for an unvoiced (or silent) frame.
fn dsp_pitch_of(frame: [Float], sr: Int) -> Float {
    let n: Int = native_list_len(frame)
    let lag_min: Int = sr / 400
    let lag_max: Int = sr / 70
    if lag_max > n - 1 {
        lag_max = n - 1
    }
    if lag_max <= lag_min {
        return 0.0
    }
    let r0: Float = 0.0
    let i: Int = 0
    while i < n {
        let v: Float = native_list_get(frame, i)
        r0 = r0 + v * v
        i = i + 1
    }
    if r0 < 0.00001 {
        return 0.0
    }
    let best_lag: Int = 0
    let best: Float = 0.0
    let lag: Int = lag_min
    while lag <= lag_max {
        let c: Float = 0.0
        let j: Int = lag
        while j < n {
            let a: Float = native_list_get(frame, j)
            let b: Float = native_list_get(frame, j - lag)
            c = c + a * b
            j = j + 1
        }
        if c > best {
            best = c
            best_lag = lag
        }
        lag = lag + 1
    }
    if best_lag > 0 {
        let ratio: Float = best / r0
        if ratio > 0.30 {
            let srf: Float = int_to_float(sr)
            return srf / int_to_float(best_lag)
        }
    }
    return 0.0
}

// A copy of x[pos..pos+n) — the El stand-in for Swift's Array(x[a..<b]).
fn dsp_slice(x: [Float], pos: Int, n: Int) -> [Float] {
    let out: [Float] = native_list_empty()
    let total: Int = native_list_len(x)
    let i: Int = 0
    while i < n {
        if pos + i < total {
            let v: Float = native_list_get(x, pos + i)
            out = native_list_append(out, v)
        }
        i = i + 1
    }
    return out
}

// median == sorted()[count/2]. There is no list-set in El, so instead of
// sorting we find the value whose rank bracket contains index count/2 —
// numerically identical to the Swift expression, without a mutable buffer.
fn dsp_median(v: [Float]) -> Float {
    let n: Int = native_list_len(v)
    if n == 0 {
        return 0.0
    }
    let target: Int = n / 2
    let i: Int = 0
    while i < n {
        let x: Float = native_list_get(v, i)
        let less: Int = 0
        let eq: Int = 0
        let j: Int = 0
        while j < n {
            let y: Float = native_list_get(v, j)
            if y < x {
                less = less + 1
            }
            if y == x {
                eq = eq + 1
            }
            j = j + 1
        }
        if less <= target {
            if target < less + eq {
                return x
            }
        }
        i = i + 1
    }
    return 0.0
}

fn dsp_min_of(v: [Float]) -> Float {
    let n: Int = native_list_len(v)
    if n == 0 {
        return 0.0
    }
    let m: Float = native_list_get(v, 0)
    let i: Int = 1
    while i < n {
        let x: Float = native_list_get(v, i)
        if x < m {
            m = x
        }
        i = i + 1
    }
    return m
}

fn dsp_max_of(v: [Float]) -> Float {
    let n: Int = native_list_len(v)
    if n == 0 {
        return 0.0
    }
    let m: Float = native_list_get(v, 0)
    let i: Int = 1
    while i < n {
        let x: Float = native_list_get(v, i)
        if x > m {
            m = x
        }
        i = i + 1
    }
    return m
}

// ---------------------------------------------------------------------------
// voiceprint — the voice-signature: median F0 over voiced frames + the median
// of each formant F1..F5 (and its bandwidth), plus the F0 min/max range.
// FRAME = 400 (25 ms @ 16k), HOP = 160 (10 ms), LPC order 16.
//
// Returns PACKED: [ f0med, f0lo, f0hi, nf, f1, b1, f2, b2, ... ] where nf is
// how many formants were recovered (0..5). Empty list if unreadable.
// ---------------------------------------------------------------------------
fn dsp_voiceprint(path: String) -> [Float] {
    let out: [Float] = native_list_empty()
    let w: [Float] = dsp_read_wav(path)
    let n: Int = dsp_wav_n(w)
    let frame_len: Int = 400
    let hop: Int = 160
    let order: Int = 16
    if n <= frame_len {
        return out
    }
    let sr: Int = dsp_wav_sr(w)
    let x: [Float] = dsp_wav_pcm(w)

    let f0s: [Float] = native_list_empty()
    // five formant banks + five bandwidth banks, flat lists each
    let f1s: [Float] = native_list_empty()
    let f2s: [Float] = native_list_empty()
    let f3s: [Float] = native_list_empty()
    let f4s: [Float] = native_list_empty()
    let f5s: [Float] = native_list_empty()
    let b1s: [Float] = native_list_empty()
    let b2s: [Float] = native_list_empty()
    let b3s: [Float] = native_list_empty()
    let b4s: [Float] = native_list_empty()
    let b5s: [Float] = native_list_empty()

    let pos: Int = 0
    while pos + frame_len <= n {
        let raw: [Float] = dsp_slice(x, pos, frame_len)
        let f0: Float = dsp_pitch_of(raw, sr)
        if f0 > 0.0 {
            f0s = native_list_append(f0s, f0)
            let win: [Float] = dsp_hamming(raw)
            let r: [Float] = dsp_autocorr(win, order)
            let r0: Float = native_list_get(r, 0)
            if r0 > 0.000001 {
                let al: [Float] = dsp_levinson(r, order)
                // strip the trailing residual energy: coeffs are a[0..order]
                let a: [Float] = native_list_empty()
                let ci: Int = 0
                while ci <= order {
                    let av: Float = native_list_get(al, ci)
                    a = native_list_append(a, av)
                    ci = ci + 1
                }
                let fs: [Float] = dsp_formants(a, sr)
                let nfs: Int = native_list_len(fs) / 2
                if nfs > 0 {
                    let fv: Float = native_list_get(fs, 0)
                    let bv: Float = native_list_get(fs, 1)
                    f1s = native_list_append(f1s, fv)
                    b1s = native_list_append(b1s, bv)
                }
                if nfs > 1 {
                    let fv2: Float = native_list_get(fs, 2)
                    let bv2: Float = native_list_get(fs, 3)
                    f2s = native_list_append(f2s, fv2)
                    b2s = native_list_append(b2s, bv2)
                }
                if nfs > 2 {
                    let fv3: Float = native_list_get(fs, 4)
                    let bv3: Float = native_list_get(fs, 5)
                    f3s = native_list_append(f3s, fv3)
                    b3s = native_list_append(b3s, bv3)
                }
                if nfs > 3 {
                    let fv4: Float = native_list_get(fs, 6)
                    let bv4: Float = native_list_get(fs, 7)
                    f4s = native_list_append(f4s, fv4)
                    b4s = native_list_append(b4s, bv4)
                }
                if nfs > 4 {
                    let fv5: Float = native_list_get(fs, 8)
                    let bv5: Float = native_list_get(fs, 9)
                    f5s = native_list_append(f5s, fv5)
                    b5s = native_list_append(b5s, bv5)
                }
            }
        }
        pos = pos + hop
    }

    let f0med: Float = dsp_median(f0s)
    let f0lo: Float = dsp_min_of(f0s)
    let f0hi: Float = dsp_max_of(f0s)
    out = native_list_append(out, f0med)
    out = native_list_append(out, f0lo)
    out = native_list_append(out, f0hi)

    let nf: Int = 0
    if native_list_len(f1s) > 0 {
        nf = nf + 1
    }
    if native_list_len(f2s) > 0 {
        nf = nf + 1
    }
    if native_list_len(f3s) > 0 {
        nf = nf + 1
    }
    if native_list_len(f4s) > 0 {
        nf = nf + 1
    }
    if native_list_len(f5s) > 0 {
        nf = nf + 1
    }
    let nff: Float = int_to_float(nf)
    out = native_list_append(out, nff)
    if native_list_len(f1s) > 0 {
        let m: Float = dsp_median(f1s)
        let b: Float = dsp_median(b1s)
        out = native_list_append(out, m)
        out = native_list_append(out, b)
    }
    if native_list_len(f2s) > 0 {
        let m2: Float = dsp_median(f2s)
        let bb2: Float = dsp_median(b2s)
        out = native_list_append(out, m2)
        out = native_list_append(out, bb2)
    }
    if native_list_len(f3s) > 0 {
        let m3: Float = dsp_median(f3s)
        let bb3: Float = dsp_median(b3s)
        out = native_list_append(out, m3)
        out = native_list_append(out, bb3)
    }
    if native_list_len(f4s) > 0 {
        let m4: Float = dsp_median(f4s)
        let bb4: Float = dsp_median(b4s)
        out = native_list_append(out, m4)
        out = native_list_append(out, bb4)
    }
    if native_list_len(f5s) > 0 {
        let m5: Float = dsp_median(f5s)
        let bb5: Float = dsp_median(b5s)
        out = native_list_append(out, m5)
        out = native_list_append(out, bb5)
    }
    return out
}

// ---------------------------------------------------------------------------
// imitate — LPC analysis-resynthesis. Per frame: autocorr + Levinson, gain =
// sqrt(residual energy), excitation = an energy-normalized glottal impulse
// train at F0 for voiced frames or white noise for unvoiced, run through the
// all-pole filter using the past-output state. The whole output is normalized
// to peak 0.9 and returned as int16 samples — hand them to speech.el's
// write_wav(samples, sr, path).
// ---------------------------------------------------------------------------
fn dsp_imitate(path: String) -> [Int] {
    let res: [Int] = native_list_empty()
    let w: [Float] = dsp_read_wav(path)
    let n: Int = dsp_wav_n(w)
    let frame_len: Int = 400
    let hop: Int = 160
    let order: Int = 16
    if n <= frame_len {
        return res
    }
    let sr: Int = dsp_wav_sr(w)
    let srf: Float = int_to_float(sr)
    let x: [Float] = dsp_wav_pcm(w)

    let out: [Float] = native_list_empty()
    // past outputs, order deep
    let state: [Float] = native_list_empty()
    let si: Int = 0
    while si < order {
        state = native_list_append(state, 0.0)
        si = si + 1
    }
    let phase: Float = 0.0
    let last_f0: Float = 0.0
    let nstate: Int = 22695            // LCG for the unvoiced source
    let written: Int = 0

    let pos: Int = 0
    while pos + frame_len <= n {
        let raw: [Float] = dsp_slice(x, pos, frame_len)
        let win: [Float] = dsp_hamming(raw)
        let r: [Float] = dsp_autocorr(win, order)
        let f0: Float = dsp_pitch_of(raw, sr)
        let r0: Float = native_list_get(r, 0)
        if r0 < 0.0000001 {
            // frame skipped in periph.swift -> those output samples stay zero
            let z: Int = 0
            while z < hop {
                if written < n {
                    out = native_list_append(out, 0.0)
                    written = written + 1
                }
                z = z + 1
            }
        } else {
            let al: [Float] = dsp_levinson(r, order)
            let errv: Float = native_list_get(al, order + 1)
            let ge: Float = errv
            if ge < 0.0 {
                ge = 0.0
            }
            let gain: Float = math_sqrt(ge)
            let use_f0: Float = f0
            if f0 <= 0.0 {
                use_f0 = 0.0
                if last_f0 > 0.0 {
                    use_f0 = last_f0
                }
            }
            last_f0 = f0
            let i: Int = 0
            while i < hop {
                if written < n {
                    let e: Float = 0.0
                    if use_f0 > 0.0 {
                        phase = phase + use_f0 / srf
                        if phase >= 1.0 {
                            phase = phase - 1.0
                            e = math_sqrt(srf / use_f0)
                        }
                    } else {
                        nstate = nstate * 1103515245 + 12345
                        nstate = nstate - (nstate / 2147483648) * 2147483648
                        if nstate < 0 {
                            nstate = 0 - nstate
                        }
                        let u: Float = int_to_float(nstate) / 2147483648.0
                        e = u * 2.0 - 1.0
                    }
                    let y: Float = gain * e
                    let k: Int = 1
                    while k <= order {
                        let ak: Float = native_list_get(al, k)
                        let sk: Float = native_list_get(state, k - 1)
                        y = y - ak * sk
                        k = k + 1
                    }
                    let ns: [Float] = native_list_empty()
                    ns = native_list_append(ns, y)
                    let m: Int = 0
                    while m < order - 1 {
                        let sv: Float = native_list_get(state, m)
                        ns = native_list_append(ns, sv)
                        m = m + 1
                    }
                    state = ns
                    out = native_list_append(out, y)
                    written = written + 1
                }
                i = i + 1
            }
        }
        pos = pos + hop
    }
    // the tail past the last full frame stays silent, exactly as in periph.swift
    while written < n {
        out = native_list_append(out, 0.0)
        written = written + 1
    }

    // normalize to peak 0.9, then to int16
    let peak: Float = 0.0
    let q: Int = 0
    while q < n {
        let v: Float = native_list_get(out, q)
        let av: Float = dsp_fabs(v)
        if av > peak {
            peak = av
        }
        q = q + 1
    }
    let scale: Float = 1.0
    if peak > 0.000000001 {
        scale = 0.9 / peak
    }
    let t: Int = 0
    while t < n {
        let v2: Float = native_list_get(out, t)
        let sv2: Float = v2 * scale * 32767.0
        if sv2 > 32767.0 {
            sv2 = 32767.0
        }
        if sv2 < 0.0 - 32767.0 {
            sv2 = 0.0 - 32767.0
        }
        let iv: Int = float_to_int(sv2)
        res = native_list_append(res, iv)
        t = t + 1
    }
    return res
}