diff --git a/pyerrors/obs.py b/pyerrors/obs.py
index 3cb22f36..461e76a8 100644
--- a/pyerrors/obs.py
+++ b/pyerrors/obs.py
@@ -1332,7 +1332,7 @@ def correlate(obs_a, obs_b):
     return o
 
 
-def covariance(obs, visualize=False, **kwargs):
+def covariance(obs, visualize=False, correlation=False, **kwargs):
     """Calculates the covariance matrix of a set of observables.
 
     covariance([obs, obs])[0,1] is equal to obs.dvalue ** 2
@@ -1343,7 +1343,9 @@ def covariance(obs, visualize=False, **kwargs):
     obs : list or numpy.ndarray
         List or one dimensional array of Obs
     visualize : bool
-        Plots the corresponding normalized correlation matrix.
+        If True plots the corresponding normalized correlation matrix (default False).
+    correlation : bool
+        If True the correlation instead of the covariance is returned (default False).
     """
 
     length = len(obs)
@@ -1368,7 +1370,10 @@ def covariance(obs, visualize=False, **kwargs):
         plt.colorbar()
         plt.draw()
 
-    return cov
+    if correlation is True:
+        return corr
+    else:
+        return cov
 
 
 def _covariance_element(obs1, obs2):
diff --git a/tests/obs_test.py b/tests/obs_test.py
index c86da060..2a80cac1 100644
--- a/tests/obs_test.py
+++ b/tests/obs_test.py
@@ -246,17 +246,6 @@ def test_gamma_method_kwargs():
     assert my_obs.N_sigma['ens'] == pe.Obs.N_sigma_global
 
 
-def test_covariance_is_variance():
-    value = np.random.normal(5, 10)
-    dvalue = np.abs(np.random.normal(0, 1))
-    test_obs = pe.pseudo_Obs(value, dvalue, 't')
-    test_obs.gamma_method()
-    assert np.abs(test_obs.dvalue ** 2 - pe.covariance([test_obs, test_obs])[0, 1]) <= 10 * np.finfo(np.float64).eps
-    test_obs = test_obs + pe.pseudo_Obs(value, dvalue, 'q', 200)
-    test_obs.gamma_method()
-    assert np.abs(test_obs.dvalue ** 2 - pe.covariance([test_obs, test_obs])[0, 1]) <= 10 * np.finfo(np.float64).eps
-
-
 def test_fft():
     value = np.random.normal(5, 100)
     dvalue = np.abs(np.random.normal(0, 5))
@@ -605,6 +594,18 @@ def test_gamma_method_irregular():
     ol = [a, b]
     o = (ol[0] - ol[1]) / (ol[1])
 
+
+def test_covariance_is_variance():
+    value = np.random.normal(5, 10)
+    dvalue = np.abs(np.random.normal(0, 1))
+    test_obs = pe.pseudo_Obs(value, dvalue, 't')
+    test_obs.gamma_method()
+    assert np.isclose(test_obs.dvalue ** 2, pe.covariance([test_obs, test_obs])[0, 1])
+    test_obs = test_obs + pe.pseudo_Obs(value, dvalue, 'q', 200)
+    test_obs.gamma_method()
+    assert np.isclose(test_obs.dvalue ** 2, pe.covariance([test_obs, test_obs])[0, 1])
+
+
 def test_covariance_symmetry():
     value1 = np.random.normal(5, 10)
     dvalue1 = np.abs(np.random.normal(0, 1))
@@ -616,7 +617,7 @@ def test_covariance_symmetry():
     test_obs2.gamma_method()
     cov_ab = pe.covariance([test_obs1, test_obs2])[0, 1]
     cov_ba = pe.covariance([test_obs2, test_obs1])[0, 1]
-    assert np.abs(cov_ab - cov_ba) <= 10 * np.finfo(np.float64).eps
+    assert np.isclose(cov_ab, cov_ba)
     assert np.abs(cov_ab) < test_obs1.dvalue * test_obs2.dvalue * (1 + 10 * np.finfo(np.float64).eps)
 
     N = 100
@@ -673,6 +674,11 @@ def test_covariance_factorizing():
     assert np.isclose(pe.covariance([mt0, tt[1]])[0, 1], -pe.covariance(tt)[0, 1])
 
 
+def test_covariance_correlation():
+    test_obs = pe.pseudo_Obs(-4, 0.8, 'test', samples=784)
+    assert np.allclose(pe.covariance([test_obs, test_obs, test_obs], correlation=True), np.ones((3, 3)))
+
+
 def test_empty_obs():
     o = pe.Obs([np.random.rand(100)], ['test'])
     q = o + pe.Obs([], [])