// tdpsola - an implementation of the TD-PSOLA algorithm
// Copyright (C) 2020 Pieter Penninckx
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU Affero General Public License as
// published by the Free Software Foundation, either version 3 of the
// License, or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU Affero General Public License for more details.
//
// You should have received a copy of the GNU Affero General Public License
// along with this program.  If not, see <https://www.gnu.org/licenses/>.
//
//! Tdpsola: an implementation of the tdpsola algorithm
//! ===================================================
//!
//! TDPSOLA (Time-Domain Pitch-Synchronous Overlap and Add) is an algorithm that allows to modify
//! the pitch and the speed of a sound independently of each other.
//!
//! # Example of reconstructing audio
//! ```
//! use tdpsola::{TdpsolaSynthesis, Speed, AlternatingHann, TdpsolaAnalysis};
//! let input = vec![1.0, 2.0, 3.0, -3.0, -2.0, -1.0, 1.0, 0.0, 1.0];
//! let source_wavelength = 4.0; // Source wavelength
//! let mut alternating_hann = AlternatingHann::new(source_wavelength);
//! let mut analysis = TdpsolaAnalysis::new(&alternating_hann);
//!
//! // Pre-padding, this is used to avoid artifacts at the beginning of the reconstructed audio.
//! // If you do not apply pre-padding and discarding of the first samples, the audio will
//! // in general "fade in" (depending on the window being used) during the first `x` samples
//! // where `x` is the source wavelength.
//! let padding_length = source_wavelength as usize + 1;
//! for _ in 0..padding_length {
//!     analysis.push_sample(0.0, &mut alternating_hann);
//! }
//!
//! // The analysis
//! for sample in input.iter() {
//!     analysis.push_sample(*sample, &mut alternating_hann);
//! }
//!
//! let target_wavelength = source_wavelength;
//! let mut synthesis = TdpsolaSynthesis::new(Speed::from_f32(1.0), target_wavelength);
//!
//!
//! // Re-synthesis and check that it matches with the input.
//!
//! for (input_sample, output_sample) in input.iter().zip(synthesis.iter(&analysis).skip(padding_length)) {
//!     assert!((output_sample - input_sample).abs() < 1e-6);
//! }
//! ```
mod self_clearing_vec_deque;

use self_clearing_vec_deque::{Active, SelfClearingVecDeque};

mod offsetted_vecdeque;

use offsetted_vecdeque::OffsettedVeqDeque;
#[cfg(test)]
use std::collections::vec_deque::VecDeque;
use std::error::Error;
use std::fmt::{Display, Formatter};

/// Describes two alternating windows.
pub trait AlternatingWindow {
    /// Proceed one step.
    fn step(&mut self);
    /// Get the value of window a at the current step.
    /// By convention, this is something between 0.0 and 1.0.
    fn window_a(&self) -> f32;
    /// Get the value of window b at the current step.
    /// By convention, this is something between 0.0 and 1.0.
    fn window_b(&self) -> f32;
    /// Indicates in which phase the windows are.
    /// If `a_is_ending()` is true, this means that
    /// window a is past its peak and is decreasing and that window b is increasing.
    /// If `a_is_ending()` returns true, this means that
    /// window_a is increasing and that window b is past its peak and is decreasing.
    fn a_is_ending(&self) -> bool;
}

/// Two alternating Hann windows.
pub struct AlternatingHann {
    s: f32,
    c: f32,
    s_square: f32,
    c_square: f32,
    s_delta: f32,
    c_delta: f32,
}

impl AlternatingHann {
    /// Create a new `AlternatingHann` with the given impulse wavelength.
    ///
    /// # Note
    /// `a_is_ending()` called on a newly created `AlternatingHann` will return `true`.
    /// # Note
    /// Aliasing occurs if the wavelength is < 1.0.
    pub fn new(impulse_wavelength: f32) -> Self {
        let window_wavelength = impulse_wavelength * 4.0;
        let delta = 2.0 * std::f32::consts::PI / window_wavelength;
        Self {
            s: 0.0,
            c: 1.0,
            s_square: 0.0,
            c_square: 1.0,
            s_delta: delta.sin(),
            c_delta: delta.cos(),
        }
    }
}

