Main Content

Allan 분산을 사용하여 관성 센서 잡음 분석

스크립트 열기

이 예제에서는 Allan 분산을 사용하여 MEMS 자이로스코프의 잡음 파라미터를 확인하는 방법을 보여줍니다. 이러한 파라미터는 시뮬레이션에서 자이로스코프를 모델링하는 데 사용할 수 있습니다. 자이로스코프 측정은 다음으로 모델링됩니다.

$$\Omega(t) = \Omega_{Ideal}(t) + Bias_N(t) + Bias_B(t) + Bias_K(t)$$

세 개의 잡음 파라미터 N(각도 무작위 행보), K(속도 무작위 행보), B(편향 불안정성)는 정상 상태의 자이로스코프에서 기록된 데이터를 사용하여 추정됩니다.

배경

Allan 분산은 원래 정밀 발진기의 주파수 안정성을 측정하기 위해 David W. Allan이 개발했습니다. 또한 정상 자이로스코프 측정값에 존재하는 다양한 잡음 소스를 식별하는 데에도 사용할 수 있습니다. 샘플 시간이 $\tau_{0}$인 자이로스코프에서 데이터의 L개 샘플을 고려합니다. 지속 시간이 $\tau_{0}$, $2\tau_{0}$, ..., $m\tau_{0}, (m < (L-1)/2)$인 데이터 군집들을 형성하고 각 군집에서 해당 군집 길이 전체에 포함된 데이터 점들의 합의 평균을 구합니다. Allan 분산은 군집 시간의 함수로서 데이터 군집 평균의 2-샘플 분산으로 정의됩니다. 이 예제에서는 겹치는 Allan 분산 추정량을 사용합니다. 즉, 계산되는 군집이 겹쳐 있음을 의미합니다. 이러한 추정량은 L 값이 클 경우 겹치지 않는 추정량보다 더 나은 성능을 보입니다.

Allan 분산 계산

Allan 분산은 다음과 같이 계산됩니다.

샘플 주기 $\tau_{0}$으로 L개의 정상 자이로스코프 샘플을 기록합니다. $\Omega$를 기록된 샘플이라고 가정합니다.

% Load logged data from one axis of a three-axis gyroscope. This recording
% was done over a six hour period with a 100 Hz sampling rate.
load('LoggedSingleAxisGyroscope', 'omega', 'Fs')
t0 = 1/Fs;

각 샘플에 대해, 다음과 같이 출력 각도 $\theta$를 계산합니다.

$$\theta(t) = \int^{t}\Omega(t')dt'$$

이산 샘플의 경우, 누적합에 $\tau_{0}$을 곱한 값을 사용할 수 있습니다.

theta = cumsum(omega, 1)*t0;

다음으로, Allan 분산을 계산합니다.

$$\sigma^2(\tau) =&#10;\frac{1}{2\tau^2}<(\theta_{k+2m}-2\theta_{k+m}+\theta_{k})^2>$$

여기서 $\tau = m\tau_{0}$이고 $<>$은 앙상블 평균입니다.

앙상블 평균은 다음과 같이 전개할 수 있습니다.

$$\sigma^2(\tau) =&#10;\frac{1}{2\tau^2(L-2m)}\sum_{k=1}^{L-2m}(\theta_{k+2m} - 2\theta_{k+m}&#10;+ \theta_{k})^2$$

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));

마지막으로, Allan 편차 $\sigma(t) = \sqrt{\sigma^2(t)}$는 자이로스코프 잡음 파라미터를 확인하는 데 사용됩니다.

adev = sqrt(avar);

figure
loglog(tau, adev)
title('Allan Deviation')
xlabel('\tau');
ylabel('\sigma(\tau)')
grid on
axis equal

allanvar 함수를 사용하여 Allan 분산을 계산할 수도 있습니다.

[avarFromFunc, tauFromFunc] = allanvar(omega, m, Fs);
adevFromFunc = sqrt(avarFromFunc);

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

잡음 파라미터 식별

자이로스코프에 대한 잡음 파라미터를 구하려면, 원래 데이터 세트 $\Omega$에서 잡음 파라미터의 Allan 분산과 양측 PSD(파워 스펙트럼 밀도) 간의 다음 관계를 사용합니다. 관계는 다음과 같습니다.

$$\sigma^2(\tau) = 4\int_{0}^{\infty}S_\Omega(f)&#10;\frac{sin^4(\pi f\tau)}{(\pi f\tau)^2}df$$

위의 방정식에서 Allan 분산은 전달 함수 $sin^4(x)/(x)^2$을 갖는 필터를 통과할 때 자이로스코프의 총 잡음 파워에 비례합니다. 이 전달 함수는 군집을 생성하고 이를 조작하기 위해 수행한 연산에서 도출됩니다.

