diff --git a/docs/reference/index.rst b/docs/reference/index.rst
index 9a1cd7b5..4229bc52 100644
--- a/docs/reference/index.rst
+++ b/docs/reference/index.rst
@@ -29,3 +29,8 @@ API reference
    noise.gauss_markov
    noise.random_walk
    noise.white_noise
+   simulate.BeatDOF
+   simulate.ChirpDOF
+   simulate.ConstantDOF
+   simulate.IMUSimulator
+   simulate.SineDOF
diff --git a/docs/user_guide/quickstart.rst b/docs/user_guide/quickstart.rst
index 8e21368f..e6ca1811 100644
--- a/docs/user_guide/quickstart.rst
+++ b/docs/user_guide/quickstart.rst
@@ -26,17 +26,28 @@ data from an IMU sensor, and ideally position and heading data from other aiding
 sensors. If you do not have access to such data, you can generate synthetic
 measurements using the code provided here.
 
-Using the ``benchmark`` module, you can generate synthetic 3D motion data with ``smsfusion``.
-For example, you can generate beating signals representing position, velocity and
-attitude (PVA) degrees of freedom using :func:`~smsfusion.benchmark.benchmark_full_pva_beat_202311A`:
+Using the :class:`~smsfusion.simulate.IMUSimulator` class in the ``simulate`` module,
+you can generate synthetic 3D motion data and corresponding IMU accelerometer and
+gyroscope measurements. For example, you can simulate beating motion:
+
 
 .. code-block:: python
 
-    from smsfusion.benchmark import benchmark_full_pva_beat_202311A
+    from smsfusion.simulate import BeatDOF, IMUSimulator
+
 
+    sim = IMUSimulator(
+        pos_x=BeatDOF(0.5, phase=0.0, phase_degrees=True),
+        pos_y=BeatDOF(0.5, phase=45.0, phase_degrees=True),
+        pos_z=BeatDOF(0.5, phase=90.0, phase_degrees=True),
+        alpha=BeatDOF(0.1, phase=135.0, phase_degrees=True),
+        beta=BeatDOF(0.1, phase=180.0, phase_degrees=True),
+        gamma=BeatDOF(0.1, phase=225.0, phase_degrees=True),
+        degrees=False,
+    )
 
     fs = 10.24  # sampling rate in Hz
-    t, pos, vel, euler, acc, gyro = benchmark_full_pva_beat_202311A(fs)
+    t, pos, vel, euler, acc, gyro = sim(fs, 10_000, degrees=False)
     head = euler[:, 2]
 
 Note that the generated position signals are in meters (m), velocity signals are in meters
diff --git a/src/smsfusion/__init__.py b/src/smsfusion/__init__.py
index fb68cdbb..9e1cd436 100644
--- a/src/smsfusion/__init__.py
+++ b/src/smsfusion/__init__.py
@@ -1,4 +1,4 @@
-from . import benchmark, calibrate, constants, noise
+from . import benchmark, calibrate, constants, noise, simulate
 from ._ins import AHRS, VRU, AidedINS, FixedNED, StrapdownINS, gravity
 from ._smoothing import FixedIntervalSmoother
 from ._transforms import quaternion_from_euler
@@ -16,4 +16,5 @@
     "StrapdownINS",
     "VRU",
     "quaternion_from_euler",
+    "simulate",
 ]
