# Simulating sensor readings based on an existing Alan Variance Analysis

My goal is to model an accelerometer and a gyroscope based on real hardware. I understand deterministic errors like bias and scaling but the different types of non-deterministic errors are difficult.

After a bit of research, I found a tutorial provided in the matlab documentation.

• $$\sigma^2(\tau) = \frac{1}{2\tau^2(n-2m)}\sum_{k=1}^{n-2m}(\theta_{k+2m} - 2\theta_{k+m} + \theta_{k})^2$$
• $$\sigma(t) = \sqrt{\sigma^2(t)}$$

Since I have no access to matlab, I implemented the function using python:

def calculate_avar(theta, t0, max_num_m):
n = theta.size # number of samples
max_m = 2**int(np.log2(n/2)) # maximum cluster size
m = np.logspace(np.log10(1), np.log10(max_m), max_num_m) # cluster sizes
m = np.ceil(m).astype(int) # m must be an integer.
m = np.unique(m) # Remove duplicates.

tau = m*t0
result = np.empty_like(m)
for i in range(m.size):
result[i] = np.sum((theta[2*m[i]:n] - 2*theta[m[i]:n-m[i]] + theta[:n-2*m[i]])**2)
result = result / (2*tau**2 * (n - 2*m))

return tau, result


Using this function I calculated:

• Noise density/Angle Random Walk/Velocity random walk (slope = $$-\frac{1}{2}$$)
• Bias (in)stability (slope = 0)
• (Rate) Random walk (slope = $$+\frac{1}{2}$$)
# find y intersection given a particular slope
logtau = np.log10(tau)

# noise density N
slope = -0.5
b, _ = get_y_intersect(slope, tau, adev)
logN = slope*np.log10(1) + b
N = 10**logN

# rate random walk R
slope = 0.5
b, _ = get_y_intersect(slope, tau, adev)
logK = slope*np.log10(3) + b
K = 10**logK

# bias (in)stability B
slope = 0
b, i = get_y_intersect(slope, tau, adev)
scfB = np.sqrt(2*np.log(2)/np.pi)
logB = b - np.log10(scfB)
B = 10**logB


The problem is that I don't have access to Matlab and don't understand how to simulate the resulting parameters N, K and B. I've seen various names for these parameters (pink noise, white noise Brownian noise, second-order random walk, etc...) which is even more confusing to me. Everything I've found so far was either too shallow or requires in-depth knowledge about signal processing.

fs = 100 # sample rate
max_num_m = 1000 # number of clusters
n = 12*60*60*fs # number of samples (12hours of 100Hz samples)
t0 = 1/fs # sample time
theta = np.zeros(n, dtype=np.float64) # ideal samples

##### I can't figure out this part #####
theta += np.random.normal(0, N, n) # noise density !?
theta += np.random.normal(0, B, n).cumsum() # bias (in)stability !?
# theta += ??? # rate random walk
########################################

tau, avar = calculate_avar(theta, t0, max_num_m)

## plotting everything
fig, ax = plt.subplots(1, 1, figsize=(16,9))
ax.grid(True)
ax.set_title("Allan Deviation")
ax.axis("equal")
ax.set_xlabel("$$\\tau$$")
ax.set_ylabel("$$\\sigma(\\tau)$$")
fig.show()


I'm not sure if the last bit of code is correct (probably not) but I'm kind of stuck here... can someone help?

• You could get SciLab, Octave, or other MatLab clones. Jun 17 '20 at 22:21