229 lines
7.9 KiB
Dart
229 lines
7.9 KiB
Dart
import 'dart:math' as math;
|
||
import 'dart:typed_data';
|
||
|
||
/// Stateless DSP utility functions used by the modem and voice codec.
|
||
abstract final class AudioMath {
|
||
// ── Window functions ──────────────────────────────────────────────
|
||
|
||
/// Hamming window of length [n].
|
||
static Float64List hammingWindow(int n) {
|
||
final w = Float64List(n);
|
||
final factor = 2.0 * math.pi / (n - 1);
|
||
for (int i = 0; i < n; i++) {
|
||
w[i] = 0.54 - 0.46 * math.cos(i * factor);
|
||
}
|
||
return w;
|
||
}
|
||
|
||
/// Hann (Hanning) window of length [n].
|
||
static Float64List hannWindow(int n) {
|
||
final w = Float64List(n);
|
||
final factor = 2.0 * math.pi / (n - 1);
|
||
for (int i = 0; i < n; i++) {
|
||
w[i] = 0.5 * (1.0 - math.cos(i * factor));
|
||
}
|
||
|
||
return w;
|
||
}
|
||
|
||
// ── Signal statistics ─────────────────────────────────────────────
|
||
|
||
/// Root-mean-square of [samples] (all values).
|
||
static double rms(Float64List samples) {
|
||
double sum = 0.0;
|
||
for (final s in samples) { sum += s * s; }
|
||
return math.sqrt(sum / samples.length);
|
||
}
|
||
|
||
/// Short-time energy (sum of squares over [samples]).
|
||
static double energy(Float64List samples) {
|
||
double e = 0.0;
|
||
for (final s in samples) { e += s * s; }
|
||
return e;
|
||
}
|
||
|
||
/// Zero-crossing rate of [samples] (crossings per sample).
|
||
static double zeroCrossingRate(Float64List samples) {
|
||
int crossings = 0;
|
||
for (int i = 1; i < samples.length; i++) {
|
||
if ((samples[i] >= 0) != (samples[i - 1] >= 0)) crossings++;
|
||
}
|
||
return crossings / samples.length;
|
||
}
|
||
|
||
// ── Int16 ↔ Float64 conversion ────────────────────────────────────
|
||
|
||
/// Convert raw 16-bit PCM [bytes] (little-endian Int16) to normalised
|
||
/// Float64 in [-1, +1].
|
||
static Float64List int16BytesToFloat(Uint8List bytes) {
|
||
final count = bytes.length >> 1;
|
||
final out = Float64List(count);
|
||
final view = ByteData.sublistView(bytes);
|
||
for (int i = 0; i < count; i++) {
|
||
out[i] = view.getInt16(i * 2, Endian.little) / 32768.0;
|
||
}
|
||
return out;
|
||
}
|
||
|
||
/// Convert normalised Float64 [-1, +1] to Int16 LE bytes.
|
||
static Uint8List floatToInt16Bytes(Float64List samples) {
|
||
final out = Uint8List(samples.length * 2);
|
||
final view = ByteData.sublistView(out);
|
||
for (int i = 0; i < samples.length; i++) {
|
||
final clamped = samples[i].clamp(-1.0, 1.0);
|
||
view.setInt16(i * 2, (clamped * 32767).round(), Endian.little);
|
||
}
|
||
return out;
|
||
}
|
||
|
||
// ── Filter ────────────────────────────────────────────────────────
|
||
|
||
/// Apply a simple first-order high-pass (pre-emphasis) filter in-place.
|
||
/// y[n] = x[n] − alpha·x[n−1]
|
||
static void preEmphasis(Float64List samples, double alpha,
|
||
[double prevSample = 0.0]) {
|
||
double prev = prevSample;
|
||
for (int i = 0; i < samples.length; i++) {
|
||
final cur = samples[i];
|
||
samples[i] = cur - alpha * prev;
|
||
prev = cur;
|
||
}
|
||
}
|
||
|
||
/// Apply de-emphasis (inverse of pre-emphasis) in-place.
|
||
/// y[n] = x[n] + alpha·y[n−1]
|
||
static void deEmphasis(Float64List samples, double alpha,
|
||
[double prevOut = 0.0]) {
|
||
double prev = prevOut;
|
||
for (int i = 0; i < samples.length; i++) {
|
||
final cur = samples[i] + alpha * prev;
|
||
samples[i] = cur;
|
||
prev = cur;
|
||
}
|
||
}
|
||
|
||
// ── Autocorrelation ───────────────────────────────────────────────
|
||
|
||
/// Compute biased autocorrelation of [frame] at lags 0..maxLag (inclusive).
|
||
/// Result is a Float64List of length maxLag+1.
|
||
static Float64List autocorrelation(Float64List frame, int maxLag) {
|
||
final n = frame.length;
|
||
final r = Float64List(maxLag + 1);
|
||
for (int lag = 0; lag <= maxLag; lag++) {
|
||
double sum = 0.0;
|
||
for (int i = 0; i < n - lag; i++) {
|
||
sum += frame[i] * frame[i + lag];
|
||
}
|
||
r[lag] = sum;
|
||
}
|
||
return r;
|
||
}
|
||
|
||
// ── AMDF pitch detector ───────────────────────────────────────────
|
||
|
||
/// Average Magnitude Difference Function.
|
||
/// Returns the pitch period in samples (within [minP, maxP]) or 0 if
|
||
/// the frame appears unvoiced (no clear minimum in AMDF).
|
||
static int amdfPitch(Float64List frame, int minP, int maxP) {
|
||
final n = frame.length;
|
||
double minAmdf = double.infinity;
|
||
int bestTau = 0;
|
||
|
||
for (int tau = minP; tau <= maxP; tau++) {
|
||
double sum = 0.0;
|
||
final len = n - tau;
|
||
for (int i = 0; i < len; i++) {
|
||
sum += (frame[i] - frame[i + tau]).abs();
|
||
}
|
||
final amdf = sum / len;
|
||
if (amdf < minAmdf) {
|
||
minAmdf = amdf;
|
||
bestTau = tau;
|
||
}
|
||
}
|
||
|
||
// If minimum AMDF is more than 60 % of mean magnitude, frame is unvoiced.
|
||
double meanMag = 0.0;
|
||
for (final s in frame) { meanMag += s.abs(); }
|
||
meanMag /= n;
|
||
|
||
if (meanMag < 1e-6 || minAmdf > 0.6 * meanMag) return 0;
|
||
return bestTau;
|
||
}
|
||
|
||
// ── Goertzel single-tone DFT ──────────────────────────────────────
|
||
|
||
/// Compute the squared magnitude of DFT bin for frequency [freqHz] over
|
||
/// [samples]. More efficient than full FFT for a single frequency.
|
||
///
|
||
/// Phase: Q0, Q1, Q2 Goertzel state variables.
|
||
/// coeff = 2·cos(2π·f/fs)
|
||
static double goertzelMagnitudeSq(
|
||
Float64List samples, double freqHz, int sampleRate) {
|
||
final coeff = 2.0 * math.cos(2.0 * math.pi * freqHz / sampleRate);
|
||
double q0 = 0.0, q1 = 0.0, q2 = 0.0;
|
||
for (final s in samples) {
|
||
q0 = coeff * q1 - q2 + s;
|
||
q2 = q1;
|
||
q1 = q0;
|
||
}
|
||
return q1 * q1 + q2 * q2 - q1 * q2 * coeff;
|
||
}
|
||
|
||
// ── CRC-16 (CCITT / X.25 polynomial 0x1021) ──────────────────────
|
||
|
||
static const int _crcPoly = 0x1021;
|
||
|
||
/// Compute CRC-16/CCITT over [data]. Initial value 0xFFFF.
|
||
static int crc16(List<int> data) {
|
||
int crc = 0xFFFF;
|
||
for (final byte in data) {
|
||
crc ^= (byte & 0xFF) << 8;
|
||
for (int i = 0; i < 8; i++) {
|
||
if ((crc & 0x8000) != 0) {
|
||
crc = ((crc << 1) ^ _crcPoly) & 0xFFFF;
|
||
} else {
|
||
crc = (crc << 1) & 0xFFFF;
|
||
}
|
||
}
|
||
}
|
||
return crc;
|
||
}
|
||
|
||
// ── Clamp / normalise ─────────────────────────────────────────────
|
||
|
||
/// Soft-clip [x] to [-1, 1] using tanh to avoid hard saturation.
|
||
static double softClip(double x) => math.max(-1.0, math.min(1.0, x));
|
||
|
||
/// Normalise [samples] to peak amplitude 1.0 (in-place).
|
||
static void normalise(Float64List samples) {
|
||
double peak = 0.0;
|
||
for (final s in samples) {
|
||
final a = s.abs();
|
||
if (a > peak) peak = a;
|
||
}
|
||
if (peak < 1e-10) return;
|
||
final scale = 1.0 / peak;
|
||
for (int i = 0; i < samples.length; i++) { samples[i] *= scale; }
|
||
}
|
||
|
||
// ── Pseudo-random noise (for unvoiced excitation) ─────────────────
|
||
|
||
static final _rng = math.Random();
|
||
|
||
/// Return a Float64List of [n] samples of white Gaussian noise
|
||
/// approximated by Box-Muller transform, amplitude ≈ 0.3 RMS.
|
||
static Float64List whiteNoise(int n, [double amplitude = 0.3]) {
|
||
final out = Float64List(n);
|
||
for (int i = 0; i < n; i += 2) {
|
||
final u1 = _rng.nextDouble();
|
||
final u2 = _rng.nextDouble();
|
||
final r = math.sqrt(-2.0 * math.log(u1.clamp(1e-12, 1.0)));
|
||
final theta = 2.0 * math.pi * u2;
|
||
out[i] = amplitude * r * math.cos(theta);
|
||
if (i + 1 < n) out[i + 1] = amplitude * r * math.sin(theta);
|
||
}
|
||
return out;
|
||
}
|
||
}
|