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

sva_samples/measure/
formants.rs
