动态Allan方差，MATLAB程序

% DAVAR Dynamic Allan Variance -- BETA VERSION!! % S = DAVAR(x,t,N,tau0,flag) computes and represents the Dynamic Allan Variance. % The function actually computes the Dynamic Allan Deviation. % The signal x is truncated to N samples at every time instant specified % by the vector t (t is the center of the window N, that must be % odd), and the square root of the Dynamic Allan Variance is % computed. When truncating at the beginning or at the end of % the signal a zero padding may be performed (the program generates a warning % whenever a zero padding is performed). % % INPUT % x --> Signal to be analyzed. x must be a column vector. % t --> Vector of time indexes at which the DAVAR is computed % (it must follow the Matlab syntax for indexes). % N --> Length of the window. N must be odd. % tau0 --> Time step on the tau axis. % flag --> If flag=0 no output is produced, otherwise a mesh plot % of the Dynamic Allan Deviation is generated. % % OUTPUT % S --> A Ntau-by-Nt matrix representing the Dynamic Allan Deviation, % where Ntau=N/3 (rounded toward zero) and Nt is the length of t. % % Example 1: Compute the DAVAR of a white noise random process. % x=randn(1000,1); % S=DAVAR(x,[1:10:1000],201,1,1); % % Example 2: Compute the DAVAR of a nonstationary white noise process. % The noise has a variance that changes with time. % x=[randn(400,1); 2*randn(200,1); randn(400,1)]; % S=DAVAR(x,[51:1:950],101,1,1); % Avoid zero padding. % Lorenzo Galleani -- July 15, 2005. % Politecnico di Torino, Torino, Italy. % This program has been developed in collaboration with % Patrizia Tavella, IEN Galileo Ferraris, Torino, Italy. function [S] = DAVAR(x,t,N,tau0,flag) Nx=length(x); Nt=length(t); if ~mod(N,2) error('N must be odd'); end L1=t(1)-1-(N-1)/2; if L1<0 x=[zeros(-L1,1);x]; disp('Warning: zero padding at the beginning of x'); else L1=0; end L2=t(Nt)+(N-1)/2-Nx; if L2>0 x=[x;zeros(L2,1)]; disp('Warning: zero padding at the end of x'); else L2=0; end Ntau=floor(N/3); S=zeros(Ntau,Nt); for n=1:Nt xN=x(-L1+t(n)-(N-1)/2:-L1+t(n)+(N-1)/2); S(:,n)=allan(xN,tau0); end if flag & Nt>1 tau=tau0*[1:Ntau]; figure,mesh(t,tau,S); set(gca,'YDir','reverse','YScale','log','ZScale','log'); axis('tight'); xlabel('t'),ylabel('\tau'),zlabel('\sigma(t,\tau)'); hidden off; view(53,32); end