diff --git a/README.md b/README.md
index 60b7938..17cb821 100644
--- a/README.md
+++ b/README.md
@@ -44,6 +44,14 @@ The following [infinite impulse response](https://en.wikipedia.org/wiki/Infinite
 
 [Numerically controlled oscillator](https://en.wikipedia.org/wiki/Numerically_controlled_oscillator) implemented using a [phasor](https://en.wikipedia.org/wiki/Phasor) and [complex multiplication](https://en.wikipedia.org/wiki/Complex_number#Multiplication) instead of a [lookup table](https://en.wikipedia.org/wiki/Lookup_table).
 
+### [fmd.hh](fmd.hh)
+
+[Frequency modulation](https://en.wikipedia.org/wiki/Frequency_modulation) [demodulation](https://en.wikipedia.org/wiki/Demodulation) with and without [atan2](https://en.wikipedia.org/wiki/Atan2).
+
+### [atan2.hh](atan2.hh)
+
+[atan](https://en.wikipedia.org/wiki/Inverse_trigonometric_functions) and [atan2](https://en.wikipedia.org/wiki/Atan2).
+
 ### [const.hh](const.hh)
 
 Some constants we need
diff --git a/atan2.hh b/atan2.hh
new file mode 100644
index 0000000..a4d1ee1
--- /dev/null
+++ b/atan2.hh
@@ -0,0 +1,74 @@
+/*
+atan and atan2 approximations
+
+Constants below lifted from the Cephes Mathematical Library:
+https://www.netlib.org/cephes/cmath.tgz
+
+Copyright 2019 Ahmet Inan <inan@aicodix.de>
+*/
+
+#pragma once
+
+#include "const.hh"
+
+namespace DSP {
+
+template <typename TYPE>
+TYPE atan_kernel(TYPE x)
+{
+	static const TYPE
+		P0 = -8.750608600031904122785E-1,
+		P1 = -1.615753718733365076637E1,
+		P2 = -7.500855792314704667340E1,
+		P3 = -1.228866684490136173410E2,
+		P4 = -6.485021904942025371773E1,
+		Q0 = 2.485846490142306297962E1,
+		Q1 = 1.650270098316988542046E2,
+		Q2 = 4.328810604912902668951E2,
+		Q3 = 4.853903996359136964868E2,
+		Q4 = 1.945506571482613964425E2;
+	TYPE xx = x * x;
+	return x + x * xx *
+		(xx*(xx*(xx*(xx*P0+P1)+P2)+P3)+P4) /
+		(xx*(xx*(xx*(xx*(xx+Q0)+Q1)+Q2)+Q3)+Q4);
+}
+
+template <typename TYPE>
+TYPE atan_approx(TYPE x)
+{
+	static const TYPE
+		MOREBITS = 6.123233995736765886130E-17;
+	if (x <= TYPE(0.66))
+		return atan_kernel(x);
+	if (x > TYPE(2.41421356237309504880))
+		return Const<TYPE>::HalfPi() + MOREBITS -
+			atan_kernel(TYPE(1) / x);
+	return Const<TYPE>::FourthPi() + TYPE(0.5) * MOREBITS +
+		atan_kernel((x - TYPE(1)) / (x + TYPE(1)));
+}
+
+template <typename TYPE>
+TYPE atan(TYPE x)
+{
+	return x >= TYPE(0) ? atan_approx(x) : -atan_approx(-x);
+}
+
+template <typename TYPE>
+TYPE atan2(TYPE y, TYPE x)
+{
+	if (y == TYPE(0))
+		return x >= TYPE(0) ?
+			copysign(TYPE(0), y) :
+			copysign(Const<TYPE>::Pi(), y);
+	if (x == TYPE(0))
+		return copysign(Const<TYPE>::HalfPi(), y);
+	TYPE phase = atan(y / x);
+	if (x < TYPE(0))
+		return phase > TYPE(0) ?
+			phase - Const<TYPE>::Pi() :
+			phase + Const<TYPE>::Pi();
+	return phase;
+}
+
+}
+
diff --git a/fmd.hh b/fmd.hh
new file mode 100644
index 0000000..900d7b5
--- /dev/null
+++ b/fmd.hh
@@ -0,0 +1,135 @@
+/*
+FM Demodulation
+
+Copyright 2019 Ahmet Inan <inan@aicodix.de>
+*/
+
+#pragma once
+
+namespace DSP {
+
+template <typename TYPE>
+class FMD1
+{
+	typedef TYPE complex_type;
+	typedef typename complex_type::value_type value_type;
+	value_type prev;
+	value_type scale;
+public:
+	constexpr FMD1() : prev(0), scale(0)
+	{
+	}
+	void bandwidth(value_type bw)
+	{
+		scale = value_type(1) / (bw * DSP::Const<value_type>::Pi());
+	}
+	value_type operator()(complex_type input)
+	{
+		value_type phase = arg(input);
+		value_type delta = phase - prev;
+		prev = phase;
+		return scale * delta;
+	}
+};
+
+template <typename TYPE>
+class FMD2
+{
+	typedef TYPE complex_type;
+	typedef typename complex_type::value_type value_type;
+	complex_type prev;
+	value_type scale;
+public:
+	constexpr FMD2() : prev(1, 0), scale(0)
+	{
+	}
+	void bandwidth(value_type bw)
+	{
+		scale = value_type(1) / (bw * DSP::Const<value_type>::Pi());
+	}
+	value_type operator()(complex_type input)
+	{
+		value_type phase = arg(input / prev);
+		prev = input;
+		return scale * phase;
+	}
+};
+
+template <typename TYPE>
+class FMD3
+{
+	typedef TYPE complex_type;
+	typedef typename complex_type::value_type value_type;
+	complex_type prev;
+	value_type scale;
+public:
+	constexpr FMD3() : prev(1, 0), scale(0)
+	{
+	}
+	void bandwidth(value_type bw)
+	{
+		scale = value_type(1) / (bw * DSP::Const<value_type>::Pi());
+	}
+	value_type operator()(complex_type input)
+	{
+		value_type phase = arg(input * conj(prev));
+		prev = input;
+		return scale * phase;
+	}
+};
+
+template <typename TYPE>
+class FMD4
+{
+	typedef TYPE complex_type;
+	typedef typename complex_type::value_type value_type;
+	complex_type prev;
+	value_type scale;
+public:
+	constexpr FMD4() : prev(1, 0), scale(0)
+	{
+	}
+	void bandwidth(value_type bw)
+	{
+		scale = value_type(1) / (bw * DSP::Const<value_type>::Pi());
+	}
+	value_type operator()(complex_type input)
+	{
+		value_type phase =
+			prev.real() * input.imag() -
+			prev.imag() * input.real();
+		phase /= norm(input);
+		prev = input;
+		return scale * phase;
+	}
+};
+
+template <typename TYPE>
+class FMD5
+{
+	typedef TYPE complex_type;
+	typedef typename complex_type::value_type value_type;
+	complex_type prev1, prev2;
+	value_type scale;
+public:
+	constexpr FMD5() : prev1(1, 0), prev2(1, 0), scale(0)
+	{
+	}
+	void bandwidth(value_type bw)
+	{
+		scale = value_type(1) / (bw * DSP::Const<value_type>::TwoPi());
+	}
+	value_type operator()(complex_type input)
+	{
+		complex_type diff = input - prev2;
+		value_type phase =
+			prev1.real() * diff.imag() -
+			prev1.imag() * diff.real();
+		phase /= norm(prev1);
+		prev2 = prev1; prev1 = input;
+		return scale * phase;
+	}
+};
+
+}
+