%Allan方差分析
clear;
clc
load TG400_5.dat;
c=TG400_5/11399*3600;   %轉(zhuǎn)換成角速率，標(biāo)度因數(shù)11300
tt=0.3;       %數(shù)據(jù)輸出周期
NA=length(c);
NN=NA/7;
for M=1:NN;    %
    K=NA/M;
    for i=1:K
       b_aver(i)=mean(c(M*(i-1)+1:M*i)); 
    end
    si=0;
    for j=1:K-1
       si=si+(b_aver(j+1)-b_aver(j))^2;
    end
    sig(M)=si/(K-1)/2;
end
figure(2);
loglog((1:NN)*tt,sig(1:NN));
%grid on;%打開分格線
tau=(1:NN)*tt;
M=2;%階數(shù)
C=nihe(tau,sig,M);
C
Q=sqrt(abs(C(M-1)/3))*(1000000/180/3600);
N=sqrt(abs(C(M)/1))/60; 
Bs=sqrt(abs(C(M+1))/(log(2)*2/pi));
K=sqrt(abs(C(M+2)*3))*60; 
R=sqrt(abs(C(M+3)*2))*3600; 
Q
N
Bs
K
R
X=tau';Y=sig';
YY=zeros(length(X),1);
for i=-M:M
YY=C(M+i+1).*X.^(i)+YY;
end
YY=sqrt(YY);sig=sqrt(sig);
loglog(tau,sig);%grid on;
hold on;loglog(X,YY,'r');xlabel('相關(guān)時間對數(shù)\tau');ylabel('Allan方差對數(shù)\sigma');title('Allan Analysis');
figure(1);
plot(1:NA,c);
xlabel('采樣點');ylabel('輸出信號');title('輸出信號波形');