diff --git a/src/processor_imu.cpp b/src/processor_imu.cpp
index 0231fa66c13356b290b687610c5641c852b8307f..0e56a1a4b28e9f16f0779a922e9977bbaaf5b0fd 100644
--- a/src/processor_imu.cpp
+++ b/src/processor_imu.cpp
@@ -108,6 +108,129 @@ bool ProcessorIMU::voteForKeyFrame()
     return false;
 }
 
+Motion ProcessorIMU::interpolate(const Motion& _motion_ref, Motion& _motion_second, TimeStamp& _ts)
+{
+    /* Note: See extensive documentation in ProcessorMotion::interpolate().
+     *
+     * Interpolate between motion_ref and motion, as in:
+     *
+     *    motion_ref ------ ts_ ------ motion_sec
+     *                    return
+     *
+     * and return the motion at the given time_stamp ts_.
+     *
+     * DATA:
+     * Data receives linear interpolation
+     *    a_ret = (ts - t_ref) / dt * a_sec
+     *    w_ret = (ts - t_ref) / dt * w_sec
+     *
+     * DELTA:
+     * The delta's position and velocity receive linear interpolation:
+     *    p_ret = (ts - t_ref) / dt * (p - p_ref)
+     *    v_ret = (ts - t_ref) / dt * (v - v_ref)
+     *
+     * The delta's quaternion receives a slerp interpolation
+     *    q_ret = q_ref.slerp( (ts - t_ref) / dt , q);
+     *
+     * DELTA_INTEGR:
+     * The interpolated delta integral is just the reference integral plus the interpolated delta
+     *
+     * The second integral does not change
+     *
+     * Covariances receive linear interpolation
+     *    Q_ret = (ts - t_ref) / dt * Q_sec
+     */
+    // resolve out-of-bounds time stamp as if the time stamp was exactly on the bounds
+    if (_ts <= _motion_ref.ts_)    // behave as if _ts == _motion_ref.ts_
+    {
+        // return null motion. Second stays the same.
+        Motion motion_int     ( _motion_ref );
+        motion_int.data_      . setZero();
+        motion_int.data_cov_  . setZero();
+        motion_int.delta_     = deltaZero();
+        motion_int.delta_cov_ . setZero();
+        return motion_int;
+    }
+    if (_motion_second.ts_ <= _ts)    // behave as if _ts == _motion_second.ts_
+    {
+        // return original second motion. Second motion becomes null motion
+        Motion motion_int         ( _motion_second );
+        _motion_second.data_      . setZero();
+        _motion_second.data_cov_  . setZero();
+        _motion_second.delta_     = deltaZero();
+        _motion_second.delta_cov_ . setZero();
+        return motion_int;
+    }
+    assert(_motion_ref.ts_ <= _ts && "Interpolation time stamp out of bounds");
+    assert(_ts <= _motion_second.ts_ && "Interpolation time stamp out of bounds");
+
+    assert(_motion_ref.delta_.size() == delta_size_ && "Wrong delta size");
+    assert(_motion_ref.delta_cov_.cols() == delta_cov_size_ && "Wrong delta cov size");
+    assert(_motion_ref.delta_cov_.rows() == delta_cov_size_ && "Wrong delta cov size");
+    assert(_motion_ref.delta_integr_.size() == delta_size_ && "Wrong delta size");
+    assert(_motion_ref.delta_integr_cov_.cols() == delta_cov_size_ && "Wrong delta cov size");
+    assert(_motion_ref.delta_integr_cov_.rows() == delta_cov_size_ && "Wrong delta cov size");
+
+    assert(_motion_second.delta_.size() == delta_size_ && "Wrong delta size");
+    assert(_motion_second.delta_cov_.cols() == delta_cov_size_ && "Wrong delta cov size");
+    assert(_motion_second.delta_cov_.rows() == delta_cov_size_ && "Wrong delta cov size");
+    assert(_motion_second.delta_integr_.size() == delta_size_ && "Wrong delta size");
+    assert(_motion_second.delta_integr_cov_.cols() == delta_cov_size_ && "Wrong delta cov size");
+    assert(_motion_second.delta_integr_cov_.rows() == delta_cov_size_ && "Wrong delta cov size");
+
+    // reference
+    TimeStamp t_ref = _motion_ref.ts_;
+
+    // second
+    TimeStamp t_sec = _motion_second.ts_;
+    Map<VectorXs>    motion_sec_dp  (_motion_second.delta_.data() + 0, 3);
+    Map<Quaternions> motion_sec_dq  (_motion_second.delta_.data() + 3   );
+    Map<VectorXs>    motion_sec_dv  (_motion_second.delta_.data() + 7, 3);
+
+    // interpolated
+    Motion motion_int = motionZero(_ts);
+
+    // Jacobians for covariance propagation
+    MatrixXs J_ref(delta_cov_size_, delta_cov_size_);
+    MatrixXs J_int(delta_cov_size_, delta_cov_size_);
+
+    // interpolation factor
+    Scalar dt  = _ts - t_ref;
+    Scalar tau = dt / (t_sec - t_ref); // interpolation factor (0 to 1)
+
+    // interpolate data
+    motion_int.data_     = tau * _motion_second.data_;
+    motion_int.data_cov_ = tau * _motion_second.data_cov_;
+
+    // interpolate delta
+    motion_int.ts_        = _ts;
+    Map<VectorXs>    motion_int_dp  (motion_int.delta_.data() + 0, 3);
+    Map<Quaternions> motion_int_dq  (motion_int.delta_.data() + 3   );
+    Map<VectorXs>    motion_int_dv  (motion_int.delta_.data() + 7, 3);
+    motion_int_dp         = tau * motion_sec_dp; // FIXME: delta_p not correctly interpolated
+    motion_int_dq         = Quaternions::Identity().slerp(tau, motion_sec_dq);
+    motion_int_dv         = tau * motion_sec_dv;
+    motion_int.delta_cov_ = tau * _motion_second.delta_cov_;
+
+    // integrate
+    deltaPlusDelta(_motion_ref.delta_integr_, motion_int.delta_, dt, motion_int.delta_integr_, J_ref, J_int);
+    motion_int.delta_integr_cov_ = J_ref * _motion_ref.delta_integr_cov_ * J_ref.transpose()
+                                 + J_int * _motion_second.delta_cov_     * J_int.transpose();
+
+    // update second delta ( in place update )
+    _motion_second.data_      *= (1 - tau);
+    _motion_second.data_cov_  *= (1 - tau);
+    motion_sec_dp = motion_int_dq.conjugate() * ((1 - tau) * motion_sec_dp);
+    motion_sec_dq = motion_int_dq.conjugate() * motion_sec_dq;
+    motion_sec_dv = motion_int_dq.conjugate() * ((1 - tau) * motion_sec_dv);
+    _motion_second.delta_cov_ *= (1 - tau); // easy interpolation // TODO check for correctness
+    //Dp            = Dp; // trivial, just leave the code commented
+    //Dq            = Dq; // trivial, just leave the code commented
+    //_motion.delta_integr_cov_ = _motion.delta_integr_cov_; // trivial, just leave the code commented
+
+    return motion_int;
+}
+
 } // namespace wolf
 
 
