jumpNoiseAnalysis.m

Go to the documentation of this file.

% =========================================================================== %
%> @brief performs a noise analysis of data of a JUMP sensor.
%>
%> this function performs an Allan Deviation analysis of the data of a JUMP
%> sensor. The data is expected to be in the format of the JUMP sensor data
%> structure defined in @ref jumpReadData.m. The function an interface between
%> the data and the @ref allanDeviationAnalysis.m function.
%>
%> @param data (struct): the data of a JUMP sensor in the format of the JUMP
%> sensor data structure defined in @ref jumpReadData.m.
%> @param options (struct): options for the analysis
%> - toAnalyze (string): the fields of the data to analyze (default: ["acc","gyr","hig","mag","prs","tmp"])
%> - plot (logical): if true, the function will plot the Allan Deviation of the data (default: false)
%>
%> @retval analysis (struct): a struct with the results of the analysis. The
%> struct has the same fields as the input data, and each field has the
%> following subfields:
%> - tau (double): the tau values of the Allan Deviation
%> - adev (double): the Allan Deviation values
%> - N (double): the Power Spectral Density
%> - K (double): the K value of the Allan Deviation
%> - B (double): the B value of the Allan Deviation
%> - tau_B (double): the tau value of the B value
%> - note (string): a note about the units of the results
%>
%> @file jumpNoiseAnalysis.m
%>
%> @copyright see the file @ref LICENSE in the root directory of the repository
% =========================================================================== %
function analysis = jumpNoiseAnalysis(data, options)
 
% input validation
arguments
    data (1,1) struct
    options.toAnalyze (1,:) string  = ["acc","gyr","hig","mag","prs","tmp"];
    options.plot      (1,:) logical = false;
end
 
% check data
assert(jumpCheckData(data), "Data structure is not valid.");
 
% define struct
analysis = struct();
 
if options.plot
    figure("Name", sprintf("Allan Deviation Analysis of Sensor %s", char(data.name)), "NumberTitle", "off");
end
 
% loop over toAnalyze
for i = 1:length(options.toAnalyze)
    % get data
    dataToAnalyze = data.tt.(options.toAnalyze(i));
    
    % get fs
    if strcmp(options.toAnalyze(i), "acc") || strcmp(options.toAnalyze(i), "gyr")
        fs = data.mpu_rate;
    elseif strcmp(options.toAnalyze(i), "mag")
        
        this_idx = all(~isnan(data.tt.mag),2);
        % get fs
        this_times = data.tt.t(this_idx);
        fs = round(1 / seconds(median(diff(this_times))));
        
        if fs >= 20
            % no timing correction was applied, therefore use 10 Hz
            fs = 10;
        end
        
        % downsampe because the data is stored at mpu_rate, but the actual rate is 10 Hz
        dataToAnalyze = resample(dataToAnalyze(this_idx,:), this_times, fs, "pchip");
        
    elseif strcmp(options.toAnalyze(i), "hig")
        fs = data.adxl_rate;
    elseif strcmp(options.toAnalyze(i), "prs") || strcmp(options.toAnalyze(i), "tmp")
        
        this_idx = all(~isnan(data.tt.(options.toAnalyze(i))),2);
        % get fs
        this_times = data.tt.t(this_idx);
        fs = round(1 / seconds(median(diff(this_times))));
        
        if fs >= data.pres_rate * 2
            % no timing correction was applied, therefore use 1 Hz
            fs = data.pres_rate;
        end
        
        % downsampe because the data is stored at mpu_rate, but the actual rate is pres_rate
        dataToAnalyze = resample(dataToAnalyze(this_idx,:), this_times, fs, "pchip");
    end
    
    [tau, adev, N, K, B, tau_B] = allanDeviationAnalysis(dataToAnalyze, fs);
    
    % store in analysis
    analysis.(options.toAnalyze(i)).tau  = tau;
    analysis.(options.toAnalyze(i)).adev = adev;
    analysis.(options.toAnalyze(i)).N    = N;
    analysis.(options.toAnalyze(i)).K    = K;
    analysis.(options.toAnalyze(i)).B    = B;
    analysis.(options.toAnalyze(i)).tau_B = tau_B;
    analysis.(options.toAnalyze(i)).note = "The units are related to the input data. Check units before applying them to a model. Also, the Power Spectral Density (N) is outputtet according to single sided spectrum.";
    
    
    % plot
    if options.plot
        subplot(2,ceil(length(options.toAnalyze)/2),i);
        
        loglog(tau, adev, '-r', 1, N, 'o', 3, K, 'x', tau_B, 0.664*B, '+');
        title(options.toAnalyze(i));
        xlabel('tau [s]');
        ylabel('Allan Deviation');
        grid on;
        
    end
end
 
 
 
end