From 749668d830da1f34cba6e38a9be01958ed651080 Mon Sep 17 00:00:00 2001
From: Ahmet Inan <inan@aicodix.de>
Date: Sun, 7 Oct 2018 19:00:57 +0200
Subject: [PATCH] added probability density function of the normal distribution

---
 README.md | 5 ++++-
 const.hh  | 2 ++
 utils.hh  | 8 ++++++++
 3 files changed, 14 insertions(+), 1 deletion(-)

diff --git a/README.md b/README.md
index c3ad47c..cb4627f 100644
--- a/README.md
+++ b/README.md
@@ -52,7 +52,10 @@ Mixed-radix [decimation-in-time](https://en.wikipedia.org/wiki/Cooley%E2%80%93Tu
 
 ### [utils.hh](utils.hh)
 
-Some everyday helpers, like the [signum function](https://en.wikipedia.org/wiki/Sign_function) or the [lerp function](https://en.wikipedia.org/wiki/Linear_interpolation).
+Some everyday helpers:
+* [signum function](https://en.wikipedia.org/wiki/Sign_function)
+* [lerp function](https://en.wikipedia.org/wiki/Linear_interpolation)
+* [probability density function](https://en.wikipedia.org/wiki/Probability_density_function) of the [normal distribution](https://en.wikipedia.org/wiki/Normal_distribution)
 
 ### [resampler.hh](resampler.hh)
 
diff --git a/const.hh b/const.hh
index 953806d..ab4acdc 100644
--- a/const.hh
+++ b/const.hh
@@ -17,6 +17,8 @@ struct Const {
 	static constexpr T Pi()		{ return   3.14159265358979323846; } // 4*a(1)
 	static constexpr T TwoPi()	{ return   6.28318530717958647693; } // 8*a(1)
 	static constexpr T FourPi()	{ return  12.56637061435917295385; } // 16*a(1)
+	static constexpr T SqrtPi()	{ return   1.77245385090551602730; } // sqrt(4*a(1))
+	static constexpr T SqrtTwoPi()	{ return   2.50662827463100050242; } // sqrt(8*a(1))
 };
 
 }
diff --git a/utils.hh b/utils.hh
index b23d265..ff987b3 100644
--- a/utils.hh
+++ b/utils.hh
@@ -7,6 +7,8 @@ Copyright 2018 Ahmet Inan <inan@aicodix.de>
 #ifndef UTILS_HH
 #define UTILS_HH
 
+#include "const.hh"
+
 namespace DSP {
 
 template <typename TYPE>
@@ -21,6 +23,12 @@ AB lerp(X x, AB a, AB b)
 	return (X(1) - x) * a + x * b;
 }
 
+template <typename TYPE>
+TYPE normal_pdf(TYPE x, TYPE m, TYPE s)
+{
+	return std::exp(-std::pow((x - m) / s, TYPE(2)) / TYPE(2)) / (Const<TYPE>::SqrtTwoPi() * s);
+}
+
 }
 
 #endif