sva_samples/measure/
formants.rs

1// Concern: resolves each frame into the resonances an all-pole model puts under it | Non-concern: the discrete spectrum (spectrum.rs), naming a peak (pitch.rs) | IO: (&[f64], sample rate) -> frames
2
3use crate::measure::spectrum::db;
4
5pub const MAX_ORDER: usize = 64;
6
7/// Nearer than this to DC or Nyquist a pole is the overall tilt, not a resonance.
8const EDGE_HZ: f64 = 40.0;
9
10#[derive(Clone, Copy, Debug, PartialEq)]
11pub struct Formant {
12    pub hz: f64,
13    pub bandwidth_hz: f64,
14    pub db: f64,
15}
16
17/// `residual` and `energy` are mean squares, and neither derives from the other.
18#[derive(Clone, Debug, PartialEq)]
19pub struct FormantFrame {
20    pub t_secs: f64,
21    pub order: usize,
22    pub energy: f64,
23    pub residual: f64,
24    pub formants: Vec<Formant>,
25}
26
27/// Two poles per kilohertz of bandwidth plus two for the tilt — 46 at 44.1 kHz, far more than a
28/// voice wants and about right for a full-band instrument.
29pub fn default_order(sample_rate: f64) -> usize {
30    (2.0 + sample_rate / 1000.0).clamp(2.0, MAX_ORDER as f64) as usize
31}
32
33pub fn track(
34    samples: &[f64],
35    sample_rate: f64,
36    start_secs: f64,
37    frame_secs: f64,
38    order: usize,
39    max_formants: usize,
40) -> Vec<FormantFrame> {
41    let stride = ((frame_secs * sample_rate).round() as usize).max(1);
42    samples
43        .chunks(stride)
44        .enumerate()
45        .map(|(n, chunk)| FormantFrame {
46            t_secs: start_secs + (n * stride) as f64 / sample_rate,
47            ..analyze(chunk, sample_rate, order, max_formants)
48        })
49        .collect()
50}
51
52/// Hamming, not Hann: a taper reaching zero throws the frame's edges away. Never
53/// pre-emphasised; `x - 0.97*x(t - 1sp)` is a composition's own to write.
54pub fn analyze(frame: &[f64], sample_rate: f64, order: usize, max_formants: usize) -> FormantFrame {
55    let order = order.clamp(2, MAX_ORDER).min(frame.len().saturating_sub(1));
56    let windowed: Vec<f64> = frame
57        .iter()
58        .enumerate()
59        .map(|(n, &x)| {
60            let turn = 2.0 * std::f64::consts::PI * n as f64 / frame.len() as f64;
61            let w = 0.54 - 0.46 * turn.cos();
62            x * w
63        })
64        .collect();
65    let energy = match frame.is_empty() {
66        true => 0.0,
67        false => windowed.iter().map(|x| x * x).sum::<f64>() / frame.len() as f64,
68    };
69    if order < 2 || energy <= 0.0 {
70        return FormantFrame {
71            t_secs: 0.0,
72            order,
73            energy,
74            residual: energy,
75            formants: Vec::new(),
76        };
77    }
78
79    let r = autocorrelation(&windowed, order);
80    let (a, error) = levinson(&r, order);
81    let residual = error / frame.len() as f64;
82    let gain = residual.max(0.0).sqrt();
83
84    let mut formants: Vec<Formant> = roots(&a)
85        .into_iter()
86        .filter(|z| z.im > 0.0)
87        .filter_map(|z| {
88            let hz = z.im.atan2(z.re) * sample_rate / (2.0 * std::f64::consts::PI);
89            let radius = z.abs();
90            if !(EDGE_HZ..sample_rate / 2.0 - EDGE_HZ).contains(&hz) || radius >= 1.0 {
91                return None;
92            }
93            Some(Formant {
94                hz,
95                bandwidth_hz: -radius.ln() * sample_rate / std::f64::consts::PI,
96                db: db(gain / response(&a, hz, sample_rate)),
97            })
98        })
99        .collect();
100    formants.sort_by(|x, y| x.hz.total_cmp(&y.hz));
101    formants.truncate(max_formants);
102
103    FormantFrame {
104        t_secs: 0.0,
105        order,
106        energy,
107        residual,
108        formants,
109    }
110}
111
112fn autocorrelation(x: &[f64], lags: usize) -> Vec<f64> {
113    (0..=lags)
114        .map(|k| x.iter().skip(k).zip(x).map(|(a, b)| a * b).sum())
115        .collect()
116}
117
118/// Levinson-Durbin: the order-`p` predictor from `p` autocorrelations, no FFT and no matrix.
119fn levinson(r: &[f64], order: usize) -> (Vec<f64>, f64) {
120    let mut a = vec![0.0; order + 1];
121    a[0] = 1.0;
122    let mut error = r[0];
123
124    for i in 1..=order {
125        if error <= 0.0 {
126            return (a, error.max(0.0));
127        }
128        let acc: f64 = r[i] + (1..i).map(|j| a[j] * r[i - j]).sum::<f64>();
129        let k = -acc / error;
130        let held: Vec<f64> = a[1..i].to_vec();
131        for (j, prior) in held.iter().enumerate() {
132            a[j + 1] = prior + k * held[i - 2 - j];
133        }
134        a[i] = k;
135        error *= 1.0 - k * k;
136    }
137    (a, error.max(0.0))
138}
139
140fn response(a: &[f64], hz: f64, sample_rate: f64) -> f64 {
141    let w = 2.0 * std::f64::consts::PI * hz / sample_rate;
142    let (mut re, mut im) = (0.0, 0.0);
143    for (k, coefficient) in a.iter().enumerate() {
144        re += coefficient * (w * k as f64).cos();
145        im -= coefficient * (w * k as f64).sin();
146    }
147    (re * re + im * im).sqrt()
148}
149
150#[derive(Clone, Copy, Debug)]
151struct C {
152    re: f64,
153    im: f64,
154}
155
156impl C {
157    fn mul(self, o: C) -> C {
158        C {
159            re: self.re * o.re - self.im * o.im,
160            im: self.re * o.im + self.im * o.re,
161        }
162    }
163
164    fn sub(self, o: C) -> C {
165        C {
166            re: self.re - o.re,
167            im: self.im - o.im,
168        }
169    }
170
171    fn div(self, o: C) -> C {
172        let d = o.re * o.re + o.im * o.im;
173        match d == 0.0 {
174            true => C { re: 0.0, im: 0.0 },
175            false => C {
176                re: (self.re * o.re + self.im * o.im) / d,
177                im: (self.im * o.re - self.re * o.im) / d,
178            },
179        }
180    }
181
182    fn abs(self) -> f64 {
183        self.re.hypot(self.im)
184    }
185}
186
187/// Durand-Kerner: every root moves at once, so no deflation loses the closely spaced pairs an
188/// envelope is made of.
189fn roots(a: &[f64]) -> Vec<C> {
190    const ROUNDS: usize = 500;
191    const SETTLED: f64 = 1e-14;
192
193    let degree = a.len() - 1;
194    let seed = C { re: 0.4, im: 0.9 };
195    let mut z: Vec<C> = Vec::with_capacity(degree);
196    let mut power = C { re: 1.0, im: 0.0 };
197    for _ in 0..degree {
198        z.push(power);
199        power = power.mul(seed);
200    }
201
202    for _ in 0..ROUNDS {
203        let mut moved: f64 = 0.0;
204        for i in 0..degree {
205            let mut denominator = C { re: 1.0, im: 0.0 };
206            for j in 0..degree {
207                if j != i {
208                    denominator = denominator.mul(z[i].sub(z[j]));
209                }
210            }
211            let step = evaluate(a, z[i]).div(denominator);
212            z[i] = z[i].sub(step);
213            moved = moved.max(step.abs());
214        }
215        if moved < SETTLED {
216            break;
217        }
218    }
219    z
220}
221
222fn evaluate(a: &[f64], at: C) -> C {
223    let mut acc = C { re: 0.0, im: 0.0 };
224    for coefficient in a {
225        acc = acc.mul(at);
226        acc.re += coefficient;
227    }
228    acc
229}
sva_samples/measure/formants.rs

sva_samples/measure/
formants.rs
