code: def integrate_gyro_quaternion(gyro_ts, gyro_data):
      
      #NB: Quaternion q = [a, n1, n2, n3], scalar first
      q_list = np.zeros((gyro_ts.shape[0], 4)) # Nx4 quaternion list
      q_list[0,:] = np.array([1, 0, 0, 0]) # Initial rotation (no rotation)
      
      # Iterate over all (except first)
      for i in range(1, gyro_ts.size):
          w = gyro_data[i]
          dt = gyro_ts[i] - gyro_ts[i - 1]
          qprev = q_list[i - 1]
          
          A = np.array([[0,    -w[0],  -w[1],  -w[2]],
                       [w[0],  0,      w[2],  -w[1]],
                       [w[1], -w[2],   0,      w[0]],
                       [w[2],  w[1],  -w[0],   0]])
          qnew = (np.eye(4) + (dt/2.0) * A).dot(qprev)
          qnorm = np.sqrt(np.sum(qnew ** 2))
          qnew /= qnorm
          q_list[i] = qnew
           
      return q_list
docs: Integrate angular velocities to rotations
      
      Parameters
      ---------------
      gyro_ts : ndarray
              Timestamps
      gyro_data : (3, N) ndarray
              Angular velocity measurements
      
      Returns
      ---------------
      rotations : (4, N) ndarray
              Rotation sequence as unit quaternions (first element scalar)