impl AlternatingWindow for AlternatingHann {
    #[inline]
    fn step(&mut self) {
        let s = self.s;
        // We apply some basic trigeometry here:
        // sin(t + delta) = cos(t) * sin(delta) + sin(t) * cos(delta).
        self.s = self.c * self.s_delta + self.s * self.c_delta;
        // Similar reasoning here.
        self.c = self.c * self.c_delta - s * self.s_delta;
        self.s_square = self.s * self.s;
        self.c_square = self.c * self.c;
        // Apply a correction factor to avoid drifting in amplitude due to floating point inaccuracies.
        // See https://dsp.stackexchange.com/a/1087 for more info.
        let factor = (3.0 - self.s_square - self.c_square) / 2.0;
        self.s *= factor;
        self.c *= factor;
        // For performance reasons, we do not update self.s_square and self.c_square here.
    }

    #[inline(always)]
    fn window_a(&self) -> f32 {
        self.c_square
    }

    #[inline(always)]
    fn window_b(&self) -> f32 {
        self.s_square
    }

    #[inline(always)]
    fn a_is_ending(&self) -> bool {
        self.s.is_sign_positive() == self.c.is_sign_positive()
    }
}

#[test]
fn a_is_ending_works() {
    let wavelength = 20;
    let mut oscilator = AlternatingHann::new((wavelength / 4) as f32);
    for i in 0..wavelength {
        if i < 5 {
            assert_eq!(true, oscilator.a_is_ending());
        } else if i > 5 && i < 10 {
            assert_eq!(false, oscilator.a_is_ending());
        } else if i > 10 && i < 15 {
            assert_eq!(true, oscilator.a_is_ending());
        } else if i > 15 {
            assert_eq!(false, oscilator.a_is_ending());
        }
        oscilator.step();
    }
}

#[test]
fn alternating_hann_window_works() {
    let wavelength = 200;
    let mut oscilator = AlternatingHann::new(wavelength as f32);
    for i in 0..wavelength {
        let expected_phase = (i as f32) * 2.0 * std::f32::consts::PI / ((4 * wavelength) as f32);
        assert!((expected_phase.sin() - oscilator.s).abs() < 1e-6);
        assert!((expected_phase.cos() - oscilator.c).abs() < 1e-6);
        oscilator.step();
    }
}

/// Two alternating hat windows.
/// This should be more performant than the `AlternatingHann`.
pub struct AlternatingHat {
    step: f32,
    current_value: f32,
}

impl AlternatingHat {
    pub fn new(wavelength: f32) -> Self {
        Self {
            step: 1.0 / wavelength,
            current_value: 0.0,
        }
    }
}

impl AlternatingWindow for AlternatingHat {
    fn step(&mut self) {
        self.current_value += self.step;
        if self.current_value < 0.0 {
            self.current_value = -self.current_value;
            self.step = -self.step;
        }
        if self.current_value > 1.0 {
            self.current_value = 2.0 - self.current_value;
            self.step = -self.step;
        }
    }

    fn window_a(&self) -> f32 {
        self.current_value
    }

    fn window_b(&self) -> f32 {
        1.0 - self.current_value
    }

    fn a_is_ending(&self) -> bool {
        self.step < 0.0
    }
}

#[test]
fn alternating_hat_window_works() {
    let wavelength = 4;
    let mut oscilator = AlternatingHat::new(wavelength as f32);
    let expected_a = vec![0.0, 0.25, 0.5, 0.75, 1.0, 0.75, 0.5, 0.25, 0.0, 0.25];
    let expected_b = vec![1.0, 0.75, 0.5, 0.25, 0.0, 0.25, 0.5, 0.75, 1.0, 0.75];

    for i in 0..10 {
        assert_eq!(oscilator.window_a(), expected_a[i]);
        assert_eq!(oscilator.window_b(), expected_b[i]);
        oscilator.step();
    }
}

// Documentation about the internal representation:
// `window_a_samples` contains the input, pointwise multiplied by window a.
// `window_b_samples` contains the input, pointwise multiplied by window b.
// `markers` indicates where `a_is_ending` switches
pub struct TdpsolaAnalysis {
    window_a_samples: OffsettedVeqDeque<f32>,
    window_b_samples: OffsettedVeqDeque<f32>,
    markers: OffsettedVeqDeque<u64>,
    current_a_is_ending: bool,
}

