diff --git a/README.md b/README.md
index 6b4a9e2..60b7938 100644
--- a/README.md
+++ b/README.md
@@ -40,6 +40,10 @@ The following [infinite impulse response](https://en.wikipedia.org/wiki/Infinite
 * [second-order Butterworth low pass filter](https://en.wikipedia.org/wiki/Butterworth_filter)
 * [2n-order Butterworth cascade of second-order low pass filters](https://en.wikipedia.org/wiki/Butterworth_filter)
 
+### [phasor.hh](phasor.hh)
+
+[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).
+
 ### [const.hh](const.hh)
 
 Some constants we need
diff --git a/phasor.hh b/phasor.hh
new file mode 100644
index 0000000..7f14c70
--- /dev/null
+++ b/phasor.hh
@@ -0,0 +1,38 @@
+/*
+Numerically controlled oscillator
+
+Copyright 2019 Ahmet Inan <inan@aicodix.de>
+*/
+
+#include "unit_circle.hh"
+
+#pragma once
+
+namespace DSP {
+
+template <typename TYPE>
+class Phasor
+{
+	typedef TYPE complex_type;
+	typedef typename complex_type::value_type value_type;
+	complex_type prev, delta;
+public:
+	constexpr Phasor() : prev(1, 0), delta(1, 0)
+	{
+	}
+	void omega(int n, int N)
+	{
+		delta = complex_type(
+			UnitCircle<value_type>::cos(n, N),
+			UnitCircle<value_type>::sin(n, N));
+	}
+	complex_type operator()()
+	{
+		complex_type tmp = prev;
+		prev *= delta;
+		prev /= abs(prev);
+		return tmp;
+	}
+};
+
+}