diff --git a/ohmpi/hardware_system.py b/ohmpi/hardware_system.py
index 488d0cbd4b03e1775d70e36c79741d49f0cb7667..acb90f83ce353b3b45ad689573cfa5fade73181f 100644
--- a/ohmpi/hardware_system.py
+++ b/ohmpi/hardware_system.py
@@ -299,6 +299,34 @@ class OhmPiHardware:
         else:
             return np.nan
 
+    def last_vmn(self, delay=0.):
+        v = np.where((self.readings[:, 0] >= delay) & (self.readings[:, 2] != 0))[0]
+        if len(v) > 1:
+            return np.mean(self.readings[v, 2] * (self.readings[v, 4] - self.sp))
+        else:
+            return np.nan
+
+    def last_vmn_dev(self, delay=0.):  # TODO: should compute std per stack because this does not account for SP...
+        v = np.where((self.readings[:, 0] >= delay) & (self.readings[:, 2] != 0))[0]
+        if len(v) > 1:
+            return 100. * np.std(self.readings[v, 2] * (self.readings[v, 4])) / self.last_vmn(delay=delay)
+        else:
+            return np.nan
+
+    def last_iab(self, delay=0.):
+        v = np.where((self.readings[:, 0] >= delay) & (self.readings[:, 2] != 0))[0]
+        if len(v) > 1:
+            return np.mean(self.readings[v, 3])
+        else:
+            return np.nan
+
+    def last_iab_dev(self, delay=0.):
+        v = np.where((self.readings[:, 0] >= delay) & (self.readings[:, 2] != 0))[0]
+        if len(v) > 1:
+            return 100. * np.std(self.readings[v, 3]) / / self.last_iab(delay=delay)
+        else:
+            return np.nan
+
     @property
     def sp(self):  # TODO: allow for different strategies for computing sp (i.e. when sp drift is not linear)
         if self.readings.shape == (0,) or len(self.readings[self.readings[:, 2] == 1, :]) < 1 or \
@@ -362,6 +390,7 @@ class OhmPiHardware:
                          iab_max=None, vmn_max=None, vmn_min=voltage_min, polarities=(1, -1), delay=0.05,
                          p_max=None, diff_vab_lim=2.5, n_steps=4):
         # TODO: Optimise how to pass iab_max, vab_max, vmn_min
+        # TODO: Update docstring
         """Estimates best Tx voltage based on different strategies.
         At first a half-cycle is made for a short duration with a fixed
         known voltage. This gives us Iab and Rab. We also measure Vmn.
diff --git a/ohmpi/ohmpi.py b/ohmpi/ohmpi.py
index 1fdbbc822bc9b265c121f7115a53b56f4f07677f..f10a5f07d2922d3a3299f2a0dadb9697a6b94a19 100644
--- a/ohmpi/ohmpi.py
+++ b/ohmpi/ohmpi.py
@@ -536,10 +536,10 @@ class OhmPi(object):
             else:
                 delay = injection_duration * 2/3  # TODO: check if this is ok and if last point is not taken the end of injection
             x = np.where((self._hw.readings[:, 0] >= delay) & (self._hw.readings[:, 2] != 0))
-            Vmn = np.mean(self._hw.readings[x, 2] * (self._hw.readings[x, 4] - self._hw.sp))
-            Vmn_std = 100. * np.std(self._hw.readings[x, 2] * (self._hw.readings[x, 4])) # - self._hw.sp))
-            I = np.mean(self._hw.readings[x, 3])
-            I_std = 100. * np.std(self._hw.readings[x, 3])
+            Vmn = self._hw.last_vmn(delay=delay)
+            Vmn_std = self._hw.last_vmn_dev(delay=delay)
+            I =  self._hw.last_iab(delay=delay)
+            I_std =  self._hw.last_iab_dev(delay=delay)
             R = self._hw.last_resistance(delay=delay)
             R_std = self._hw.last_dev(delay=delay)
             d = {