Fix phase response in the Crossover plugin

Didn't have time to do this, now I do. This nudges the phases from the
lower bands to match the higher bands, making the frequency response
magnitudes sum to unity again.

plugins/crossover/src/crossover/iir.rs

@@ -22,6 +22,8 @@ use std::simd::f32x2;
 
 use crate::NUM_BANDS;
 
+const NEUTRAL_Q: f32 = std::f32::consts::FRAC_1_SQRT_2;
+
 #[derive(Debug)]
 pub struct IirCrossover {
     /// The kind of crossover to use. `.update_filters()` must be called after changing this.
@@ -30,6 +32,8 @@ pub struct IirCrossover {
     /// The crossovers. Depending on the number of bands argument passed to `.process()` one to four
     /// of these may be used.
     crossovers: [Crossover; NUM_BANDS - 1],
+    /// Used to compensate the earlier bands for the phase shift introduced in the higher bands.
+    all_passes: AllPassCascade,
 }
 
 /// The type of IIR crossover to use.
@@ -52,6 +56,26 @@ struct Crossover {
     hp_filters: [Biquad<f32x2>; 2],
 }
 
+/// The crossover is super simple and feeds the low-passed result to the next band output while
+/// using the high-passed version as the input for the next band. Because the higher bands will thus
+/// have had more filters applied to them, the lower bands need to have their phase response
+/// adjusted to match the higher bands. So for the LR24 crossovers, low-passed band `n` will get a
+/// second order all-pass for the frequencies corresponding to crossovers `n + 1..NUM_CROSSOVERS`
+/// applied to it.
+#[derive(Debug, Default)]
+struct AllPassCascade {
+    /// The aforementioned all-pass filters. This is indexed by `[crossover_idx][0..num_bands -
+    /// crossover_index - 1]`. Ergo, if there are three crossovers, then the low-pass section from
+    /// the first crossover needs to have `[0][0]` and `[0][1]` applied to it. The last band doesn't
+    /// need any compensation, hence the `NUM_BANDS - 2`. The outer array is equal to the number of
+    /// crossovers. It will never contain any filters, but this makes the code a bit nicer by
+    /// needing an explicit check for this.
+    ap_filters: [[Biquad<f32x2>; NUM_BANDS - 2]; NUM_BANDS - 1],
+
+    /// The number of activate bands. Only coefficients for used bands are computed in `ap_filters`.
+    num_bands: usize,
+}
+
 impl IirCrossover {
     /// Create a new multiband crossover processor. All filters will be configured to pass audio
     /// through as it. `.update()` needs to be called first to set up the filters, and `.reset()`
@@ -60,6 +84,7 @@ impl IirCrossover {
         Self {
             mode,
             crossovers: Default::default(),
+            all_passes: Default::default(),
         }
     }
 
@@ -81,14 +106,19 @@ impl IirCrossover {
         let mut samples: f32x2 = unsafe { main_io.to_simd_unchecked() };
         match self.mode {
             IirCrossoverType::LinkwitzRiley24 => {
-                for (crossover, band_channel_samples) in self
+                for (crossover_idx, (crossover, band_channel_samples)) in self
                     .crossovers
                     .iter_mut()
                     .zip(band_outputs.iter_mut())
                     .take(num_bands as usize - 1)
+                    .enumerate()
                 {
                     let (lp_samples, hp_samples) = crossover.process_lr24(samples);
 
+                    // The low-pass result needs to have the same phase shift applied to it that
+                    // higher bands would get
+                    let lp_samples = self.all_passes.compensate_lr24(lp_samples, crossover_idx);
+
                     unsafe { band_channel_samples.from_simd_unchecked(lp_samples) };
                     samples = hp_samples;
                 }
@@ -106,31 +136,34 @@ impl IirCrossover {
         &mut self,
         sample_rate: f32,
         num_bands: usize,
-        mut frequencies: [f32; NUM_BANDS - 1],
+        frequencies: [f32; NUM_BANDS - 1],
     ) {
-        // Make sure the frequencies are monotonic by pushing bands down when they are too close to
+        // TODO: Can probably get rid of the /2 now
-        // the next band
+        // // Make sure the frequencies are monotonic by pushing bands down when they are too close to
-        for frequency_idx in (1..num_bands - 1).rev() {
+        // // the next band
-            if frequencies[frequency_idx - 1] > frequencies[frequency_idx] / 2.0 {
+        // for frequency_idx in (1..num_bands - 1).rev() {
-                frequencies[frequency_idx - 1] = frequencies[frequency_idx] / 2.0;
+        //     if frequencies[frequency_idx - 1] > frequencies[frequency_idx] / 2.0 {
-            }
+        //         frequencies[frequency_idx - 1] = frequencies[frequency_idx] / 2.0;
-        }
+        //     }
+        // }
 
         match self.mode {
             IirCrossoverType::LinkwitzRiley24 => {
-                const Q: f32 = std::f32::consts::FRAC_1_SQRT_2;
                 for (crossover, frequency) in self
                     .crossovers
                     .iter_mut()
                     .zip(frequencies)
                     .take(num_bands - 1)
                 {
-                    let lp_coefs = BiquadCoefficients::lowpass(sample_rate, frequency, Q);
+                    let lp_coefs = BiquadCoefficients::lowpass(sample_rate, frequency, NEUTRAL_Q);
-                    let hp_coefs = BiquadCoefficients::highpass(sample_rate, frequency, Q);
+                    let hp_coefs = BiquadCoefficients::highpass(sample_rate, frequency, NEUTRAL_Q);
                     crossover.update_coefficients(lp_coefs, hp_coefs);
                 }
             }
         }
+
+        self.all_passes
+            .update_coefficients(sample_rate, num_bands, &frequencies);
     }
 
     /// Reset the internal filter state for all crossovers.
@@ -138,6 +171,8 @@ impl IirCrossover {
         for crossover in &mut self.crossovers {
             crossover.reset();
         }
+
+        self.all_passes.reset();
     }
 }
 
@@ -183,6 +218,69 @@ impl Crossover {
     }
 }
 
+impl AllPassCascade {
+    /// Compensate lower bands for the additional phase shift introduced in higher bands when using
+    /// LR24 filters to split those bands.
+    pub fn compensate_lr24(&mut self, lp_samples: f32x2, band_idx: usize) -> f32x2 {
+        // The all-pass filters are set up based on the crossover that produced the low-passed
+        // samples
+        let crossover_idx = band_idx;
+
+        // The idea here is that if `band_idx == 0`, and `self.num_bands == 3`, then there are two
+        // crossovers, and `lp_samples` only needs to be filtered by `self.ap_filters[0][0]`. If
+        // `self.num_bands` were 4 then it would additionally also be filtered by
+        // `self.ap_filters[0][1]`.
+        let mut compensated = lp_samples;
+        for filter in &mut self.ap_filters[crossover_idx][..self.num_bands - band_idx - 2] {
+            compensated = filter.process(compensated)
+        }
+
+        compensated
+    }
+
+    /// Update the coefficients for all filters in the cascade. For every active band, this adds up
+    /// to `num_bands - band_idx - 1` filters. The filter state of course cannot be shared between
+    /// bands, but the coefficients along the matrix's diagonals are identical.
+    pub fn update_coefficients(
+        &mut self,
+        sample_rate: f32,
+        num_bands: usize,
+        frequencies: &[f32; NUM_BANDS - 1],
+    ) {
+        self.num_bands = num_bands;
+
+        // All output bands go through the first filter, so we don't compensate for that. `band_idx`
+        // starts at 1
+        for (crossover_idx, crossover_frequency) in
+            frequencies.iter().enumerate().take(num_bands - 1).skip(1)
+        {
+            let ap_coefs =
+                BiquadCoefficients::allpass(sample_rate, *crossover_frequency, NEUTRAL_Q);
+
+            // This sets the coefficients in a diagonal pattern. If `crossover_idx == 2`, then this
+            // will set the coefficients for these filters:
+            // ```
+            // [_, x, ...] // Crossover 1 filters
+            // [x, ...]    // Crossover 2 filters
+            // ...
+            // ```
+            for target_crossover_idx in 0..crossover_idx {
+                self.ap_filters[target_crossover_idx][crossover_idx - target_crossover_idx - 1]
+                    .coefficients = ap_coefs;
+            }
+        }
+    }
+
+    /// Reset the internal filter state.
+    pub fn reset(&mut self) {
+        for filters in &mut self.ap_filters {
+            for filter in filters.iter_mut() {
+                filter.reset();
+            }
+        }
+    }
+}
+
 /// A simple biquad filter with functions for generating coefficients for second order low-pass and
 /// high-pass filters. Since these filters have 3 dB of attenuation at the center frequency, we'll
 /// two of them in series to get 6 dB of attenutation at the crossover point for the LR24
@@ -318,6 +416,30 @@ impl<T: SimdType> BiquadCoefficients<T> {
 
         Self::from_f32s(BiquadCoefficients { b0, b1, b2, a1, a2 })
     }
+
+    /// Compute the coefficients for an all-pass filter.
+    ///
+    /// Based on <http://shepazu.github.io/Audio-EQ-Cookbook/audio-eq-cookbook.html>.
+    pub fn allpass(sample_rate: f32, frequency: f32, q: f32) -> Self {
+        nih_debug_assert!(sample_rate > 0.0);
+        nih_debug_assert!(frequency > 0.0);
+        nih_debug_assert!(frequency < sample_rate / 2.0);
+        nih_debug_assert!(q > 0.0);
+
+        let omega0 = consts::TAU * (frequency / sample_rate);
+        let cos_omega0 = omega0.cos();
+        let alpha = omega0.sin() / (2.0 * q);
+
+        // We'll prenormalize everything with a0
+        let a0 = 1.0 + alpha;
+        let b0 = (1.0 - alpha) / a0;
+        let b1 = (-2.0 * cos_omega0) / a0;
+        let b2 = (1.0 + alpha) / a0;
+        let a1 = (-2.0 * cos_omega0) / a0;
+        let a2 = (1.0 - alpha) / a0;
+
+        Self::from_f32s(BiquadCoefficients { b0, b1, b2, a1, a2 })
+    }
 }
 
 impl SimdType for f32 {