From 330d6d13592ee770c668f6068b01390aa2415d0b Mon Sep 17 00:00:00 2001 From: Robbert van der Helm Date: Thu, 2 Jun 2022 14:30:48 +0200 Subject: [PATCH] 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 | 146 +++++++++++++++++++++++-- 1 file changed, 134 insertions(+), 12 deletions(-) diff --git a/plugins/crossover/src/crossover/iir.rs b/plugins/crossover/src/crossover/iir.rs index 70c5e964..8c22570d 100644 --- a/plugins/crossover/src/crossover/iir.rs +++ b/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; 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; 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 - // the next band - for frequency_idx in (1..num_bands - 1).rev() { - if frequencies[frequency_idx - 1] > frequencies[frequency_idx] / 2.0 { - frequencies[frequency_idx - 1] = frequencies[frequency_idx] / 2.0; - } - } + // TODO: Can probably get rid of the /2 now + // // Make sure the frequencies are monotonic by pushing bands down when they are too close to + // // the next band + // for frequency_idx in (1..num_bands - 1).rev() { + // 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 hp_coefs = BiquadCoefficients::highpass(sample_rate, frequency, Q); + let lp_coefs = BiquadCoefficients::lowpass(sample_rate, frequency, NEUTRAL_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 BiquadCoefficients { Self::from_f32s(BiquadCoefficients { b0, b1, b2, a1, a2 }) } + + /// Compute the coefficients for an all-pass filter. + /// + /// Based on . + 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 {