Difference between gazebo and rotors IMU noise model

[ghost]

Hi! I read in the GaussianNoiseModel::ApplyImpl(double _in, double _dt) function that it's implementation is based on the one available in Rotors. By comparing both of them, I notice that a "-1" had been removed in the gazebo model's calculation of sigma_b_d :

• gazebo : double sigma_b_d = sqrt(-sigma_b * sigma_b * tau / 2 * expm1(-2 * _dt / tau));
• rotors : double sigma_b_a_d = sqrt(-sigma_b_a * sigma_b_a * tau_a / 2.0 * (exp(-2.0 * dt / tau_a) - 1.0));

Is there a particular reason behind this change? I'd be interested to know.

Additionally, the gazebo noise model adds an extra sign randomization step to the bias in :

// With equal probability, we pick a negative bias (by convention,
  // rateBiasMean should be positive, though it would work fine if
  // negative).
  if (math::Rand::DblUniform() < 0.5)
    this->dataPtr->bias = -this->dataPtr->bias;

Which I found no equivalent for in the rotors code and would be interested in knowing the motivation for.

Lastly, it is said that "More information about the parameters and their derivation can be found at https://github.com/ethz-asl/kalibr/wiki/IMU-Noise-Model". However, I find no mention of the dynamicBiasStdDev (sigma_b) and dynamicBiasCorrTime (tau) parameters, or how to obtain them. Could you share which IMU calibration method you used to obtain these? I found some values in the files for the gazebo-fuel X4 drone model but they seem quite platform-specific.

For context, I'm trying to match my gazebo simulation to tests I'm running outdoors with an IMU calibrated by Kalibr, to see how my filter performs in sim vs in real life. This is why I'm trying to understand the differences between both implementations, and decide whether to write a gazebo plugin to ensure they match.

Thank you!

[mjcarroll]

Is there a particular reason behind this change? I'd be interested to know.

For this change, we make use of std::expm1

Computes the e (Euler's number, 2.7182818...) raised to the given power num, minus 1.0. This function is more accurate than the expression std::exp(num) - 1.0 if num is close to zero.

Additionally, the gazebo noise model adds an extra sign randomization step to the bias in :

This predates the inclusion of the rotors noise model, but my assumptions is that because the bias is drawn from a gaussian distribution specified by mean/stddev, we want to make sure that it can additionally be drawn as a positive or negative bias.

In the case that bias mean was 1, bias stddev is 0.1, there is also the chance that the bias may be centered around -1, which is handled by the sign flip.

[ghost]

Is there a particular reason behind this change? I'd be interested to know.

For this change, we make use of std::expm1

Computes the e (Euler's number, 2.7182818...) raised to the given power num, minus 1.0. This function is more accurate than the expression std::exp(num) - 1.0 if num is close to zero.

Additionally, the gazebo noise model adds an extra sign randomization step to the bias in :

This predates the inclusion of the rotors noise model, but my assumptions is that because the bias is drawn from a gaussian distribution specified by mean/stddev, we want to make sure that it can additionally be drawn as a positive or negative bias.

In the case that bias mean was 1, bias stddev is 0.1, there is also the chance that the bias may be centered around -1, which is handled by the sign flip.

Ah I see that makes alot more sense, thank you for your lightning-fast reply!

[mjcarroll]

Lastly, it is said that "More information about the parameters and their derivation can be found at https://github.com/ethz-asl/kalibr/wiki/IMU-Noise-Model". However, I find no mention of the dynamicBiasStdDev (sigma_b) and dynamicBiasCorrTime (tau) parameters

Those are also copied from rotors:
https://github.com/ethz-asl/rotors_simulator/blob/513bb92da0c1a0c968bdc679dffc8fe7d77de918/rotors_gazebo_plugins/include/rotors_gazebo_plugins/gazebo_imu_plugin.h#L39-L55

computing tau can be done using the allan variance method. I have a script that does it from gazebo data here: https://github.com/gazebosim/gz-sensors/blob/gz-sensors7/examples/imu_noise/plot_samples.py and there is another ROS package that does the same thing here: https://github.com/ori-drs/allan_variance_ros

It is the x (tau) value of the minimum of the allan variance plot, like this one from the kalibr wiki page:

I believe the reason you don't see it in the kalibr parameters is because while we (and rotors) simulate it, kalibr does not make use of it.

[mjcarroll]

Ah I see that makes alot more sense, thank you for your lightning-fast reply!

Happy to help, let me know if you uncover any more oddities or concerns.

[ghost]

Ah I see that makes alot more sense, thank you for your lightning-fast reply!

Happy to help, let me know if you uncover any more oddities or concerns.

Thanks! Is dynamicBiasStdDev (sigma_b) in the gazebo IMU model the same as sigma(tau) In the gazebo script you linked?

[mjcarroll]

I believe you are looking for this:

gz-sensors/examples/imu_noise/plot_samples.py

Line 67 in a5a0c34

(tau_b, b, line_b) = estimate_bias_instability(t2, ad)

For that minimum point in the graph. tau_b is the correlation time (X value). b is the bias instability (Y value).

b is the value you are looking for.

[ghost]

I believe you are looking for this:

gz-sensors/examples/imu_noise/plot_samples.py

Line 67 in a5a0c34

(tau_b, b, line_b) = estimate_bias_instability(t2, ad)

For that minimum point in the graph. tau_b is the correlation time (X value). b is the bias instability (Y value).

b is the value you are looking for.

I see, one last question :

• dynamicBiasStdDev = sigma_b = b = bias instability = y-value at the minima of a which curve(s)?
• dynamicBiasCorrTime = tau = tau_b = correlation time = x-value at the minima of which curve(s)?

Thanks you for your responsivity!

[mjcarroll]

dynamicBiasStdDev = sigma_b = b = bias instability = y-value at the minima of a which curve(s)?
dynamicBiasCorrTime = tau = tau_b = correlation time = x-value at the minima of which curve(s)?

Correct. If you are planning to reproduce this on a real, physical IMU, make sure that you capture a long enough sequence of data (a few hours) to be sure that you can get the full graph.

[ghost]

dynamicBiasStdDev = sigma_b = b = bias instability = y-value at the minima of a which curve(s)?
dynamicBiasCorrTime = tau = tau_b = correlation time = x-value at the minima of which curve(s)?

Correct. If you are planning to reproduce this on a real, physical IMU, make sure that you capture a long enough sequence of data (a few hours) to be sure that you can get the full graph.

Yep will run it overnight, could you clarify which curve I should be taking the minima (x=tau, y=b) of? That wasnt super clear to me :)

[mjcarroll]

Yep will run it overnight, could you clarify which curve I should be taking the minima (x=tau, y=b) of? That wasnt super clear to me

It comes from the allan variance curve. You will have one per channel of your IMU (3x for accel, 3x for gyro). The minimum of each of those curves will be what you are looking for for bias instability.

[ghost]

Yep will run it overnight, could you clarify which curve I should be taking the minima (x=tau, y=b) of? That wasnt super clear to me

It comes from the allan variance curve. You will have one per channel of your IMU (3x for accel, 3x for gyro). The minimum of each of those curves will be what you are looking for for bias instability.

Thank you for confirming that, have a great day!