이 전달 함수 해석을 사용할 때 필터 대역통과는 $\tau$에 따라 달라집니다. 다시 말해 필터 대역통과를 변경하여(즉, $\tau$를 변경하여) 다양한 잡음 파라미터를 식별할 수 있습니다.

각도 무작위 행보

각도 무작위 행보는 자이로스코프 출력값의 백색 잡음 스펙트럼으로 특징지을 수 있습니다. PSD는 다음으로 표현됩니다.

$$S_\Omega(f) = N^2$$

여기서

N = 각도 무작위 행보 계수

원래 PSD 방정식에 대입하고 적분을 수행하면 다음 결과를 얻게 됩니다.

$$\sigma^2(\tau) = \frac{N^2}{\tau}$$

위의 방정식은 $\sigma(\tau)$ 대 $\tau$의 로그-로그 플롯에 플로팅될 때 기울기가 -1/2인 직선입니다. N의 값은 $\tau = 1$일 때 이 선에서 직접 읽어올 수 있습니다. N의 단위는 $(rad/s)/\sqrt{Hz}$입니다.

% Find the index where the slope of the log-scaled Allan deviation is equal
% to the slope specified.
slope = -0.5;
logtau = log10(tau);
logadev = log10(adev);
dlogadev = diff(logadev) ./ diff(logtau);
[~, i] = min(abs(dlogadev - slope));

% Find the y-intercept of the line.
b = logadev(i) - slope*logtau(i);

% Determine the angle random walk coefficient from the line.
logN = slope*log(1) + b;
N = 10^logN

% Plot the results.
tauN = 1;
lineN = N ./ sqrt(tau);
figure
loglog(tau, adev, tau, lineN, '--', tauN, N, 'o')
title('Allan Deviation with Angle Random Walk')
xlabel('\tau')
ylabel('\sigma(\tau)')
legend('\sigma', '\sigma_N')
text(tauN, N, 'N')
grid on
axis equal

N =

    0.0126

속도 무작위 행보

속도 무작위 행보는 자이로스코프 출력값의 적색 잡음(브라운 잡음) 스펙트럼으로 특징지을 수 있습니다. PSD는 다음으로 표현됩니다.

$$S_\Omega(f) = (\frac{K}{2\pi})^2\frac{1}{f^2}$$

여기서

K = 속도 무작위 행보 계수

원래 PSD 방정식에 대입하고 적분을 수행하면 다음 결과를 얻게 됩니다.

$$\sigma^2(\tau) = \frac{K^2\tau}{3}$$

위의 방정식은 $\sigma(\tau)$ 대 $\tau$의 로그-로그 플롯에 플로팅될 때 기울기가 1/2인 직선입니다. K의 값은 $\tau = 3$일 때 이 선에서 직접 읽어올 수 있습니다. K의 단위는 $(rad/s)\sqrt{Hz}$입니다.

% Find the index where the slope of the log-scaled Allan deviation is equal
% to the slope specified.
slope = 0.5;
logtau = log10(tau);
logadev = log10(adev);
dlogadev = diff(logadev) ./ diff(logtau);
[~, i] = min(abs(dlogadev - slope));

% Find the y-intercept of the line.
b = logadev(i) - slope*logtau(i);

% Determine the rate random walk coefficient from the line.
logK = slope*log10(3) + b;
K = 10^logK

% Plot the results.
tauK = 3;
lineK = K .* sqrt(tau/3);
figure
loglog(tau, adev, tau, lineK, '--', tauK, K, 'o')
title('Allan Deviation with Rate Random Walk')
xlabel('\tau')
ylabel('\sigma(\tau)')
legend('\sigma', '\sigma_K')
text(tauK, K, 'K')
grid on
axis equal

K =

   9.0679e-05

편향 불안정성

편향 불안정성은 자이로스코프 출력값의 핑크 잡음(플리커 잡음) 스펙트럼으로 특징지을 수 있습니다. PSD는 다음으로 표현됩니다.

