Simulating sensor readings based on an existing Alan Variance Analysis

Question

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
def get_y_intersect(slope, tau, adev):
    logtau = np.log10(tau)
    logadev = np.log10(adev)
    dlogadev = np.diff(logadev) / np.diff(logtau)
    i = np.argmin(np.abs(dlogadev - slope))

    return logadev[i] - slope*logtau[i], i

# 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)
adev = np.sqrt(avar)

## plotting everything
fig, ax = plt.subplots(1, 1, figsize=(16,9))
ax.loglog(tau, adev, "-")
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?

Answer 1

From the same documentation you linked,

The three noise parameters N (angle random walk), K (rate random walk), and B (bias instability) are estimated using data logged from a stationary gyroscope.

You shouldn't need to simulate these; they're typically terms provided to you on datasheets or, if you have a physical device, you can calculate them from the test device.