allanDeviationAnalysis.m

Go to the documentation of this file.

% =========================================================================== %
%> @brief allanDeviationAnalysis calculates the Allan deviation and estimates
%> the noise parameters of a given measurement.
%>
%> @param meas - the measurement data (Mxu vector) [unit: any = X, e.g. m/s^2]
%> @param fs - the sampling frequency of the measurement (scalar) [unit: Hz]
%> @param plt - a flag to plot the Allan deviation (logical)
%>
%> @retval tau - the time intervals used for Allan deviation calculation (Mxu vector) [unit: s]
%> @retval adev - the Allan deviation values corresponding to each tau (Mxu vector) [unit: X, e.g. m/s^2]
%> @retval N - the noise density factor (1xu vector) [unit: X/sqrt(Hz), e.g. (m/s^2)/sqrt(Hz)]
%> @retval K - the random walk coefficient (1xu vector) [unit: X * sqrt(Hz), e.g. (m/s^2) * sqrt(Hz)]
%> @retval B - the bias instability coefficient (1xu vector) [unit: X, e.g. m/s^2]
%> @retval tau_B - the time interval where the bias instability coefficient is estimated (1xu vector) [unit: s]
%>
%> The implementation is based on the following references:
%> - https://ch.mathworks.com/help/nav/ug/inertial-sensor-noise-analysis-using-allan-variance.html
%> - IEEE Standard Specification Format Guide and Test Procedure for Single-Axis Laser Gyros, IEEE Std 647-2006
%>
%> @file allanDeviationAnalysis.m
%>
%> @copyright see the file @ref LICENSE in the root directory of the repository
% =========================================================================== %
function [tau, adev, N, K, B, tau_B] = allanDeviationAnalysis(meas, fs)
 
arguments
    meas (:,:) double {mustBeNumeric}
    fs (1,1) double {mustBeNumeric}
end
 
% calculate Allan variance
[avar, tau] = allanvar(meas, 'octave', fs);
 
% remove last element of tau and avar (because we have very long measurements)
tau  = tau(1:end-1,:);
avar = avar(1:end-1,:);
 
% calculate Allan deviation
adev = sqrt(avar);
 
u = width(meas);
N = nan(1,u);
K = nan(1,u);
B = nan(1,u);
tau_B = nan(1,u);
 
% find Noise Parameters
logtau   = log10(tau);
logadev  = log10(adev);
dlogadev = diff(logadev) ./ diff(logtau);
 
for i = 1:u
    
    % first, find the global minimum of the Allan deviation
    [~, idx_min] = min(adev(1:end-1,i));
    
    % find Noise Density N (must be left of the minimum of the Allan deviation curve)
    slope     = -0.5;
    [~, idxN] = min(abs(dlogadev(1:idx_min,i) - slope));
    b         = logadev(idxN,i) - slope * logtau(idxN); % find y-intercept
    logN      = slope * log10(1) + b;
    N(i)      = 10^logN; % Noise Density Factor N [unit: X/sqrt(Hz), e.g. (m/s^2)/sqrt(Hz)]
    
    % find Random Walk Coefficient K (must be right of the minimum of the Allan deviation curve)
    slope     = 0.5;
    [~, idxK] = min(abs(dlogadev(idx_min:end,i) - slope));
    idxK      = idxK + idx_min - 1; % correct index
    b         = logadev(idxK,i) - slope * logtau(idxK); % find y-intercept
    logK      = slope * log10(3) + b;
    K(i)      = 10^logK; % Random Walk Coefficient K [unit: X * sqrt(Hz), e.g. (m/s^2) * sqrt(Hz)]
    
    % find bias instability Coefficient B (must be between idxN and idxK)
    slope     = 0;
    [~, idxB] = min(abs(dlogadev(idxN:idxK,i) - slope));
    idxB      = idxB + idxN - 1; % correct index
    b         = logadev(idxB,i) - slope * logtau(idxB); % find y-intercept
    scfB      = sqrt(2*log(2)/pi);
    logB      = b - log10(scfB);
    B(i)      = 10^logB; % bias instability Coefficient B [unit: X, e.g. m/s^2]
    
    % find tau_B
    tau_B(i) = tau(idxB);
    
end
end