Add custom operations to Point and simplify algorithm

rir_generator/src/lib.rs

@@ -3,6 +3,7 @@ extern crate itertools;
 
 use std::f64::consts::PI;
 use std::f64::EPSILON;
+use std::ops::Div;
 use thiserror::Error;
 
 #[derive(Debug, Clone)]
@@ -14,15 +15,17 @@ pub enum MicrophoneType {
     Omnidirectional,
 }
 
-impl From<MicrophoneType> for f64 {
-    fn from(val: MicrophoneType) -> f64 {
-        match val {
-            MicrophoneType::Bidirectional => 0.0,
-            MicrophoneType::Hypercardioid => 0.25,
-            MicrophoneType::Cardioid => 0.5,
-            MicrophoneType::Subcardioid => 0.75,
-            MicrophoneType::Omnidirectional => 1.0,
-        }
+impl MicrophoneType {
+    fn adjust_gain(&self, gain: f64) -> f64 {
+        let rho = match self {
+            Self::Bidirectional => 0.0,
+            Self::Hypercardioid => 0.25,
+            Self::Cardioid => 0.5,
+            Self::Subcardioid => 0.75,
+            Self::Omnidirectional => 1.0,
+        };
+
+        rho + (1.0 - rho) * gain
     }
 }
 
@@ -75,6 +78,25 @@ impl From<Vec<f64>> for Position {
     }
 }
 
+impl Position {
+    pub fn radius(&self) -> f64 {
+        (self.x.powi(2) + self.y.powi(2) + self.z.powi(2)).sqrt()
+    }
+}
+
+impl Div<f64> for &Position {
+    // The division of rational numbers is a closed operation.
+    type Output = Position;
+
+    fn div(self, rhs: f64) -> Self::Output {
+        Position {
+            x: self.x / rhs,
+            y: self.y / rhs,
+            z: self.z / rhs,
+        }
+    }
+}
+
 #[derive(Debug, Clone)]
 pub struct Betas {
     pub x: [f64; 2],
@@ -174,17 +196,17 @@ impl FloatSinc for f64 {
     }
 }
 
-fn sim_microphone(p: &Position, a: &Angle, t: &MicrophoneType) -> f64 {
-    match t {
+fn sim_microphone(p: Position, recv: &Receiver) -> f64 {
+    match recv.microphone_type {
         MicrophoneType::Omnidirectional => 1.0,
         _ => {
-            let rho: f64 = t.clone().into();
-            let vartheta = (p.z / (p.x.powi(2) + p.y.powi(2) + p.z.powi(2)).sqrt()).acos();
-            let varphi = p.y.atan2(p.x);
+            let (vartheta, varphi) = ((p.z / p.radius()).acos(), p.y.atan2(p.x));
+            let Angle { theta, phi } = recv.angle;
 
-            let gain = (PI / 2.0 - a.theta).sin() * (vartheta).sin() * (a.phi - varphi).cos()
-                + (PI / 2.0 - a.theta).cos() * (vartheta).cos();
-            rho + (1.0 - rho) * gain
+            let gain = (PI / 2.0 - theta).sin() * vartheta.sin() * (phi - varphi).cos()
+                + (PI / 2.0 - theta).cos() * vartheta.cos();
+
+            recv.microphone_type.adjust_gain(gain)
         }
     }
 }
@@ -201,37 +223,16 @@ pub fn compute_rir(
     n_order: FilterOrder,
     enable_highpass_filter: bool,
 ) -> ndarray::Array2<f64> {
+    const FC: f64 = 0.5;
     // Temporary variables and constants (image-method)
-    let fc = 0.5; // The normalized cut-off frequency equals (fs/2) / fs = 0.5
     let tw = (2.0 * (0.004 * fs).round()) as usize; // The width of the low-pass FIR equals 8 ms
     let cts = c / fs;
     let mut imp = ndarray::Array2::<f64>::zeros((n_samples, receivers.len()));
 
-    let source = Position {
-        x: source.x / cts,
-        y: source.y / cts,
-        z: source.z / cts,
-    };
-    let room = Room {
-        x: room.x / cts,
-        y: room.y / cts,
-        z: room.z / cts,
-    };
-
-    for (
-        Receiver {
-            position: receiver,
-            microphone_type: mtype,
-            angle,
-        },
-        mut imp,
-    ) in receivers.iter().zip(imp.axis_iter_mut(ndarray::Axis(1)))
-    {
-        let receiver = Position {
-            x: receiver.x / cts,
-            y: receiver.y / cts,
-            z: receiver.z / cts,
-        };
+    let (source, room) = (source / cts, room / cts);
+
+    for (receiver, mut imp) in receivers.iter().zip(imp.axis_iter_mut(ndarray::Axis(1))) {
+        let scaled_receiver = &receiver.position / cts;
 
         let n1 = (n_samples as f64 / (2.0 * room.x)).ceil() as i64;
         let n2 = (n_samples as f64 / (2.0 * room.y)).ceil() as i64;
@@ -246,18 +247,18 @@ pub fn compute_rir(
             };
 
             let rp_plus_rm = Position {
-                x: (1 - 2 * q) as f64 * source.x - receiver.x + rm.x,
-                y: (1 - 2 * j) as f64 * source.y - receiver.y + rm.y,
-                z: (1 - 2 * k) as f64 * source.z - receiver.z + rm.z,
+                x: (1 - 2 * q) as f64 * source.x - scaled_receiver.x + rm.x,
+                y: (1 - 2 * j) as f64 * source.y - scaled_receiver.y + rm.y,
+                z: (1 - 2 * k) as f64 * source.z - scaled_receiver.z + rm.z,
             };
             let refl = [
                 beta.x[0].powi((mx - q).abs() as i32) * beta.x[1].powi((mx).abs() as i32),
                 beta.y[0].powi((my - j).abs() as i32) * beta.y[1].powi((my).abs() as i32),
                 beta.z[0].powi((mz - k).abs() as i32) * beta.z[1].powi((mz).abs() as i32),
             ];
 
-            let dist = (rp_plus_rm.x.powi(2) + rp_plus_rm.y.powi(2) + rp_plus_rm.z.powi(2)).sqrt();
-            let fdist = (dist).floor();
+            let dist = rp_plus_rm.radius();
+            let fdist = dist.floor();
 
             let run_block = match n_order {
                 FilterOrder::Any => true,
@@ -267,7 +268,7 @@ pub fn compute_rir(
             };
 
             if run_block && (fdist as usize) < n_samples {
-                let gain = sim_microphone(&rp_plus_rm, angle, mtype) * refl[0] * refl[1] * refl[2]
+                let gain = sim_microphone(rp_plus_rm, receiver) * refl[0] * refl[1] * refl[2]
                     / (4.0 * PI * dist * cts);
 
                 let mut lpi = ndarray::Array::zeros(tw);
@@ -277,8 +278,8 @@ pub fn compute_rir(
                     *lpi = 0.5
                         * (1.0 + (2.0 * PI * t / (tw as f64)).cos())
                         * 2.0
-                        * fc
-                        * (2.0 * PI * fc * t).sinc();
+                        * FC
+                        * (2.0 * PI * FC * t).sinc();
                 }
 
                 let start_position = (fdist - (tw as f64 / 2.0) + 1.0) as usize;