surge-biquad 0.2.12-alpha.0

/*!
  | For testing, here are some suggested tests:
  |
  | 1. Test that the coefficients are set
  | correctly for a known cutoff frequency and
  | quality factor. This can be done by comparing
  | the coefficients generated by the filter to
  | the expected coefficients for a known input
  | signal.
  |
  | 2. Test the frequency response of the filter
  | by passing in sine waves at various
  | frequencies and comparing the output to the
  | expected values. This can be done using an FFT
  | to compute the magnitude and phase of the
  | output signal, and comparing to the expected
  | frequency response of a low-pass filter.
  |
  */

crate::ix!();

/// ### Testing
/// 
/// To test the `BiquadCoeffLP` implementation, we
/// can write unit tests that check whether the
/// coefficients are calculated correctly for
/// various input parameters. We can also write
/// integration tests that check whether the
/// filter behaves as expected when processing
/// input signals.
/// 
/// Here's an example of a unit test that checks
/// whether the coefficients are calculated
/// correctly for a given cutoff frequency and
/// Q factor:
/// 
#[test]
fn test_coeff_lp() {

    let mut filter = BiquadFilter::default();

    filter.coeff_lp(0.1, 0.5);

    let expected = (
        1.1800682076597505,
        -1.690035837119334,
        0.7436025846804989,
        0.12968589617062476,
        0.25937179234124953,
        0.12968589617062476,
    );

    let actual = filter.get_coef();

    assert_approx_eq!(expected.0, actual.0, 1e-12);
    assert_approx_eq!(expected.1, actual.1, 1e-12);
    assert_approx_eq!(expected.2, actual.2, 1e-12);
    assert_approx_eq!(expected.3, actual.3, 1e-12);
    assert_approx_eq!(expected.4, actual.4, 1e-12);
    assert_approx_eq!(expected.5, actual.5, 1e-12);
}

/// This test generates a range of sine wave
/// signals at different frequencies, filters each
/// signal with the low-pass filter, and
/// calculates the frequency response of the
/// filtered signal using a Fourier transform. The
/// expected frequency response is calculated
/// based on the filter's cutoff frequency and
/// Q value, and the measured frequency response
/// is compared to the expected response at each
/// frequency. The test fails if the relative
/// error between the measured and expected
///
 #[test]