diff --git a/src/smsfusion/benchmark/__init__.py b/src/smsfusion/benchmark/__init__.py
index 40e4ab88..848a4b83 100644
--- a/src/smsfusion/benchmark/__init__.py
+++ b/src/smsfusion/benchmark/__init__.py
@@ -1,7 +1,6 @@
 from ._benchmark import (
     BeatSignal,
     ChirpSignal,
-    _benchmark_helper,
     benchmark_full_pva_beat_202311A,
     benchmark_full_pva_chirp_202311A,
     benchmark_pure_attitude_beat_202311A,
@@ -9,7 +8,6 @@
 )
 
 __all__ = [
-    "_benchmark_helper",
     "BeatSignal",
     "benchmark_full_pva_beat_202311A",
     "benchmark_full_pva_chirp_202311A",
diff --git a/src/smsfusion/benchmark/_benchmark.py b/src/smsfusion/benchmark/_benchmark.py
index 1c66446d..d8dff047 100644
--- a/src/smsfusion/benchmark/_benchmark.py
+++ b/src/smsfusion/benchmark/_benchmark.py
@@ -1,12 +1,12 @@
 from __future__ import annotations
 
 import abc
+from warnings import warn
 
 import numpy as np
 from numpy.typing import ArrayLike, NDArray
 
-from smsfusion._ins import gravity
-from smsfusion._transforms import _inv_angular_matrix_from_euler, _rot_matrix_from_euler
+from smsfusion.simulate import BeatDOF, ChirpDOF, IMUSimulator
 
 
 class _Signal(abc.ABC):
@@ -86,47 +86,29 @@ class BeatSignal(_Signal):
     """
 
     def __init__(self, f_main: float, f_beat: float, freq_hz: bool = True) -> None:
-        self._f_main = f_main
-        self._f_beat = f_beat
-
-        if freq_hz:
-            self._f_main = self._f_main * 2.0 * np.pi
-            self._f_beat = self._f_beat * 2.0 * np.pi
+        warn("`BeatSignal`` is deprecated, use ``simulate.BeatDOF`` instead.")
+        self._f_main = 2.0 * np.pi * f_main if freq_hz else f_main
+        self._f_beat = f_beat * 2.0 * np.pi if freq_hz else f_beat
 
     def _y(self, t: NDArray[np.float64], phase: float) -> NDArray[np.float64]:
         """
         Generate a unit amplitude beating signal.
         """
-        main = np.cos(self._f_main * t + phase)
-        beat = np.sin(self._f_beat / 2.0 * t)
-        y = beat * main
+        y = BeatDOF(1.0, self._f_main, self._f_beat, phase=phase)._y(t)
         return y  # type: ignore[no-any-return]
 
     def _dydt(self, t: NDArray[np.float64], phase: float) -> NDArray[np.float64]:
         """
         Generate the first time derivative of a unit amplitue beating signal.
         """
-        main = np.cos(self._f_main * t + phase)
-        beat = np.sin(self._f_beat / 2.0 * t)
-        dmain = -self._f_main * np.sin(self._f_main * t + phase)
-        dbeat = self._f_beat / 2.0 * np.cos(self._f_beat / 2.0 * t)
-
-        dydt = dbeat * main + beat * dmain
+        dydt = BeatDOF(1.0, self._f_main, self._f_beat, phase=phase)._dydt(t)
         return dydt  # type: ignore[no-any-return]
 
     def _d2ydt2(self, t: NDArray[np.float64], phase: float) -> NDArray[np.float64]:
         """
         Generate the second time derivative of a unit amplitue beating signal.
         """
-        main = np.cos(self._f_main * t + phase)
-        beat = np.sin(self._f_beat / 2.0 * t)
-        dmain = -self._f_main * np.sin(self._f_main * t + phase)
-        dbeat = self._f_beat / 2.0 * np.cos(self._f_beat / 2.0 * t)
-        d2main = -((self._f_main) ** 2) * np.cos(self._f_main * t + phase)
-        d2beat = -((self._f_beat / 2.0) ** 2) * np.sin(self._f_beat / 2.0 * t)
-
-        d2ydt2 = dbeat * dmain + d2beat * main + beat * d2main + dbeat * dmain
-
+        d2ydt2 = BeatDOF(1.0, self._f_main, self._f_beat, phase=phase)._d2ydt2(t)
         return d2ydt2  # type: ignore[no-any-return]
 
 
@@ -157,19 +139,15 @@ class ChirpSignal(_Signal):
     """
 
     def __init__(self, f_max: float, f_os: float, freq_hz: bool = True) -> None:
-        self._f_max = f_max
-        self._f_os = f_os
-
-        if freq_hz:
-            self._f_max *= 2.0 * np.pi
-            self._f_os *= 2.0 * np.pi
+        warn("`ChirpSignal` is deprecated, use `simulate.ChirpDOF` instead.")
+        self._f_max = 2.0 * np.pi * f_max if freq_hz else f_max
+        self._f_os = 2.0 * np.pi * f_os if freq_hz else f_os
 
     def _y(self, t: NDArray[np.float64], phase: float) -> NDArray[np.float64]:
         """
         Generate a unit amplitude chirp signal with oscillating frequency.
         """
-        phi = 2.0 * self._f_max / self._f_os * np.sin(self._f_os / 2.0 * t)
-        y = np.sin(phi + phase)
+        y = ChirpDOF(1.0, self._f_max, self._f_os, phase=phase)._y(t)
         return y  # type: ignore[no-any-return]
 
     def _dydt(self, t: NDArray[np.float64], phase: float) -> NDArray[np.float64]:
@@ -177,9 +155,7 @@ def _dydt(self, t: NDArray[np.float64], phase: float) -> NDArray[np.float64]:
         Generate the time derivative of a unit amplitude chirp signal with
         oscillating frequency.
         """
-        phi = 2.0 * self._f_max / self._f_os * np.sin(self._f_os / 2.0 * t)
-        dphi = self._f_max * np.cos(self._f_os / 2.0 * t)
-        dydt = dphi * np.cos(phi + phase)
+        dydt = ChirpDOF(1.0, self._f_max, self._f_os, phase=phase)._dydt(t)
         return dydt  # type: ignore[no-any-return]
 
     def _d2ydt2(self, t: NDArray[np.float64], phase: float) -> NDArray[np.float64]:
@@ -187,161 +163,10 @@ def _d2ydt2(self, t: NDArray[np.float64], phase: float) -> NDArray[np.float64]:
         Generate the second time derivative of a unit amplitude chirp signal with
         oscillating frequency.
         """
-        phi = 2.0 * self._f_max / self._f_os * np.sin(self._f_os / 2.0 * t)
-        dphi = self._f_max * np.cos(self._f_os / 2.0 * t)
-        d2phi = -self._f_max * self._f_os / 2.0 * np.sin(self._f_os / 2.0 * t)
-        d2ydt2 = -(dphi**2) * np.sin(phi + phase) + d2phi * np.cos(phi + phase)
+        d2ydt2 = ChirpDOF(1.0, self._f_max, self._f_os, phase=phase)._d2ydt2(t)
         return d2ydt2  # type: ignore[no-any-return]
 
 
-def _benchmark_helper(
-    duration: float,
-    amplitude: tuple[float, float, float, float, float, float],
-    mean: tuple[float, float, float, float, float, float],
-    phase: tuple[float, float, float, float, float, float],
-    signal_family: _Signal,
-    fs: float,
-) -> tuple[
-    NDArray[np.float64],
-    NDArray[np.float64],
-    NDArray[np.float64],
-    NDArray[np.float64],
-    NDArray[np.float64],
-    NDArray[np.float64],
-]:
-    """
-    Generate benchmark scenarios for INS sensor fusion performance evaluation.
-
-    The general helper function works by first generating position and attitude
-    based on the provided input, and then derives the remaining variables. The
-    final output consists of:
-
-    * Time [s]
-    * Position [m]
-    * Velocity [m/s]
-    * Attitude [rad] (Euler angles, see Notes)
-    * Accelerations [m/s**2] in body frame (corresponding to accelerometer measurements)
-    * Angular rates [rad/s] in body frame (corresponding to gyroscope measurements)
-
-    Position, velocity, and attitude are relative to the *NED* frame and expressed
-    in the *NED* frame. The position origin is at (0., 0., 0.).
-
-    Parameters
-    ----------
-    duration : float
-        Duration of the generated signals in seconds.
-    amplitude : tuple of (float, float, float, float, float, float)
-        Amplitude of the signals generated by the generator defined in ``signal_family``.
-        The values represent the amplitude of postions in x, y, and z directions in
-        meters, and Euler angles as roll, pitch, and yaw in radians (see Notes).
-        The order is given as [pos_x, pos_y, pos_z, roll, pitch, yaw].
-    mean : tuple of (float, float, float, float, float, float)
-        Mean values of the generated signals. Follows the same conventions as
-        ``amplitude``.
-    phase : tuple of (float, float, float, float, float, float)
-        Phase values in radians passed on to the signal generator in
-        ``signal_family``. Otherwise, follows the same conventions as
-        ``amplitude``.
-    signal_family : {:clas:`smsfusion.benchmark.BeatSignal`, :clas:`smsfusion.benchmark.BeatSignal`}
-        Instance of a class used to generate the signals.
-    fs : float
-        Sampling rate of the outputs in hertz.
-
-    Returns
-    -------
-    t : numpy.ndarray, shape (N,)
-        Time in seconds.
-    position : numpy.ndarray, shape (N, 3)
-        Position [m] of the body.
-    velocity : numpy.ndarray, shape (N, 3)
-        Velocity [m/s] of the body.
-    euler : numpy.ndarray, shape (N, 3)
-        Attitude of the body in Euler angles [rad], see Notes.
-    acc : numpy.ndarray, shape (N, 3)
-        Accelerations [m/s**2] in body frame (corresponding to accelerometer measurements).
-    gyro : numpy.ndarray, shape (N, 3)
-        Angular rates [rad/s] in body frame (corresponding to gyroscope measurements).
-
-    Notes
-    -----
-    The length of the signals are determined by ``np.arange(0.0, duration, 1.0 / fs)``.
-
-    The Euler angles describe how to transition from the 'NED' frame to the 'body'
-    frame through three consecutive intrinsic and passive rotations in the ZYX order:
-
-    #. A rotation by an angle gamma (often called yaw) about the z-axis.
-    #. A subsequent rotation by an angle beta (often called pitch) about the y-axis.
-    #. A final rotation by an angle alpha (often called roll) about the x-axis.
-
-    This sequence of rotations is used to describe the orientation of the 'body' frame
-    relative to the 'NED' frame in 3D space.
-
-    Intrinsic rotations mean that the rotations are with respect to the changing
-    coordinate system; as one rotation is applied, the next is about the axis of
-    the newly rotated system.
-
-    Passive rotations mean that the frame itself is rotating, not the object
-    within the frame.
-    """
-    t = np.arange(0.0, duration, 1.0 / fs)
-
-    amp_x, amp_y, amp_z, amp_roll, amp_pitch, amp_yaw = amplitude
-    mean_x, mean_y, mean_z, mean_roll, mean_pitch, mean_yaw = mean
-    phase_x, phase_y, phase_z, phase_roll, phase_pitch, phase_yaw = phase
-
-    pos_x_, vel_x_, acc_x_ = signal_family(t, phase_x, phase_degrees=False)
-    pos_x = mean_x + amp_x * pos_x_
-    vel_x = amp_x * vel_x_
-    acc_x = amp_x * acc_x_
-
-    pos_y_, vel_y_, acc_y_ = signal_family(t, phase_y, phase_degrees=False)
-    pos_y = mean_y + amp_y * pos_y_
-    vel_y = amp_y * vel_y_
-    acc_y = amp_y * acc_y_
-
-    pos_z_, vel_z_, acc_z_ = signal_family(t, phase_z, phase_degrees=False)
-    pos_z = mean_z + amp_z * pos_z_
-    vel_z = amp_z * vel_z_
-    acc_z = amp_z * acc_z_
-
-    roll_, droll_, _ = signal_family(t, phase_roll, phase_degrees=False)
-    roll = mean_roll + amp_roll * roll_
-    droll = amp_roll * droll_
-
-    pitch_, dpitch_, _ = signal_family(t, phase_pitch, phase_degrees=False)
-    pitch = mean_pitch + amp_pitch * pitch_
-    dpitch = amp_pitch * dpitch_
-
-    yaw_, dyaw_, _ = signal_family(t, phase_yaw, phase_degrees=False)
-    yaw = mean_yaw + amp_yaw * yaw_
-    dyaw = amp_yaw * dyaw_
-
-    position = np.column_stack((pos_x, pos_y, pos_z))
-    velocity = np.column_stack((vel_x, vel_y, vel_z))
-    acceleration = np.column_stack((acc_x, acc_y, acc_z))
-    attitude = np.column_stack((roll, pitch, yaw))
-    dattitude = np.column_stack((droll, dpitch, dyaw))
-
-    accelerometer = []
-    gyroscope = []
-    for euler_i, deuler_i, acc_i in zip(attitude, dattitude, acceleration):
-        gyroscope.append(_inv_angular_matrix_from_euler(euler_i).dot(deuler_i))
-        accelerometer.append(
-            _rot_matrix_from_euler(euler_i).T.dot(
-                acc_i + np.array([0.0, 0.0, -gravity()])
-            )
-        )
-
-    return (
-        t,
-        position,
-        velocity,
-        attitude,
-        np.asarray(accelerometer),
-        np.asarray(gyroscope),
-    )
-
-
 def benchmark_pure_attitude_beat_202311A(
     fs: float = 10.24,
 ) -> tuple[
@@ -358,7 +183,7 @@ def benchmark_pure_attitude_beat_202311A(
     The generated signals are "beating" with maximum amplitudes corresponding to
     5 degrees. The main signal frequency is 0.1 Hz while the "beating" frequency
     is 0.01 Hz. The phases for roll, pitch, and yaw are 0.0, 45.0, and 90.0 degrees
-    respectively. See :class:`smsfusion.benchmark.BeatSignal` for details.
+    respectively. See :class:`smsfusion.simulate.BeatDOF` for details.
 
     Parameters
     ----------
@@ -396,17 +221,20 @@ def benchmark_pure_attitude_beat_202311A(
     within the frame.
     """
     duration = 1800.0  # 30 minutes
-    amplitude = np.radians(np.array([0.0, 0.0, 0.0, 5.0, 5.0, 5.0]))
-    mean = np.radians(np.array([0.0, 0.0, 0.0, 0.0, 0.0, 0.0]))
-    phase = np.radians(np.array([0.0, 0.0, 0.0, 0.0, 45.0, 90.0]))
 
     f_main = 0.1
     f_beat = 0.01
 
-    t, _, _, euler, acc, gyro = _benchmark_helper(
-        duration, amplitude, mean, phase, BeatSignal(f_main, f_beat), fs
-    )
-    return t, euler, acc, gyro
+    amp = np.radians(5.0)
+    alpha = BeatDOF(amp, f_main, f_beat, freq_hz=True, phase=0.0)
+    beta = BeatDOF(amp, f_main, f_beat, freq_hz=True, phase=45.0, phase_degrees=True)
+    gamma = BeatDOF(amp, f_main, f_beat, freq_hz=True, phase=90.0, phase_degrees=True)
+    sim = IMUSimulator(alpha=alpha, beta=beta, gamma=gamma)
+
+    n = int(duration * fs)
+    t, _, _, euler, f, w = sim(fs, n, degrees=False)
+
+    return t, euler, f, w
 
 
 def benchmark_pure_attitude_chirp_202311A(
@@ -426,7 +254,7 @@ def benchmark_pure_attitude_chirp_202311A(
     amplitudes corresponding to 5 degrees. The signal frequency oscillates
     between 0.0 Hz and 0.25 Hz every 100 seconds. The phases for roll, pitch,
     and yaw are 0.0, 45.0, and 90.0 degrees respectively. See
-    :class:`smsfusion.benchmark.ChirpSignal` for details.
+    :class:`smsfusion.simulate.ChirpDOF` for details.
 
     Parameters
     ----------
@@ -464,17 +292,20 @@ def benchmark_pure_attitude_chirp_202311A(
     within the frame.
     """
     duration = 1800.0  # 30 minutes
-    amplitude = np.radians(np.array([0.0, 0.0, 0.0, 5.0, 5.0, 5.0]))
-    mean = np.radians(np.array([0.0, 0.0, 0.0, 0.0, 0.0, 0.0]))
-    phase = np.radians(np.array([0.0, 0.0, 0.0, 0.0, 45.0, 90.0]))
 
     f_max = 0.25
     f_os = 0.01
 
-    t, _, _, euler, acc, gyro = _benchmark_helper(
-        duration, amplitude, mean, phase, ChirpSignal(f_max, f_os), fs
-    )
-    return t, euler, acc, gyro
+    amp = np.radians(5.0)
+    alpha = ChirpDOF(amp, f_max, f_os, freq_hz=True, phase=0.0)
+    beta = ChirpDOF(amp, f_max, f_os, freq_hz=True, phase=45.0, phase_degrees=True)
+    gamma = ChirpDOF(amp, f_max, f_os, freq_hz=True, phase=90.0, phase_degrees=True)
+    sim = IMUSimulator(alpha=alpha, beta=beta, gamma=gamma)
+
+    n = int(duration * fs)
+    t, _, _, euler, f, w = sim(fs, n, degrees=False)
+
+    return t, euler, f, w
 
 
 def benchmark_full_pva_beat_202311A(
@@ -496,13 +327,13 @@ def benchmark_full_pva_beat_202311A(
 
     The generated position reference signals are characterized by:
 
-    * "Beating" signal. See ``BeatSignal`` for details.
+    * "Beating" signal. See ``smsfusion.simulate.BeatDOF`` for details.
     * Maximum amplitude of 0.5 m.
     * The phases for x-, y-, and z-axis are 0.0, 30.0, and 60.0 degrees respectively.
 
     The generated Euler reference signals are characterized by:
 
-    * "Beating" signal. See :class:`smsfusion.benchmark.BeatSignal` for details.
+    * "Beating" signal. See :class:`smsfusion.simulate.BeatDOF` for details.
     * Maximum amplitude of 5 degrees.
     * The phases for roll, pitch, and yaw are 90.0, 120.0, and 150.0 degrees respectively.
 
@@ -548,17 +379,26 @@ def benchmark_full_pva_beat_202311A(
     within the frame.
     """
     duration = 1800.0  # 30 minutes
-    amplitude = (0.5, 0.5, 0.5, np.radians(5.0), np.radians(5.0), np.radians(5.0))
-    mean = (0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
-    phase = np.radians(np.array([0.0, 30.0, 60.0, 90.0, 120.0, 150.0]))
 
     f_main = 0.1
     f_beat = 0.01
 
-    t, pos, vel, euler, acc, gyro = _benchmark_helper(
-        duration, amplitude, mean, phase, BeatSignal(f_main, f_beat), fs
+    amp_p = 0.5
+    amp_r = np.radians(5.0)
+    pos_x = BeatDOF(amp_p, f_main, f_beat, freq_hz=True, phase=0.0, phase_degrees=True)
+    pos_y = BeatDOF(amp_p, f_main, f_beat, freq_hz=True, phase=30.0, phase_degrees=True)
+    pos_z = BeatDOF(amp_p, f_main, f_beat, freq_hz=True, phase=60.0, phase_degrees=True)
+    alpha = BeatDOF(amp_r, f_main, f_beat, freq_hz=True, phase=90.0, phase_degrees=True)
+    beta = BeatDOF(amp_r, f_main, f_beat, freq_hz=True, phase=120.0, phase_degrees=True)
+    gamma = BeatDOF(
+        amp_r, f_main, f_beat, freq_hz=True, phase=150.0, phase_degrees=True
     )
-    return t, pos, vel, euler, acc, gyro
+    sim = IMUSimulator(pos_x, pos_y, pos_z, alpha, beta, gamma)
+
+    n = int(duration * fs)
+    t, pos, vel, euler, f, w = sim(fs, n, degrees=False)
+
+    return t, pos, vel, euler, f, w
 
 
 def benchmark_full_pva_chirp_202311A(
@@ -580,13 +420,13 @@ def benchmark_full_pva_chirp_202311A(
 
     The generated position reference signals are characterized by:
 
-    * "Chirp" signal. See :class:`smsfusion.benchmark.ChirpSignal` for details.
+    * "Chirp" signal. See :class:`smsfusion.simulate.ChirpDOF` for details.
     * Maximum amplitude of 0.5 m.
     * The phases for x-, y-, and z-axis are 0.0, 30.0, and 60.0 degrees respectively.
 
     The generated Euler reference signals are characterized by:
 
-    * "Chirp" signal. See :class:`smsfusion.benchmark.ChirpSignal` for details.
+    * "Chirp" signal. See :class:`smsfusion.simulate.ChirpDOF` for details.
     * Maximum amplitude of 5 degrees.
     * The phases for roll, pitch, and yaw are 90.0, 120.0, and 150.0 degrees respectively.
 
@@ -632,14 +472,21 @@ def benchmark_full_pva_chirp_202311A(
     within the frame.
     """
     duration = 1800.0  # 30 minutes
-    amplitude = (0.5, 0.5, 0.5, np.radians(5.0), np.radians(5.0), np.radians(5.0))
-    mean = (0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
-    phase = np.radians(np.array([0.0, 30.0, 60.0, 90.0, 120.0, 150.0]))
 
     f_max = 0.25
     f_os = 0.01
 
-    t, pos, vel, euler, acc, gyro = _benchmark_helper(
-        duration, amplitude, mean, phase, ChirpSignal(f_max, f_os), fs
-    )
-    return t, pos, vel, euler, acc, gyro
+    amp_p = 0.5
+    amp_r = np.radians(5.0)
+    pos_x = ChirpDOF(amp_p, f_max, f_os, freq_hz=True, phase=0.0, phase_degrees=True)
+    pos_y = ChirpDOF(amp_p, f_max, f_os, freq_hz=True, phase=30.0, phase_degrees=True)
+    pos_z = ChirpDOF(amp_p, f_max, f_os, freq_hz=True, phase=60.0, phase_degrees=True)
+    alpha = ChirpDOF(amp_r, f_max, f_os, freq_hz=True, phase=90.0, phase_degrees=True)
+    beta = ChirpDOF(amp_r, f_max, f_os, freq_hz=True, phase=120.0, phase_degrees=True)
+    gamma = ChirpDOF(amp_r, f_max, f_os, freq_hz=True, phase=150.0, phase_degrees=True)
+    sim = IMUSimulator(pos_x, pos_y, pos_z, alpha, beta, gamma)
+
+    n = int(duration * fs)
+    t, pos, vel, euler, f, w = sim(fs, n, degrees=False)
+
+    return t, pos, vel, euler, f, w
diff --git a/src/smsfusion/simulate/__init__.py b/src/smsfusion/simulate/__init__.py
new file mode 100644
index 00000000..d9b7dec5
--- /dev/null
+++ b/src/smsfusion/simulate/__init__.py
@@ -0,0 +1,3 @@
+from ._simulate import BeatDOF, ChirpDOF, ConstantDOF, IMUSimulator, SineDOF
+
+__all__ = ["IMUSimulator", "SineDOF", "ConstantDOF", "BeatDOF", "ChirpDOF"]
diff --git a/src/smsfusion/simulate/_simulate.py b/src/smsfusion/simulate/_simulate.py
new file mode 100644
index 00000000..5908e786
--- /dev/null
+++ b/src/smsfusion/simulate/_simulate.py
@@ -0,0 +1,533 @@
+from abc import ABC, abstractmethod
+
+import numpy as np
+from numpy.typing import ArrayLike, NDArray
+
+from smsfusion._transforms import _rot_matrix_from_euler
+
+
+class DOF(ABC):
+    """
+    Abstract base class for degree of freedom (DOF) signal generators.
+    """
+
+    @abstractmethod
+    def _y(self, t: NDArray[np.float64]) -> NDArray[np.float64]:
+        raise NotImplementedError("Not implemented.")
+
+    @abstractmethod
+    def _dydt(self, t: NDArray[np.float64]) -> NDArray[np.float64]:
+        raise NotImplementedError("Not implemented.")
+
+    @abstractmethod
+    def _d2ydt2(self, t: NDArray[np.float64]) -> NDArray[np.float64]:
+        raise NotImplementedError("Not implemented.")
+
+    def y(self, t: ArrayLike) -> NDArray[np.float64]:
+        """
+        Generates y(t) signal.
+
+        Parameters
+        ----------
+        t : array_like, shape (n,)
+            Time vector in seconds.
+        """
+        t = np.asarray_chkfinite(t)
+        return self._y(t)
+
+    def dydt(self, t: ArrayLike) -> NDArray[np.float64]:
+        """
+        Generates dy(t)/dt signal.
+
+        Parameters
+        ----------
+        t : array_like, shape (n,)
+            Time vector in seconds.
+        """
+        t = np.asarray_chkfinite(t)
+        return self._dydt(t)
+
+    def d2ydt2(self, t):
+        """
+        Generates d2y(t)/dt2 signal.
+
+        Parameters
+        ----------
+        t : array_like, shape (n,)
+            Time vector in seconds.
+        """
+        t = np.asarray_chkfinite(t)
+        return self._d2ydt2(t)
+
+    def __call__(
+        self, t: ArrayLike
+    ) -> tuple[NDArray[np.float64], NDArray[np.float64], NDArray[np.float64]]:
+        """
+        Generates y(t), dy(t)/dt, and d2y(t)/dt2 signals.
+
+        Parameters
+        ----------
+        t : array_like, shape (n,)
+            Time vector in seconds.
+
+        Returns
+        -------
+        y : ndarray, shape (n,)
+            DOF signal y(t).
+        dydt : ndarray, shape (n,)
+            Time derivative, dy(t)/dt, of DOF signal.
+        d2ydt2 : ndarray, shape (n,)
+            Second time derivative, d2y(t)/dt2, of DOF signal.
+        """
+        y = self._y(t)
+        dydt = self._dydt(t)
+        d2ydt2 = self._d2ydt2(t)
+
+        return y, dydt, d2ydt2
+
+
+class SineDOF(DOF):
+    """
+    1D sine wave DOF signal generator.
+
+    Defined as:
+
+        y(t) = A * sin(w * t + phi) + B
+        dy(t)/dt = A * w * cos(w * t + phi)
+        d2y(t)/dt2 = -A * w^2 * sin(w * t + phi)
+
+    where,
+
+    - A  : Amplitude of the sine wave.
+    - w  : Angular frequency of the sine wave.
+    - phi: Phase offset of the sine wave.
+    - B  : Constant offset of the sine wave.
+
+    Parameters
+    ----------
+    amp : float, default 1.0
+        Amplitude of the sine wave. Default is 1.0.
+    freq : float, default 1.0
+        Frequency of the sine wave in rad/s. Default is 1.0 rad/s.
+    freq_hz : bool, optional
+        If True, interpret `omega` as frequency in Hz. If False, interpret as angular
+        frequency in radians per second. Default is False.
+    phase : float, default 0.0
+        Phase offset of the sine wave. Default is 0.0.
+    phase_degrees : bool, optional
+        If True, interpret `phase` in degrees. If False, interpret in radians.
+        Default is False.
+    offset : float, default 0.0
+        Offset of the sine wave. Default is 0.0.
+    """
+
+    def __init__(
+        self,
+        amp: float = 1.0,
+        freq: float = 1.0,
+        freq_hz: bool = False,
+        phase: float = 0.0,
+        phase_degrees: bool = False,
+        offset: float = 0.0,
+    ) -> None:
+        self._amp = amp
+        self._w = 2.0 * np.pi * freq if freq_hz else freq
+        self._phase = np.deg2rad(phase) if phase_degrees else phase
+        self._offset = offset
+
+    def _y(self, t: NDArray[np.float64]) -> NDArray[np.float64]:
+        y = self._amp * np.sin(self._w * t + self._phase) + self._offset
+        return y
+
+    def _dydt(self, t: NDArray[np.float64]) -> NDArray[np.float64]:
+        dydt = self._amp * self._w * np.cos(self._w * t + self._phase)
+        return dydt
+
+    def _d2ydt2(self, t: NDArray[np.float64]) -> NDArray[np.float64]:
+        d2ydt2 = -self._amp * self._w**2 * np.sin(self._w * t + self._phase)
+        return d2ydt2
+
+
+class ConstantDOF(DOF):
+    """
+    1D constant DOF signal generator.
+
+    Defined as:
+
+        y(t) = C
+        dy(t)/dt = 0
+        d2y(t)/dt2 = 0
+
+    where,
+
+    - C : Constant value of the signal.
+
+    Parameters
+    ----------
+    value : float, default 0.0
+        Constant value of the signal. Default is 0.0.
+    """
+
+    def __init__(self, value: float = 0.0) -> None:
+        self._value = value
+
+    def _y(self, t: NDArray[np.float64]) -> NDArray[np.float64]:
+        return np.full_like(t, self._value)
+
+    def _dydt(self, t: NDArray[np.float64]) -> NDArray[np.float64]:
+        return np.zeros_like(t)
+
+    def _d2ydt2(self, t: NDArray[np.float64]) -> NDArray[np.float64]:
+        return np.zeros_like(t)
+
+
+class BeatDOF(DOF):
+    """
+    1D beating sinusoidal DOF signal generator.
+
+    Defined as:
+
+        y = A * sin(f_beat / 2.0 * t) * cos(f_main * t + phi) + B
+
+    where,
+
+    - A      : Amplitude of the sine waves.
+    - w_main : Angular frequency of the main sine wave.
+    - w_beat : Angular frequency of the beat sine wave.
+    - phi    : Phase offset of the main sine wave.
+    - B      : Constant offset of the beat signal.
+
+    Parameters
+    ----------
+    f_main : float
+        The main frequency of the sinusoidal signal, y(t).
+    f_beat : float
+        The beating frequency, which controls the variation in amplitude.
+    freq_hz : bool, default True.
+        Whether the frequencies, ``f_main`` and ``f_beat``, are in Hz or rad/s (default).
+    phase : float, default 0.0
+        Phase offset of the beat signal. Default is 0.0.
+    phase_degrees : bool, optional
+        If True, interpret `phase` in degrees. If False, interpret in radians.
+        Default is False.
+    offset : float, default 0.0
+        Offset of the beat signal. Default is 0.0.
+    """
+
+    def __init__(
+        self,
+        amp: float = 1.0,
+        freq_main: float = 0.1,
+        freq_beat: float = 0.01,
+        freq_hz: bool = False,
+        phase: float = 0.0,
+        phase_degrees: bool = False,
+        offset: float = 0.0,
+    ) -> None:
+        self._amp = amp
+        self._w_main = 2.0 * np.pi * freq_main if freq_hz else freq_main
+        self._w_beat = 2.0 * np.pi * freq_beat if freq_hz else freq_beat
+        self._phase = np.deg2rad(phase) if phase_degrees else phase
+        self._offset = offset
+
+    def _y(self, t: NDArray[np.float64]) -> NDArray[np.float64]:
+        amp = self._amp
+        w_main = self._w_main
+        w_beat = self._w_beat
+        phase = self._phase
+        offset = self._offset
+
+        main = np.cos(w_main * t + phase)
+        beat = np.sin(w_beat / 2.0 * t)
+        y = amp * beat * main + offset
+        return y  # type: ignore[no-any-return]
+
+    def _dydt(self, t: NDArray[np.float64]) -> NDArray[np.float64]:
+        amp = self._amp
+        w_main = self._w_main
+        w_beat = self._w_beat
+        phase = self._phase
+
+        main = np.cos(w_main * t + phase)
+        beat = np.sin(w_beat / 2.0 * t)
+        dmain = -w_main * np.sin(w_main * t + phase)
+        dbeat = w_beat / 2.0 * np.cos(w_beat / 2.0 * t)
+
+        dydt = amp * (dbeat * main + beat * dmain)
+        return dydt  # type: ignore[no-any-return]
+
+    def _d2ydt2(self, t: NDArray[np.float64]) -> NDArray[np.float64]:
+
+        amp = self._amp
+        w_main = self._w_main
+        w_beat = self._w_beat
+        phase = self._phase
+
+        main = np.cos(w_main * t + phase)
+        beat = np.sin(w_beat / 2.0 * t)
+        dmain = -w_main * np.sin(w_main * t + phase)
+        dbeat = w_beat / 2.0 * np.cos(w_beat / 2.0 * t)
+        d2main = -(w_main**2) * np.cos(w_main * t + phase)
+        d2beat = -((w_beat / 2.0) ** 2) * np.sin(w_beat / 2.0 * t)
+        d2ydt2 = amp * (dbeat * dmain + d2beat * main + beat * d2main + dbeat * dmain)
+
+        return d2ydt2  # type: ignore[no-any-return]
+
+
+class ChirpDOF(DOF):
+    """
+    1D chirp sinusoidal DOF signal generator.
+
+    This class creates a signal with a frequency that appears to vary in time.
+    The frequency oscillates between 0. and a maximum frequency at a specific
+    rate.
+
+    Defined as:
+
+        phi = 2 * f_max / f_os * sin(f_os * t)
+        y = sin(phi + phase)
+
+    where,
+
+    - f_max : Maximum frequency the signal ramps up to, before ramping down to zero.
+    - f_os  : Frequency oscillation rate.
+    - phase : Phase offset of the chirp signal.
+
+    Parameters
+    ----------
+    f_max : float
+        The maximum frequency of the chirp signal, y(t). Default is 0.25 Hz.
+    f_os : float
+        The frequency oscillation rate. Default is 0.01 Hz.
+    freq_hz : bool, default False.
+        Whether the frequencies, ``f_max`` and ``f_os``, are in Hz or rad/s (default).
+    phase : float, default 0.0
+        Phase offset of the chirp signal. Default is 0.0.
+    phase_degrees : bool, optional
+        If True, interpret `phase` in degrees. If False, interpret in radians.
+        Default is False.
+    offset : float, default 0.0
+        Offset of the chirp signal. Default is 0.0.
+    """
+
+    def __init__(
+        self,
+        amp: float = 1.0,
+        f_max: float = 0.5 * np.pi,
+        f_os: float = 0.02 * np.pi,
+        freq_hz: bool = False,
+        phase: float = 0.0,
+        phase_degrees: bool = False,
+        offset: float = 0.0,
+    ) -> None:
+        self._amp = amp
+        self._w_max = 2.0 * np.pi * f_max if freq_hz else f_max
+        self._w_os = 2.0 * np.pi * f_os if freq_hz else f_os
+        self._phase = np.deg2rad(phase) if phase_degrees else phase
+        self._offset = offset
+
+    def _y(self, t: NDArray[np.float64]) -> NDArray[np.float64]:
+        amp = self._amp
+        w_max = self._w_max
+        w_os = self._w_os
+        phase = self._phase
+        offset = self._offset
+
+        phi = 2.0 * w_max / w_os * np.sin(w_os / 2.0 * t)
+        y = amp * np.sin(phi + phase) + offset
+        return y  # type: ignore[no-any-return]
+
+    def _dydt(self, t: NDArray[np.float64]) -> NDArray[np.float64]:
+        amp = self._amp
+        w_max = self._w_max
+        w_os = self._w_os
+        phase = self._phase
+
+        phi = 2.0 * w_max / w_os * np.sin(w_os / 2.0 * t)
+        dphi = w_max * np.cos(w_os / 2.0 * t)
+        dydt = amp * dphi * np.cos(phi + phase)
+        return dydt  # type: ignore[no-any-return]
+
+    def _d2ydt2(self, t: NDArray[np.float64]) -> NDArray[np.float64]:
+        amp = self._amp
+        w_max = self._w_max
+        w_os = self._w_os
+        phase = self._phase
+
+        phi = 2.0 * w_max / w_os * np.sin(w_os / 2.0 * t)
+        dphi = w_max * np.cos(w_os / 2.0 * t)
+        d2phi = -w_max * w_os / 2.0 * np.sin(w_os / 2.0 * t)
+        d2ydt2 = -amp * (dphi**2) * np.sin(phi + phase) + amp * d2phi * np.cos(
+            phi + phase
+        )
+        return d2ydt2  # type: ignore[no-any-return]
+
+
+class IMUSimulator:
+    """
+    IMU simulator.
+
+    Parameters
+    ----------
+    pos_x : float or DOF, default 0.0
+        X position signal.
+    pos_y : float or DOF, default 0.0
+        Y position signal.
+    pos_z : float or DOF, default 0.0
+        Z position signal.
+    alpha : float or DOF, default 0.0
+        Roll signal.
+    beta : float or DOF, default 0.0
+        Pitch signal
+    gamma : float or DOF, default 0.0
+        Yaw signal
+    degrees: bool, default False
+        Whether to interpret the Euler angle signals as degrees (True) or radians (False).
+        Default is False.
+    g : float, default 9.80665
+        The gravitational acceleration. Default is 'standard gravity' of 9.80665.
+    nav_frame : str, default "NED"
+        Navigation frame. Either "NED" (North-East-Down) or "ENU" (East-North-Up).
+        Default is "NED".
+    """
+
+    def __init__(
+        self,
+        pos_x: float | DOF = 0.0,
+        pos_y: float | DOF = 0.0,
+        pos_z: float | DOF = 0.0,
+        alpha: float | DOF = 0.0,
+        beta: float | DOF = 0.0,
+        gamma: float | DOF = 0.0,
+        degrees: bool = False,
+        g: float = 9.80665,
+        nav_frame: str = "NED",
+    ) -> None:
+        self._pos_x = pos_x if isinstance(pos_x, DOF) else ConstantDOF(pos_x)
+        self._pos_y = pos_y if isinstance(pos_y, DOF) else ConstantDOF(pos_y)
+        self._pos_z = pos_z if isinstance(pos_z, DOF) else ConstantDOF(pos_z)
+        self._alpha = alpha if isinstance(alpha, DOF) else ConstantDOF(alpha)
+        self._beta = beta if isinstance(beta, DOF) else ConstantDOF(beta)
+        self._gamma = gamma if isinstance(gamma, DOF) else ConstantDOF(gamma)
+        self._degrees = degrees
+        self._nav_frame = nav_frame.lower()
+
+        if self._nav_frame == "ned":
+            self._g_n = np.array([0.0, 0.0, g])
+        elif self._nav_frame == "enu":
+            self._g_n = np.array([0.0, 0.0, -g])
+        else:
+            raise ValueError("Invalid navigation frame. Must be 'NED' or 'ENU'.")
+
+    def _specific_force_body(
+        self,
+        pos: NDArray[np.float64],
+        acc: NDArray[np.float64],
+        euler: NDArray[np.float64],
+    ) -> NDArray[np.float64]:
+        """
+        Specific force in the body frame.
+
+        Parameters
+        ----------
+        pos : ndarray, shape (n, 3)
+            Position [x, y, z]^T in meters.
+        vel : ndarray, shape (n, 3)
+            Velocity [x_dot, y_dot, z_dot]^T in meters per second.
+        acc : ndarray, shape (n, 3)
+            Acceleration [x_ddot, y_ddot, z_ddot]^T in meters per second squared.
+        euler : ndarray, shape (n, 3)
+            Euler angles [alpha, beta, gamma]^T in radians.
+        """
+        n = pos.shape[0]
+        f_b = np.zeros((n, 3))
+
+        for i in range(n):
+            R_i = _rot_matrix_from_euler(euler[i])
+            f_b[i] = R_i.T.dot(acc[i] - self._g_n)
+
+        return f_b
+
+    def _angular_velocity_body(
+        self, euler: NDArray[np.float64], euler_dot: NDArray[np.float64]
+    ) -> NDArray[np.float64]:
+        """
+        Angular velocity in the body frame.
+
+        Parameters
+        ----------
+        euler : ndarray, shape (n, 3)
+            Euler angles [alpha, beta, gamma]^T in radians.
+        euler_dot : ndarray, shape (n, 3)
+            Time derivatives of Euler angles [alpha_dot, beta_dot, gamma_dot]^T
+            in radians per second.
+        """
+        alpha, beta, _ = euler.T
+        alpha_dot, beta_dot, gamma_dot = euler_dot.T
+
+        w_x = alpha_dot - np.sin(beta) * gamma_dot
+        w_y = np.cos(alpha) * beta_dot + np.sin(alpha) * np.cos(beta) * gamma_dot
+        w_z = -np.sin(alpha) * beta_dot + np.cos(alpha) * np.cos(beta) * gamma_dot
+
+        w_b = np.column_stack([w_x, w_y, w_z])
+
+        return w_b
+
+    def __call__(self, fs: float, n: int, degrees: bool | None = None):
+        """
+        Generate a length-n gyroscope signal and corresponding Euler angles.
+
+        Parameters
+        ----------
+        fs : float
+            Sampling frequency in Hz.
+        n : int
+            Number of samples to generate.
+        degrees : bool, optional
+            Whether to return Euler angles and angular velocities in degrees and
+            degrees per second (True) or radians and radians per second (False).
+            Defaults to the value specified at initialization.
+
+        Returns
+        -------
+        t : ndarray, shape (n,)
+            Time vector in seconds.
+        euler : ndarray, shape (n, 3)
+            Simulated (ZYX) Euler angles [roll, pitch, yaw]^T.
+        w_b : ndarray, shape (n, 3)
+            Simulated angular velocities, [w_x, w_y, w_z]^T, in the body frame.
+        """
+        if degrees is None:
+            degrees = self._degrees
+
+        # Time
+        dt = 1.0 / fs
+        t = dt * np.arange(n)
+
+        # DOFs and corresponding rates and accelerations
+        pos_x, pos_x_dot, pos_x_ddot = self._pos_x(t)
+        pos_y, pos_y_dot, pos_y_ddot = self._pos_y(t)
+        pos_z, pos_z_dot, pos_z_ddot = self._pos_z(t)
+        alpha, alpha_dot, _ = self._alpha(t)
+        beta, beta_dot, _ = self._beta(t)
+        gamma, gamma_dot, _ = self._gamma(t)
+
+        pos = np.column_stack([pos_x, pos_y, pos_z])
+        vel = np.column_stack([pos_x_dot, pos_y_dot, pos_z_dot])
+        acc = np.column_stack([pos_x_ddot, pos_y_ddot, pos_z_ddot])
+        euler = np.column_stack([alpha, beta, gamma])
+        euler_dot = np.column_stack([alpha_dot, beta_dot, gamma_dot])
+
+        if self._degrees:
+            euler = np.deg2rad(euler)
+            euler_dot = np.deg2rad(euler_dot)
+
+        # IMU measurements (i.e., specific force and angular velocity in body frame)
+        f_b = self._specific_force_body(pos, acc, euler)
+        w_b = self._angular_velocity_body(euler, euler_dot)
+
+        if degrees:
+            euler = np.rad2deg(euler)
+            w_b = np.rad2deg(w_b)
+
+        return t, pos, vel, euler, f_b, w_b
diff --git a/tests/test_benchmark.py b/tests/test_benchmark.py
index d8fc1566..6ab0575d 100644
--- a/tests/test_benchmark.py
+++ b/tests/test_benchmark.py
@@ -95,211 +95,6 @@ def test_signal(self):
         assert not np.array_equal(d2ydt2, d2ydt2_phase)
 
 
-class Test__benchmark_helper:
-    def test_signal_family(self):
-        f_main = 0.1
-        f_beat = 0.01
-        signal_family = benchmark.BeatSignal(f_main, f_beat)
-
-        duration = 100.0
-        fs = 10.0
-        amplitude = (0.5, 0.5, 0.5, *np.radians((5.0, 5.0, 5.0)))
-        mean = (0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
-        phase = tuple(np.radians((0.0, 22.5, 45.0, 67.5, 90.0, 112.5)))
-        t, pos, vel, euler, acc, gyro = benchmark._benchmark_helper(
-            duration, amplitude, mean, phase, signal_family, fs
-        )
-
-        assert len(t) == int(duration * fs)
-        assert pos.shape == (len(t), 3)
-        assert vel.shape == (len(t), 3)
-        assert euler.shape == (len(t), 3)
-        assert acc.shape == (len(t), 3)
-        assert gyro.shape == (len(t), 3)
-
-        np.testing.assert_allclose(
-            pos[:, 0], amplitude[0] * signal_family(t, phase[0], phase_degrees=False)[0]
-        )
-        np.testing.assert_allclose(
-            pos[:, 1], amplitude[1] * signal_family(t, phase[1], phase_degrees=False)[0]
-        )
-        np.testing.assert_allclose(
-            pos[:, 2], amplitude[2] * signal_family(t, phase[2], phase_degrees=False)[0]
-        )
-        np.testing.assert_allclose(
-            euler[:, 0],
-            amplitude[3] * signal_family(t, phase[3], phase_degrees=False)[0],
-        )
-        np.testing.assert_allclose(
-            euler[:, 1],
-            amplitude[4] * signal_family(t, phase[4], phase_degrees=False)[0],
-        )
-        np.testing.assert_allclose(
-            euler[:, 2],
-            amplitude[5] * signal_family(t, phase[5], phase_degrees=False)[0],
-        )
-
-    def test_pure_roll(self):
-        f_main = 0.1
-        f_beat = 0.01
-        signal_family = benchmark.BeatSignal(f_main, f_beat)
-
-        duration = 100.0
-        fs = 10.0
-        amplitude = (0.0, 0.0, 0.0, 0.1, 0.0, 0.0)
-        mean = (0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
-        phase = (0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
-        t, pos, vel, euler, acc, gyro = benchmark._benchmark_helper(
-            duration, amplitude, mean, phase, signal_family, fs
-        )
-
-        assert len(t) == int(duration * fs)
-        assert pos.shape == (len(t), 3)
-        assert vel.shape == (len(t), 3)
-        assert euler.shape == (len(t), 3)
-        assert acc.shape == (len(t), 3)
-        assert gyro.shape == (len(t), 3)
-
-        np.testing.assert_allclose(acc[:, 0], np.zeros_like(t))
-        np.testing.assert_allclose(acc[:, 1], -9.80665 * np.sin(euler[:, 0]))
-        np.testing.assert_allclose(acc[:, 2], -9.80665 * np.cos(euler[:, 0]))
-
-        np.testing.assert_allclose(gyro[:, 0], np.gradient(euler[:, 0], t), atol=0.003)
-        np.testing.assert_allclose(gyro[:, 1], np.zeros_like(t))
-        np.testing.assert_allclose(gyro[:, 2], np.zeros_like(t))
-
-    def test_pure_pitch(self):
-        f_main = 0.1
-        f_beat = 0.01
-        signal_family = benchmark.BeatSignal(f_main, f_beat)
-
-        duration = 100.0
-        fs = 10.0
-        amplitude = (0.0, 0.0, 0.0, 0.0, 0.1, 0.0)
-        mean = (0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
-        phase = (0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
-        t, pos, vel, euler, acc, gyro = benchmark._benchmark_helper(
-            duration, amplitude, mean, phase, signal_family, fs
-        )
-
-        assert len(t) == int(duration * fs)
-        assert pos.shape == (len(t), 3)
-        assert vel.shape == (len(t), 3)
-        assert euler.shape == (len(t), 3)
-        assert acc.shape == (len(t), 3)
-        assert gyro.shape == (len(t), 3)
-
-        np.testing.assert_allclose(acc[:, 0], 9.80665 * np.sin(euler[:, 1]))
-        np.testing.assert_allclose(acc[:, 1], np.zeros_like(t))
-        np.testing.assert_allclose(acc[:, 2], -9.80665 * np.cos(euler[:, 1]))
-
-        np.testing.assert_allclose(gyro[:, 0], np.zeros_like(t))
-        np.testing.assert_allclose(gyro[:, 1], np.gradient(euler[:, 1], t), atol=0.003)
-        np.testing.assert_allclose(gyro[:, 2], np.zeros_like(t))
-
-    def test_pure_yaw(self):
-        f_main = 0.1
-        f_beat = 0.01
-        signal_family = benchmark.BeatSignal(f_main, f_beat)
-
-        duration = 100.0
-        fs = 10.0
-        amplitude = (0.0, 0.0, 0.0, 0.0, 0.0, 0.1)
-        mean = (0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
-        phase = (0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
-        t, pos, vel, euler, acc, gyro = benchmark._benchmark_helper(
-            duration, amplitude, mean, phase, signal_family, fs
-        )
-
-        assert len(t) == int(duration * fs)
-        assert pos.shape == (len(t), 3)
-        assert vel.shape == (len(t), 3)
-        assert euler.shape == (len(t), 3)
-        assert acc.shape == (len(t), 3)
-        assert gyro.shape == (len(t), 3)
-
-        np.testing.assert_allclose(acc[:, 0], np.zeros_like(t))
-        np.testing.assert_allclose(acc[:, 1], np.zeros_like(t))
-        np.testing.assert_allclose(acc[:, 2], -9.80665 * np.ones_like(t))
-
-        np.testing.assert_allclose(gyro[:, 0], np.zeros_like(t))
-        np.testing.assert_allclose(gyro[:, 1], np.zeros_like(t))
-        np.testing.assert_allclose(gyro[:, 2], np.gradient(euler[:, 2], t), atol=0.003)
-
-    def test_attitude_mean(self):
-        f_main = 0.1
-        f_beat = 0.01
-        signal_family = benchmark.BeatSignal(f_main, f_beat)
-
-        duration = 100.0
-        fs = 10.0
-        amplitude = (0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
-        mean = (0.0, 0.0, 0.0, 0.1, 0.2, 0.3)
-        phase = (0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
-        t, pos, vel, euler, acc, gyro = benchmark._benchmark_helper(
-            duration, amplitude, mean, phase, signal_family, fs
-        )
-
-        assert len(t) == int(duration * fs)
-        assert pos.shape == (len(t), 3)
-        assert vel.shape == (len(t), 3)
-        assert euler.shape == (len(t), 3)
-        assert acc.shape == (len(t), 3)
-        assert gyro.shape == (len(t), 3)
-
-        np.testing.assert_allclose(acc[:, 0], -9.80665 * -np.sin(0.2))
-        np.testing.assert_allclose(acc[:, 1], -9.80665 * np.sin(0.1) * np.cos(0.2))
-        np.testing.assert_allclose(acc[:, 2], -9.80665 * np.cos(0.1) * np.cos(0.2))
-
-        np.testing.assert_allclose(
-            np.sqrt(acc[:, 0] ** 2 + acc[:, 1] ** 2 + acc[:, 2] ** 2),
-            9.80665 * np.ones_like(t),
-        )
-
-        np.testing.assert_allclose(gyro[:, 0], np.zeros_like(t))
-        np.testing.assert_allclose(gyro[:, 1], np.zeros_like(t))
-        np.testing.assert_allclose(gyro[:, 2], np.zeros_like(t))
-
-    def test_pure_translation(self):
-        f_main = 0.1
-        f_beat = 0.01
-        signal_family = benchmark.BeatSignal(f_main, f_beat)
-
-        duration = 100.0
-        fs = 10.0
-        amplitude = (1.0, 2.0, 3.0, 0.0, 0.0, 0.0)
-        mean = (-1.0, -2.0, -3.0, 0.0, 0.0, 0.0)
-        phase = (0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
-        t, pos, vel, euler, acc, gyro = benchmark._benchmark_helper(
-            duration, amplitude, mean, phase, signal_family, fs
-        )
-
-        assert len(t) == int(duration * fs)
-        assert pos.shape == (len(t), 3)
-        assert vel.shape == (len(t), 3)
-        assert euler.shape == (len(t), 3)
-        assert acc.shape == (len(t), 3)
-        assert gyro.shape == (len(t), 3)
-
-        np.testing.assert_allclose(pos[:, 0].mean(), -1.0, atol=0.1)
-        np.testing.assert_allclose(pos[:, 1].mean(), -2.0, atol=0.1)
-        np.testing.assert_allclose(pos[:, 2].mean(), -3.0, atol=0.1)
-
-        np.testing.assert_allclose(vel[:, 0], np.gradient(pos[:, 0], t), atol=0.005)
-        np.testing.assert_allclose(vel[:, 1], np.gradient(pos[:, 1], t), atol=0.005)
-        np.testing.assert_allclose(vel[:, 2], np.gradient(pos[:, 2], t), atol=0.005)
-
-        np.testing.assert_allclose(acc[:, 0], np.gradient(vel[:, 0], t), atol=0.01)
-        np.testing.assert_allclose(acc[:, 1], np.gradient(vel[:, 1], t), atol=0.01)
-        np.testing.assert_allclose(
-            acc[:, 2], -9.80665 + np.gradient(vel[:, 2], t), atol=0.01
-        )
-
-        np.testing.assert_allclose(gyro[:, 0], np.zeros_like(t))
-        np.testing.assert_allclose(gyro[:, 1], np.zeros_like(t))
-        np.testing.assert_allclose(gyro[:, 2], np.zeros_like(t))
-
-
 def test_benchmark_pure_attitude_beat_202311A():
     signature_signal, _, _ = benchmark.BeatSignal(0.1, 0.01)(
         np.arange(0.0, 1800.0, 1.0 / 10.24)
@@ -312,7 +107,7 @@ def test_benchmark_pure_attitude_beat_202311A():
     assert acc.shape == (len(t), 3)
     assert gyro.shape == (len(t), 3)
 
-    np.testing.assert_array_equal(euler[:, 0], np.radians(5.0) * signature_signal)
+    np.testing.assert_allclose(euler[:, 0], np.radians(5.0) * signature_signal)
 
 
 def test_benchmark_pure_attitude_chirp_202311A():
diff --git a/tests/test_simulate.py b/tests/test_simulate.py
new file mode 100644
index 00000000..4e25ab56
--- /dev/null
+++ b/tests/test_simulate.py
@@ -0,0 +1,607 @@
+import numpy as np
+import pytest
+
+import smsfusion as sf
+from smsfusion.simulate import BeatDOF, ChirpDOF, ConstantDOF, IMUSimulator, SineDOF
+from smsfusion.simulate._simulate import DOF
+
+
+@pytest.fixture
+def t():
+    return np.linspace(0, 10, 100)
+
+
+class Test_DOF:
+    @pytest.fixture
+    def some_dof(self):
+        class SomeDOF(DOF):
+
+            def _y(self, t):
+                return np.ones_like(t)
+
+            def _dydt(self, t):
+                return 2 * np.ones_like(t)
+
+            def _d2ydt2(self, t):
+                return 3 * np.ones_like(t)
+
+        return SomeDOF()
+
+    def test_y(self, some_dof, t):
+        y = some_dof.y(t)
+        np.testing.assert_allclose(y, np.ones(100))
+
+    def test_dydt(self, some_dof, t):
+        dydt = some_dof.dydt(t)
+        np.testing.assert_allclose(dydt, 2 * np.ones(100))
+
+    def test_d2ydt2(self, some_dof, t):
+        d2ydt2 = some_dof.d2ydt2(t)
+        np.testing.assert_allclose(d2ydt2, 3 * np.ones(100))
+
+    def test__call__(self, some_dof, t):
+        y, dydt, dy2dt2 = some_dof(t)
+        np.testing.assert_allclose(y, np.ones(100))
+        np.testing.assert_allclose(dydt, 2 * np.ones(100))
+        np.testing.assert_allclose(dy2dt2, 3 * np.ones(100))
+
+
+class Test_ConstantDOF:
+    @pytest.fixture
+    def constant_dof(self):
+        return ConstantDOF(value=5.0)
+
+    def test__init__(self):
+        constant_dof = ConstantDOF(value=123.0)
+        assert isinstance(constant_dof, DOF)
+        assert constant_dof._value == 123.0
+
+    def test_y(self, constant_dof, t):
+        y = constant_dof.y(t)
+        np.testing.assert_allclose(y, 5.0 * np.ones(100))
+
+    def test_dydt(self, constant_dof, t):
+        dydt = constant_dof.dydt(t)
+        np.testing.assert_allclose(dydt, np.zeros(100))
+
+    def test_d2ydt2(self, constant_dof, t):
+        d2ydt2 = constant_dof.d2ydt2(t)
+        np.testing.assert_allclose(d2ydt2, np.zeros(100))
+
+    def test__call__(self, constant_dof, t):
+        y, dydt, dy2dt2 = constant_dof(t)
+        np.testing.assert_allclose(y, 5.0 * np.ones(100))
+        np.testing.assert_allclose(dydt, np.zeros(100))
+        np.testing.assert_allclose(dy2dt2, np.zeros(100))
+
+
+class Test_SineDOF:
+    @pytest.fixture
+    def sine_dof(self):
+        return SineDOF(2.0, 1.0)
+
+    def test__init__(self):
+        sine_dof = SineDOF(
+            amp=2.0, freq=3.0, freq_hz=True, phase=4.0, phase_degrees=True, offset=5.0
+        )
+
+        assert isinstance(sine_dof, DOF)
+        assert sine_dof._amp == 2.0
+        assert sine_dof._w == pytest.approx(2.0 * np.pi * 3.0)
+        assert sine_dof._phase == pytest.approx((np.pi / 180.0) * 4.0)
+        assert sine_dof._offset == 5.0
+
+    def test_y(self, sine_dof, t):
+        y = sine_dof.y(t)
+        expected_y = 2.0 * np.sin(1.0 * t + 0.0)
+        np.testing.assert_allclose(y, expected_y)
+
+    def test_dydt(self, sine_dof, t):
+        dydt = sine_dof.dydt(t)
+        expected_dydt = 2.0 * 1.0 * np.cos(1.0 * t + 0.0)
+        np.testing.assert_allclose(dydt, expected_dydt)
+
+    def test_d2ydt2(self, sine_dof, t):
+        d2ydt2 = sine_dof.d2ydt2(t)
+        expected_d2ydt2 = -2.0 * (1.0**2) * np.sin(1.0 * t + 0.0)
+        np.testing.assert_allclose(d2ydt2, expected_d2ydt2)
+
+    def test__call__(self, sine_dof, t):
+        y, dydt, dy2dt2 = sine_dof(t)
+        expected_y = 2.0 * np.sin(1.0 * t + 0.0)
+        expected_dydt = 2.0 * 1.0 * np.cos(1.0 * t + 0.0)
+        expected_d2ydt2 = -2.0 * (1.0**2) * np.sin(1.0 * t + 0.0)
+        np.testing.assert_allclose(y, expected_y)
+        np.testing.assert_allclose(dydt, expected_dydt)
+        np.testing.assert_allclose(dy2dt2, expected_d2ydt2)
+
+    def test_amp(self):
+        sine_dof = SineDOF(amp=3.0)
+        t = np.linspace(0, 2 * np.pi, 100)
+        y, dydt, dy2dt2 = sine_dof(t)
+        y_expect = 3.0 * np.sin(t)
+        dydt_expect = 3.0 * np.cos(t)
+        dy2dt2_expect = -3.0 * np.sin(t)
+        np.testing.assert_allclose(y, y_expect)
+        np.testing.assert_allclose(dydt, dydt_expect)
+        np.testing.assert_allclose(dy2dt2, dy2dt2_expect)
+
+    def test_freq_hz(self, t):
+        sine_dof = SineDOF(freq=0.5, freq_hz=True)
+        y, dydt, dy2dt2 = sine_dof(t)
+        y_expect = np.sin(np.pi * t)
+        dydt_expect = np.pi * np.cos(np.pi * t)
+        dy2dt2_expect = -np.pi**2 * np.sin(np.pi * t)
+        np.testing.assert_allclose(y, y_expect)
+        np.testing.assert_allclose(dydt, dydt_expect)
+        np.testing.assert_allclose(dy2dt2, dy2dt2_expect)
+
+    def test_freq_rads(self, t):
+        sine_dof = SineDOF(freq=np.pi, freq_hz=False)
+        y, dydt, dy2dt2 = sine_dof(t)
+        y_expect = np.sin(np.pi * t)
+        dydt_expect = np.pi * np.cos(np.pi * t)
+        dy2dt2_expect = -np.pi**2 * np.sin(np.pi * t)
+        np.testing.assert_allclose(y, y_expect)
+        np.testing.assert_allclose(dydt, dydt_expect)
+        np.testing.assert_allclose(dy2dt2, dy2dt2_expect)
+
+    def test_phase_degrees(self, t):
+        sine_dof = SineDOF(phase=90.0, phase_degrees=True)
+        y, dydt, dy2dt2 = sine_dof(t)
+        y_expect = np.sin(t + np.pi / 2)
+        dydt_expect = np.cos(t + np.pi / 2)
+        dy2dt2_expect = -np.sin(t + np.pi / 2)
+        np.testing.assert_allclose(y, y_expect)
+        np.testing.assert_allclose(dydt, dydt_expect)
+        np.testing.assert_allclose(dy2dt2, dy2dt2_expect)
+
+    def test_phase_radians(self, t):
+        sine_dof = SineDOF(phase=np.pi / 2, phase_degrees=False)
+        y, dydt, dy2dt2 = sine_dof(t)
+        y_expect = np.sin(t + np.pi / 2)
+        dydt_expect = np.cos(t + np.pi / 2)
+        dy2dt2_expect = -np.sin(t + np.pi / 2)
+        np.testing.assert_allclose(y, y_expect)
+        np.testing.assert_allclose(dydt, dydt_expect)
+        np.testing.assert_allclose(dy2dt2, dy2dt2_expect)
+
+    def test_offset(self, t):
+        sine_dof = SineDOF(offset=2.0)
+        y, dydt, dy2dt2 = sine_dof(t)
+        y_expect = 2.0 + np.sin(t)
+        dydt_expect = np.cos(t)
+        dy2dt2_expect = -np.sin(t)
+        np.testing.assert_allclose(y, y_expect)
+        np.testing.assert_allclose(dydt, dydt_expect)
+        np.testing.assert_allclose(dy2dt2, dy2dt2_expect)
+
+
+class Test_IMUSimulator:
+
+    @pytest.fixture
+    def sim(self):
+        pos_x = SineDOF(1.0, 1.0)
+        pos_y = SineDOF(2.0, 0.5)
+        pos_z = SineDOF(3.0, 0.1)
+        alpha = SineDOF(4.0, 1.0)
+        beta = SineDOF(5.0, 0.5)
+        gamma = SineDOF(6.0, 0.1)
+        sim = IMUSimulator(pos_x, pos_y, pos_z, alpha, beta, gamma, degrees=True)
+        return sim
+
+    def test__init__default(self):
+        sim = IMUSimulator()
+        assert isinstance(sim._pos_x, ConstantDOF)
+        assert isinstance(sim._pos_y, ConstantDOF)
+        assert isinstance(sim._pos_z, ConstantDOF)
+        assert isinstance(sim._alpha, ConstantDOF)
+        assert isinstance(sim._beta, ConstantDOF)
+        assert isinstance(sim._gamma, ConstantDOF)
+        assert sim._pos_x._value == 0.0
+        assert sim._pos_y._value == 0.0
+        assert sim._pos_z._value == 0.0
+        assert sim._alpha._value == 0.0
+        assert sim._beta._value == 0.0
+        assert sim._gamma._value == 0.0
+        assert sim._degrees is False
+        assert sim._nav_frame == "ned"
+        np.testing.assert_allclose(sim._g_n, np.array([0.0, 0.0, 9.80665]))
+
+    def test__init__float(self):
+        sim = IMUSimulator(1.0, 2.0, 3.0, 4.0, 5.0, 6.0)
+        assert isinstance(sim._pos_x, ConstantDOF)
+        assert isinstance(sim._pos_y, ConstantDOF)
+        assert isinstance(sim._pos_z, ConstantDOF)
+        assert isinstance(sim._alpha, ConstantDOF)
+        assert isinstance(sim._beta, ConstantDOF)
+        assert isinstance(sim._gamma, ConstantDOF)
+        assert sim._pos_x._value == 1.0
+        assert sim._pos_y._value == 2.0
+        assert sim._pos_z._value == 3.0
+        assert sim._alpha._value == 4.0
+        assert sim._beta._value == 5.0
+        assert sim._gamma._value == 6.0
+
+    def test__init__dof(self):
+        pos_x = SineDOF(1.0, 1.0)
+        pos_y = ConstantDOF(2.0)
+        pos_z = SineDOF(0.5, 0.5)
+        alpha = ConstantDOF(10.0)
+        beta = SineDOF(5.0, 2.0)
+        gamma = ConstantDOF(-5.0)
+
+        sim = IMUSimulator(
+            pos_x=pos_x,
+            pos_y=pos_y,
+            pos_z=pos_z,
+            alpha=alpha,
+            beta=beta,
+            gamma=gamma,
+            degrees=True,
+            g=9.84,
+            nav_frame="ENU",
+        )
+
+        assert sim._pos_x is pos_x
+        assert sim._pos_y is pos_y
+        assert sim._pos_z is pos_z
+        assert sim._alpha is alpha
+        assert sim._beta is beta
+        assert sim._gamma is gamma
+        assert sim._degrees is True
+        assert sim._nav_frame == "enu"
+        np.testing.assert_allclose(sim._g_n, np.array([0.0, 0.0, -9.84]))
+
+    def test__call__default(self, sim):
+        fs = 10.24
+        n = 100
+        t, pos, vel, euler, f, w = sim(fs, n)
+
+        np.testing.assert_allclose(t, np.arange(n) / fs)
+
+        # Position
+        assert pos.shape == (n, 3)
+        np.testing.assert_allclose(pos[:, 0], sim._pos_x.y(t))
+        np.testing.assert_allclose(pos[:, 1], sim._pos_y.y(t))
+        np.testing.assert_allclose(pos[:, 2], sim._pos_z.y(t))
+
+        # Velocity
+        assert vel.shape == (n, 3)
+        np.testing.assert_allclose(vel[:, 0], sim._pos_x.dydt(t))
+        np.testing.assert_allclose(vel[:, 1], sim._pos_y.dydt(t))
+        np.testing.assert_allclose(vel[:, 2], sim._pos_z.dydt(t))
+
+        # Euler angles
+        assert euler.shape == (n, 3)
+        np.testing.assert_allclose(euler[:, 0], sim._alpha.y(t))
+        np.testing.assert_allclose(euler[:, 1], sim._beta.y(t))
+        np.testing.assert_allclose(euler[:, 2], sim._gamma.y(t))
+
+        # Specific force
+        assert f.shape == (n, 3)
+        acc_x = sim._pos_x.d2ydt2(t)
+        acc_y = sim._pos_y.d2ydt2(t)
+        acc_z = sim._pos_z.d2ydt2(t)
+        acc_expect = np.column_stack((acc_x, acc_y, acc_z))
+        for f_i, euler_i, acc_i in zip(f, euler, acc_expect):
+            R_nb_i = sf._transforms._rot_matrix_from_euler(np.radians(euler_i))
+            f_i_expect = R_nb_i.T.dot(acc_i - sim._g_n)
+            np.testing.assert_allclose(f_i, f_i_expect)
+
+        # Angular rate
+        assert w.shape == (n, 3)
+        alpha, beta = np.radians(euler[:, 0:2]).T
+        alpha_dot = sim._alpha.dydt(t)
+        beta_dot = sim._beta.dydt(t)
+        gamma_dot = sim._gamma.dydt(t)
+        w_x = alpha_dot - np.sin(beta) * gamma_dot
+        w_y = np.cos(alpha) * beta_dot + np.sin(alpha) * np.cos(beta) * gamma_dot
+        w_z = -np.sin(alpha) * beta_dot + np.cos(alpha) * np.cos(beta) * gamma_dot
+        w_b = np.column_stack([w_x, w_y, w_z])
+        np.testing.assert_allclose(w, w_b)
+
+    def test__call__degrees(self, sim):
+        fs = 10.24
+        n = 100
+        t, pos, vel, euler, f, w = sim(fs, n, degrees=True)
+
+        np.testing.assert_allclose(t, np.arange(n) / fs)
+
+        # Position
+        assert pos.shape == (n, 3)
+        np.testing.assert_allclose(pos[:, 0], sim._pos_x.y(t))
+        np.testing.assert_allclose(pos[:, 1], sim._pos_y.y(t))
+        np.testing.assert_allclose(pos[:, 2], sim._pos_z.y(t))
+
+        # Velocity
+        assert vel.shape == (n, 3)
+        np.testing.assert_allclose(vel[:, 0], sim._pos_x.dydt(t))
+        np.testing.assert_allclose(vel[:, 1], sim._pos_y.dydt(t))
+        np.testing.assert_allclose(vel[:, 2], sim._pos_z.dydt(t))
+
+        # Euler angles
+        assert euler.shape == (n, 3)
+        np.testing.assert_allclose(euler[:, 0], sim._alpha.y(t))
+        np.testing.assert_allclose(euler[:, 1], sim._beta.y(t))
+        np.testing.assert_allclose(euler[:, 2], sim._gamma.y(t))
+
+        # Specific force
+        assert f.shape == (n, 3)
+        acc_x = sim._pos_x.d2ydt2(t)
+        acc_y = sim._pos_y.d2ydt2(t)
+        acc_z = sim._pos_z.d2ydt2(t)
+        acc_expect = np.column_stack((acc_x, acc_y, acc_z))
+        for f_i, euler_i, acc_i in zip(f, euler, acc_expect):
+            R_nb_i = sf._transforms._rot_matrix_from_euler(np.radians(euler_i))
+            f_i_expect = R_nb_i.T.dot(acc_i - sim._g_n)
+            np.testing.assert_allclose(f_i, f_i_expect)
+
+        # Angular rate
+        assert w.shape == (n, 3)
+        alpha, beta = np.radians(euler[:, 0:2]).T
+        alpha_dot = sim._alpha.dydt(t)
+        beta_dot = sim._beta.dydt(t)
+        gamma_dot = sim._gamma.dydt(t)
+        w_x = alpha_dot - np.sin(beta) * gamma_dot
+        w_y = np.cos(alpha) * beta_dot + np.sin(alpha) * np.cos(beta) * gamma_dot
+        w_z = -np.sin(alpha) * beta_dot + np.cos(alpha) * np.cos(beta) * gamma_dot
+        w_b = np.column_stack([w_x, w_y, w_z])
+        np.testing.assert_allclose(w, w_b)
+
+    def test__call__radians(self, sim):
+        fs = 10.24
+        n = 100
+        t, pos, vel, euler, f, w = sim(fs, n, degrees=False)  # radians
+
+        np.testing.assert_allclose(t, np.arange(n) / fs)
+
+        # Position
+        assert pos.shape == (n, 3)
+        np.testing.assert_allclose(pos[:, 0], sim._pos_x.y(t))
+        np.testing.assert_allclose(pos[:, 1], sim._pos_y.y(t))
+        np.testing.assert_allclose(pos[:, 2], sim._pos_z.y(t))
+
+        # Velocity
+        assert vel.shape == (n, 3)
+        np.testing.assert_allclose(vel[:, 0], sim._pos_x.dydt(t))
+        np.testing.assert_allclose(vel[:, 1], sim._pos_y.dydt(t))
+        np.testing.assert_allclose(vel[:, 2], sim._pos_z.dydt(t))
+
+        # Euler angles
+        assert euler.shape == (n, 3)
+        np.testing.assert_allclose(euler[:, 0], np.radians(sim._alpha.y(t)))
+        np.testing.assert_allclose(euler[:, 1], np.radians(sim._beta.y(t)))
+        np.testing.assert_allclose(euler[:, 2], np.radians(sim._gamma.y(t)))
+
+        # Specific force
+        assert f.shape == (n, 3)
+        acc_x = sim._pos_x.d2ydt2(t)
+        acc_y = sim._pos_y.d2ydt2(t)
+        acc_z = sim._pos_z.d2ydt2(t)
+        acc_expect = np.column_stack((acc_x, acc_y, acc_z))
+        for f_i, euler_i, acc_i in zip(f, euler, acc_expect):
+            R_nb_i = sf._transforms._rot_matrix_from_euler(euler_i)
+            f_i_expect = R_nb_i.T.dot(acc_i - sim._g_n)
+            np.testing.assert_allclose(f_i, f_i_expect)
+
+        # Angular rate
+        assert w.shape == (n, 3)
+        alpha, beta = euler[:, 0:2].T
+        alpha_dot = np.radians(sim._alpha.dydt(t))
+        beta_dot = np.radians(sim._beta.dydt(t))
+        gamma_dot = np.radians(sim._gamma.dydt(t))
+        w_x = alpha_dot - np.sin(beta) * gamma_dot
+        w_y = np.cos(alpha) * beta_dot + np.sin(alpha) * np.cos(beta) * gamma_dot
+        w_z = -np.sin(alpha) * beta_dot + np.cos(alpha) * np.cos(beta) * gamma_dot
+        w_b = np.column_stack([w_x, w_y, w_z])
+        np.testing.assert_allclose(w, w_b)
+
+
+class Test_BeatDOF:
+    @pytest.fixture
+    def beat(self):
+        dof = BeatDOF(amp=2.0, freq_main=1.0, freq_beat=0.1, freq_hz=False, offset=1.0)
+        return dof
+
+    def test__init__(self):
+        beat = BeatDOF(
+            amp=3.0,
+            freq_main=2.0,
+            freq_beat=0.2,
+            freq_hz=True,
+            phase=4.0,
+            phase_degrees=True,
+            offset=5.0,
+        )
+
+        assert isinstance(beat, DOF)
+        assert beat._amp == 3.0
+        assert beat._w_main == pytest.approx(2.0 * np.pi * 2.0)
+        assert beat._w_beat == pytest.approx(2.0 * np.pi * 0.2)
+        assert beat._phase == pytest.approx((np.pi / 180.0) * 4.0)
+        assert beat._offset == 5.0
+
+    def test__init__default(self):
+        beat_dof = BeatDOF()
+
+        assert isinstance(beat_dof, DOF)
+        assert beat_dof._amp == 1.0
+        assert beat_dof._w_main == pytest.approx(0.1)
+        assert beat_dof._w_beat == pytest.approx(0.01)
+        assert beat_dof._phase == pytest.approx(0.0)
+        assert beat_dof._offset == 0.0
+
+    def test_y(self, beat, t):
+        y = beat.y(t)
+
+        amp = beat._amp
+        w_main = beat._w_main
+        w_beat = beat._w_beat
+        phase = beat._phase
+        offset = beat._offset
+
+        main = np.cos(w_main * t + phase)
+        beat = np.sin(w_beat / 2.0 * t)
+
+        y_expect = amp * beat * main + offset
+
+        np.testing.assert_allclose(y, y_expect)
+
+    def test_dydt(self, beat, t):
+        dydt = beat.dydt(t)
+
+        amp = beat._amp
+        w_main = beat._w_main
+        w_beat = beat._w_beat
+        phase = beat._phase
+
+        main = np.cos(w_main * t + phase)
+        beat = np.sin(w_beat / 2.0 * t)
+        dmain = -w_main * np.sin(w_main * t + phase)
+        dbeat = (w_beat / 2.0) * np.cos(w_beat / 2.0 * t)
+
+        dydt_expect = amp * (dbeat * main + beat * dmain)
+
+        np.testing.assert_allclose(dydt, dydt_expect)
+
+    def test_d2ydt2(self, beat, t):
+        d2ydt2 = beat.d2ydt2(t)
+
+        amp = beat._amp
+        w_main = beat._w_main
+        w_beat = beat._w_beat
+        phase = beat._phase
+
+        main = np.cos(w_main * t + phase)
+        beat = np.sin(w_beat / 2.0 * t)
+        dmain = -w_main * np.sin(w_main * t + phase)
+        dbeat = (w_beat / 2.0) * np.cos(w_beat / 2.0 * t)
+        d2main = -(w_main**2) * np.cos(w_main * t + phase)
+        d2beat = -(w_beat**2 / 4.0) * np.sin(w_beat / 2.0 * t)
+
+        d2ydt2_expect = amp * (d2beat * main + 2.0 * dbeat * dmain + beat * d2main)
+
+        np.testing.assert_allclose(d2ydt2, d2ydt2_expect)
+
+    def test__call__(self, beat, t):
+        y, dydt, d2ydt2 = beat(t)
+
+        amp = beat._amp
+        w_main = beat._w_main
+        w_beat = beat._w_beat
+        phase = beat._phase
+        offset = beat._offset
+
+        main = np.cos(w_main * t + phase)
+        beat = np.sin(w_beat / 2.0 * t)
+        dmain = -w_main * np.sin(w_main * t + phase)
+        dbeat = (w_beat / 2.0) * np.cos(w_beat / 2.0 * t)
+        d2main = -(w_main**2) * np.cos(w_main * t + phase)
+        d2beat = -(w_beat**2 / 4.0) * np.sin(w_beat / 2.0 * t)
+
+        y_expect = amp * beat * main + offset
+        dydt_expect = amp * (dbeat * main + beat * dmain)
+        d2ydt2_expect = amp * (d2beat * main + 2.0 * dbeat * dmain + beat * d2main)
+
+        np.testing.assert_allclose(y, y_expect)
+        np.testing.assert_allclose(dydt, dydt_expect)
+        np.testing.assert_allclose(d2ydt2, d2ydt2_expect)
+
+
+class Test_ChirpDOF:
+    @pytest.fixture
+    def chirp(self):
+        chirp = ChirpDOF(amp=2.0, f_max=0.25, f_os=0.01, freq_hz=True, offset=1.0)
+        return chirp
+
+    def test__init__(self):
+        chirp_dof = ChirpDOF(
+            3.0, 2.0, 1.0, freq_hz=True, phase=4.0, phase_degrees=True, offset=5.0
+        )
+
+        assert isinstance(chirp_dof, DOF)
+        assert chirp_dof._amp == 3.0
+        assert chirp_dof._w_max == pytest.approx(2.0 * np.pi * 2.0)
+        assert chirp_dof._w_os == pytest.approx(2.0 * np.pi * 1.0)
+        assert chirp_dof._phase == pytest.approx((np.pi / 180.0) * 4.0)
+        assert chirp_dof._offset == 5.0
+
+    def test__init__default(self):
+        chirp_dof = ChirpDOF()
+
+        assert isinstance(chirp_dof, DOF)
+        assert chirp_dof._amp == 1.0
+        assert chirp_dof._w_max == pytest.approx(np.pi / 2.0)
+        assert chirp_dof._w_os == pytest.approx(2.0 * np.pi * 0.01)
+        assert chirp_dof._phase == pytest.approx(0.0)
+        assert chirp_dof._offset == 0.0
+
+    def test_y(self, chirp, t):
+        y = chirp.y(t)
+
+        amp = chirp._amp
+        w_max = chirp._w_max
+        w_os = chirp._w_os
+        phase = chirp._phase
+        offset = chirp._offset
+
+        phi = 2.0 * w_max / w_os * np.sin(w_os / 2.0 * t)
+        y_expect = amp * np.sin(phi + phase) + offset
+
+        np.testing.assert_allclose(y, y_expect)
+
+    def test_dydt(self, chirp, t):
+        dydt = chirp.dydt(t)
+
+        amp = chirp._amp
+        w_max = chirp._w_max
+        w_os = chirp._w_os
+        phase = chirp._phase
+
+        phi = 2.0 * w_max / w_os * np.sin(w_os / 2.0 * t)
+        dphi_dt = w_max * np.cos(w_os / 2.0 * t)
+
+        dydt_expect = amp * dphi_dt * np.cos(phi + phase)
+
+        np.testing.assert_allclose(dydt, dydt_expect)
+
+    def test_d2ydt2(self, chirp, t):
+        d2ydt2 = chirp.d2ydt2(t)
+
+        amp = chirp._amp
+        w_max = chirp._w_max
+        w_os = chirp._w_os
+        phase = chirp._phase
+
+        phi = 2.0 * w_max / w_os * np.sin(w_os / 2.0 * t)
+        dphi = w_max * np.cos(w_os / 2.0 * t)
+        d2phi = -w_max * w_os / 2.0 * np.sin(w_os / 2.0 * t)
+        d2ydt2_expect = -amp * (dphi**2) * np.sin(phi + phase) + amp * d2phi * np.cos(
+            phi + phase
+        )
+
+        np.testing.assert_allclose(d2ydt2, d2ydt2_expect)
+
+    def test__call__(self, chirp, t):
+        y, dydt, d2ydt2 = chirp(t)
+
+        amp = chirp._amp
+        w_max = chirp._w_max
+        w_os = chirp._w_os
+        phase = chirp._phase
+        offset = chirp._offset
+
+        phi = 2.0 * w_max / w_os * np.sin(w_os / 2.0 * t)
+        dphi = w_max * np.cos(w_os / 2.0 * t)
+        d2phi = -w_max * w_os / 2.0 * np.sin(w_os / 2.0 * t)
+
+        y_expect = amp * np.sin(phi + phase) + offset
+        dydt_expect = amp * dphi * np.cos(phi + phase)
+        d2ydt2_expect = -amp * (dphi**2) * np.sin(phi + phase) + amp * d2phi * np.cos(
+            phi + phase
+        )
+
+        np.testing.assert_allclose(y, y_expect)
+        np.testing.assert_allclose(dydt, dydt_expect)
+        np.testing.assert_allclose(d2ydt2, d2ydt2_expect)