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: (planes, 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
33/// A fit of `order` poles reads `order + 1` lags.
34#[derive(Clone, Copy, Debug, PartialEq)]
35pub enum UnderOrder {
36    Frame {
37        frame_secs: f64,
38        held: usize,
39        order: usize,
40        least_secs: f64,
41    },
42    Target {
43        held: usize,
44        order: usize,
45        least_secs: f64,
46    },
47}
48
49pub fn track(
50    planes: &[&[f64]],
51    sample_rate: f64,
52    start_secs: f64,
53    frame_secs: f64,
54    order: usize,
55    max_formants: usize,
56) -> Result<Vec<FormantFrame>, UnderOrder> {
57    let stride = ((frame_secs * sample_rate).round() as usize).max(1);
58    let len = crate::measure::held(planes);
59    let order = order.clamp(2, MAX_ORDER);
60    let least_secs = ((order + 1) as f64 / sample_rate * 1e4).ceil() / 1e4;
61    if stride <= order {
62        return Err(UnderOrder::Frame {
63            frame_secs,
64            held: stride,
65            order,
66            least_secs,
67        });
68    }
69    if len <= order {
70        return Err(UnderOrder::Target {
71            held: len,
72            order,
73            least_secs,
74        });
75    }
76    Ok((0..len)
77        .step_by(stride)
78        .enumerate()
79        .map(|(n, from)| {
80            let chunk: Vec<&[f64]> = planes
81                .iter()
82                .map(|p| &p[from..len.min(from + stride)])
83                .collect();
84            FormantFrame {
85                t_secs: start_secs + (n * stride) as f64 / sample_rate,
86                ..analyze(&chunk, sample_rate, order, max_formants)
87            }
88        })
89        .collect())
90}
91
92/// Hamming, not Hann: a taper reaching zero throws the frame's edges away. Never
93/// pre-emphasised; `x - 0.97*x(t - 1sp)` is a composition's own to write. Zero past a short
94/// frame's end.
95pub fn analyze(
96    planes: &[&[f64]],
97    sample_rate: f64,
98    order: usize,
99    max_formants: usize,
100) -> FormantFrame {
101    let len = crate::measure::held(planes);
102    let order = order.clamp(2, MAX_ORDER);
103    let windowed: Vec<Vec<f64>> = planes
104        .iter()
105        .map(|frame| {
106            frame
107                .iter()
108                .enumerate()
109                .map(|(n, &x)| {
110                    let turn = 2.0 * std::f64::consts::PI * n as f64 / len as f64;
111                    let w = 0.54 - 0.46 * turn.cos();
112                    x * w
113                })
114                .collect()
115        })
116        .collect();
117    let energy = match len {
118        0 => 0.0,
119        _ => {
120            let each = windowed
121                .iter()
122                .map(|w| w.iter().map(|x| x * x).sum::<f64>());
123            each.sum::<f64>() / (len * planes.len()) as f64
124        }
125    };
126    if energy <= 0.0 {
127        return FormantFrame {
128            t_secs: 0.0,
129            order,
130            energy,
131            residual: energy,
132            formants: Vec::new(),
133        };
134    }
135
136    let r = crate::measure::averaged(windowed.iter().map(|w| autocorrelation(w, order)).collect());
137    let (a, error) = levinson(&r, order);
138    let residual = error / len as f64;
139    let gain = residual.max(0.0).sqrt();
140
141    let mut formants: Vec<Formant> = roots(&a)
142        .into_iter()
143        .filter(|z| z.im > 0.0)
144        .filter_map(|z| {
145            let hz = z.im.atan2(z.re) * sample_rate / (2.0 * std::f64::consts::PI);
146            let radius = z.abs();
147            if !(EDGE_HZ..sample_rate / 2.0 - EDGE_HZ).contains(&hz) || radius >= 1.0 {
148                return None;
149            }
150            Some(Formant {
151                hz,
152                bandwidth_hz: -radius.ln() * sample_rate / std::f64::consts::PI,
153                db: db(gain / response(&a, hz, sample_rate)),
154            })
155        })
156        .collect();
157    formants.sort_by(|x, y| x.hz.total_cmp(&y.hz));
158    formants.truncate(max_formants);
159
160    FormantFrame {
161        t_secs: 0.0,
162        order,
163        energy,
164        residual,
165        formants,
166    }
167}
168
169fn autocorrelation(x: &[f64], lags: usize) -> Vec<f64> {
170    (0..=lags)
171        .map(|k| x.iter().skip(k).zip(x).map(|(a, b)| a * b).sum())
172        .collect()
173}
174
175/// Levinson-Durbin: the order-`p` predictor from `p` autocorrelations, no FFT and no matrix.
176fn levinson(r: &[f64], order: usize) -> (Vec<f64>, f64) {
177    let mut a = vec![0.0; order + 1];
178    a[0] = 1.0;
179    let mut error = r[0];
180
181    for i in 1..=order {
182        if error <= 0.0 {
183            return (a, error.max(0.0));
184        }
185        let acc: f64 = r[i] + (1..i).map(|j| a[j] * r[i - j]).sum::<f64>();
186        let k = -acc / error;
187        let held: Vec<f64> = a[1..i].to_vec();
188        for (j, prior) in held.iter().enumerate() {
189            a[j + 1] = prior + k * held[i - 2 - j];
190        }
191        a[i] = k;
192        error *= 1.0 - k * k;
193    }
194    (a, error.max(0.0))
195}
196
197fn response(a: &[f64], hz: f64, sample_rate: f64) -> f64 {
198    let w = 2.0 * std::f64::consts::PI * hz / sample_rate;
199    let (mut re, mut im) = (0.0, 0.0);
200    for (k, coefficient) in a.iter().enumerate() {
201        re += coefficient * (w * k as f64).cos();
202        im -= coefficient * (w * k as f64).sin();
203    }
204    (re * re + im * im).sqrt()
205}
206
207#[derive(Clone, Copy, Debug)]
208struct C {
209    re: f64,
210    im: f64,
211}
212
213impl C {
214    fn mul(self, o: C) -> C {
215        C {
216            re: self.re * o.re - self.im * o.im,
217            im: self.re * o.im + self.im * o.re,
218        }
219    }
220
221    fn sub(self, o: C) -> C {
222        C {
223            re: self.re - o.re,
224            im: self.im - o.im,
225        }
226    }
227
228    fn div(self, o: C) -> C {
229        let d = o.re * o.re + o.im * o.im;
230        match d == 0.0 {
231            true => C { re: 0.0, im: 0.0 },
232            false => C {
233                re: (self.re * o.re + self.im * o.im) / d,
234                im: (self.im * o.re - self.re * o.im) / d,
235            },
236        }
237    }
238
239    fn abs(self) -> f64 {
240        self.re.hypot(self.im)
241    }
242}
243
244/// Durand-Kerner: every root moves at once, so no deflation loses the closely spaced pairs an
245/// envelope is made of.
246fn roots(a: &[f64]) -> Vec<C> {
247    const ROUNDS: usize = 500;
248    const SETTLED: f64 = 1e-14;
249
250    let degree = a.len() - 1;
251    let seed = C { re: 0.4, im: 0.9 };
252    let mut z: Vec<C> = Vec::with_capacity(degree);
253    let mut power = C { re: 1.0, im: 0.0 };
254    for _ in 0..degree {
255        z.push(power);
256        power = power.mul(seed);
257    }
258
259    for _ in 0..ROUNDS {
260        let mut moved: f64 = 0.0;
261        for i in 0..degree {
262            let mut denominator = C { re: 1.0, im: 0.0 };
263            for j in 0..degree {
264                if j != i {
265                    denominator = denominator.mul(z[i].sub(z[j]));
266                }
267            }
268            let step = evaluate(a, z[i]).div(denominator);
269            z[i] = z[i].sub(step);
270            moved = moved.max(step.abs());
271        }
272        if moved < SETTLED {
273            break;
274        }
275    }
276    z
277}
278
279fn evaluate(a: &[f64], at: C) -> C {
280    let mut acc = C { re: 0.0, im: 0.0 };
281    for coefficient in a {
282        acc = acc.mul(at);
283        acc.re += coefficient;
284    }
285    acc
286}
sva_samples/measure/formants.rs

sva_samples/measure/
formants.rs
