sound_to_patch/dsp/
analytic.rs

1//! The analytic signal: envelope, phase and instantaneous frequency, plus
2//! the band-pass heterodyne that moves one partial to baseband.
3//!
4//! The analytic signal means something only for one narrowband component
5//! (the Bedrosian condition). Band-pass around one partial first; the phase
6//! of a full-band stem is meaningless.
7
8use std::error::Error;
9use std::f64::consts::{PI, TAU};
10use std::fmt;
11
12use biquad::{Biquad, Coefficients, DirectForm2Transposed, ToHertz, Type};
13use realfft::num_complex::Complex64;
14use rustfft::FftPlanner;
15
16use crate::audio::SampleRate;
17use crate::cast;
18use crate::units::{Hz, Normalized};
19
20/// The default low-pass for an attack envelope. A 5 Hz smoother smears any
21/// attack shorter than tens of milliseconds; 100 to 200 Hz keeps it.
22pub const ATTACK_LOW_PASS_HZ: Hz = Hz::from_const(150.0);
23
24/// Q values of the two sections of a fourth-order Butterworth low-pass:
25/// `1 / (2 cos(pi / 8))` and `1 / (2 cos(3 pi / 8))`.
26const BUTTERWORTH_4: [f64; 2] = [0.541_196_100_146_197, 1.306_562_964_876_376_7];
27
28/// The analytic signal `z = x + j H{x}` of a real signal.
29#[derive(Clone, Debug, PartialEq)]
30pub struct AnalyticSignal {
31    z: Vec<Complex64>,
32    rate: Hz,
33}
34
35impl AnalyticSignal {
36    /// Computes the analytic signal by FFT: the negative bins are zeroed,
37    /// the positive bins doubled, and the DC bin and, for an even length,
38    /// the Nyquist bin keep gain 1.
39    ///
40    /// The FFT wraps circularly, so `pad` zeros go before and after the
41    /// signal and are cut after the transform. The cut ends still carry an
42    /// error that decays as `1 / d` with the distance `d` from an end,
43    /// about `fs / (2 pi^2 f d)` for a sine of frequency `f`: 0.1 % needs
44    /// some 4000 samples at 123 Hz and 8 kHz. Mirroring the signal instead
45    /// of zeros is worse, because a reversed partial runs backwards in phase.
46    ///
47    /// It is meant for note-sized segments. The transform holds the padded
48    /// signal as complex values, 16 bytes a sample, and an FFT of a length
49    /// that is not a power of two needs several times that in scratch: a
50    /// whole 10-minute stem at 48 kHz peaks near 2 GB. Heterodyne a long
51    /// signal instead.
52    #[must_use]
53    pub fn new(samples: &[f64], rate: SampleRate, pad: usize) -> Self {
54        let mut buffer: Vec<Complex64> = std::iter::repeat_n(0.0, pad)
55            .chain(samples.iter().copied())
56            .chain(std::iter::repeat_n(0.0, pad))
57            .map(|re| Complex64::new(re, 0.0))
58            .collect();
59        let len = buffer.len();
60        if len > 0 {
61            let mut planner = FftPlanner::<f64>::new();
62            planner.plan_fft_forward(len).process(&mut buffer);
63            let half = len / 2;
64            for (bin, value) in buffer.iter_mut().enumerate() {
65                let gain = if bin == 0 || (len.is_multiple_of(2) && bin == half) {
66                    1.0
67                } else if bin <= (len - 1) / 2 {
68                    2.0
69                } else {
70                    0.0
71                };
72                *value *= gain / cast::usize_to_f64(len);
73            }
74            planner.plan_fft_inverse(len).process(&mut buffer);
75        }
76        buffer.truncate(pad + samples.len());
77        buffer.drain(..pad);
78        Self {
79            z: buffer,
80            rate: rate.hz(),
81        }
82    }
83
84    /// Returns the complex samples.
85    #[must_use]
86    pub fn samples(&self) -> &[Complex64] {
87        &self.z
88    }
89
90    /// Returns the envelope `|z|`.
91    #[must_use]
92    pub fn envelope(&self) -> Vec<f64> {
93        self.z.iter().map(|value| value.norm()).collect()
94    }
95
96    /// Returns the unwrapped phase `arg z` in radians.
97    #[must_use]
98    pub fn phase(&self) -> Vec<f64> {
99        unwrap_phase(&self.z)
100    }
101
102    /// Returns the instantaneous frequency in hertz between consecutive
103    /// samples, `f[n] = (fs / 2 pi) arg(z[n+1] conj(z[n]))`. Avoids
104    /// differentiating a noisy phase. One value fewer than samples. The
105    /// per-sample track stays bare `f64` hertz: a value can be negative,
106    /// which [`Hz`] rejects.
107    #[must_use]
108    pub fn instantaneous_frequency(&self) -> Vec<f64> {
109        instantaneous_frequency(&self.z, self.rate)
110    }
111}
112
113fn unwrap_phase(z: &[Complex64]) -> Vec<f64> {
114    let mut phase = Vec::with_capacity(z.len());
115    let mut previous: Option<(f64, f64)> = None;
116    for value in z {
117        let raw = value.arg();
118        let unwrapped = match previous {
119            None => raw,
120            Some((last_raw, last)) => {
121                let mut step = raw - last_raw;
122                step -= TAU * ((step + PI) / TAU).floor();
123                last + step
124            }
125        };
126        previous = Some((raw, unwrapped));
127        phase.push(unwrapped);
128    }
129    phase
130}
131
132fn instantaneous_frequency(z: &[Complex64], rate: Hz) -> Vec<f64> {
133    z.windows(2)
134        .map(|pair| (pair[1] * pair[0].conj()).arg() * rate.get() / TAU)
135        .collect()
136}
137
138/// The modulation depth of an envelope, `(max - min) / (max + min)`.
139/// Returns `None` for an empty or all-zero envelope.
140#[must_use]
141pub fn am_depth(envelope: &[f64]) -> Option<Normalized> {
142    let max = envelope.iter().copied().reduce(f64::max)?;
143    let min = envelope.iter().copied().reduce(f64::min)?;
144    if max + min <= 0.0 {
145        return None;
146    }
147    Normalized::new((max - min) / (max + min)).ok()
148}
149
150/// Settings of a band-pass heterodyne.
151#[derive(Clone, Copy, Debug, PartialEq)]
152pub struct HeterodyneParams {
153    centre: Hz,
154    cutoff: Hz,
155    decimation: usize,
156}
157
158impl HeterodyneParams {
159    /// Shifts `centre` to 0 Hz and keeps `cutoff` hertz on either side. The
160    /// decimation factor defaults to the largest that keeps four samples per
161    /// cutoff period: `floor(fs / (4 cutoff))`, at least 1.
162    #[must_use]
163    pub fn new(centre: Hz, cutoff: Hz) -> Self {
164        Self {
165            centre,
166            cutoff,
167            decimation: 0,
168        }
169    }
170
171    /// Sets the decimation factor; 0 restores the default.
172    #[must_use]
173    pub fn with_decimation(self, decimation: usize) -> Self {
174        Self { decimation, ..self }
175    }
176
177    fn decimation_for(self, rate: SampleRate) -> usize {
178        if self.decimation > 0 {
179            return self.decimation;
180        }
181        cast::floor_to_usize(rate.hz().get() / (4.0 * self.cutoff.get()))
182            .unwrap_or(1)
183            .max(1)
184    }
185}
186
187/// One partial moved to baseband: complex samples at a reduced rate.
188#[derive(Clone, Debug, PartialEq)]
189pub struct Baseband {
190    z: Vec<Complex64>,
191    rate: Hz,
192}
193
194impl Baseband {
195    /// Returns the complex baseband samples.
196    #[must_use]
197    pub fn samples(&self) -> &[Complex64] {
198        &self.z
199    }
200
201    /// Returns the baseband sample rate, `fs / decimation`.
202    #[must_use]
203    pub fn rate(&self) -> Hz {
204        self.rate
205    }
206
207    /// Returns the envelope `|z|`.
208    #[must_use]
209    pub fn envelope(&self) -> Vec<f64> {
210        self.z.iter().map(|value| value.norm()).collect()
211    }
212
213    /// Returns the instantaneous frequency offset from the centre in hertz;
214    /// negative below the centre, so the per-sample track is bare `f64`.
215    #[must_use]
216    pub fn instantaneous_frequency(&self) -> Vec<f64> {
217        instantaneous_frequency(&self.z, self.rate)
218    }
219}
220
221/// Moves the band around `centre` to baseband, so a tone at `f_c + d`
222/// comes out at `d`.
223///
224/// It multiplies the signal by `2 exp(-j 2 pi f_c t)`, low-passes the
225/// product with a fourth-order Butterworth at `cutoff`, and keeps every
226/// `decimation`-th sample. A partial `A cos(2 pi f t + phi)` comes out as
227/// `A exp(j (2 pi (f - f_c) t + phi))`, so the envelope reads `A`. The
228/// first samples carry the filter's start-up transient.
229///
230/// # Errors
231///
232/// Returns a [`FilterError`] when the cutoff is zero or not below the
233/// Nyquist frequency.
234pub fn heterodyne(
235    samples: &[f64],
236    rate: SampleRate,
237    params: HeterodyneParams,
238) -> Result<Baseband, FilterError> {
239    let mut re = LowPass::new(rate, params.cutoff)?;
240    let mut im = LowPass::new(rate, params.cutoff)?;
241    let decimation = params.decimation_for(rate);
242    let step = TAU * params.centre.get() / rate.hz().get();
243    let z = samples
244        .iter()
245        .enumerate()
246        .map(|(n, &x)| {
247            // Reduce the phase per sample so it stays exact on long stems.
248            let phase = (step * cast::usize_to_f64(n)) % TAU;
249            // The factor 2 restores the amplitude that the mixing halves.
250            let mixed = 2.0 * x;
251            Complex64::new(re.run(mixed * phase.cos()), im.run(-mixed * phase.sin()))
252        })
253        .step_by(decimation)
254        .collect();
255    Ok(Baseband {
256        z,
257        rate: Hz::clamped(rate.hz().get() / cast::usize_to_f64(decimation)),
258    })
259}
260
261/// Low-passes a signal forward and then backward with a fourth-order
262/// Butterworth, so the result has no delay: an attack stays where it is.
263/// For attack envelopes use [`ATTACK_LOW_PASS_HZ`].
264///
265/// The two passes square the magnitude response, so the gain at `cutoff` is
266/// -6 dB, not -3 dB, and the slope is 48 dB per octave. As `filtfilt` does,
267/// each end is extended by an odd reflection, `2 x[0] - x[k]`, of three cutoff
268/// periods, and each pass starts in the steady state of its first input. A
269/// constant or a ramp then passes without a start-up or an end transient,
270/// so a release at the end of a stem keeps its level.
271///
272/// # Errors
273///
274/// Returns a [`FilterError`] when the cutoff is zero or not below the
275/// Nyquist frequency.
276pub fn low_pass_zero_phase(
277    samples: &[f64],
278    rate: SampleRate,
279    cutoff: Hz,
280) -> Result<Vec<f64>, FilterError> {
281    let pad = cast::ceil_to_usize(3.0 * rate.hz().get() / cutoff.get()).unwrap_or(0);
282    zero_phase(samples, pad, Reflection::Odd, || LowPass::new(rate, cutoff))
283}
284
285/// How [`zero_phase`] extends the ends of a signal.
286#[derive(Clone, Copy)]
287enum Reflection {
288    /// `2 x[0] - x[k]`: keeps the value and the slope, for a smooth signal.
289    Odd,
290    /// `x[k]`: keeps the value only. For a signal with a fast component that
291    /// the filter removes, an odd reflection would turn that component's end
292    /// value into a step the filter passes.
293    Even,
294}
295
296/// Runs a filter forward and backward over `samples` extended at each end
297/// by a reflection of up to `pad` samples.
298fn zero_phase(
299    samples: &[f64],
300    pad: usize,
301    reflection: Reflection,
302    filter: impl Fn() -> Result<LowPass, FilterError>,
303) -> Result<Vec<f64>, FilterError> {
304    let (mut forward, mut backward) = (filter()?, filter()?);
305    let (Some(&first), Some(&last)) = (samples.first(), samples.last()) else {
306        return Ok(Vec::new());
307    };
308    let pad = pad.min(samples.len() - 1);
309    let mirror = |end: f64, x: f64| match reflection {
310        Reflection::Odd => 2.0 * end - x,
311        Reflection::Even => x,
312    };
313    let mut out: Vec<f64> = samples[1..=pad]
314        .iter()
315        .rev()
316        .map(|&x| mirror(first, x))
317        .chain(samples.iter().copied())
318        .chain(
319            samples[samples.len() - 1 - pad..samples.len() - 1]
320                .iter()
321                .rev()
322                .map(|&x| mirror(last, x)),
323        )
324        .collect();
325    for pass in [&mut forward, &mut backward] {
326        if let Some(&start) = out.first() {
327            pass.settle(start);
328        }
329        for value in &mut out {
330            *value = pass.run(*value);
331        }
332        out.reverse();
333    }
334    out.truncate(out.len() - pad);
335    out.drain(..pad);
336    Ok(out)
337}
338
339/// Band-passes a signal with a fourth-order Butterworth high-pass at `lo`
340/// and low-pass at `hi`, each run forward and backward, so the band has no
341/// delay.
342///
343/// Both edges get the ends of [`low_pass_zero_phase`]: an odd reflection of
344/// three periods of `lo` and steady-state starts.
345///
346/// # Errors
347///
348/// Returns a [`FilterError`] when an edge is zero or not below the Nyquist
349/// frequency.
350pub fn band_pass_zero_phase(
351    samples: &[f64],
352    rate: SampleRate,
353    lo: Hz,
354    hi: Hz,
355) -> Result<Vec<f64>, FilterError> {
356    let pad = cast::ceil_to_usize(3.0 * rate.hz().get() / lo.get()).unwrap_or(0);
357    let high = zero_phase(samples, pad, Reflection::Odd, || {
358        LowPass::of(Type::HighPass, rate, lo)
359    })?;
360    zero_phase(&high, pad, Reflection::Odd, || LowPass::new(rate, hi))
361}
362
363/// The envelope of the band from `lo` to `hi`.
364///
365/// The signal is mixed down by `2 exp(-j 2 pi f_c t)` at the centre `f_c`,
366/// low-passed forward and backward at half the bandwidth, and the magnitude
367/// taken. A partial of amplitude `A` inside the band reads `A`.
368///
369/// Each low-pass extends the ends by an even reflection of three periods of
370/// its cutoff, so the ends carry no transient, and memory stays linear in
371/// the length, unlike the FFT of [`AnalyticSignal`].
372///
373/// # Errors
374///
375/// Returns a [`FilterError`] when `lo` is not above 0 Hz, `hi` not below
376/// the Nyquist frequency, or `hi` not above `lo`.
377pub fn band_magnitude(
378    samples: &[f64],
379    rate: SampleRate,
380    lo: Hz,
381    hi: Hz,
382) -> Result<Vec<f64>, FilterError> {
383    let nyquist = rate.nyquist().get();
384    if lo.get() <= 0.0 || lo.get() >= nyquist {
385        return Err(FilterError { cutoff: lo, rate });
386    }
387    if hi.get() <= lo.get() || hi.get() >= nyquist {
388        return Err(FilterError { cutoff: hi, rate });
389    }
390    let centre = f64::midpoint(lo.get(), hi.get());
391    let half = Hz::clamped(0.5 * (hi.get() - lo.get()));
392    let step = TAU * centre / rate.hz().get();
393    let (re, im): (Vec<f64>, Vec<f64>) = samples
394        .iter()
395        .enumerate()
396        .map(|(n, &x)| {
397            // Reduce the phase per sample so it stays exact on long stems.
398            let (sin, cos) = ((step * cast::usize_to_f64(n)) % TAU).sin_cos();
399            (2.0 * x * cos, -2.0 * x * sin)
400        })
401        .unzip();
402    let pad = cast::ceil_to_usize(3.0 * rate.hz().get() / half.get()).unwrap_or(0);
403    let filter = || LowPass::new(rate, half);
404    let re = zero_phase(&re, pad, Reflection::Even, filter)?;
405    let im = zero_phase(&im, pad, Reflection::Even, filter)?;
406    Ok(re.iter().zip(&im).map(|(a, b)| a.hypot(*b)).collect())
407}
408
409/// A fourth-order Butterworth low-pass or high-pass as two biquad sections.
410struct LowPass([(DirectForm2Transposed<f64>, Coefficients<f64>); 2]);
411
412impl LowPass {
413    fn new(rate: SampleRate, cutoff: Hz) -> Result<Self, FilterError> {
414        Self::of(Type::LowPass, rate, cutoff)
415    }
416
417    fn of(kind: Type<f64>, rate: SampleRate, cutoff: Hz) -> Result<Self, FilterError> {
418        let fail = || FilterError { cutoff, rate };
419        if cutoff.get() <= 0.0 || cutoff.get() >= rate.nyquist().get() {
420            return Err(fail());
421        }
422        let fs = rate.hz().get().hz();
423        let f0 = cutoff.get().hz();
424        let section = |q| {
425            Coefficients::<f64>::from_params(kind, fs, f0, q)
426                .map(|c| (DirectForm2Transposed::<f64>::new(c), c))
427                .map_err(|_| fail())
428        };
429        Ok(Self([
430            section(BUTTERWORTH_4[0])?,
431            section(BUTTERWORTH_4[1])?,
432        ]))
433    }
434
435    /// Sets the state that a constant input `x` holds forever, as
436    /// `lfilter_zi` does.
437    fn settle(&mut self, x: f64) {
438        let mut input = x;
439        for (section, c) in &mut self.0 {
440            let output = input * (c.b0 + c.b1 + c.b2) / (1.0 + c.a1 + c.a2);
441            section.s2 = c.b2 * input - c.a2 * output;
442            section.s1 = section.s2 + c.b1 * input - c.a1 * output;
443            input = output;
444        }
445    }
446
447    fn run(&mut self, x: f64) -> f64 {
448        let [(first, _), (second, _)] = &mut self.0;
449        second.run(first.run(x))
450    }
451}
452
453/// A filter cutoff that the sample rate cannot carry.
454#[derive(Clone, Copy, Debug, PartialEq)]
455pub struct FilterError {
456    cutoff: Hz,
457    rate: SampleRate,
458}
459
460impl fmt::Display for FilterError {
461    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
462        write!(
463            f,
464            "a filter edge at {} cannot run at {}: it must lie above 0 Hz and below {}",
465            self.cutoff,
466            self.rate,
467            self.rate.nyquist()
468        )
469    }
470}
471
472impl Error for FilterError {}
473
474#[cfg(test)]
475mod tests {
476    use super::*;
477
478    #[test]
479    fn the_attack_low_pass_constant_passes_the_unit_check() {
480        assert_eq!(ATTACK_LOW_PASS_HZ, Hz::new(150.0).unwrap());
481    }
482
483    fn rate() -> SampleRate {
484        SampleRate::new(8_000).unwrap()
485    }
486
487    #[test]
488    fn the_analytic_signal_of_a_cosine_is_a_unit_phasor() {
489        // An integer number of cycles needs no padding.
490        let x: Vec<f64> = (0..800)
491            .map(|n| (TAU * 100.0 * f64::from(n) / 8000.0).cos())
492            .collect();
493        let z = AnalyticSignal::new(&x, rate(), 0);
494        assert!(z.envelope().iter().all(|a| (a - 1.0).abs() < 1e-12));
495        assert!(
496            z.instantaneous_frequency()
497                .iter()
498                .all(|f| (f - 100.0).abs() < 1e-9)
499        );
500        let phase = z.phase();
501        let expected = TAU * 100.0 * 799.0 / 8000.0;
502        assert!((phase[799] - expected).abs() < 1e-9, "{}", phase[799]);
503    }
504
505    #[test]
506    fn am_depth_of_a_flat_or_empty_envelope() {
507        assert_eq!(am_depth(&[2.0, 2.0]).unwrap().get(), 0.0);
508        assert!((am_depth(&[0.5, 1.5]).unwrap().get() - 0.5).abs() < 1e-15);
509        assert!(am_depth(&[]).is_none());
510        assert!(am_depth(&[0.0, 0.0]).is_none());
511    }
512
513    #[test]
514    fn a_cutoff_at_nyquist_fails_with_the_limit() {
515        let error = low_pass_zero_phase(&[0.0], rate(), Hz::new(4000.0).unwrap()).unwrap_err();
516        assert_eq!(
517            error.to_string(),
518            "a filter edge at 4000 Hz cannot run at 8000 Hz: it must lie above 0 Hz and below 4000 Hz"
519        );
520    }
521
522    #[test]
523    fn the_zero_phase_low_pass_passes_dc_and_keeps_its_length() {
524        let out = low_pass_zero_phase(&vec![1.0; 4000], rate(), Hz::new(150.0).unwrap()).unwrap();
525        assert_eq!(out.len(), 4000);
526        assert!((out[2000] - 1.0).abs() < 1e-9, "{}", out[2000]);
527    }
528
529    #[test]
530    fn the_zero_phase_low_pass_keeps_a_constant_and_a_ramp_at_both_ends() {
531        let cutoff = Hz::new(150.0).unwrap();
532        let ones = low_pass_zero_phase(&vec![1.0; 4000], rate(), cutoff).unwrap();
533        assert!(ones.iter().all(|v| (v - 1.0).abs() < 1e-9), "{ones:?}");
534        let ramp: Vec<f64> = (0..4000).map(|n| f64::from(n) / 4000.0).collect();
535        let out = low_pass_zero_phase(&ramp, rate(), cutoff).unwrap();
536        for (x, y) in ramp.iter().zip(&out) {
537            assert!((x - y).abs() < 1e-5, "{x} {y}");
538        }
539    }
540
541    #[test]
542    fn the_zero_phase_low_pass_takes_one_or_no_sample() {
543        let cutoff = Hz::new(150.0).unwrap();
544        assert_eq!(
545            low_pass_zero_phase(&[], rate(), cutoff).unwrap(),
546            Vec::<f64>::new()
547        );
548        let one = low_pass_zero_phase(&[0.5], rate(), cutoff).unwrap();
549        assert!((one[0] - 0.5).abs() < 1e-12);
550    }
551
552    #[test]
553    fn default_decimation_keeps_four_samples_per_cutoff_period() {
554        let params = HeterodyneParams::new(Hz::new(1000.0).unwrap(), Hz::new(10.0).unwrap());
555        assert_eq!(params.decimation_for(rate()), 200);
556        assert_eq!(params.with_decimation(7).decimation_for(rate()), 7);
557    }
558}
sound_to_patch/dsp/analytic.rs

sound_to_patch/dsp/
analytic.rs
