spectral.kt

src/jvmMain/kotlin/afield/audio/analysis/spectral.kt · 914 lines

@file:afield.runtime.completion.CompletionCategory("Math")

package afield.audio.analysis

import afield.data.collections.get
import afield.data.collections.invoke
import afield.data.collections.times
import afield.data.graph.main
import afield.document.context.Why
import afield.document.graph.minusAssign
import afield.document.graph.plusAssign
import afield.document.properties.compareTo
import afield.document.workspace.file
import afield.geometry.mesh.normalized
import afield.geometry.mesh.scaled
import afield.graphics.render.rangeTo
import afield.graphics.vulkan.buffers.add
import afield.io.osc.main
import afield.math.autodiff.compareTo
import afield.math.autodiff.div
import afield.math.autodiff.log2
import afield.math.autodiff.mean
import afield.math.autodiff.minus
import afield.math.autodiff.plus
import afield.math.autodiff.sign
import afield.math.autodiff.sum
import afield.math.autodiff.times
import afield.math.control.main
import afield.math.linalg.invoke
import afield.math.linalg.remAssign
import afield.math.linalg.set
import afield.math.linalg.times
import afield.math.linalg.x
import afield.math.statistics.main
import afield.math.statistics.noise
import afield.math.units.div
import afield.math.units.minus
import afield.math.units.plus
import afield.math.units.times
import afield.media.image.normalized
import afield.media.image.pad
import afield.simulation.agents.step
import afield.simulation.automata.main
import afield.simulation.automata.step
import afield.study.ui.numbers

import afield.audio.synth.*
import afield.data.collections.m
import afield.data.collections.s
import afield.geometry.mesh.centroid
import afield.graphics.render.filled
import afield.math.autodiff.atan2
import afield.math.autodiff.cos
import afield.math.autodiff.exp
import afield.math.autodiff.hypot
import afield.math.autodiff.ln
import afield.math.autodiff.log10
import afield.math.autodiff.max
import afield.math.autodiff.min
import afield.math.autodiff.norm
import afield.math.autodiff.pow
import afield.math.autodiff.sin
import afield.math.autodiff.sqrt
import afield.math.functions.atan2
import afield.math.functions.cos
import afield.math.functions.sin
import afield.math.linalg.ar
import afield.math.linalg.br
import afield.math.linalg.d
import afield.math.linalg.f
import afield.math.linalg.i
import afield.math.linalg.l
import afield.math.statistics.FloatFilter
import afield.runtime.execution.start
import kotlin.math.*

Frequency-domain / spectral analysis family for the stdlib DSP toolkit.

This is the frequency-domain companion to timeseries.kt (which is entirely time-domain): a pure-Kotlin FFT, analysis windows, a Spectrum and a Spectrogram (STFT), classic spectral features (centroid / spread / rolloff / flux / flatness / RMS / energy), and streaming FloatFilter wrappers so those features drop straight into the time-domain filter chain in timeseries.kt (compose with then).

Everything is plain Double/Float math — no external DSP dependency — so it works headless in a box, a test, or a docs capture. The FFT is an iterative Cooley–Tukey radix-2 core with a Bluestein (chirp-z) fallback so arbitrary (non-power-of-two) lengths transform exactly.

FFT core — iterative Cooley–Tukey radix-2 with a Bluestein (chirp-z) fallback.

Conventions: forward Xk = Σ_n xn · exp(-2πi·kn/N) inverse xn = (1/N) Σ_k Xk · exp(+2πi·kn/N) The radix-2 kernel is the classic public-domain in-place decimation-in-time transform; the Bluestein routine reduces an arbitrary-length DFT to a power-of-two convolution.

object FFT

The transform itself: an in-place, self-contained complex FFT that accepts any length, not just a power of two.

A power-of-two length goes through the fast iterative radix-2 kernel; anything else is handled by a Bluestein (chirp-z) reduction to a power-of-two convolution, which is correct but allocates several arrays of the next power of two above 2n+1 and is markedly slower. So a length of 1024 is cheap and a length of 1000 is not — pad to a power of two when you can.

The sign convention is the usual one: forward carries exp(-2πi·kn/N) and no scaling, inverse carries exp(+2πi·kn/N) and the 1/N. Most callers want Spectrum and the stft family rather than this directly.

object FFT {

fun isPowerOfTwo(n: Int): Boolean = n > 0 && (n and (n - 1)) == 0

True for 1, 2, 4, 8, … — and false for zero and negatives, so it is safe on unvalidated input.

fun nextPowerOfTwo(n: Int): Int

Smallest power of two >= n (>= 1).

    fun nextPowerOfTwo(n: Int): Int {
        var m = 1
        while (m < n) m = m shl 1
        return m
    }

fun transform(re: DoubleArray, im: DoubleArray)

Forward FFT in place. re/im must be the same length; any length >= 1 is accepted (power-of-two goes through the radix-2 kernel, otherwise Bluestein).