impl TdpsolaAnalysis {
    /// Create a new `TdpsolaAnalysis`.
    /// The `window` is used to initialise the `TdpsolaAnalysis` and is assumed to be the first
    /// window used when calling the `push_sample` method.
    pub fn new<W: AlternatingWindow>(window: &W) -> Self {
        // TODO: specify capacity.
        let window_a_samples = OffsettedVeqDeque::with_capacity(0);
        let window_b_samples = OffsettedVeqDeque::with_capacity(0);
        let mut markers = OffsettedVeqDeque::with_capacity(0);
        markers.push(0);
        Self {
            window_a_samples,
            window_b_samples,
            markers,
            current_a_is_ending: window.a_is_ending(),
        }
    }

    fn start_new_impulse(&mut self) {
        self.markers.push(self.window_a_samples.len());
        self.current_a_is_ending = !self.current_a_is_ending;
    }

    /// Analyse a new sample.
    ///
    /// Different structs, even of different types can be used in subsequent calls to this method.
    /// It is assumed that the window changes fluidly between subsequent calls to `push_sample`.
    ///
    /// In order to avoid edge effects, it is advised to add one window of pre-padding.
    pub fn push_sample<W: AlternatingWindow>(&mut self, sample: f32, window: &mut W) {
        self.window_a_samples.push(sample * window.window_a());
        self.window_b_samples.push(sample * window.window_b());
        window.step();
        if window.a_is_ending() != self.current_a_is_ending {
            self.start_new_impulse();
        }
    }

    fn get_overlap(&self, time: Time) -> AlternatingOverlap {
        let approximate_sample_index = time.floor();
        let closest_marker = self.markers.binary_closest(&approximate_sample_index);
        AlternatingOverlap {
            window_a: closest_marker % 2 == 1,
            overlap: Overlap {
                is_active: true,
                start_index: self
                    .markers
                    .get(closest_marker)
                    .cloned()
                    .expect("`binary_closest` should only return a valid index"),
                end_marker_index: closest_marker + 2,
            },
        }
    }

    fn get_first_overlap() -> AlternatingOverlap {
        let closest_marker = 0;
        AlternatingOverlap {
            window_a: closest_marker % 2 == 1,
            overlap: Overlap {
                is_active: true,
                start_index: 0,
                end_marker_index: closest_marker + 2,
            },
        }
    }
}

#[test]
fn get_overlap_works() {
    struct TestWindow {
        index: usize,
    }

    impl AlternatingWindow for TestWindow {
        fn step(&mut self) {
            self.index += 1;
        }

        fn window_a(&self) -> f32 {
            1.0
        }

        fn window_b(&self) -> f32 {
            2.0
        }

        fn a_is_ending(&self) -> bool {
            match self.index % 4 {
                0 => true,
                1 => true,
                2 => false,
                3 => false,
                _ => unreachable!(),
            }
        }
    }
    let mut window = TestWindow { index: 0 };
    let mut analysis = TdpsolaAnalysis::new(&window);
    analysis.push_sample(1.0, &mut window);
    analysis.push_sample(2.0, &mut window);
    analysis.push_sample(3.0, &mut window);
    analysis.push_sample(4.0, &mut window);
    analysis.push_sample(5.0, &mut window);
    analysis.push_sample(6.0, &mut window);
    analysis.push_sample(7.0, &mut window);
    analysis.push_sample(8.0, &mut window);
    let expected = AlternatingOverlap {
        window_a: false,
        overlap: Overlap {
            start_index: 0,
            is_active: true,
            end_marker_index: 2,
        },
    };
    assert_eq!(analysis.get_overlap(Time::from_u64(0)), expected);
    assert_eq!(analysis.get_overlap(Time::from_u64(1)), expected);
    let expected = AlternatingOverlap {
        window_a: true,
        overlap: Overlap {
            start_index: 2,
            is_active: true,
            end_marker_index: 3,
        },
    };
    assert_eq!(analysis.get_overlap(Time::from_u64(2)), expected);
    assert_eq!(analysis.get_overlap(Time::from_u64(3)), expected);
    let expected = AlternatingOverlap {
        window_a: false,
        overlap: Overlap {
            start_index: 4,
            is_active: true,
            end_marker_index: 4,
        },
    };
    assert_eq!(analysis.get_overlap(Time::from_u64(4)), expected);
    assert_eq!(analysis.get_overlap(Time::from_u64(5)), expected);
    let expected = AlternatingOverlap {
        window_a: true,
        overlap: Overlap {
            start_index: 6,
            is_active: true,
            end_marker_index: 5,
        },
    };
    assert_eq!(analysis.get_overlap(Time::from_u64(6)), expected);
    assert_eq!(analysis.get_overlap(Time::from_u64(7)), expected);
    let expected = AlternatingOverlap {
        window_a: false,
        overlap: Overlap {
            start_index: 8,
            is_active: true,
            end_marker_index: 6,
        },
    };
    assert_eq!(analysis.get_overlap(Time::from_u64(8)), expected);
    assert_eq!(analysis.get_overlap(Time::from_u64(9)), expected);
}

