Why doesn't the noise density of the Allan deviation plot match the noise density of imuSensor?

My test is as follows:

1. Define a imuSensor simulating only `NoiseDensity`

2. Generate gyro readings for a stationary IMU

3. Plot the Allan deviation using `adev`

4. Read the noise density as the deviation at tau=1

Unfortunately, the noise density specified in step 1 does not match the plot value read in step 4. The code below implements the above test for a specified noise density 0.1 but the Allan deviation indicates a noise density of 0.07. What's going on?

fs = 100;
totalNumSamples = 1E5;

traj = kinematicTrajectory('SampleRate',fs);

% Groundtruth accel and gyro 

accBody = zeros(totalNumSamples,3);
angVelBody = zeros(totalNumSamples,3);

% Groundtruth trajectory in NED 

[~,orientationNED,~,accNED,angVelNED] = traj(accBody,angVelBody);

% Define IMU model and its parameters 

IMU = imuSensor('accel-gyro','SampleRate',fs);

IMU.Gyroscope = gyroparams('NoiseDensity', 0.1);    
   
% Simulate the gyro readings

[~, gyroReadings] = IMU(accNED,angVelNED,orientationNED);
gyroReadings = gyroReadings(:, 1);   % only take the x-axis 

% Allan plot as instructed in 
% https://uk.mathworks.com/help/nav/ug/inertial-sensor-noise-analysis-using-allan-variance.html

t0 = 1/fs;
theta = cumsum(gyroReadings, 1)*t0;

maxNumM = 100;
L = size(theta, 1);
maxM = 2.^floor(log2(L/2));
m = logspace(log10(1), log10(maxM), maxNumM).';
m = ceil(m); % m must be an integer.
m = unique(m); % Remove duplicates.

tau = m*t0;

avar = zeros(numel(m), 1);
for i = 1:numel(m)
    mi = m(i);
    avar(i,:) = sum( ...
        (theta(1+2*mi:L) - 2*theta(1+mi:L-mi) + theta(1:L-2*mi)).^2, 1);
end
avar = avar ./ (2*tau.^2 .* (L - 2*m));

[avarFromFunc, tauFromFunc] = allanvar(gyroReadings, m, fs);
adevFromFunc = sqrt(avarFromFunc);
adev = sqrt(avar);

% plot it 
figure
loglog(tau, adev, tauFromFunc, adevFromFunc);
title('Allan Deviations')
xlabel('\tau')
ylabel('\sigma(\tau)')
legend('Manual Calculation', 'allanvar Function')
grid on

Answer 1

[Paul]

Hi Tsz Hey Lam,

The example in the link was for angle random walk (ARW), but the example in the question is using white noise.

rng(100);
fs = 100;
totalNumSamples = 1E5;
traj = kinematicTrajectory('SampleRate',fs);
% Groundtruth accel and gyro 
accBody = zeros(totalNumSamples,3);
angVelBody = zeros(totalNumSamples,3);
% Groundtruth trajectory in NED 
[~,orientationNED,~,accNED,angVelNED] = traj(accBody,angVelBody);
% Define IMU model and its parameters 
IMU = imuSensor('accel-gyro','SampleRate',fs);
IMU.Gyroscope = gyroparams('NoiseDensity', 0.1);    
   
% Simulate the gyro readings
[~, gyroReadings] = IMU(accNED,angVelNED,orientationNED);
gyroReadings = gyroReadings(:, 1);   % only take the x-axis 
% Allan plot as instructed in 
% https://uk.mathworks.com/help/nav/ug/inertial-sensor-noise-analysis-using-allan-variance.html
t0 = 1/fs;
theta = cumsum(gyroReadings, 1)*t0;
maxNumM = 100;
L = size(theta, 1);
maxM = 2.^floor(log2(L/2));
m = logspace(log10(1), log10(maxM), maxNumM).';
m = ceil(m); % m must be an integer.
m = unique(m); % Remove duplicates.
tau = m*t0;
avar = zeros(numel(m), 1);
for i = 1:numel(m)
    mi = m(i);
    avar(i,:) = sum( ...
        (theta(1+2*mi:L) - 2*theta(1+mi:L-mi) + theta(1:L-2*mi)).^2, 1);
end
avar = avar ./ (2*tau.^2 .* (L - 2*m));
[avarFromFunc, tauFromFunc] = allanvar(gyroReadings, m, fs);
adevFromFunc = sqrt(avarFromFunc);
adev = sqrt(avar);
% plot it 
figure
loglog(tau, adev, tauFromFunc, adevFromFunc);
title('Allan Deviations')
xlabel('\tau')
ylabel('\sigma(\tau)')
legend('Manual Calculation', 'allanvar Function')
grid on

According to the doc page imusensor, the white noise model uses standard white noise scaled by

IMU.Gyroscope.NoiseDensity(1)*sqrt(fs/2)

ans = 0.7071

which should be the standard deviation of our gyro readings

std(gyroReadings)

ans = 0.7083

When tau = 1, I think that the Allan deviation should be the standard devation of the gyroReadings scaled by sqrt(fs). However, I am not 100% certain that's the case, so check into what the actual relationship should be.

hold on
yline(std(gyroReadings)/sqrt(fs))
xline(1)

Answer 2

[Tsz Hey Lam]

This MATLAB article uses ARW and noise density interchangeably. The article refers to the parameter N as "N (angle random walk)" and states that "The value of N can be read directly off of this line at "; subsequently using `N` as `NoiseDensity` for gyroparams. We believe that this is incorrect and that NoiseDensity should be read as the value multiplied by .

We have tested this using imuSensor and allanvar for a range of NoiseDensity values and can reliably recover those values as the value multiplied by .

The article concludes with a plot and caption stating, "The model measurements contain slightly less noise since the quantization and temperature-related parameters are not set using gyroparams". We believe this explanation is mistaken and the reduced noise is accounted for by the failure to scale by .

Our plot below shows hardware gyroscope data plotted alongside simulated gyroscope data. The simulated data was generated using the White Noise Drift model. The hardware value is 0.007. This plot confirms that accurate simulation requires that noise density is read as the value multiplied by .

We are demonstrating a significant error in the MATLAB article. Please can someone from MATLAB respond?

[Tsz Hey Lam]

Update - I have not had any reponses to my previous comment but the MATLAB article has now been updated to correct for the error. The plot at the end of the article now shows the hardware data and simulated data overlapping.

This is because gyroparams now includes NoiseType specified as single-sided. Doing so changes the expression in the noise model from to, thereby addressing the scaling described in my previous comment.