fn test_lp_filter_frequency_response() {

    // Create a filter with a cutoff frequency of 500 Hz and Q of 0.7
    let mut filter = BiquadFilter::new();
    filter.coeff_lp_with_bw(2.0 * PI * 500.0, 0.7);

    // Set up test parameters
    let num_samples = 10000;
    let sample_rate = 44100.0;
    let duration = num_samples as f64 / sample_rate;
    let f_start = 20.0;
    let f_stop = 20000.0;
    let num_freqs = 100;
    let log_spacing = true;

    // Generate test signals at a range of frequencies
    let mut signals = Vec::new();
    let freqs = if log_spacing {
        logspace(f_start, f_stop, num_freqs)
    } else {
        linspace(f_start, f_stop, num_freqs)
    };
    for f in freqs {
        let mut signal = Vec::new();
        for n in 0..num_samples {
            let t = n as f64 / sample_rate;
            signal.push(0.5 * (2.0 * PI * f * t).sin());
        }
        signals.push(signal);
    }

    // Filter the test signals and calculate the frequency response
    let mut outputs = Vec::new();

    let mut freq_responses = Vec::new();

    for signal in signals {
        let output = filter.filter(&signal);
        outputs.push(output.clone());
        let freq_response = fourier_transform::magnitude_spectrum(&output, sample_rate);
        freq_responses.push(freq_response);
    }

    // Calculate expected frequency response
    let expected_response = if log_spacing {

        logspace(f_start, f_stop, num_freqs)
            .iter()
            .map(|&f| if f <= 500.0 {
                1.0
            } else {
                1.0 / (1.0 + ((f / 500.0) / 0.7).powi(4))
            })
        .collect()

    } else {

        linspace(f_start, f_stop, num_freqs)
            .iter()
            .map(|&f| if f <= 500.0 {
                1.0
            } else {
                1.0 / (1.0 + ((f / 500.0) / 0.7).powi(4))
            })
        .collect()
    };

    // Compare output to expected frequency response
    for i in 0..num_freqs {

        let freq          = freqs[i];
        let freq_response = &freq_responses[i];
        let expected      = expected_response[i];
        let measured      = freq_response[0]; // assume signal is pure sine wave, take magnitude of first frequency bin
        let error         = (measured - expected).abs() / expected;

        assert!(error < 0.1, "Error too large at frequency {:.1} Hz: expected {:.4}, measured {:.4}", freq, expected, measured);
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_coeff_lp() {

        // Test case 1: omega > PI
        let mut filter = BiquadFilter::new();

        filter.coeff_lp(3.5, 0.5);

        let expected = [1.0, 0.0, 0.0, 1.0, 0.0, 0.0];

        assert_eq!(filter.get_coef(), expected);

        // Test case 2: omega < PI
        let mut filter = BiquadFilter::new();

        filter.coeff_lp(1.0, 0.7);

        let expected = [
            1.4289375376768227,
            -0.5486291845786517,
            0.4161249249193546,
            0.7915624891192163,
            1.5831249782384326,
            0.7915624891192163,
        ];

        assert_approx_eq!(filter.get_coef(), expected);
    }

    #[test]
    fn test_coeff_lp_with_bw() {

        // Test case 1: bandwidth = 0.5
        let mut filter = BiquadFilter::new();

        filter.coeff_lp_with_bw(1.0, 0.5);

        let expected = [
            1.2405145438263382,
            -0.28010248420823334,
            0.24051454382633833,
            0.6202572719131692,
            1.2405145438263382,
            0.6202572719131692,
        ];

        assert_approx_eq!(filter.get_coef(), expected);

        // Test case 2: bandwidth = 1.0
        let mut filter = BiquadFilter::new();

        filter.coeff_lp_with_bw(2.5, 1.0);

        let expected = [
            1.1037195231547624,
            -0.42750521737628485,
            0.1869475083210603,
            0.09347375416053014,
            0.1869475083210603,
            0.09347375416053014,
        ];

        assert_approx_eq!(filter.get_coef(), expected);
    }
}

/// This tests the filter with a sine wave input
/// signal, generating an expected output signal
/// based on the filter coefficients, and then
/// comparing the actual output signal from the
/// filter to the expected output signal.
///
#[test]
fn test_filter_sine_wave() {

    // Set up filter coefficients
    let mut filter = BiquadFilter::new();

    let f0          = 1000.0;
    let q           = 0.707;
    let sample_rate = 44100.0;
    let omega       = 2.0 * PI * f0 / sample_rate;

    filter.coeff_lp(omega, q);

    // Generate input signal
    let duration     = 1.0;
    let amplitude    = 0.5;
    let sample_count = (duration * sample_rate) as usize;

    let mut input_signal = vec![0.0; sample_count];

    let t = (0..sample_count).map(|i| i as f64 / sample_rate);

    let sine_wave = t.map(|t| amplitude * (2.0 * PI * f0 * t).sin());

    input_signal.copy_from_slice(&sine_wave.collect::<Vec<f64>>());

    // Filter input signal
    let output_signal = filter.process(&input_signal);

    // Test output signal against expected output
    let mut expected_output_signal = vec![0.0; sample_count];

    let b0 = 0.5 * (1.0 - (2.0 * f0 / sample_rate * PI).cos());
    let b1 = 1.0 - (2.0 * f0 / sample_rate * PI).cos();
    let b2 = 0.5 * (1.0 - (2.0 * f0 / sample_rate * PI).cos());
    let a0 = 1.0 + (q * omega).sin();
    let a1 = -2.0 * (2.0 * f0 / sample_rate * PI).cos();
    let a2 = 1.0 - (q * omega).sin();

    for i in 0..sample_count {

        if i == 0 {
            expected_output_signal[i] = b0 * input_signal[i] + b1 * 0.0 + b2 * 0.0;

        } else if i == 1 {

            expected_output_signal[i] = b0 * input_signal[i] + b1 * input_signal[i-1] + b2 * 0.0
                - a1 * expected_output_signal[i-1];

        } else {

            expected_output_signal[i] = b0 * input_signal[i] + b1 * input_signal[i-1] + b2 * input_signal[i-2]
                - a1 * expected_output_signal[i-1] - a2 * expected_output_signal[i-2];
        }
    }

    assert!(output_signal.iter().zip(expected_output_signal.iter())
        .all(|(a, b)| (a-b).abs() < 1e-12));
}