sound_to_patch/estimate/pitch/
line.rs

1//! The sub-band Hilbert line: the instantaneous frequency of the
2//! fundamental of a low note, one value per sample.
3//!
4//! An f0 tracker averages over its frame, 93 ms at 22,050 Hz, so it smears
5//! a pitch drop that settles in 20 ms. The line keeps the fundamental
6//! alone with a band-pass, takes the frequency of its analytic signal
7//! sample by sample and removes the phase spikes with a short median.
8
9use crate::audio::SampleRate;
10use crate::cast;
11use crate::dsp::{AnalyticSignal, FilterError, band_pass_zero_phase};
12use crate::units::{Hz, Seconds};
13
14/// The low edge of the band.
15pub(crate) const LINE_LO: Hz = Hz::from_const(20.0);
16/// The low edge of the band for the refit of a fast fall.
17///
18/// Where a tone below 40 Hz starts from silence or rises, the 20 Hz high-pass raises
19/// the line for about 10 ms, which fits as a fast fall; an octave below
20/// the band it does not.
21pub(crate) const LINE_LO_FAST: Hz = Hz::from_const(10.0);
22/// The high edge of the band.
23pub(crate) const LINE_HI: Hz = Hz::from_const(200.0);
24/// The length of the median filter, about 1.4 ms: 31 samples at 22,050 Hz,
25/// 63 at 44,100 Hz and 67 at 48,000 Hz.
26pub(crate) const LINE_MEDIAN: Seconds = Seconds::from_const(0.0014);
27
28/// The convention of the line, as `layers.json` names it.
29pub const LINE_CONVENTION: &str = "fourth-order Butterworth high-pass at 20 Hz and low-pass at \
30200 Hz, each forward and backward; the frequency of the analytic signal between consecutive \
31samples, arg(z[n+1] conj(z[n])) fs / 2 pi; a centred median of round(1.4 ms fs) samples, made \
32odd, shortened at the ends; the refit of a drop uses a high-pass at 10 Hz";
33
34/// The instantaneous frequency of the 20 to 200 Hz band of a segment, with
35/// the envelope of that band.
36#[derive(Clone, Debug, PartialEq)]
37pub struct SubBandLine {
38    start: Seconds,
39    rate: SampleRate,
40    hz: Vec<f64>,
41    envelope: Vec<f64>,
42}
43
44impl SubBandLine {
45    /// Builds the line of `segment`, whose first sample lies `start` into
46    /// the stem.
47    ///
48    /// The analytic signal is meant for note-sized segments. Each cut end
49    /// carries an error that decays with the distance from it, so a caller
50    /// keeps a margin of silence or signal around the span it reads.
51    ///
52    /// # Errors
53    ///
54    /// Returns a [`FilterError`] when the rate puts 200 Hz at or above the
55    /// Nyquist frequency.
56    pub fn of(segment: &[f64], rate: SampleRate, start: Seconds) -> Result<Self, FilterError> {
57        Self::with_low_edge(segment, rate, start, LINE_LO)
58    }
59
60    /// Builds the line as [`Self::of`] does, with the high-pass at `lo`.
61    ///
62    /// # Errors
63    ///
64    /// Returns a [`FilterError`] when `lo` is not above 0 Hz or the rate
65    /// puts 200 Hz at or above the Nyquist frequency.
66    pub fn with_low_edge(
67        segment: &[f64],
68        rate: SampleRate,
69        start: Seconds,
70        lo: Hz,
71    ) -> Result<Self, FilterError> {
72        Self::with_edges(segment, rate, start, lo, LINE_HI)
73    }
74
75    /// Builds the line as [`Self::of`] does, with the high-pass at `lo` and
76    /// the low-pass at `hi`.
77    ///
78    /// # Errors
79    ///
80    /// Returns a [`FilterError`] when an edge is not above 0 Hz or the rate
81    /// puts `hi` at or above the Nyquist frequency.
82    pub fn with_edges(
83        segment: &[f64],
84        rate: SampleRate,
85        start: Seconds,
86        lo: Hz,
87        hi: Hz,
88    ) -> Result<Self, FilterError> {
89        let band = band_pass_zero_phase(segment, rate, lo, hi)?;
90        // Padding as long as the segment keeps the circular wrap of the FFT
91        // off the segment.
92        let analytic = AnalyticSignal::new(&band, rate, band.len());
93        let hz = median_filter(&analytic.instantaneous_frequency(), median_len(rate));
94        Ok(Self {
95            start,
96            rate,
97            hz,
98            envelope: analytic.envelope(),
99        })
100    }
101
102    /// Returns the frequency between samples `n` and `n + 1`, in hertz. A
103    /// value can be negative where the band holds no signal.
104    #[must_use]
105    pub fn hz(&self) -> &[f64] {
106        &self.hz
107    }
108
109    /// Returns the rate of the segment.
110    #[must_use]
111    pub fn sample_rate(&self) -> SampleRate {
112        self.rate
113    }
114
115    /// Returns the envelope of the band, one value per sample.
116    #[must_use]
117    pub fn envelope(&self) -> &[f64] {
118        &self.envelope
119    }
120
121    /// Returns the time of frequency value `n`: halfway between samples `n`
122    /// and `n + 1`.
123    #[must_use]
124    pub fn time_of(&self, n: usize) -> Seconds {
125        Seconds::clamped(self.start.get() + (cast::usize_to_f64(n) + 0.5) / self.rate.hz().get())
126    }
127
128    /// Returns the time of sample `n`.
129    #[must_use]
130    pub fn sample_time(&self, n: usize) -> Seconds {
131        Seconds::clamped(self.start.get() + cast::usize_to_f64(n) / self.rate.hz().get())
132    }
133
134    /// Returns the index of the first sample at or after `time`.
135    #[must_use]
136    pub fn index_at(&self, time: Seconds) -> usize {
137        let offset = (time.get() - self.start.get()) * self.rate.hz().get();
138        cast::ceil_to_usize(offset.max(0.0)).unwrap_or(0)
139    }
140}
141
142/// The median length at `rate`: `round(1.4 ms fs)`, plus one when even.
143#[must_use]
144pub(crate) fn median_len(rate: SampleRate) -> usize {
145    let len = cast::round_to_usize(LINE_MEDIAN.get() * rate.hz().get()).unwrap_or(1);
146    len.max(1) | 1
147}
148
149/// A centred running median of odd length `len`. Near an end the window
150/// shrinks to the samples that exist, so the ends are not pulled towards a
151/// padding value.
152#[must_use]
153pub(crate) fn median_filter(x: &[f64], len: usize) -> Vec<f64> {
154    let half = len / 2;
155    let mut window = Vec::with_capacity(len);
156    (0..x.len())
157        .map(|n| {
158            window.clear();
159            window.extend_from_slice(&x[n.saturating_sub(half)..(n + half + 1).min(x.len())]);
160            window.sort_by(f64::total_cmp);
161            let mid = window.len() / 2;
162            if window.len() % 2 == 1 {
163                window[mid]
164            } else {
165                f64::midpoint(window[mid - 1], window[mid])
166            }
167        })
168        .collect()
169}
170
171#[cfg(test)]
172mod tests {
173    use super::*;
174
175    #[test]
176    fn the_convention_states_the_constants() {
177        use super::super::{num, says};
178        says(
179            LINE_CONVENTION,
180            &[
181                format!(
182                    "high-pass at {} Hz and low-pass at {} Hz",
183                    num(LINE_LO.get()),
184                    num(LINE_HI.get())
185                ),
186                format!("round({} ms fs)", num(1000.0 * LINE_MEDIAN.get())),
187                format!("uses a high-pass at {} Hz", num(LINE_LO_FAST.get())),
188            ],
189        );
190    }
191
192    #[test]
193    fn the_median_length_is_odd_and_near_1_4_ms() {
194        for (rate, len) in [(22_050, 31), (44_100, 63), (48_000, 67), (8_000, 11)] {
195            assert_eq!(median_len(SampleRate::new(rate).unwrap()), len, "{rate} Hz");
196        }
197    }
198
199    #[test]
200    fn the_median_removes_a_spike_and_keeps_the_ends() {
201        let x = [1.0, 1.0, 9.0, 1.0, 1.0, 2.0, 3.0];
202        assert_eq!(median_filter(&x, 3), [1.0, 1.0, 1.0, 1.0, 1.0, 2.0, 2.5]);
203        assert_eq!(median_filter(&[], 31), Vec::<f64>::new());
204        assert_eq!(median_filter(&[4.0], 31), [4.0]);
205    }
206
207    #[test]
208    fn a_steady_low_sine_reads_its_frequency() {
209        let rate = SampleRate::new(22_050).unwrap();
210        let x: Vec<f64> = (0..22_050)
211            .map(|n| (std::f64::consts::TAU * 41.2 * cast::usize_to_f64(n) / 22_050.0).sin())
212            .collect();
213        let line = SubBandLine::of(&x, rate, Seconds::new(1.0).unwrap()).unwrap();
214        assert_eq!(line.hz().len(), x.len() - 1);
215        // Away from the cut ends the frequency ripples at twice the tone,
216        // by an amount that falls as one over the distance from an end.
217        let middle = &line.hz()[5_000..17_000];
218        let mean = middle.iter().sum::<f64>() / cast::usize_to_f64(middle.len());
219        assert!((mean - 41.2).abs() < 0.01, "{mean}");
220        assert!(middle.iter().all(|f| (f - 41.2).abs() < 0.5));
221        assert!((line.time_of(0).get() - (1.0 + 0.5 / 22_050.0)).abs() < 1e-12);
222        assert_eq!(line.index_at(Seconds::new(1.5).unwrap()), 11_025);
223        assert!((line.envelope()[11_025] - 1.0).abs() < 0.01);
224    }
225}
sound_to_patch/estimate/pitch/line.rs

sound_to_patch/estimate/pitch/
line.rs