#[test]
fn push_sample_works() {
    struct AlternatingWindowMock {
        index: usize,
    }

    #[cfg(test)]
    impl AlternatingWindow for AlternatingWindowMock {
        fn step(&mut self) {
            self.index += 1;
        }

        fn window_a(&self) -> f32 {
            [4.0, 2.0, 1.0, 3.0][self.index % 4]
        }

        fn window_b(&self) -> f32 {
            [1.0, 3.0, 5.0, 2.0][self.index % 4]
        }

        fn a_is_ending(&self) -> bool {
            self.index % 4 == 0 || self.index % 4 == 1
        }
    }

    let mut window = AlternatingWindowMock { index: 0 };
    let mut ir = TdpsolaAnalysis::new(&window);
    ir.push_sample(1.0, &mut window);
    ir.push_sample(2.0, &mut window);
    ir.push_sample(3.0, &mut window);
    ir.push_sample(4.0, &mut window);
    ir.push_sample(5.0, &mut window);
    ir.push_sample(6.0, &mut window);
    ir.push_sample(7.0, &mut window);
    ir.push_sample(8.0, &mut window);
    // a: 4 2 1 3
    // b: 1 3 5 2
    let expected_a = vec![
        1.0 * 4.0,
        2.0 * 2.0,
        3.0 * 1.0,
        4.0 * 3.0,
        5.0 * 4.0,
        6.0 * 2.0,
        7.0 * 1.0,
        8.0 * 3.0,
    ];
    let expected_b = vec![
        1.0 * 1.0,
        2.0 * 3.0,
        3.0 * 5.0,
        4.0 * 2.0,
        5.0 * 1.0,
        6.0 * 3.0,
        7.0 * 5.0,
        8.0 * 2.0,
    ];
    assert_eq!(ir.window_a_samples.len(), expected_a.len() as u64);
    assert_eq!(ir.window_b_samples.len(), expected_b.len() as u64);
    for i in 0..expected_a.len() {
        assert_eq!(expected_a[i], *ir.window_a_samples.get(i as u64).unwrap());
    }
    for i in 0..expected_b.len() {
        assert_eq!(expected_b[i], *ir.window_b_samples.get(i as u64).unwrap());
    }
    let expected_markers = vec![0, 2, 4, 6, 8];
    assert_eq!(ir.markers.len(), expected_markers.len() as u64);
    for i in 0..expected_markers.len() {
        assert_eq!(*ir.markers.get(i as u64).unwrap(), expected_markers[i]);
    }
}

const SPEED_PRECISION: u32 = 16;

/// The playback speed.
#[derive(Clone, Copy, PartialEq, Eq, PartialOrd, Ord)]
pub struct Speed {
    // 32 bits fixed point number with 16 bits before the binary fixed point
    inner: u32,
}

impl Speed {
    /// Create a new `Speed`.
    /// Speeds `> 1` represent fast-forward.
    /// Speeds `< 1` represent slower speeds.
    /// Panics if the speed is `< 0`.
    pub fn from_f32(speed: f32) -> Self {
        assert!(speed >= 0.0);
        Self {
            inner: ((speed as f64) * ((1 << SPEED_PRECISION) as f64)) as u32,
        }
    }
}

const TIME_PRECISION: u64 = 16;

#[derive(Clone, Copy, PartialEq, Eq, PartialOrd, Ord)]
struct Time {
    // inner: 64 bits fixed point integer with TIME_PRECISION bits after the binary fixed point
    inner: u64,
}

impl Time {
    fn from_u64(time: u64) -> Self {
        Self {
            inner: time << TIME_PRECISION,
        }
    }
    fn increment(&mut self, speed: Speed) {
        self.inner += speed.inner as u64;
    }
    fn floor(self) -> u64 {
        self.inner >> TIME_PRECISION
    }
}

