Анализ шума инерционного датчика с использованием Отклонения Allan

Открыть скрипт

Этот пример показывает, как использовать отклонение Аллана для определения шумовых параметров гироскопа MEMS. Эти параметры могут использоваться, чтобы смоделировать гироскоп в симуляции. Измерение гироскопа моделируется как:

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

Три параметра шума N (случайная прогулка под углом), K (случайное блуждание по скорости) и B (нестабильность смещения) оцениваются с помощью данных, записанных со стационарного гироскопа.

Фон

Аллан отклонением был первоначально разработан Дэвидом У. Алланом, чтобы измерить стабильность частоты прецизионных генераторов. Он также может использоваться, чтобы идентифицировать различные источники шума, присутствующие в стационарных измерениях гироскопа. Рассмотрим L выборок данных из гироскопа со шаг расчета. $\tau_{0}$Формируйте кластеры данных длительности,,$\tau_{0}$..$2\tau_{0}$. и$m\tau_{0}, (m < (L-1)/2)$ получайте средние значения суммы точек данных, содержащихся в каждом кластере, по длине кластера. Отклонение Allan определяется как двухвыборочное отклонение средних значений кластера данных как функция времени кластера. Этот пример использует перекрывающуюся оценку отклонения Аллана. Это означает, что вычисленные кластеры перекрываются. Оценка работает лучше, чем непересекающиеся оценки для больших значений L.

Расчет отклонений Аллана

Отклонение Allan вычисляется следующим образом:

Логгирование выборок стационарного гироскопа L с периодом дискретизации. $\tau_{0}$Позвольте$\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;

Затем вычислите отклонение Аллана:

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

где$\tau = m\tau_{0}$ и$<&#62;$ является ли ансамбль среднего значения.

Средний ансамбль может быть расширен до:

$$\sigma^2(\tau) =&#xA;\frac{1}{2\tau^2(L-2m)}\sum_{k=1}^{L-2m}(\theta_{k+2m} - 2\theta_{k+m}&#xA;+ \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));

Наконец, отклонение Аллана$\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 функция.

[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

Идентификация параметра шума

Чтобы получить параметры шума для гироскопа, используйте следующее соотношение между отклонением Аллана и двусторонней спектральной плотностью степени (PSD) параметров шума в исходном наборе данных. $\Omega$Отношения следующие:

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

Из приведенного выше уравнения отклонение Аллана пропорционально общей степени гироскопа при прохождении через фильтр с передаточной функцией от. $sin^4(x)/(x)^2$Эта передаточная функция возникает из-за операций, выполненных для создания и работы с кластерами.

Используя эту интерпретацию передаточной функции, полосу пропускания фильтра зависит от. $\tau$Это означает, что различные параметры шума могут быть идентифицированы путем изменения полосы пропускания фильтра или изменения.$\tau$

Случайная прогулка по углу

Случайная прогулка под углом характеризуется белым шумовым спектром выходного сигнала гироскопа. PSD представлен:

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

где

N = угол случайного коэффициента обхода

Подстановка в исходное уравнение PSD и выполнение интегрирования приводит к:

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

Приведенное выше уравнение является линией с наклоном -1/2 при построении на логарифмическом графике$\sigma(\tau)$ от. $\tau$Значение N может считываться непосредственно вне этой линии.$\tau = 1$

% 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}$$

Приведенное выше уравнение является линией с наклоном 1/2 при построении на логарифмическом графике$\sigma(\tau)$ от. $\tau$Значение K может считываться непосредственно вне этой линии.$\tau = 3$

% 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}&#xA;&#38; : f \leq f_0\\0 &#38; : f &#62; f_0\end{array}\right.$$

где

B = коэффициент нестабильности смещения

$f_0$ = частота отключения

Подстановка в исходное уравнение PSD и выполнение интегрирования приводит к:

$$\sigma^2(\tau) = \frac{2B^2}{\pi}[\ln{2} + \\&#xA;-\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}$$

Приведенное выше уравнение является линией с наклоном 0 при построении на логарифмическом графике$\sigma(\tau)$ от. $\tau$Значение B может считываться непосредственно вне этой линии с масштабированием.$\sqrt{\frac{2\ln{2}}{\pi}} \approx 0.664$

% 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

Теперь, когда все параметры шума были вычислены, постройте график отклонения Аллана со всеми линиями, используемыми для количественной оценки параметров.

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', '\sigma_N', '\sigma_K', '\sigma_B')
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.
    gyro = gyroparams('NoiseDensity', N, 'RandomWalk', K, ...
        'BiasInstability', B);
    omegaSim = helperAllanVarianceExample(L, Fs, gyro);
else
    load('SimulatedSingleAxisGyroscope', 'omegaSim')
end

Вычислите моделируемое отклонение Аллана и сравните его с записанными данными.

[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 генерирует измерения с аналогичным отклонением Аллана к записанным данным. Измерения модели содержат немного меньше шума, поскольку параметры квантования и температуры не заданы с помощью gyroparams. Модель гироскопа может использоваться, чтобы сгенерировать измерения с помощью движений, которые не легко захватываются с помощью оборудования.

Ссылки

  • IEEE Std. 647-2006 Руководство по формату стандартных спецификаций IEEE и процедура тестирования для одноосных лазерных гироскопов

Документация по Sensor Fusion and Tracking Toolbox

Поддержка

Для просмотра документации необходимо авторизоваться на сайте