最近在学习计算IMU的allan方差时查找到一些资料，现在分享出来：

最先看到的参考文档是2004年的国标《光纤陀螺测试方法》里面的公式，对其中的公式多少有些费解，后面又找来了其它参考文献：analysis-and-modeling-of-inertial-sensors-using-allan-variance，最后反复查资料确认，国标上的公式是有问题的。

1.关于allan公式对比：

感兴趣的同学可以在附件上面下载两篇文章进行对比。

下图是国标给出的公式：

下图是参考论文里面的公式：

通过上面的公式对比，你会发现两份文档里面的公式定义是不相同的。通过查找其他参考资料核实，实际上国标的定义不太对。

考虑到在计算时如果以角度为观测量，那么ALLAN方差的计算公式如下所示：

2.双对数曲线

在计算出allan方差的时候我们需要画出双对数曲线，关于双LOG曲线可以理解如下：

   对数坐标的刻度位置是以10^x为等距标定的（所谓等距刻度就是类似普通坐标系中两1,2,3的位置），如图横坐标最左端的10^1和最右端的10^2间的距离为一个等距刻度。下一个等距刻度是10^3的位置，而对数坐标系中等距刻度上的标注值为真数，实际上就是将普通坐标系中的（1,2,3....）位置不变而改标为(10^1,10^2,10^3...)，然后数据点按照自己的值标在图中对应位置，就可以了。坐标系改变，不会改变数据点，只会改变曲线的样子。

3.源码参考：

关于参考源码，最先在篇博客上看到了源码，博客链接：https://blog.csdn.net/u012325601/article/details/60882949

该博客所对应的ALLAN方差的matlab源码地址为：https://github.com/XinLiGH/GyroAllan

allan方差计算的源码如下所示：

function [T,sigma] = allan(omega,fs,pts)

[N,M] = size(omega);             % figure out how big the output data set is

n = 2.^(0:floor(log2(N/2)))';    % determine largest bin size

maxN = n(end);

endLogInc = log10(maxN);

m = unique(ceil(logspace(0,endLogInc,pts)))';    % create log spaced vector average factor

t0 = 1/fs;                                       % t0 = sample interval

T = m*t0;                                        % T = length of time for each cluster

theta = cumsum(omega)/fs;       % integration of samples over time to obtain output angle θ

sigma2 = zeros(length(T),M);    % array of dimensions (cluster periods) X (#variables)

for i=1:length(m)               % loop over the various cluster sizes

    for k=1:N-2*m(i)            % implements the summation in the AV equation

        sigma2(i,:) = sigma2(i,:) + (theta(k+2*m(i),:) - 2*theta(k+m(i),:) + theta(k,:)).^2;

    end

end

sigma2 = sigma2./repmat((2*T.^2.*(N-2*m)),1,M);

sigma = sqrt(sigma2);

关于allan方差的拟合如下所示：

function C=nihe(tau,sig,M)

X=tau';Y=sig';

B=zeros(1,2*M+1);

F=zeros(length(X),2*M+1);

for i=1:2*M+1

    kk=i-M-1;

    F(:,i)=X.^kk;

end

A=F'*F;

B=F'*Y;

C=A\B;

主函数如下所示：

clc;    %清空命令行窗口

clear;  %清空工作区

data = dlmread('data.dat');         %从文本中读取数据，单位：deg/s，速率：100Hz

data = data(1:720000, 3:5)*3600;    %截取两个小时的数据，把 deg/s 转为 deg/h

[A, B] = allan(data, 100, 100);     %求Allan标准差，用100个点来描述

loglog(A, B, 'o');                  %画双对数坐标图

xlabel('time:sec');                 %添加x轴标签

ylabel('Sigma:deg/h');              %添加y轴标签

legend('X axis','Y axis','Z axis'); %添加标注

grid on;                            %添加网格线

hold on;                            %使图像不被覆盖

C(1, :) = nihe(sqrt(A'), B(:, 1)', 2)';   %拟合

C(2, :) = nihe(sqrt(A'), B(:, 2)', 2)';

C(3, :) = nihe(sqrt(A'), B(:, 3)', 2)';

Q = abs(C(:, 1)) / sqrt(3);         %量化噪声，单位：arcsec

N = abs(C(:, 2)) / 60;            %角度随机游走，单位：deg/h^0.5

Bs = abs(C(:, 3)) / 0.6643;        %零偏不稳定性，单位：deg/h

K = abs(C(:, 4)) * sqrt(3) * 60;%角速率游走，单位：deg/h/h^0.5

R = abs(C(:, 5)) * sqrt(2) * 3600;%速率斜坡，单位：deg/h/h

fprintf('量化噪声      X轴：%f Y轴：%f Z轴：%f  单位：arcsec\n', Q(1), Q(2), Q(3));

fprintf('角度随机游走  X轴：%f Y轴：%f Z轴：%f  单位：deg/h^0.5\n', N(1), N(2), N(3));

fprintf('零偏不稳定性  X轴：%f Y轴：%f Z轴：%f  单位：deg/h\n', Bs(1), Bs(2), Bs(3));

fprintf('角速率游走    X轴：%f Y轴：%f Z轴：%f  单位：deg/h/h^0.5\n', K(1), K(2), K(3));

fprintf('速率斜坡      X轴：%f Y轴：%f Z轴：%f  单位：deg/h/h\n', R(1), R(2), R(3));

D(:, 1) = C(1, 1)*sqrt(A).^(-2) + C(1, 2)*sqrt(A).^(-1) + C(1, 3)*sqrt(A).^(0) + C(1, 4)*sqrt(A).^(1) + C(1, 5)*sqrt(A).^(2);    %生成拟合函数

D(:, 2) = C(2, 1)*sqrt(A).^(-2) + C(2, 2)*sqrt(A).^(-1) + C(2, 3)*sqrt(A).^(0) + C(2, 4)*sqrt(A).^(1) + C(2, 5)*sqrt(A).^(2);

D(:, 3) = C(3, 1)*sqrt(A).^(-2) + C(3, 2)*sqrt(A).^(-1) + C(3, 3)*sqrt(A).^(0) + C(3, 4)*sqrt(A).^(1) + C(3, 5)*sqrt(A).^(2);

loglog(A, D);   %画双对数坐标图

4.利用gyo数据算出的Allan方差图片如下：

5.参考文献下载地址如下，如有侵权行为，请联系作者删除附件，liangwang_nx@163.com

链接：https://pan.baidu.com/s/1TTzZywJb1tdaVWuI-XMclw 密码：l7p4

关于allan方差拟合后的系数与零偏不稳定性以及角速度随机游走等在下一篇讨论。