#[derive(PartialEq, Eq, Debug)]
struct Overlap {
    start_index: u64,
    end_marker_index: u64,
    is_active: bool,
}

#[derive(Debug, PartialEq, Eq)]
pub enum SamplingError {
    NotEnoughSamplesAvailable,
}

impl Display for SamplingError {
    fn fmt(&self, f: &mut Formatter) -> std::fmt::Result {
        write!(f, "Not enough samples available")
    }
}

impl Error for SamplingError {
    fn source(&self) -> Option<&(dyn Error + 'static)> {
        None
    }
}

impl Overlap {
    /// Return `Ok(None)` when no future samples can be expected.
    /// Return an error if not enough samples available yet.
    fn try_get_sample(
        &self,
        source: &OffsettedVeqDeque<f32>,
    ) -> Result<Option<f32>, SamplingError> {
        if !self.is_active {
            return Ok(None);
        }
        let result = source
            .get(self.start_index)
            .ok_or(SamplingError::NotEnoughSamplesAvailable)?;
        Ok(Some(*result))
    }

    // Return `true` if the overlap is still active and `false` otherwise.
    fn step(&mut self, markers: &OffsettedVeqDeque<u64>) -> bool {
        self.start_index += 1;
        // Performance note: we're always doing a lookup here, may be worth caching
        // the result.
        if let Some(end_index) = markers.get(self.end_marker_index).cloned() {
            if self.start_index >= end_index {
                self.is_active = false;
            }
        }
        self.is_active
    }
}

#[test]
fn overlap_try_get_sample_works() {
    let partial_source = OffsettedVeqDeque {
        inner: VecDeque::from(vec![1.0]),
        offset: 2,
    };
    let source = OffsettedVeqDeque {
        inner: VecDeque::from(vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0]),
        offset: 2,
    };
    let markers = OffsettedVeqDeque {
        inner: VecDeque::from(vec![0, 4, 8]),
        // Corresponds to ??, 3.0, 7.0
        offset: 1,
    };
    let mut overlap = Overlap {
        start_index: 2,
        is_active: true,
        end_marker_index: 2, // Corresponds to 4
    };
    assert_eq!(Ok(Some(1.0)), overlap.try_get_sample(&source));
    assert_eq!(Ok(Some(1.0)), overlap.try_get_sample(&partial_source));
    assert!(overlap.step(&markers));
    assert_eq!(Ok(Some(2.0)), overlap.try_get_sample(&source));
    assert_eq!(
        Err(SamplingError::NotEnoughSamplesAvailable),
        overlap.try_get_sample(&partial_source)
    );
    assert!(!overlap.step(&markers));
    assert_eq!(Ok(None), overlap.try_get_sample(&source));
    assert!(!overlap.step(&markers));
    assert_eq!(Ok(None), overlap.try_get_sample(&source));
}

impl Active for Overlap {
    fn is_active(&self) -> bool {
        self.is_active
    }
}

#[derive(PartialEq, Eq, Debug)]
struct AlternatingOverlap {
    window_a: bool,
    overlap: Overlap,
}

impl AlternatingOverlap {
    /// Return `Ok(None)` when the overlap is no longer active.
    fn try_get_sample(&self, analysis: &TdpsolaAnalysis) -> Result<Option<f32>, SamplingError> {
        if self.window_a {
            self.overlap.try_get_sample(&analysis.window_a_samples)
        } else {
            self.overlap.try_get_sample(&analysis.window_b_samples)
        }
    }

    /// Return `true` when still active and `false` otherwise.
    fn step(&mut self, analysis: &TdpsolaAnalysis) -> bool {
        self.overlap.step(&analysis.markers)
    }
}

impl Active for AlternatingOverlap {
    fn is_active(&self) -> bool {
        self.overlap.is_active()
    }
}

/// Struct used to re-synthesize a `TdpsolaSynthesis` into audio again.
pub struct TdpsolaSynthesis {
    speed: Speed,
    source_time: Time,
    windows: SelfClearingVecDeque<AlternatingOverlap>,
    wavelength: f32,
    remaining_wavelength: f32,
}

pub struct SynthesisIter<'s, 'a> {
    synthesis: &'s mut TdpsolaSynthesis,
    analysis: &'a TdpsolaAnalysis,
}

impl<'s, 'a> Iterator for SynthesisIter<'s, 'a> {
    type Item = f32;

