allon方差计算方法

admin

[mw_shl_code=applescript,true]function [Allan,T]=allan_cal(N,dt)
%% 计算Allan方差 %%%%%%%%%%%%%%%%%%%%%
n=length(N);              %%%%%  N:样本（浓度）；n：样本数量
for j=2:length(N)         %%%%%  j:分组数（从2组--n组）
    A(j)=0;               %%%%%  每组的Allan方差初值为0
    for k=1:j           %%%%%  每组平均值的计算
        y(k)=0;           %%%%%  每组平均值初值为0
        for p=1:1:fix(n/j)                           %%%%%  每组的样本数量为 fix(n/j)
            y(k)=y(k)+N(p+fix(n/j)*(k-1));           %%%%%  每组加和的计算
        end
        y(k)=y(k)/fix(n/j);                          %%%%%  每组平均值的计算
        if k>1
            A(j)=A(j)+1/(j-1)/2*(y(k)-y(k-1))^2;     %%%%%  每组Allan方差的计算
        end
    end
end

m=2:length(N);                  %%%%%  分组数
T1=fix(n./m)*dt;                %%%%%  每组的样本数量为fix(n./m)，间隔率为dt，则积分时间为T1；

for j=1:length(A)-1
    Allan1(j)=sqrt(A(j+1));      %%%%%  每个积分时间所对应的Allan方差
end

dT=diff(T1);
ndT=find(dT~=0);
T=T1(ndT);
Allan=Allan1(ndT);

subplot(211)
loglog(T,Allan);                %%%%%  画对数坐标图
subplot(212)
plot(N)
   
%% 验证浓度分布是否满足Gauss分布 %%%%%%%%%%%%%%%%%%%%%
n1=0;  n2=0;  n3=0;  n4=0;  n5=0;
dN=(max(N)-min(N))/5;
for j=1:length(N)
if N(j)<min(N)+dN
    n1=n1+1;
end
if N(j)>=min(N)+dN && N(j)<min(N)+2*dN
    n2=n2+1;
end
if N(j)>=min(N)+2*dN && N(j)<min(N)+3*dN
    n3=n3+1;
end
if N(j)>=min(N)+3*dN && N(j)<min(N)+4*dN
    n4=n4+1;
end
if N(j)>=min(N)+4*dN && N(j)<max(N)
    n5=n5+1;
end
end
nn=[n1 n2 n3 n4 n5];
x=[min(N)+0.5*dN  min(N)+1.5*dN  min(N)+2.5*dN  min(N)+3.5*dN  min(N)+4.5*dN ];
bar(x,nn);
end[/mw_shl_code]