    fun transform(re: DoubleArray, im: DoubleArray) {
        require(re.size == im.size) { "real/imag arrays must match length" }
        val n = re.size
        if (n == 0) return
        if (isPowerOfTwo(n)) transformRadix2(re, im) else transformBluestein(re, im)
    }

fun inverseTransform(re: DoubleArray, im: DoubleArray)

Inverse FFT in place, normalized (divides by N). By the standard identity the inverse is the forward transform with real/imag swapped, then scaled by 1/N.

    fun inverseTransform(re: DoubleArray, im: DoubleArray) {
        val n = re.size
        if (n == 0) return
        transform(im, re) // swap -> unnormalized inverse
        val inv = 1.0 / n
        for (i in 0 until n) {
            re[i] *= inv
            im[i] *= inv
        }
    }

fun transformRadix2(re: DoubleArray, im: DoubleArray)

In-place radix-2 decimation-in-time forward transform; requires power-of-two length.

    fun transformRadix2(re: DoubleArray, im: DoubleArray) {
        val n = re.size
        require(isPowerOfTwo(n)) { "radix-2 requires a power-of-two length, got $n" }
        if (n == 1) return
        val levels = Integer.numberOfTrailingZeros(n) // = log2(n)

        // Trig lookup tables: cos/sin of 2πk/n for k in 0 until n/2.
        val half = n / 2
        val cosTable = DoubleArray(half)
        val sinTable = DoubleArray(half)
        for (i in 0 until half) {
            val angle = 2.0 * PI * i / n
            cosTable[i] = cos(angle)
            sinTable[i] = sin(angle)
        }

        // Bit-reversed addressing.
        for (i in 0 until n) {
            val j = Integer.reverse(i) ushr (32 - levels)
            if (j > i) {
                var t = re[i]; re[i] = re[j]; re[j] = t
                t = im[i]; im[i] = im[j]; im[j] = t
            }
        }

        // Butterflies.
        var size = 2
        while (size <= n) {
            val halfsize = size / 2
            val tablestep = n / size
            var i = 0
            while (i < n) {
                var j = i
                var k = 0
                while (j < i + halfsize) {
                    val l = j + halfsize
                    val tpre = re[l] * cosTable[k] + im[l] * sinTable[k]
                    val tpim = -re[l] * sinTable[k] + im[l] * cosTable[k]
                    re[l] = re[j] - tpre
                    im[l] = im[j] - tpim
                    re[j] += tpre
                    im[j] += tpim
                    j++
                    k += tablestep
                }
                i += size
            }
            if (size == n) break
            size = size shl 1
        }
    }

fun transformBluestein(re: DoubleArray, im: DoubleArray)

In-place arbitrary-length forward transform via Bluestein's chirp-z algorithm.

    fun transformBluestein(re: DoubleArray, im: DoubleArray) {
        val n = re.size
        require(re.size == im.size)
        // Convolution length: a power of two >= 2n+1 so the cyclic convolution is acyclic.
        val m = nextPowerOfTwo(2 * n + 1)

        // Chirp trig tables: angle = π·(i²  mod 2n)/n.
        val cosTable = DoubleArray(n)
        val sinTable = DoubleArray(n)
        for (i in 0 until n) {
            val j = ((i.toLong() * i) % (2L * n)).toInt()
            val angle = PI * j / n
            cosTable[i] = cos(angle)
            sinTable[i] = sin(angle)
        }

        val areal = DoubleArray(m)
        val aimag = DoubleArray(m)
        for (i in 0 until n) {
            areal[i] = re[i] * cosTable[i] + im[i] * sinTable[i]
            aimag[i] = -re[i] * sinTable[i] + im[i] * cosTable[i]
        }

        val breal = DoubleArray(m)
        val bimag = DoubleArray(m)
        breal[0] = cosTable[0]
        bimag[0] = sinTable[0]
        for (i in 1 until n) {
            breal[i] = cosTable[i]; breal[m - i] = cosTable[i]
            bimag[i] = sinTable[i]; bimag[m - i] = sinTable[i]
        }

        val creal = DoubleArray(m)
        val cimag = DoubleArray(m)
        convolveComplex(areal, aimag, breal, bimag, creal, cimag)

        for (i in 0 until n) {
            re[i] = creal[i] * cosTable[i] + cimag[i] * sinTable[i]
            im[i] = -creal[i] * sinTable[i] + cimag[i] * cosTable[i]
        }
    }

Cyclic complex convolution of length-m power-of-two vectors, via radix-2 FFT.

