1
0
Fork 0

Parameterize FirCoefficients over the kernel size

This commit is contained in:
Robbert van der Helm 2022-06-07 15:32:29 +02:00
parent b32cd27e8c
commit 7e3dfe904d
2 changed files with 27 additions and 29 deletions

View file

@ -200,7 +200,7 @@ impl FirCrossover {
for coef in fir_hp_coefs.0.iter_mut() {
*coef = -*coef;
}
fir_hp_coefs.0[FILTER_SIZE / 2] += f32x2::splat(1.0);
fir_hp_coefs.0[FILTER_SIZE / 2] += 1.0;
self.band_filters[num_bands - 1].coefficients = fir_hp_coefs;
}

View file

@ -26,10 +26,7 @@ use crate::crossover::iir::biquad::{Biquad, BiquadCoefficients};
pub struct FirFilter {
/// The coefficients for this filter. The filters for both channels should be equivalent, this
/// just avoids broadcasts in the filter process.
///
/// TODO: Profile to see if storing this as f32x2 rather than f32s plus splatting makes any
/// difference in performance at all
pub coefficients: FirCoefficients,
pub coefficients: FirCoefficients<FILTER_SIZE>,
/// A ring buffer storing the last `FILTER_SIZE - 1` samples. The capacity is `FILTER_SIZE`
/// rounded up to the next power of two.
@ -40,11 +37,11 @@ pub struct FirFilter {
delay_buffer_next_idx: usize,
}
/// Coefficients for an FIR filter. This struct includes ways to design the filter. Parameterized
/// over `f32x2` only for the time being since that's what we need here.
/// Coefficients for a (linear-phase) FIR filter. This struct includes ways to design the filter.
/// `T` is the sample type and `N` is the number of taps/coefficients and should be odd for linear-phase filters.
#[repr(transparent)]
#[derive(Debug, Clone)]
pub struct FirCoefficients(pub [f32x2; FILTER_SIZE]);
pub struct FirCoefficients<const N: usize>(pub [f32; N]);
impl Default for FirFilter {
fn default() -> Self {
@ -56,12 +53,12 @@ impl Default for FirFilter {
}
}
impl Default for FirCoefficients {
impl<const N: usize> Default for FirCoefficients<N> {
fn default() -> Self {
// Initialize this to a delay with the same amount of latency as we'd introduce with our
// linear-phase filters
let mut coefficients = [f32x2::default(); FILTER_SIZE];
coefficients[FILTER_SIZE / 2] = f32x2::splat(1.0);
let mut coefficients = [0.0; N];
coefficients[N / 2] = 1.0;
Self(coefficients)
}
@ -73,7 +70,7 @@ impl FirFilter {
// TODO: Replace direct convolution with FFT convolution, would make the implementation much
// more complex though because of the multi output part
let coefficients = &self.coefficients.0;
let mut result = coefficients[0] * samples;
let mut result = f32x2::splat(coefficients[0]) * samples;
// Now multiply `self.coefficients[1..]` with the delay buffer starting at
// `self.delay_buffer_next_idx - 1`, wrapping around to the end when that is reached
@ -89,7 +86,7 @@ impl FirFilter {
.rev(),
) {
// `result += coefficient * sample`, but with explicit FMA
result = coefficient.mul_add(*delayed_sample, result);
result = f32x2::splat(*coefficient).mul_add(*delayed_sample, result);
}
let after_wraparound_begin_idx =
@ -100,7 +97,7 @@ impl FirFilter {
.iter()
.rev(),
) {
result = coefficient.mul_add(*delayed_sample, result);
result = f32x2::splat(*coefficient).mul_add(*delayed_sample, result);
}
// And finally write the samples to the delay buffer for the enxt sample
@ -117,7 +114,7 @@ impl FirFilter {
}
}
impl FirCoefficients {
impl<const N: usize> FirCoefficients<N> {
/// A somewhat crude but very functional and relatively fast way create linear phase FIR
/// **low-pass** filter that matches the frequency response of a fourth order biquad low-pass
/// filter. As in, this matches the frequency response magnitudes of applying those biquads to a
@ -150,23 +147,24 @@ impl FirCoefficients {
///
/// The corresponding high-pass filter can be computed through spectral inversion.
pub fn design_fourth_order_linear_phase_low_pass_from_biquad(
biquad_coefs: BiquadCoefficients<f32x2>,
biquad_coefs: BiquadCoefficients<f32>,
) -> Self {
const CENTER_IDX: usize = FILTER_SIZE / 2;
// Ruest doesn't allow you to define this as a constant
let center_idx = N / 2;
// We'll start with an impulse (at exactly half of this odd sized buffer)...
let mut impulse_response = [f32x2::default(); FILTER_SIZE];
impulse_response[CENTER_IDX] = f32x2::splat(1.0);
let mut impulse_response = [0.0; N];
impulse_response[center_idx] = 1.0;
// ...and filter that in both directions
let mut biquad = Biquad::default();
biquad.coefficients = biquad_coefs;
for sample in impulse_response.iter_mut().skip(CENTER_IDX - 1) {
for sample in impulse_response.iter_mut().skip(center_idx - 1) {
*sample = biquad.process(*sample);
}
biquad.reset();
for sample in impulse_response.iter_mut().skip(CENTER_IDX - 1).rev() {
for sample in impulse_response.iter_mut().skip(center_idx - 1).rev() {
*sample = biquad.process(*sample);
}
@ -176,19 +174,19 @@ impl FirCoefficients {
// Adopted from `nih_plug::util::window`. We only end up applying the right half of the
// window, starting at the top of the window.
let blackman_scale_1 = (2.0 * f32::consts::PI) / (impulse_response.len() - 1) as f32;
let blackman_scale_1 = (2.0 * f32::consts::PI) / (N - 1) as f32;
let blackman_scale_2 = blackman_scale_1 * 2.0;
for (sample_idx, sample) in impulse_response.iter_mut().enumerate().skip(CENTER_IDX - 1) {
for (sample_idx, sample) in impulse_response.iter_mut().enumerate().skip(center_idx - 1) {
let cos_1 = (blackman_scale_1 * sample_idx as f32).cos();
let cos_2 = (blackman_scale_2 * sample_idx as f32).cos();
*sample *= f32x2::splat(0.42 - (0.5 * cos_1) + (0.08 * cos_2));
*sample *= 0.42 - (0.5 * cos_1) + (0.08 * cos_2);
}
// Since this final filter will be symmetrical around `impulse_response[CENTER_IDX]`, we
// can simply normalize based on that fact:
let would_be_impulse_response_sum =
(impulse_response.iter().skip(CENTER_IDX).sum::<f32x2>() * f32x2::splat(2.0))
- impulse_response[CENTER_IDX];
let would_be_impulse_response_sum = (impulse_response.iter().skip(center_idx).sum::<f32>()
* 2.0)
- impulse_response[center_idx];
let would_be_impulse_response_recip = would_be_impulse_response_sum.recip();
for sample in &mut impulse_response {
*sample *= would_be_impulse_response_recip;
@ -196,8 +194,8 @@ impl FirCoefficients {
// And finally we can simply copy the right half of the filter kernel to the left half
// around the `CENTER_IDX`.
for source_idx in CENTER_IDX + 1..impulse_response.len() {
let target_idx = CENTER_IDX - (source_idx - CENTER_IDX);
for source_idx in center_idx + 1..N {
let target_idx = center_idx - (source_idx - center_idx);
impulse_response[target_idx] = impulse_response[source_idx];
}