diff --git a/src/processor_imu.h b/src/processor_imu.h
index be263aaa86c3afd8c25ba0f62fb7832705c40d9f..55109586c593c607c7d8fcafc699f58f05a5306d 100644
--- a/src/processor_imu.h
+++ b/src/processor_imu.h
@@ -406,102 +406,6 @@ inline Eigen::VectorXs ProcessorIMU::deltaZero() const
     return (Eigen::VectorXs(10) << 0,0,0,  0,0,0,1,  0,0,0 ).finished(); // p, q, v
 }
 
-inline Motion ProcessorIMU::interpolate(const Motion& _motion_ref,
-                                        Motion& _motion_second,
-                                        TimeStamp& _ts)
-{
-    // resolve out-of-bounds time stamp as if the time stamp was exactly on the bounds
-    if (_ts <= _motion_ref.ts_)     // behave as if _ts == _motion_ref.ts_
-    { // return null motion. Second stays the same.
-        Motion motion_int(_motion_ref);
-        motion_int.delta_ = deltaZero();
-        motion_int.delta_cov_.setZero();
-        return motion_int;
-    }
-    if (_motion_second.ts_ <= _ts)  // behave as if _ts == _motion_second.ts_
-    { // return original second motion. Second motion becomes null motion
-        Motion motion_int(_motion_second);
-        _motion_second.delta_ = deltaZero();
-        _motion_second.delta_cov_.setZero();
-        return motion_int;
-    }
-
-    assert(_motion_ref.ts_ <= _ts && "Interpolation time stamp out of bounds");
-    assert(_ts <= _motion_second.ts_ && "Interpolation time stamp out of bounds");
-
-    assert(_motion_ref.delta_.size() == delta_size_ && "Wrong delta size");
-    assert(_motion_ref.delta_cov_.cols() == delta_cov_size_ && "Wrong delta cov size");
-    assert(_motion_ref.delta_cov_.rows() == delta_cov_size_ && "Wrong delta cov size");
-    assert(_motion_ref.delta_integr_.size() == delta_size_ && "Wrong delta size");
-    //assert(_motion_ref.delta_integr_cov_.cols() == delta_cov_size_ && "Wrong delta cov size");
-    //assert(_motion_ref.delta_integr_cov_.rows() == delta_cov_size_ && "Wrong delta cov size");
-    assert(_motion_second.delta_.size() == delta_size_ && "Wrong delta size");
-    assert(_motion_second.delta_cov_.cols() == delta_cov_size_ && "Wrong delta cov size");
-    assert(_motion_second.delta_cov_.rows() == delta_cov_size_ && "Wrong delta cov size");
-    assert(_motion_second.delta_integr_.size() == delta_size_ && "Wrong delta size");
-    //assert(_motion_second.delta_integr_cov_.cols() == delta_cov_size_ && "Wrong delta cov size");
-    //assert(_motion_second.delta_integr_cov_.rows() == delta_cov_size_ && "Wrong delta cov size");
-
-    // Interpolate between motion_ref and motion, as in:
-    //
-    // motion_ref ------ ts_ ------ motion
-    //                 return
-    //
-    // and return the value at the given time_stamp ts_.
-    //
-    // The position and velocity receive linear interpolation:
-    //    p_ret = (ts - t_ref) / dt * (p - p_ref)
-    //    v_ret = (ts - t_ref) / dt * (v - v_ref)
-    //
-    // the quaternion receives a slerp interpolation
-    //    q_ret = q_ref.slerp( (ts - t_ref) / dt , q);
-    //
-    // See extensive documentation in ProcessorMotion::interpolate().
-
-    // reference
-    TimeStamp           t_ref       = _motion_ref.ts_;
-
-    // final
-    TimeStamp           t_sec       = _motion_second.ts_;
-    Map<VectorXs>       dp_sec(_motion_second.delta_.data(), 3);
-    Map<Quaternions>    dq_sec(_motion_second.delta_.data() + 3);
-    Map<VectorXs>       dv_sec(_motion_second.delta_.data() + 7, 3);
-
-    // interpolated
-    Motion              motion_int  = motionZero(_ts);
-    Map<VectorXs>       dp_int(motion_int.delta_.data(), 3);
-    Map<Quaternions>    dq_int(motion_int.delta_.data() + 3);
-    Map<VectorXs>       dv_int(motion_int.delta_.data() + 7, 3);
-
-    // Jacobians for covariance propagation
-    MatrixXs            J_ref(delta_cov_size_, delta_cov_size_);
-    MatrixXs            J_int(delta_cov_size_, delta_cov_size_);
-
-    // interpolate delta
-    Scalar     tau  = (_ts - t_ref) / (t_sec - t_ref); // interpolation factor (0 to 1)
-    motion_int.ts_  = _ts;
-    dp_int          = tau * dp_sec;
-    dq_int          = Quaternions::Identity().slerp(tau, dq_sec);
-    dv_int          = tau * dv_sec;
-    deltaPlusDelta(_motion_ref.delta_integr_, motion_int.delta_, (t_sec - t_ref), motion_int.delta_integr_, J_ref, J_int);
-
-    // interpolate covariances
-    motion_int.delta_cov_           = tau * _motion_second.delta_cov_;
-    //motion_int.delta_integr_cov_    = J_ref * _motion_ref.delta_integr_cov_ * J_ref.transpose() + J_int * _motion_second.delta_cov_ * J_int.transpose();
-
-    // update second delta ( in place update )
-    dp_sec          = dq_int.conjugate() * ((1 - tau) * dp_sec);
-    dq_sec          = dq_int.conjugate() * dq_sec;
-    dv_sec          = dq_int.conjugate() * ((1 - tau) * dv_sec);
-    _motion_second.delta_cov_ = (1 - tau) * _motion_second.delta_cov_; // easy interpolation // TODO check for correctness
-    //Dp            = Dp; // trivial, just leave the code commented
-    //Dq            = Dq; // trivial, just leave the code commented
-    //_motion.delta_integr_cov_ = _motion.delta_integr_cov_; // trivial, just leave the code commented
-
-    return motion_int;
-}
-
-
 inline void ProcessorIMU::resetDerived()
 {
     // Cast a pointer to origin IMU frame
diff --git a/src/test/gtest_processor_imu.cpp b/src/test/gtest_processor_imu.cpp
index 0532a63e8839ca362a6b60e23794af77d81bd02a..f50e2a4e97173b385907ef7f741f38a8fd2d45df 100644
--- a/src/test/gtest_processor_imu.cpp
+++ b/src/test/gtest_processor_imu.cpp
@@ -206,6 +206,47 @@ TEST(ProcessorIMU, voteForKeyFrame)
 }
 
 //replace TEST by TEST_F if SetUp() needed