    private fun convolveComplex(
        xre: DoubleArray, xim: DoubleArray,
        yre: DoubleArray, yim: DoubleArray,
        outRe: DoubleArray, outIm: DoubleArray,
    ) {
        val n = xre.size
        val ar = xre.copyOf(); val ai = xim.copyOf()
        val br = yre.copyOf(); val bi = yim.copyOf()
        transformRadix2(ar, ai)
        transformRadix2(br, bi)
        for (i in 0 until n) {
            val tr = ar[i] * br[i] - ai[i] * bi[i]
            ai[i] = ai[i] * br[i] + ar[i] * bi[i]
            ar[i] = tr
        }
        inverseTransform(ar, ai)
        for (i in 0 until n) {
            outRe[i] = ar[i]
            outIm[i] = ai[i]
        }
    }
}

Analysis windows.

enum class WindowType { RECTANGULAR, HANN, HAMMING, BLACKMAN }

Which analysis window to taper a frame with before transforming it — the trade between how sharply you can separate two nearby frequencies and how much a strong one smears into its neighbours.

RECTANGULAR is no window at all: the narrowest main lobe, and the worst leakage, so a loud tone spatters energy right across the spectrum. HANN is the sane default. HAMMING is similar with a slightly narrower main lobe and a lower first sidelobe (about -47 dB against Hann's -31 dB), paid for with a far-field skirt that decays more slowly. BLACKMAN suppresses leakage hardest at the cost of the widest lobe — reach for it when you need to see something quiet next to something loud.

object Windows

The window functions themselves, as plain DoubleArrays you multiply a frame by.

All but rectangular, which has nothing to taper, take a periodic flag that is false by default, meaning symmetric — the right choice for filter design and for looking at a single frame. Pass true for the periodic form, which is the correct one for STFT analysis where frames overlap and are expected to sum smoothly; the difference is one sample in the divisor and it matters more than it looks.

All of them return an all-ones array for n <= 1 rather than dividing by zero.

object Windows {

fun make(type: WindowType, n: Int, periodic: Boolean = false): DoubleArray

Build a length-n window of the given type. Symmetric (periodic=false) by default.

    fun make(type: WindowType, n: Int, periodic: Boolean = false): DoubleArray = when (type) {
        WindowType.RECTANGULAR -> rectangular(n)
        WindowType.HANN -> hann(n, periodic)
        WindowType.HAMMING -> hamming(n, periodic)
        WindowType.BLACKMAN -> blackman(n, periodic)
    }

fun rectangular(n: Int): DoubleArray = DoubleArray(n) { 1.0 }

All ones — i.e. no window. The periodic distinction does not apply, so there is no such parameter.

fun hann(n: Int, periodic: Boolean = false): DoubleArray

The raised-cosine window: 0.5 - 0.5·cos(2πi/d), zero at both ends, and the default worth reaching for. Symmetric unless periodic is true.

    fun hann(n: Int, periodic: Boolean = false): DoubleArray {
        if (n <= 1) return DoubleArray(n) { 1.0 }
        val denom = if (periodic) n.toDouble() else (n - 1).toDouble()
        return DoubleArray(n) { 0.5 - 0.5 * cos(2.0 * PI * it / denom) }
    }

fun hamming(n: Int, periodic: Boolean = false): DoubleArray

0.54 - 0.46·cos(2πi/d) — like hann but it does not reach zero at the ends, which is what buys its slightly better resolution and slightly worse far-field leakage.

    fun hamming(n: Int, periodic: Boolean = false): DoubleArray {
        if (n <= 1) return DoubleArray(n) { 1.0 }
        val denom = if (periodic) n.toDouble() else (n - 1).toDouble()
        return DoubleArray(n) { 0.54 - 0.46 * cos(2.0 * PI * it / denom) }
    }

fun blackman(n: Int, periodic: Boolean = false): DoubleArray

The three-term 0.42 - 0.5·cos(w) + 0.08·cos(2w) window: the strongest leakage suppression here, and the widest main lobe — about 1.5 times hann's null-to-null width, so two frequencies Hann would just separate blur into one here.

    fun blackman(n: Int, periodic: Boolean = false): DoubleArray {
        if (n <= 1) return DoubleArray(n) { 1.0 }
        val denom = if (periodic) n.toDouble() else (n - 1).toDouble()
        return DoubleArray(n) {
            val w = 2.0 * PI * it / denom
            0.42 - 0.5 * cos(w) + 0.08 * cos(2.0 * w)
        }
    }

fun coherentGain(window: DoubleArray): Double

Coherent gain (mean) of a window — used to compensate amplitude when a window is applied.

    fun coherentGain(window: DoubleArray): Double {
        if (window.isEmpty()) return 1.0
        var s = 0.0
        for (w in window) s += w
        return s / window.size
    }
}

Spectrum — the result of one FFT, with sample-rate-aware helpers.

class Spectrum

A single complex spectrum plus the sampleRate it was measured at. re/im are the full length-N transform; for a real input the useful ("one-sided") part is bins 0..N/2.

class Spectrum(
    val re: DoubleArray,
    val im: DoubleArray,
    val sampleRate: Double,
) {

val size: Int get() = re.size

val binCount: Int get() = size / 2 + 1

Number of one-sided bins (DC .. Nyquist), i.e. N/2 + 1 for even N.

fun binFrequency(k: Int): Double = k * sampleRate / size

Frequency (Hz) of bin k.

fun frequencies(): DoubleArray = DoubleArray(binCount) { binFrequency(it) }

Frequencies (Hz) of the one-sided bins 0..N/2.

fun magnitude(): DoubleArray = DoubleArray(size) { hypot(re[it], im[it]) }

Full-length magnitude |Xk|.

fun magnitudeOneSided(): DoubleArray = DoubleArray(binCount) { hypot(re[it], im[it]) }

One-sided magnitude (bins 0..N/2).

fun powerOneSided(): DoubleArray = DoubleArray(binCount) { re[it] * re[it] + im[it] * im[it] }

One-sided power |Xk|².

fun magnitudeDbOneSided(floorDb: Double = -120.0, ref: Double = 1.0): DoubleArray

One-sided magnitude in decibels, floored at floorDb.

    fun magnitudeDbOneSided(floorDb: Double = -120.0, ref: Double = 1.0): DoubleArray {
        val mags = magnitudeOneSided()
        return DoubleArray(mags.size) {
            val db = 20.0 * log10((mags[it] / ref).coerceAtLeast(1e-12))
            max(db, floorDb)
        }
    }

fun phase(): DoubleArray = DoubleArray(size) { atan2(im[it], re[it]) }

Phase (radians) per bin.

fun peakBin(): Int

Bin index of the peak one-sided magnitude.

    fun peakBin(): Int {
        val mags = magnitudeOneSided()
        var best = 0
        for (i in mags.indices) if (mags[i] > mags[best]) best = i
        return best
    }

fun peakFrequency(): Double = binFrequency(peakBin())

Frequency (Hz) of the dominant one-sided component.

fun centroid(): Double = SpectralFeatures.centroid(magnitudeOneSided(), sampleRate)

Feature shortcuts (see SpectralFeatures).

fun spread(): Double = SpectralFeatures.spread(magnitudeOneSided(), sampleRate)

fun rolloff(rolloffFraction: Double = 0.85): Double

    fun rolloff(rolloffFraction: Double = 0.85): Double =
        SpectralFeatures.rolloff(magnitudeOneSided(), sampleRate, rolloffFraction)

fun flatness(): Double = SpectralFeatures.flatness(magnitudeOneSided())

Forward / inverse convenience over real signals.

fun fft(signal: DoubleArray, sampleRate: Double = signal.size.toDouble()): Spectrum

Forward FFT of a real signal. sampleRate defaults to the signal length (bins == cycles/frame).

fun fft(signal: DoubleArray, sampleRate: Double = signal.size.toDouble()): Spectrum {
    val re = signal.copyOf()
    val im = DoubleArray(signal.size)
    FFT.transform(re, im)
    return Spectrum(re, im, sampleRate)
}

fun fft(signal: FloatArray, sampleRate: Double = signal.size.toDouble()): Spectrum

fun fft(signal: FloatArray, sampleRate: Double = signal.size.toDouble()): Spectrum =
    fft(DoubleArray(signal.size) { signal[it].toDouble() }, sampleRate)

fun fft(signal: List<Float>, sampleRate: Double = signal.size.toDouble()): Spectrum

fun fft(signal: List<Float>, sampleRate: Double = signal.size.toDouble()): Spectrum =
    fft(DoubleArray(signal.size) { signal[it].toDouble() }, sampleRate)

fun fft(signal: DoubleArray, window: WindowType, sampleRate: Double = signal.size.toDouble()): Spectrum

Windowed forward FFT: multiplies by window first (amplitude compensated by coherent gain).

fun fft(signal: DoubleArray, window: WindowType, sampleRate: Double = signal.size.toDouble()): Spectrum {
    val w = Windows.make(window, signal.size)
    val gain = Windows.coherentGain(w).coerceAtLeast(1e-9)
    val re = DoubleArray(signal.size) { signal[it] * w[it] / gain }
    val im = DoubleArray(signal.size)
    FFT.transform(re, im)
    return Spectrum(re, im, sampleRate)
}

fun ifft(spectrum: Spectrum): DoubleArray

Inverse FFT; returns the real part of the reconstructed signal.

fun ifft(spectrum: Spectrum): DoubleArray {
    val re = spectrum.re.copyOf()
    val im = spectrum.im.copyOf()
    FFT.inverseTransform(re, im)
    return re
}

fun ifftComplex(spectrum: Spectrum): Pair<DoubleArray, DoubleArray>

Inverse FFT keeping both real and imaginary parts.

fun ifftComplex(spectrum: Spectrum): Pair<DoubleArray, DoubleArray> {
    val re = spectrum.re.copyOf()
    val im = spectrum.im.copyOf()
    FFT.inverseTransform(re, im)
    return re to im
}

Spectrogram (STFT).

class Spectrogram

The result of a short-time Fourier transform: a stack of per-frame one-sided magnitude spectra. magnitudesframe``, bin in 0..fftSize/2. Complex frames are kept too so an inverse-STFT overlap-add reconstruction is possible (istft).

class Spectrogram(
    /** Complex STFT frames (full length-fftSize), kept for reconstruction. */
    val framesRe: Array<DoubleArray>,
    val framesIm: Array<DoubleArray>,
    val sampleRate: Double,
    val fftSize: Int,
    val hopSize: Int,
    val window: WindowType,
) {

val numFrames: Int get() = framesRe.size

val numBins: Int get() = fftSize / 2 + 1

fun binFrequency(k: Int): Double = k * sampleRate / fftSize

fun frameTime(i: Int): Double = (i * hopSize + fftSize / 2.0) / sampleRate

Center time (seconds) of frame i.

fun magnitude(i: Int): DoubleArray

One-sided magnitude of frame i.

    fun magnitude(i: Int): DoubleArray {
        val re = framesRe[i]; val im = framesIm[i]
        return DoubleArray(numBins) { hypot(re[it], im[it]) }
    }

fun magnitudes(): Array<DoubleArray> = Array(numFrames) { magnitude(it) }

All frames as one-sided magnitudes: frame.

fun maxMagnitude(): Double

Largest magnitude anywhere in the spectrogram (for normalizing a heatmap).

    fun maxMagnitude(): Double {
        var m = 0.0
        for (f in 0 until numFrames) {
            val re = framesRe[f]; val im = framesIm[f]
            for (b in 0 until numBins) m = max(m, hypot(re[b], im[b]))
        }
        return m
    }

fun centroidTrace(): DoubleArray

Per-frame spectral centroid (Hz).

    fun centroidTrace(): DoubleArray =
        DoubleArray(numFrames) { SpectralFeatures.centroid(magnitude(it), sampleRate) }

fun rolloffTrace(rolloffFraction: Double = 0.85): DoubleArray

Per-frame spectral rolloff (Hz).

    fun rolloffTrace(rolloffFraction: Double = 0.85): DoubleArray =
        DoubleArray(numFrames) { SpectralFeatures.rolloff(magnitude(it), sampleRate, rolloffFraction) }

fun fluxTrace(): DoubleArray

Per-frame spectral flux (positive change vs previous frame).

    fun fluxTrace(): DoubleArray {
        val out = DoubleArray(numFrames)
        var prev: DoubleArray? = null
        for (i in 0 until numFrames) {
            val cur = magnitude(i)
            val p = prev
            out[i] = if (p == null) 0.0 else SpectralFeatures.flux(p, cur)
            prev = cur
        }
        return out
    }
}

fun stft

Short-time Fourier transform. Frames of fftSize samples are windowed by window and advanced by hopSize; the tail is zero-padded so the final partial frame is kept.

fun stft(
    signal: DoubleArray,
    fftSize: Int = 1024,
    hopSize: Int = fftSize / 4,
    window: WindowType = WindowType.HANN,
    sampleRate: Double = signal.size.toDouble(),
): Spectrogram {
    require(fftSize > 0) { "fftSize must be positive" }
    require(hopSize > 0) { "hopSize must be positive" }
    val w = Windows.make(window, fftSize)
    val numFrames = if (signal.isEmpty()) 0 else max(1, (signal.size + hopSize - 1) / hopSize)
    val framesRe = Array(numFrames) { DoubleArray(fftSize) }
    val framesIm = Array(numFrames) { DoubleArray(fftSize) }
    for (f in 0 until numFrames) {
        val start = f * hopSize
        val re = framesRe[f]
        val im = framesIm[f]
        for (i in 0 until fftSize) {
            val idx = start + i
            re[i] = if (idx < signal.size) signal[idx] * w[i] else 0.0
        }
        FFT.transform(re, im)
    }
    return Spectrogram(framesRe, framesIm, sampleRate, fftSize, hopSize, window)
}

fun stft

fun stft(
    signal: FloatArray,
    fftSize: Int = 1024,
    hopSize: Int = fftSize / 4,
    window: WindowType = WindowType.HANN,
    sampleRate: Double = signal.size.toDouble(),
): Spectrogram = stft(DoubleArray(signal.size) { signal[it].toDouble() }, fftSize, hopSize, window, sampleRate)

fun stft

fun stft(
    signal: List<Float>,
    fftSize: Int = 1024,
    hopSize: Int = fftSize / 4,
    window: WindowType = WindowType.HANN,
    sampleRate: Double = signal.size.toDouble(),
): Spectrogram = stft(DoubleArray(signal.size) { signal[it].toDouble() }, fftSize, hopSize, window, sampleRate)

fun istft(spec: Spectrogram): DoubleArray

Inverse STFT — weighted overlap-add reconstruction of a real signal from a Spectrogram. Each frame is inverse-transformed, re-windowed (synthesis window == analysis window) and overlap-added; the running window-energy is divided out so a Hann-windowed COLA STFT reconstructs faithfully.

fun istft(spec: Spectrogram): DoubleArray {
    val fftSize = spec.fftSize
    val hop = spec.hopSize
    val w = Windows.make(spec.window, fftSize)
    val length = if (spec.numFrames == 0) 0 else (spec.numFrames - 1) * hop + fftSize
    val out = DoubleArray(length)
    val norm = DoubleArray(length)
    for (f in 0 until spec.numFrames) {
        val re = spec.framesRe[f].copyOf()
        val im = spec.framesIm[f].copyOf()
        FFT.inverseTransform(re, im)
        val start = f * hop
        for (i in 0 until fftSize) {
            out[start + i] += re[i] * w[i]
            norm[start + i] += w[i] * w[i]
        }
    }
    for (i in out.indices) if (norm[i] > 1e-9) out[i] /= norm[i]
    return out
}

Spectral features — over a one-sided magnitude spectrum unless noted.

object SpectralFeatures

Single numbers that describe the shape of a spectrum — brightness, spread, where the energy stops — for driving something from a sound without looking at every bin.

Most of it takes a one-sided magnitude spectrum (as from Spectrum.magnitudeOneSided) plus the sample rate, and returns Hz where the answer is a frequency — but not all: flux and flatness need no sample rate, and rms and energy take a time-domain buffer instead of a spectrum. The bin width is derived as sampleRate / (2·(n-1)), so passing a two-sided spectrum by mistake silently halves every frequency it reports rather than failing. An empty spectrum gives 0.0 throughout instead of NaN.

object SpectralFeatures {

fun centroid(magnitude: DoubleArray, sampleRate: Double): Double

Spectral centroid (Hz): magnitude-weighted mean frequency.

    fun centroid(magnitude: DoubleArray, sampleRate: Double): Double {
        val n = magnitude.size
        if (n == 0) return 0.0
        val binHz = sampleRate / (2.0 * (n - 1).coerceAtLeast(1))
        var num = 0.0
        var den = 0.0
        for (k in 0 until n) {
            num += k * binHz * magnitude[k]
            den += magnitude[k]
        }
        return if (den <= 0.0) 0.0 else num / den
    }

fun spread(magnitude: DoubleArray, sampleRate: Double): Double

Spectral spread / bandwidth (Hz): magnitude-weighted std-dev of frequency about the centroid.

    fun spread(magnitude: DoubleArray, sampleRate: Double): Double {
        val n = magnitude.size
        if (n == 0) return 0.0
        val binHz = sampleRate / (2.0 * (n - 1).coerceAtLeast(1))
        val c = centroid(magnitude, sampleRate)
        var num = 0.0
        var den = 0.0
        for (k in 0 until n) {
            val f = k * binHz
            num += (f - c) * (f - c) * magnitude[k]
            den += magnitude[k]
        }
        return if (den <= 0.0) 0.0 else sqrt(num / den)
    }

fun rolloff(magnitude: DoubleArray, sampleRate: Double, rolloffFraction: Double = 0.85): Double

Spectral rolloff (Hz): the frequency below which rolloffFraction of the total spectral energy (Σ magnitude) lies. Classic timbre / brightness descriptor.

    fun rolloff(magnitude: DoubleArray, sampleRate: Double, rolloffFraction: Double = 0.85): Double {
        val n = magnitude.size
        if (n == 0) return 0.0
        val binHz = sampleRate / (2.0 * (n - 1).coerceAtLeast(1))
        var total = 0.0
        for (m in magnitude) total += m
        if (total <= 0.0) return 0.0
        val target = total * rolloffFraction.coerceIn(0.0, 1.0)
        var acc = 0.0
        for (k in 0 until n) {
            acc += magnitude[k]
            if (acc >= target) return k * binHz
        }
        return (n - 1) * binHz
    }

fun flux(previous: DoubleArray, current: DoubleArray): Double

Spectral flux: the (positive, L2) change between two consecutive magnitude spectra. A standard onset / novelty descriptor. Only positive differences are summed (half-wave rectified).

    fun flux(previous: DoubleArray, current: DoubleArray): Double {
        val n = min(previous.size, current.size)
        var acc = 0.0
        for (k in 0 until n) {
            val d = current[k] - previous[k]
            if (d > 0.0) acc += d * d
        }
        return sqrt(acc)
    }

fun flatness(magnitude: DoubleArray): Double

Spectral flatness (Wiener entropy): geometric mean / arithmetic mean of the power spectrum. ~1 for noise-like (flat) spectra, ->0 for tonal (peaky) spectra.

    fun flatness(magnitude: DoubleArray): Double {
        val n = magnitude.size
        if (n == 0) return 0.0
        var logSum = 0.0
        var arithSum = 0.0
        var count = 0
        for (m in magnitude) {
            val p = m * m
            if (p > 0.0) {
                logSum += ln(p)
                arithSum += p
                count++
            }
        }
        if (count == 0 || arithSum <= 0.0) return 0.0
        val geo = exp(logSum / count)
        val arith = arithSum / count
        return (geo / arith).coerceIn(0.0, 1.0)
    }

fun rms(signal: DoubleArray): Double

Root-mean-square level of a time-domain buffer.

    fun rms(signal: DoubleArray): Double {
        if (signal.isEmpty()) return 0.0
        var s = 0.0
        for (v in signal) s += v * v
        return sqrt(s / signal.size)
    }

fun rms(signal: FloatArray): Double

Root-mean-square level of a time-domain buffer — loudness, roughly, and the number to drive a meter from. 0.0 for an empty buffer.

A FloatArray convenience over the DoubleArray overload; the accumulation is in Double regardless, so a long buffer does not lose precision.

    fun rms(signal: FloatArray): Double {
        if (signal.isEmpty()) return 0.0
        var s = 0.0
        for (v in signal) s += v.toDouble() * v
        return sqrt(s / signal.size)
    }

fun energy(signal: DoubleArray): Double

Total energy Σ x² of a time-domain buffer.

    fun energy(signal: DoubleArray): Double {
        var s = 0.0
        for (v in signal) s += v * v
        return s
    }
}

Mel scale — a perceptual log-frequency warping + triangular filter bank (bonus).

object Mel

The mel scale and its triangular filter bank — frequency re-spaced the way hearing actually works, so that a step of one band means about as much low down as high up.

Why bother: a linear FFT spends most of its bins on the top few octaves, where we hear the least detail. Collapsing to mel bands gives a compact vector that changes when a sound changes audibly, which is what you want for driving visuals or comparing two sounds. Use filterBank on a one-sided magnitude spectrum; hzToMel and melToHz are the scale itself, in the standard 2595·log10(1 + hz/700) form.

object Mel {

fun hzToMel(hz: Double): Double = 2595.0 * log10(1.0 + hz / 700.0)

Hz to mels, 2595·log10(1 + hz/700). Monotonic and 0 at 0 Hz.

Defined down to but not including -700 Hz: exactly -700 gives -infinity and anything below it NaN, which is the log's domain rather than a guard this adds.

fun melToHz(mel: Double): Double = 700.0 * (10.0.pow(mel / 2595.0) - 1.0)

Mels back to Hz — the exact inverse of hzToMel.

fun filterBank

Apply a triangular mel filter bank to a one-sided magnitude spectrum, returning numBands mel-band energies spanning minHz..maxHz.

    fun filterBank(
        magnitude: DoubleArray,
        sampleRate: Double,
        numBands: Int = 26,
        minHz: Double = 0.0,
        maxHz: Double = sampleRate / 2.0,
    ): DoubleArray {
        val nBins = magnitude.size
        if (nBins == 0 || numBands <= 0) return DoubleArray(numBands.coerceAtLeast(0))
        val binHz = sampleRate / (2.0 * (nBins - 1).coerceAtLeast(1))
        val melMin = hzToMel(minHz)
        val melMax = hzToMel(maxHz)
        // numBands+2 equally-spaced mel points -> triangular band edges.
        val points = DoubleArray(numBands + 2) { melToHz(melMin + (melMax - melMin) * it / (numBands + 1)) }
        val out = DoubleArray(numBands)
        for (b in 0 until numBands) {
            val lo = points[b]; val mid = points[b + 1]; val hi = points[b + 2]
            var acc = 0.0
            for (k in 0 until nBins) {
                val f = k * binHz
                val weight = when {
                    f < lo || f > hi -> 0.0
                    f <= mid -> if (mid > lo) (f - lo) / (mid - lo) else 0.0
                    else -> if (hi > mid) (hi - f) / (hi - mid) else 0.0
                }
                acc += weight * magnitude[k]
            }
            out[b] = acc
        }
        return out
    }
}

Streaming spectral FloatFilters — one Float in per sample, one Float out, so they compose with the time-domain chain in timeseries.kt (someFloatFilter then SpectralCentroidFilter(...)).

Each maintains a ring buffer of the last fftSize samples; every hopSize samples it recomputes its feature and holds that value between hops (sample-and-hold), so the output is a smooth per-sample stream suitable for driving a control signal or an FLine.

abstract class SpectralFrameFilter

Base class: buffers samples, fires onFrame every hop with the current windowed frame.

abstract class SpectralFrameFilter(
    val fftSize: Int = 1024,
    val hopSize: Int = fftSize / 4,
    val sampleRate: Float = fftSize.toFloat(),
    val window: WindowType = WindowType.HANN,
) : FloatFilter {
    private val ring = DoubleArray(fftSize)
    private var writePos = 0
    private var filled = 0
    private var sinceHop = 0
    private val win = Windows.make(window, fftSize)

protected var held: Float = 0f

protected abstract fun onFrame(spectrum: Spectrum): Float

Compute the held output value from a freshly windowed one-sided magnitude spectrum.

    private fun snapshotSpectrum(): Spectrum {
        val re = DoubleArray(fftSize)
        val im = DoubleArray(fftSize)
        // Unroll the ring buffer oldest-first, applying the analysis window.
        for (i in 0 until fftSize) {
            val idx = (writePos + i) % fftSize
            re[i] = ring[idx] * win[i]
        }
        FFT.transform(re, im)
        return Spectrum(re, im, sampleRate.toDouble())
    }

override fun invoke(value: Float): Float

    override fun invoke(value: Float): Float {
        ring[writePos] = value.toDouble()
        writePos = (writePos + 1) % fftSize
        if (filled < fftSize) filled++
        sinceHop++
        if (filled >= fftSize && sinceHop >= hopSize) {
            sinceHop = 0
            held = onFrame(snapshotSpectrum())
        }
        return held
    }

fun reset()

    fun reset() {
        java.util.Arrays.fill(ring, 0.0)
        writePos = 0; filled = 0; sinceHop = 0; held = 0f
    }
}

class SpectralCentroidFilter

Streaming spectral centroid (Hz), held between hops.

class SpectralCentroidFilter(
    fftSize: Int = 1024,
    hopSize: Int = fftSize / 4,
    sampleRate: Float = fftSize.toFloat(),
    window: WindowType = WindowType.HANN,
) : SpectralFrameFilter(fftSize, hopSize, sampleRate, window) {

override fun onFrame(spectrum: Spectrum): Float = spectrum.centroid().toFloat()

class SpectralRolloffFilter

Streaming spectral rolloff (Hz), held between hops.

class SpectralRolloffFilter(
    fftSize: Int = 1024,
    hopSize: Int = fftSize / 4,
    sampleRate: Float = fftSize.toFloat(),
    window: WindowType = WindowType.HANN,
    var rolloffFraction: Double = 0.85,
) : SpectralFrameFilter(fftSize, hopSize, sampleRate, window) {

override fun onFrame(spectrum: Spectrum): Float = spectrum.rolloff(rolloffFraction).toFloat()

class SpectralFluxFilter

Streaming spectral flux (onset novelty), held between hops.

class SpectralFluxFilter(
    fftSize: Int = 1024,
    hopSize: Int = fftSize / 4,
    sampleRate: Float = fftSize.toFloat(),
    window: WindowType = WindowType.HANN,
) : SpectralFrameFilter(fftSize, hopSize, sampleRate, window) {
    private var prev: DoubleArray? = null

override fun onFrame(spectrum: Spectrum): Float

    override fun onFrame(spectrum: Spectrum): Float {
        val mag = spectrum.magnitudeOneSided()
        val p = prev
        val flux = if (p == null) 0.0 else SpectralFeatures.flux(p, mag)
        prev = mag
        return flux.toFloat()
    }
}

class RmsFilter(val windowSize: Int = 256) : FloatFilter

Streaming windowed RMS: a rolling root-mean-square over the last windowSize samples. Cheap O(1)-per-sample running sum-of-squares — the time-domain energy companion to the spectral filters.

class RmsFilter(val windowSize: Int = 256) : FloatFilter {
    private val ring = FloatArray(windowSize)
    private var pos = 0
    private var filled = 0
    private var sumSq = 0.0

override fun invoke(value: Float): Float

    override fun invoke(value: Float): Float {
        val old = ring[pos]
        sumSq += value.toDouble() * value - old.toDouble() * old
        ring[pos] = value
        pos = (pos + 1) % windowSize
        if (filled < windowSize) filled++
        if (sumSq < 0.0) sumSq = 0.0 // guard fp drift
        return sqrt(sumSq / filled).toFloat()
    }

fun reset()

    fun reset() {
        java.util.Arrays.fill(ring, 0f); pos = 0; filled = 0; sumSq = 0.0
    }
}