deviation: Implement peak_signal_to_noise_ratio

src/deviation.rs

@@ -52,6 +52,10 @@ where
     fn root_mean_sq_dev(&self, other: &ArrayBase<S, D>) -> A
     where
         A: AddAssign + Clone + FromPrimitive + Signed + Float;
+
+    fn peak_signal_to_noise_ratio(&self, other: &ArrayBase<S, D>, maxv: A) -> A
+    where
+        A: AddAssign + Clone + FromPrimitive + Signed + Float;
 }
 
 impl<A, S, D> DeviationExt<A, S, D> for ArrayBase<S, D>
@@ -176,6 +180,14 @@ where
     {
         self.mean_sq_dev(other).sqrt()
     }
+
+    fn peak_signal_to_noise_ratio(&self, other: &ArrayBase<S, D>, maxv: A) -> A
+    where
+        A: AddAssign + Clone + FromPrimitive + Signed + Float,
+    {
+        let ten = A::from(10.).unwrap();
+        ten * Float::log10((maxv * maxv) / self.mean_sq_dev(&other))
+    }
 }
 
 #[cfg(test)]
@@ -327,4 +339,19 @@ mod tests {
         assert_abs_diff_eq!(a.root_mean_sq_dev(&b), 10.0.sqrt(), epsilon = f64::EPSILON);
         assert_abs_diff_eq!(b.root_mean_sq_dev(&a), 10.0.sqrt(), epsilon = f64::EPSILON);
     }
+
+    #[test]
+    fn test_peak_signal_to_noise_ratio() {
+        let a = array![1., 1.];
+        assert!(a.peak_signal_to_noise_ratio(&a, 1.).is_infinite());
+
+        let a = array![1., 2., 3., 4., 5., 6., 7.];
+        let b = array![1., 3., 3., 4., 6., 7., 8.];
+        let maxv = 2.;
+
+        let expected = maxv * 10. * Float::log10(maxv) - 10. * Float::log10(a.mean_sq_dev(&b));
+        let actual = a.peak_signal_to_noise_ratio(&b, maxv);
+
+        assert_abs_diff_eq!(actual, expected, epsilon = f64::EPSILON * 8.);
+    }
 }