$$S_{\Omega}(f) = \left\{\begin{array}{lr}(\frac{B^2}{2\pi})\frac{1}{f}&#10;&#38; : f \leq f_0\\0 &#38; : f > f_0\end{array}\right.$$

여기서

B = 편향 불안정성 계수

$f_0$ = 차단 주파수

원래 PSD 방정식에 대입하고 적분을 수행하면 다음 결과를 얻게 됩니다.

$$\sigma^2(\tau) = \frac{2B^2}{\pi}[\ln{2} + \\&#10;-\frac{sin^3x}{2x^2}(sinx + 4xcosx) + Ci(2x) - Ci(4x)]$$

여기서

$x = \pi f_0\tau$

Ci = 코사인 적분 함수

$\tau$가 차단 주파수의 역수보다 훨씬 긴 경우 PSD 방정식은 다음과 같습니다.

$$\sigma^2(\tau) = \frac{2B^2}{\pi}\ln{2}$$

위의 방정식은 $\sigma(\tau)$ 대 $\tau$의 로그-로그 플롯에 플로팅될 때 기울기가 0인 직선입니다. B의 값은 스케일링 $\sqrt{\frac{2\ln{2}}{\pi}} \approx 0.664$를 사용해 이 선에서 직접 읽어올 수 있습니다. B의 단위는 $rad/s$입니다.

% Find the index where the slope of the log-scaled Allan deviation is equal
% to the slope specified.
slope = 0;
logtau = log10(tau);
logadev = log10(adev);
dlogadev = diff(logadev) ./ diff(logtau);
[~, i] = min(abs(dlogadev - slope));

% Find the y-intercept of the line.
b = logadev(i) - slope*logtau(i);

% Determine the bias instability coefficient from the line.
scfB = sqrt(2*log(2)/pi);
logB = b - log10(scfB);
B = 10^logB

% Plot the results.
tauB = tau(i);
lineB = B * scfB * ones(size(tau));
figure
loglog(tau, adev, tau, lineB, '--', tauB, scfB*B, 'o')
title('Allan Deviation with Bias Instability')
xlabel('\tau')
ylabel('\sigma(\tau)')
legend('\sigma', '\sigma_B')
text(tauB, scfB*B, '0.664B')
grid on
axis equal

B =

    0.0020

모든 잡음 파라미터를 계산했으므로 파라미터 정량화에 사용된 모든 선으로 Allan 편차를 플로팅합니다.

tauParams = [tauN, tauK, tauB];
params = [N, K, scfB*B];
figure
loglog(tau, adev, tau, [lineN, lineK, lineB], '--', ...
    tauParams, params, 'o')
title('Allan Deviation with Noise Parameters')
xlabel('\tau')
ylabel('\sigma(\tau)')
legend('$\sigma (rad/s)$', '$\sigma_N ((rad/s)/\sqrt{Hz})$', ...
    '$\sigma_K ((rad/s)\sqrt{Hz})$', '$\sigma_B (rad/s)$', 'Interpreter', 'latex')
text(tauParams, params, {'N', 'K', '0.664B'})
grid on
axis equal

자이로스코프 시뮬레이션

imuSensor 객체를 사용하여 위에서 식별한 잡음 파라미터로 자이로스코프 측정값을 시뮬레이션합니다.

% Simulating the gyroscope measurements takes some time. To avoid this, the
% measurements were generated and saved to a MAT-file. By default, this
% example uses the MAT-file. To generate the measurements instead, change
% this logical variable to true.
generateSimulatedData = false;

if generateSimulatedData
    % Set the gyroscope parameters to the noise parameters determined
    % above. Use single-sided noise to match the parameter scaling and 500
    % poles for the bias instability model to match the hardware output.
    numpoles = 500;
    gyro = gyroparams('NoiseDensity', N, 'RandomWalk', K, ...
        'BiasInstability', B, 'NoiseType', 'single-sided', ...
        'BiasInstabilityCoefficients', fractalcoef(numpoles));
    omegaSim = helperAllanVarianceExample(L, Fs, gyro);
else
    load('SimulatedSingleAxisGyroscope', 'omegaSim')
end

시뮬레이션된 Allan 편차를 계산하고 이를 기록된 데이터와 비교합니다.

[avarSim, tauSim] = allanvar(omegaSim, 'octave', Fs);
adevSim = sqrt(avarSim);
adevSim = mean(adevSim, 2); % Use the mean of the simulations.

figure
loglog(tau, adev, tauSim, adevSim, '--')
title('Allan Deviation of HW and Simulation')
xlabel('\tau');
ylabel('\sigma(\tau)')
legend('HW', 'SIM')
grid on
axis equal

플롯은 imuSensor에서 만들어진 자이로스코프 모델이, 기록된 데이터와 유사한 Allan 편차를 가진 측정값을 생성하고 있음을 보여줍니다. 양자화 파라미터와 온도 관련 파라미터가 gyroparams를 사용하여 설정되지 않았기 때문에 모델 측정값은 약간 적은 잡음을 포함합니다. 자이로스코프 모델은 하드웨어로 쉽게 캡처할 수 없는 움직임을 사용해 측정값을 생성하는 데 사용할 수 있습니다.

참고 문헌

  • IEEE Std. 647-2006 IEEE Standard Specification Format Guide and Test Procedure for Single-Axis Laser Gyros