Add selectable RTA filter characteristics

src/rta.rs

@@ -1,8 +1,7 @@
 //! Realtime-analyzer filter primitives.
 //!
-//! Fractional-octave bands are designed as complete Butterworth bandpasses:
-//! an analog low-pass prototype is transformed to a bandpass, pre-warped and
-//! mapped with the bilinear transform, then emitted as distinct SOS sections.
+//! Provides both the original repeated-RBJ Phoenix characteristic and complete
+//! Butterworth fractional-octave bandpasses for controlled comparison.
 
 pub const RTW_THIRD_OCTAVE_CENTERS: &[f32] = &[
     20.0, 25.0, 31.5, 40.0, 50.0, 63.0, 80.0, 100.0, 125.0, 160.0, 200.0, 250.0, 315.0, 400.0,
@@ -155,6 +154,35 @@ pub fn design_fractional_octave_band(
     sections
 }
 
+/// Reproduces the original Phoenix analyzer characteristic: one constant-peak
+/// RBJ bandpass is repeated for every second order of the requested cascade.
+/// It is retained for comparison with existing Phoenix/RTW captures and is
+/// not presented as a filter implementation prescribed by IEC 61260.
+pub fn design_legacy_repeated_band(
+    center_hz: f32,
+    sample_rate: u32,
+    bands_per_octave: usize,
+    total_order: usize,
+) -> Vec<BiquadCoeffs> {
+    let fs = f64::from(sample_rate.max(8_000));
+    let center = f64::from(center_hz).clamp(0.01, fs * 0.499_9);
+    let bpo = bands_per_octave.max(1) as f64;
+    let upper = 2.0f64.powf(1.0 / (2.0 * bpo));
+    let lower = 2.0f64.powf(-1.0 / (2.0 * bpo));
+    let q = 1.0 / (upper - lower);
+    let omega = 2.0 * std::f64::consts::PI * center / fs;
+    let alpha = omega.sin() / (2.0 * q);
+    let a0 = 1.0 + alpha;
+    let section = BiquadCoeffs {
+        b0: alpha / a0,
+        b1: 0.0,
+        b2: -alpha / a0,
+        a1: -2.0 * omega.cos() / a0,
+        a2: (1.0 - alpha) / a0,
+    };
+    vec![section; (total_order.clamp(2, 8) / 2).max(1)]
+}
+
 pub fn cascade_magnitude(sections: &[BiquadCoeffs], freq_hz: f32, sample_rate: u32) -> f64 {
     let omega = 2.0 * std::f64::consts::PI * f64::from(freq_hz) / f64::from(sample_rate.max(8_000));
     let z1 = Complex::new(omega.cos(), -omega.sin());
@@ -396,6 +424,27 @@ mod tests {
         assert!(((lower * upper).sqrt() - center).abs() < 0.001);
     }
 
+    #[test]
+    fn legacy_sixth_order_repeats_three_identical_sections() {
+        let sections = design_legacy_repeated_band(1_000.0, 48_000, 6, 6);
+        assert_eq!(sections.len(), 3);
+        assert_eq!(sections[0].b0, sections[1].b0);
+        assert_eq!(sections[1].a2, sections[2].a2);
+        assert!(cascade_db(&sections, 1_000.0, 48_000).abs() < 0.001);
+    }
+
+    #[test]
+    fn legacy_and_butterworth_characteristics_remain_distinct() {
+        let factor = 2.0f32.powf(1.0 / 12.0);
+        let legacy = design_legacy_repeated_band(1_000.0, 48_000, 6, 6);
+        let butterworth =
+            design_fractional_octave_band(1_000.0, 1_000.0 / factor, 1_000.0 * factor, 48_000, 3);
+        let legacy_lower = cascade_db(&legacy, 1_000.0 / factor, 48_000);
+        let butterworth_lower = cascade_db(&butterworth, 1_000.0 / factor, 48_000);
+        assert!(legacy_lower < -8.5);
+        assert!((butterworth_lower + 3.0103).abs() < 0.08);
+    }
+
     #[test]
     fn neighboring_third_octave_centers_are_suppressed() {
         let factor = 2.0f32.powf(1.0 / 6.0);