scripts/dwm1000.py

   1 #!/usr/bin/env python
   2 # -*- coding: iso-8859-15 -*-
   3
   4 VISULAIZE = False
   5
   6 import threading
   7 import struct
   8 import rospy
   9 import tf
  10 import numpy as np
  11 from math import *
  12 from datetime import datetime
  13 from i2c import i2c
  14 from time import sleep
  15 from std_msgs.msg import Float32
  16 from nav_msgs.msg import Odometry
  17 from wild_thumper.srv import DWM1000Center, DWM1000CenterResponse
  18 if VISULAIZE:
  19         import matplotlib.pyplot as plt
  20
  21 class simple_kalman:
  22         def __init__(self, x_est, P_est, Q, R):
  23                 self.x_est = x_est # Systemzustand
  24                 self.P_est = P_est # Fehlerkovarianz
  25                 self.Q = Q # Systemrauschen
  26                 self.R = R # Varianz des Messfehlers
  27
  28         def run(self, y):
  29                 # Korrektur mit der Messung
  30                 # (1) Berechnung der Kalman Verstärkung
  31                 K = self.P_est/(self.R + self.P_est)
  32                 # (2) Korrektur der Schätzung mit der Messung y
  33                 x = self.x_est + K*(y - self.x_est)
  34                 # (3) Korrektur der Fehlerkovarianzmatrix
  35                 P = (1-K)*self.P_est
  36                 #
  37                 # Prädiktion
  38                 # (1) Prädiziere den Systemzustand
  39                 self.x_est = x
  40                 # (2) Präzidiere die Fehlerkovarianzmatrix
  41                 self.P_est = P + self.Q
  42
  43                 return x
  44
  45         def set_measure_cov(self, R):
  46                 if R > 0:
  47                         self.R = R
  48
  49 class DW1000(threading.Thread):
  50         def __init__(self, name, addr, offset):
  51                 threading.Thread.__init__(self)
  52                 self.setDaemon(1)
  53                 self.dist = 0
  54                 self.offset = offset
  55                 self.addr = addr
  56                 self.name = name
  57                 self.last_update = datetime.min
  58
  59                 self.pub = rospy.Publisher(name, Float32, queue_size=16)
  60
  61                 self.start()
  62
  63         def get_value(self):
  64                 dev = i2c(self.addr)
  65                 ret = struct.unpack("f", dev.read(4))
  66                 dev.close()
  67                 return ret[0]
  68
  69         def distance(self):
  70                 return self.dist
  71
  72         # Returns each distance only if current
  73         def distance_valid(self):
  74                 if (datetime.now() - self.last_update).seconds < 1:
  75                         return self.dist
  76                 return None
  77
  78         def run(self):
  79                 last_val = 10
  80                 while not rospy.is_shutdown():
  81                         val = self.get_value()
  82                         if abs(val - last_val)  > 10:
  83                                 rospy.logwarn("Ignoring values too far apart %s: %.2f - %.2f", self.name, val, last_val)
  84                         elif not isnan(val):
  85                                 self.dist = val + self.offset
  86                                 self.last_update = datetime.now()
  87                                 self.pub.publish(self.distance())
  88                                 last_val = val
  89                         sleep(0.1)
  90
  91 class Position:
  92         def __init__(self):
  93                 # Varianz des Messfehlers
  94                 Rx = 0.2
  95                 Ry = 0.05
  96                 # Fehlerkovarianz
  97                 P_est_x = 0.02
  98                 P_est_y = 0.01
  99                 # Systemrauschen
 100                 Q = 0.002
 101                 self.filter_x = simple_kalman(1.0, P_est_x, Q, Rx)
 102                 self.filter_y = simple_kalman(0.0, P_est_y, Q, Ry)
 103                 self.speed_x = 0
 104                 self.speed_y = 0
 105                 self.speed_z = 0
 106                 self.last_time = rospy.Time.now()
 107                 rospy.Subscriber("/odom_combined", Odometry, self.odomReceived)
 108
 109         def odomReceived(self, msg):
 110                 self.speed_x = msg.twist.twist.linear.x
 111                 self.speed_y = msg.twist.twist.linear.y
 112                 self.speed_z = msg.twist.twist.angular.z
 113
 114         def filter(self, x, y):
 115                 # Correct estimation with speed
 116                 current_time = rospy.Time.now()
 117                 dt = (current_time - self.last_time).to_sec()
 118                 # Subtract vehicle speed
 119                 pos = np.array([self.filter_x.x_est, self.filter_y.x_est])
 120                 # translation
 121                 pos -= np.array([self.speed_x*dt, self.speed_y*dt])
 122                 # rotation
 123                 rot = np.array([[np.cos(self.speed_z*dt), -np.sin(self.speed_z*dt)],
 124                                 [np.sin(self.speed_z*dt),  np.cos(self.speed_z*dt)]])
 125                 pos = np.dot(pos, rot)
 126                 # copy back
 127                 self.filter_x.x_est = pos[0]
 128                 self.filter_y.x_est = pos[1]
 129
 130                 # run kalman if new measurements are valid
 131                 if x != None and y != None:
 132                         x = self.filter_x.run(x)
 133                         y = self.filter_y.run(y)
 134
 135                         # Update covariance
 136                         dist = np.linalg.norm([x, y])
 137                         self.filter_x.set_measure_cov(np.polyval([0.017795,  -0.021832, 0.010968], dist))
 138                         self.filter_y.set_measure_cov(np.polyval([0.0060314, -0.013387, 0.0065049], dist))
 139                 else:
 140                         x = self.filter_x.x_est
 141                         y = self.filter_y.x_est
 142
 143                 self.last_time = current_time
 144                 return x,y
 145
 146
 147 def handle_center_call(req):
 148         diff = dwleft.distance_valid() - dwright.distance_valid()
 149         dwleft.offset -= diff/2
 150         dwright.offset += diff/2
 151         rospy.loginfo("Centering to %.2f %.2f", dwleft.offset, dwright.offset)
 152         return DWM1000CenterResponse()
 153
 154 if __name__ == "__main__":
 155         rospy.init_node('DWM1000', log_level=rospy.DEBUG)
 156         dwleft  = DW1000("uwb_dist_left",  0xc2, +0.00)
 157         dwright = DW1000("uwb_dist_right", 0xc0, -0.00)
 158         dist_l_r = 0.285 # Distance between both DWM1000
 159         rate = rospy.Rate(10)
 160         pos = Position()
 161         tf_broadcaster = tf.broadcaster.TransformBroadcaster()
 162         rospy.Service('/DWM1000/center', DWM1000Center, handle_center_call)
 163
 164         while not rospy.is_shutdown() and dwleft.is_alive() and dwright.is_alive():
 165                 dist_left = dwleft.distance_valid()
 166                 dist_right = dwright.distance_valid()
 167                 if dist_left == None or dist_right == None:
 168                         rospy.logerr_throttle(10, "no valid sensor update")
 169                         # run kalman prediction only
 170                         pos.filter(None, None)
 171                 else:
 172                         dir = "left" if (dist_left < dist_right) else "right"
 173
 174                         diff = abs(dist_left - dist_right)
 175                         if diff >= dist_l_r:
 176                                 # difference to high, correct to maximum
 177                                 off = diff - dist_l_r + 0.01
 178                                 if dist_left > dist_right:
 179                                         dist_left -= off/2
 180                                         dist_right += off/2
 181                                 else:
 182                                         dist_left += off/2
 183                                         dist_right -= off/2
 184                         rospy.logdebug("%.2f %.2f %.2f %.2f %s", dwleft.distance(), dwright.distance(), dist_left, dist_right, dir)
 185
 186                         a_r = (-dist_right**2 + dist_left**2 - dist_l_r**2) / (-2*dist_l_r)
 187                         x = dist_l_r/2 - a_r
 188                         t = dist_right**2 - a_r**2
 189                         if t >= 0:
 190                                 y = sqrt(t)
 191                                 rospy.logdebug("x=%.2f, y=%.2f", x, y)
 192                                 # Rotate 90 deg
 193                                 x, y = (y, -x)
 194
 195                                 x, y = pos.filter(x, y)
 196                                 tf_broadcaster.sendTransform((x, y, 0.0), (0, 0, 0, 1), rospy.Time.now(), "uwb_beacon", "base_footprint")
 197
 198                                 if VISULAIZE:
 199                                         circle_left = plt.Circle((-dist_l_r/2, 0), dwleft.distance, color='red', fill=False)
 200                                         circle_right = plt.Circle((dist_l_r/2, 0), dwright.distance, color='green', fill=False)
 201                                         plt.gca().add_patch(circle_left)
 202                                         plt.gca().add_patch(circle_right)
 203                                         plt.grid(True)
 204                                         plt.axis('scaled')
 205                                         plt.show()
 206                         else:
 207                                 # No current position, still need up update kalman prediction
 208                                 pos.filter(None, None)
 209
 210                 rate.sleep()