diff --git a/GPy/kern/_src/rbf.py b/GPy/kern/_src/rbf.py
index ce769879..01cfdc16 100644
--- a/GPy/kern/_src/rbf.py
+++ b/GPy/kern/_src/rbf.py
@@ -37,6 +37,10 @@ class RBF(Stationary):
     def dK_dr(self, r):
         return -r*self.K_of_r(r)
 
+    def spectrum(self, omega):
+        assert self.input_dim == 1 #TODO: higher dim spectra?
+        return self.variance*np.sqrt(2*np.pi)*self.lengthscale*np.exp(-self.lengthscale*2*omega**2/2)
+
     #---------------------------------------#
     #             PSI statistics            #
     #---------------------------------------#
diff --git a/GPy/kern/_src/stationary.py b/GPy/kern/_src/stationary.py
index ef8900f9..5500bb88 100644
--- a/GPy/kern/_src/stationary.py
+++ b/GPy/kern/_src/stationary.py
@@ -90,13 +90,14 @@ class Stationary(Kern):
             Xsq = np.sum(np.square(X),1)
             r2 = -2.*tdot(X) + (Xsq[:,None] + Xsq[None,:])
             util.diag.view(r2)[:,]= 0. # force diagnoal to be zero: sometime numerically a little negative
+            r2 = np.clip(r2, 0, np.inf)
             return np.sqrt(r2)
         else:
             #X2, = self._slice_X(X2)
             X1sq = np.sum(np.square(X),1)
             X2sq = np.sum(np.square(X2),1)
             r2 = -2.*np.dot(X, X2.T) + X1sq[:,None] + X2sq[None,:]
-            r2[r2<0] = 0. # A bit hacky
+            r2 = np.clip(r2, 0, np.inf)
             return np.sqrt(r2)
 
     @Cache_this(limit=5, ignore_args=())