    fn next(&mut self) -> Option<Self::Item> {
        match self.synthesis.try_get_sample(self.analysis) {
            Ok(sample) => {
                self.synthesis.step(self.analysis);
                Some(sample)
            }
            Err(_) => None,
        }
    }
}

impl TdpsolaSynthesis {
    /// Create a new `TdpsolaSynthesis`.
    ///
    /// `initial_speed` is the initial speed used to synthesize the audio.
    /// `wavelength` is the wavelength of the synthesized audio.
    pub fn new(initial_speed: Speed, wavelength: f32) -> Self {
        let mut windows = SelfClearingVecDeque::with_capacity(0); // TODO: specify capacity.
        windows.push(TdpsolaAnalysis::get_first_overlap());
        Self {
            speed: initial_speed,
            source_time: Time::from_u64(0),
            windows,
            wavelength,
            remaining_wavelength: wavelength,
        }
    }

    /// Set the new speed to be used in reconstruction.
    pub fn set_speed(&mut self, new_speed: Speed) {
        self.speed = new_speed;
    }

    pub fn set_wavelength(&mut self, new_wavelength: f32) {
        self.remaining_wavelength *= new_wavelength / self.wavelength;
        self.wavelength = new_wavelength;
    }

    /// Get the current re-synthesized sample.
    ///
    /// # Note
    /// The first few samples (the first `n` samples where `n` equals the wavelength,
    /// rounded up) are different, preventing "perfect" reconstruction. This can be
    /// compensated for by pre-padding the input with `n` zero samples and discarding
    /// the first `n` samples in the reconstructed data.
    pub fn try_get_sample(&self, analysis: &TdpsolaAnalysis) -> Result<f32, SamplingError> {
        let mut sample = 0.0;
        for window in self.windows.iter() {
            sample += window.try_get_sample(analysis)?.unwrap_or(0.0)
        }
        Ok(sample)
    }

    pub fn step(&mut self, analysis: &TdpsolaAnalysis) {
        self.source_time.increment(self.speed);
        self.remaining_wavelength -= 1.0;
        for window in self.windows.iter_mut() {
            window.step(analysis);
        }
        if self.remaining_wavelength <= 0.0 {
            self.remaining_wavelength += self.wavelength;
            self.windows.push(analysis.get_overlap(self.source_time));
        }
    }

    pub fn iter<'s, 'a>(&'s mut self, analysis: &'a TdpsolaAnalysis) -> SynthesisIter<'s, 'a> {
        SynthesisIter {
            synthesis: self,
            analysis,
        }
    }
}

#[test]
fn sample_works() {
    let a = vec![2.0, -2.0, 4.0, -4.0, 8.0, -8.0, 16.0, -16.0];
    let b = vec![1.0, -1.0, 3.0, -3.0, 5.0, -5.0, 7.0, -7.0];
    let analysis = TdpsolaAnalysis {
        window_a_samples: OffsettedVeqDeque {
            offset: 0,
            inner: VecDeque::from(a.clone()),
        },
        window_b_samples: OffsettedVeqDeque {
            offset: 0,
            inner: VecDeque::from(b.clone()),
        },
        markers: OffsettedVeqDeque {
            offset: 0,
            inner: VecDeque::from(vec![0, 2, 5]),
        },
        current_a_is_ending: true, // shouldn't matter here.
    };
    let mut synthesis = TdpsolaSynthesis::new(Speed::from_f32(0.5), 2.0);
    assert_eq!(synthesis.try_get_sample(&analysis), Ok(b[0]));
    synthesis.step(&analysis);
    assert_eq!(synthesis.try_get_sample(&analysis), Ok(b[1]));
    synthesis.step(&analysis);
    assert_eq!(synthesis.try_get_sample(&analysis), Ok(b[2] + b[0]));
    synthesis.step(&analysis);
    assert_eq!(synthesis.try_get_sample(&analysis), Ok(b[3] + b[1]));
    synthesis.step(&analysis);
    assert_eq!(synthesis.try_get_sample(&analysis), Ok(b[4] + b[2] + a[2]));
    synthesis.step(&analysis);
    assert_eq!(synthesis.try_get_sample(&analysis), Ok(b[3] + a[3]));
    synthesis.step(&analysis);
    assert_eq!(synthesis.try_get_sample(&analysis), Ok(b[4] + a[4] + a[2]));
    synthesis.step(&analysis);
}