+TEST_F(ProcessorIMUt, interpolate)
+{
+    using namespace wolf;
+
+    t.set(0);
+    x0 << 0,0,0,  0,0,0,1,  0,0,0,  0,0,0,  0,0,0; // Try some non-zero biases
+
+    problem->getProcessorMotionPtr()->setOrigin(x0, t);
+
+    data << 2, 0, 0, 0, 0, 0; // only acc_x, but measure gravity!
+
+    // make two steps with half data, then simulate it was only one step
+    Motion mot_ref = problem->getProcessorMotionPtr()->getMotion();
+    cap_imu_ptr->setData(data/2);
+    cap_imu_ptr->setTimeStamp(0.05);
+    sensor_ptr->process(cap_imu_ptr);
+    Motion mot_int_gt = problem->getProcessorMotionPtr()->getMotion();
+    cap_imu_ptr->setTimeStamp(0.1);
+    sensor_ptr->process(cap_imu_ptr);
+    Motion mot_final = problem->getProcessorMotionPtr()->getMotion();
+    mot_final.data_ *= 2;
+    mot_final.delta_ = mot_final.delta_integr_;
+    Motion mot_sec = mot_final;
+
+class P : wolf::ProcessorIMU
+{
+    public:
+        Motion interp(const Motion& ref, Motion& sec, TimeStamp ts)
+        {
+            return ProcessorIMU::interpolate(ref, sec, ts);
+        }
+} imu;
+
+Motion mot_int = imu.interp(mot_ref, mot_sec, TimeStamp(0.05));
+
+ASSERT_MATRIX_APPROX(mot_int.data_,  mot_int_gt.data_, 1e-6);
+//ASSERT_MATRIX_APPROX(mot_int.delta_, mot_int_gt.delta_, 1e-6); // FIXME: delta_p not correctly interpolated
+
+
+}
+
 TEST_F(ProcessorIMUt, acc_x)
 {
     t.set(0); // clock in 0,1 ms ticks