diff --git a/README.md b/README.md
index 67da554..d0ca5ae 100644
--- a/README.md
+++ b/README.md
@@ -183,6 +183,10 @@ Very useful for data interpolation.
 
 Implemented [Simple linear regression](https://en.wikipedia.org/wiki/Simple_linear_regression) for [Regression analysis](https://en.wikipedia.org/wiki/Regression_analysis) of data.
 
+### [theil_sen.hh](theil_sen.hh)
+
+The [Theil-Sen estimator](https://en.wikipedia.org/wiki/Theil%E2%80%93Sen_estimator) is a [robust](https://en.wikipedia.org/wiki/Robust_statistics) [line fitting](https://en.wikipedia.org/wiki/Line_fitting) algorithm.
+
 ### [complex.hh](complex.hh)
 
 Faster alternative (no Inf/NaN handling) to the std::complex implementation.
diff --git a/theil_sen.hh b/theil_sen.hh
new file mode 100644
index 0000000..c565153
--- /dev/null
+++ b/theil_sen.hh
@@ -0,0 +1,56 @@
+/*
+Theil–Sen estimator
+
+Copyright 2021 Ahmet Inan <inan@aicodix.de>
+*/
+
+#pragma once
+
+#include <algorithm>
+
+namespace DSP {
+
+template <typename TYPE, int LEN_MAX>
+class TheilSenEstimator
+{
+	static const int size_ = ((LEN_MAX-1)*LEN_MAX)/2;
+	TYPE temp_[size_];
+	TYPE xint_, yint_, slope_;
+public:
+	TheilSenEstimator() : xint_(0), yint_(0), slope_(0) {}
+	void compute(TYPE *x, TYPE *y, int LEN)
+	{
+		int count = 0;
+		for (int i = 0; count < size_ && i < LEN; ++i)
+			for (int j = i+1; count < size_ && j < LEN; ++j)
+				if (x[j] != x[i])
+					temp_[count++] = (y[j] - y[i]) / (x[j] - x[i]);
+		std::nth_element(temp_, temp_+count/2, temp_+count);
+		slope_ = temp_[count/2];
+		count = 0;
+		for (int i = 0; count < size_ && i < LEN; ++i)
+			temp_[count++] = y[i] - slope_ * x[i];
+		std::nth_element(temp_, temp_+count/2, temp_+count);
+		yint_ = temp_[count/2];
+		xint_ = - yint_ / slope_;
+	}
+	TYPE xint()
+	{
+		return xint_;
+	}
+	TYPE slope()
+	{
+		return slope_;
+	}
+	TYPE yint()
+	{
+		return yint_;
+	}
+	TYPE operator () (TYPE x)
+	{
+		return yint_ + slope_ * x;
+	}
+};
+
+}
+