Skip to content
Snippets Groups Projects
Commit 3a853f85 authored by Neumeier, Sebastian's avatar Neumeier, Sebastian
Browse files

Merge branch 'development' into 'main' 

Refactor the code from Manuel Rexer

See merge request sebastian.neumeier/frequency-domain-errors-from-systematic-sensor-errors-package-matlab!1
parents 5c12adcf 9d0ec8cf
Branches main
Tags v0.1.0
1 merge request!1Refactor the code from Manuel Rexer
Show whitespace changes
Inline Side-by-side
Showing
with 453 additions and 186 deletions
.gitignore 0 → 100644
+ 2
− 0
View file @ 3a853f85
# Ignore autogenerated matlab files
**.asv
\ No newline at end of file
LICENSE 0 → 100644
+ 21
− 0
View file @ 3a853f85
MIT License
Copyright (c) 2024 Rexer, Manuel; Neumeier, Sebastian
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
README.md
+ 51
− 10
View file @ 3a853f85
# Frequency Domain Errors From Systematic Sensor Errors package matlab
## Introduction
This is the repository mentioned in the "Propagation of Systematic Sensor Errors into the Frequency Domain - A MATLAB Software Framework" paper
that provides a small matlab package with one main function with which you are able to calculate the error in the Frequency Domain (FFT) by the given
This is the repository mentioned in the "Propagation of Systematic Sensor Errors into the Frequency Domain - A MATLAB Software Framework" paper (see https://dx.doi.org/10.2139/ssrn.4452038 )
that provides a small matlab package with one main function with which you are able to calculate the error in the Frequency Domain (calculated with a FFT) by the given
arguments of one harmonic oscillation input signal with at least one complete cycle in the format of a one dimensional array y(t) and the systematic uncertainty,
which can be obtained from the data sheet or the calibration protocol of the sensor. <br>
......@@ -10,16 +10,57 @@ Please Note: <br>
2. The systematic uncertainty is a collective term that can be calculated from other errors. For example the a) bias error b) sensitivity error, c) linearity
error and d) hysteresis error. <br>
The framework uses a sophisticated sensor error model in combination with Monte Carlo simulations to propagate the
modelled errors into the frequency domain. The framework intends to enable the user to easily determine the statistical errors
and to transfer them together with the given systematic errors into the frequency domain. Therefore, the software framework
The package uses a sophisticated sensor error model in combination with a easy mathematical correlation expression, obtained with Monte Carlo simulations, to propagate the
modelled errors into the frequency domain. The package intends to enable the user to easily determine the statistical errors
and to transfer them together with the given systematic errors into the frequency domain. Therefore, the software package
represents a significant contribution in the realm of uncertainty propagation, especially in the field of system identification. <br>
## Current maintainers:
## Associated Publications
**The associated paper:** <br>
Rexer, Manuel and Pelz, Peter F. and Kuhr, Maximilian M.G., Uncertainty Propagation of Systematic Sensor Errors into the Frequency Domain. Available at SSRN: https://ssrn.com/abstract=4452038 or http://dx.doi.org/10.2139/ssrn.4452038 <br>
**The publication of this software package on Zenodo to obtain a global persistent Identifier (persitent ID):** <br>
TODO:
## Getting started
The only API function this package provides is the `calculateFrequencyDomainERRfromSystematicSensorERR`(...) function
```matlab
result = ...
calculateFrequencyDomainERRfromSystematicSensorERR(...
value_vector, ...
time_vector, ...
sampletime, ...
'bias', 0.01, ...
'slope', 0.02, ...
'lin', 0.03, ...
'hys', 0.04);
```
Where <br>
1. `value_vector` must be (1-BY-N) vector with your data,
2. `time_vector` must be (1-BY-N) vector with the corresponding time values,
3. `sampletime` a scalar (1-BY-1) real number of type double that is the step time between the single values of the time_vector and
4. `'bias'` `'slope'` `'lin'` `'hys'` **optional arguments** for the different systematic uncertainties "bias", "slope", "linearity" and "hysteresis" which can be usually found in the data sheet or the calibration protocol of the sensor you have used to record the data signal.
The "linearity" and "hysteresis" systematic uncertainties are not always given by the datasheets and, if left out in the function call, assumed as 0. <br>
Returns: <br>
`result_struct` that contains the input data, interim results, the absolute systematic uncertainty in the frequency domain and the phase systematic uncertainty for every frequency gathered from the FFT. <br>
**TODO: explaination of the data structure**
## Contributing
If you find errors in this software please open a issue on the corresponding RWTH-GitLab repository.
If you encounter severe problems or nobody answers on your issue (after 2-3 weeks) you can also try to contact the current maintainers through their emails, or better contact the secretariat of the Chair of Fluid Systems (Institut für Fluidsystemtechnik) through
the contact person mentioned on https://www.fst.tu-darmstadt.de/fachgebiet/mitarbeiter/index.en.jsp . <br>
Thank you! <br>
## Current maintainers (Latest state October 2023):
Manuel Rexer, Main Author of the paper. manuel.rexer[at]tu-darmstadt.de <br>
Sebastian Neumeier, Research Aide (german: HiWi) with the responsibility to refactor the code and maintain the repository. sebastian.neumeier[at]stud.tu-darmstadt.de
Sebastian Neumeier, Research Aide (german: SHK) with the responsibility to refactor the code and maintain the repository. sebastian.neumeier[at]stud.tu-darmstadt.de
## Improvement Suggestions:
1. Explaination of the returned data structure.
2. If there should be advances in the use of RDF in calculation packages/software we should come back to the returned data structure and implement it there.
3. Maybe add a real world example with a plot of a data signal as a small 'how-to guide' on how this package could be used correctly.
## TODOs:
TODO: Code example how to use the function.
TODO: Welche Lizenz?
\ No newline at end of file
addStatunc.m deleted 100644 → 0
+ 0
− 6
View file @ 5c12adcf
function data = addStatunc(data,sigma)
%UNTITLED2 Summary of this function goes here
% Detailed explanation goes here
data.unc.stat.sigma=sigma;
end
\ No newline at end of file
addSysunc.m deleted 100644 → 0
+ 0
− 25
View file @ 5c12adcf
function [data]=addSysunc(data,varargin)
%UNTITLED Summary of this function goes here
% Detailed explanation goes here
%% Input adaption
defaultBias=0;
defaultSlope=0;
defaultLin=0;
defaultHys=0;
p = inputParser;
validScalarPosNum = @(x) isnumeric(x) && (x > 0);
addOptional(p,'bias',defaultBias,validScalarPosNum);
addOptional(p,'slope',defaultSlope,validScalarPosNum);
addOptional(p,'lin',defaultLin,validScalarPosNum);
addOptional(p,'hys',defaultHys,validScalarPosNum);
parse(p,varargin{:})
bias = p.Results.bias;
slope = p.Results.slope;
lin = p.Results.lin;
hys = p.Results.hys;
%% generate uncertainty data structure
data.unc.sys.bias = bias;
data.unc.sys.slope = slope;
data.unc.sys.lin = lin;
data.unc.sys.hys = hys;
end
\ No newline at end of file
cropData.m deleted 100644 → 0
+ 0
− 13
View file @ 5c12adcf
function [croptime,cropvalue,nperiods] = cropData(data,excitationfreq,nperiods)
%UNTITLED6 Summary of this function goes here
% Detailed explanation goes here
y = data.value;
t = data.time.value;
nperiods=floor(nperiods);
Ndft = round(( nperiods * (1/excitationfreq) ) / data.time.sampletime);
% Cutting nperiods from timeseries
cropvalue = y(1:Ndft);
croptime= t(1:Ndft);
end
\ No newline at end of file
getMainfreq.m deleted 100644 → 0
+ 0
− 9
View file @ 5c12adcf
function freq = getMainfreq(data);
%UNTITLED3 Summary of this function goes here
% Detailed explanation goes here
y=data.value;
t=data.time.value;
[c2,~] = fit(t',y','fourier1');
freq = c2.w/2/pi;
end
\ No newline at end of file
getNumberofcycles.m deleted 100644 → 0
+ 0
− 12
View file @ 5c12adcf
function numberofcycles = getNumberofcycles(data,f_curr)
%UNTITLED4 Summary of this function goes here
% Detailed explanation goes here
t=data.time.value;
numberofcycles = (t(end) - t(1)) / (1/f_curr);
if numberofcycles<1
error('There is no complete cycle. The data cannot be analyzed')
end
end
\ No newline at end of file
getSampleData.m deleted 100644 → 0
+ 0
− 13
View file @ 5c12adcf
function [data] = getSampleData(amp,freq,nCycles,sampleTime)
%getSampleData Summary of this function goes here
% Detailed explanation goes here
data.time.sampletime=sampleTime;
data.time.unit="Seconds";
data.numberpoints= round(nCycles/sampleTime/freq);
data.time.value=[0:sampleTime:data.numberpoints*sampleTime];
data.numberpoints=data.numberpoints+1;
data.value=amp*sin(2*pi*freq*data.time.value);
data.unit="";
end
\ No newline at end of file
getStatisticoscilation.m deleted 100644 → 0
+ 0
− 18
View file @ 5c12adcf
function [data_mean, data_std] = getStatisticoscilation (value, numberOscilations)
% Umbau des Vektors in mean und std
ind_start = 1;
pointsOscilation = round(length(value)/numberOscilations)-1;
mat=[];
for ii=1:numberOscilations
ind_end=ind_start+pointsOscilation;
mat(ii,:)=value(ind_start:ind_end);
ind_start=ind_end+1;
end
for jj=1:pointsOscilation+1
data_mean(jj) = mean(mat(:,jj));
data_std(jj) = std(mat(:, jj));
end
end
\ No newline at end of file
getSysunc.m deleted 100644 → 0
+ 0
− 8
View file @ 5c12adcf
function [uncsysabs uncsysphase] = getSysunc(sysunc,Y)
%UNTITLED9 Summary of this function goes here
% Detailed explanation goes here
uncsysabs = sysunc.bias+sysunc.slope+0.608*sysunc.lin
uncsysphase = atan(0.681* sysunc.hys ./abs(Y));
end
\ No newline at end of file
propagateSystematicUncertainityToFrequencyDomain_package/propagateSysUncToFreqDomain.m 0 → 100644
+ 135
− 0
View file @ 3a853f85
function result_struct = propagateSysUncToFreqDomain(value_vector, time_vector, sampletime, varargin)
% PROPAGATESYSUNCTOFREQDOMAIN
%
% result_struct = propagateSysUncToFreqDomain( ...
% value_vector,
% time_vector,
% sampletime,
% varargin)
%
% -- value_vector needs to be a 1-by-N vector containing the values of
% the signal.
%
% -- time_vector needs to be a 1-by-N vector containing the time values
% of the signal.
%
% -- sampletime needs to be a scalar natural number (1-by-1) of type
% integer that describes the time between steps of the time_vector.
%
% -- varargin (optional parameters)
% -- 'bias' scalar (1-by-1) real number of type double that
% represents the *bias* systematic uncertainty and can be obtained
% from the data sheet of the sensor that was used for the
% measurement or from calibration protocols.
%
% -- 'slope' scalar (1-by-1) real number of type double that
% represents the *slope* systematic uncertainty and can be obtained
% from the data sheet of the sensor that was used for the
% measurement or from calibration protocols.
%
% -- 'lin' scalar (1-by-1) real number of type double that
% represents the *linearity* systematic uncertainty and can be
% *sometimes (not always given)* obtained from the data sheet of the
% sensor that was used for the measurement or from calibration
% protocols.
% Use 0 if there isn't a given value for your sensor.
%
% -- 'hys' scalar (1-by-1) real number of type double that
% represents the *hysteresis* systematic uncertainty and can be
% *sometimes (not always given)* obtained from the data sheet of the
% sensor that was used for the measurement or from calibration
% protocols.
% Use 0 if there isn't a given value for your sensor.
%
%
% Returns:
% -- result_struct A struct that contains the input data, interim
% results, the absolute systematic uncertainty in the frequency
% domain and the phase systematic uncertainty for every frequency
% gathered from the FFT.
%
% -- TODO: explanation of the data structure
%
%
%% Input adaption
defaultBias = 0;
defaultSlope = 0;
defaultLin = 0;
defaultHys = 0;
p = inputParser;
validScalarPosNum = @(x) isnumeric(x) && (x > 0);
addOptional(p, 'bias', defaultBias, validScalarPosNum);
addOptional(p, 'slope', defaultSlope, validScalarPosNum);
addOptional(p, 'lin', defaultLin, validScalarPosNum);
addOptional(p, 'hys', defaultHys, validScalarPosNum);
parse(p, varargin{:})
%% Generate the uncertainty data structure with the given input arguments
data.uncertainty.systematic.bias = p.Results.bias;
data.uncertainty.systematic.slope = p.Results.slope;
data.uncertainty.systematic.lin = p.Results.lin;
data.uncertainty.systematic.hys = p.Results.hys;
%% Pre Process the data
% Analyse main frequency
result_struct.excitation_frequency = getMainFreq(value_vector, time_vector);
% Get number of cycles
result_struct.number_of_cycles = getNumberOfCycles(time_vector, result_struct.excitation_frequency);
% Crop the number of cycles
[result_struct.valuecrop, ...
result_struct.timecrop, ...
result_struct.number_of_cycles] = cropData( ...
value_vector, ...
time_vector, ...
sampletime, ...
result_struct.excitation_frequency, ...
result_struct.number_of_cycles);
% Statistical analysis
[result_struct.value_mean_vector, ...
result_struct.value_standartdeviation_vector] ...
= getTheStatisticOscillation( ...
result_struct.valuecrop, ...
result_struct.number_of_cycles);
%% FFT with uncertainty
% FFT of mean values
result_struct.FFT.N_data_points_in_FFT = length(result_struct.value_mean_vector);
result_struct.FFT.mean = fft(result_struct.value_mean_vector);
% Calculate the statistical uncertainty
% TODO: Maybe also write a test for the statistical uncertainty
temp_term1 = result_struct.value_standartdeviation_vector * tinv(1 - (0.05/2), result_struct.number_of_cycles-1);
temp_term2 = (sqrt(result_struct.number_of_cycles))^2 * result_struct.FFT.N_data_points_in_FFT;
result_struct.FFT.uncertainty.statistical = sqrt(sum(temp_term1 / temp_term2));
clearvars temp_term1 temp_term2
% Scale coefficients of the FFT
result_struct.FFT.mean = 2 / result_struct.FFT.N_data_points_in_FFT * result_struct.FFT.mean ;
result_struct.FFT.mean(1) = result_struct.FFT.mean(1) / 2;
% Get frequency vector of the FFT
temp_term1 = (0:(result_struct.FFT.N_data_points_in_FFT / 2));
temp_term2 = (result_struct.FFT.N_data_points_in_FFT * sampletime);
result_struct.FFT.frequencies = temp_term1 / temp_term2;
clearvars temp_term1 temp_term2
% Propagate the systematic uncertainty into the frequency domain
[result_struct.FFT.uncertainty.systematic.absolute, ...
result_struct.FFT.uncertainty.systematic.phase] = ...
calculateSysUncInFreqDomain(...
data.uncertainty.systematic.bias, ...
data.uncertainty.systematic.slope, ...
data.uncertainty.systematic.lin, ...
data.uncertainty.systematic.hys, ...
result_struct.FFT.mean);
end
propagateSystematicUncertainityToFrequencyDomain_package/underlying_functions/calculateSysUncInFreqDomain.m 0 → 100644
+ 63
− 0
View file @ 3a853f85
function [systematic_uncertainty_absolute, ...
systematic_uncertainty_phase] = ...
calculateSysUncInFreqDomain( ...
sysunc_bias, ...
sysunc_slope, ...
sysunc_linearity, ...
sysunc_hysteresis, ...
FFT_Mean)
% CALCULATE SYSUNCINFREQDOMAIN Calculate the systematic uncertainty in the
% frequency domain with the mathematical expression given by the paper
% mentioned in the README.md of this repository/package.
%
% [systematic_uncertainty_absolute, ...
% systematic_uncertainty_phase] = ...
% calculateSysUncInFreqDomain( ...
% sysunc_bias, ...
% sysunc_slope, ...
% sysunc_linearity, ...
% sysunc_hysteresis, ...
% FFT_Mean)
%
% -- sysunc_bias scalar (1-by-1) real number of type double that
% represents the *bias* systematic uncertainty and can be obtained from
% the data sheet of the sensor that was used for the measurement or from
% calibration protocols.
%
% -- sysunc_slope scalar (1-by-1) real number of type double that
% represents the *slope* systematic uncertainty and can be obtained from
% the data sheet of the sensor that was used for the measurement or from
% calibration protocols.
%
% -- sysunc_linearity scalar (1-by-1) real number of type double that
% represents the *linearity* systematic uncertainty and can be *sometimes
% (not always given)* obtained from the data sheet of the sensor that was
% used for the measurement or from calibration protocols.
% Use 0 if there isn't a given value for your sensor.
% -- sysunc_hysteresis scalar (1-by-1) real number of type double that
% represents the *hysteresis* systematic uncertainty and can be
% *sometimes (not always given)* obtained from the data sheet of the
% sensor that was used for the measurement or from calibration protocols.
% Use 0 if there isn't a given value for your sensor.
% -- FFT_Mean FFT of the mean value vector of the "statistic
% oscillation".
%
%
% Returns:
% -- systematic_uncertainty_absolute scalar (1-by-1) real number of type
% 'double' that represents the absolute systematic uncertainty of the
% measurement signal with the given sensor in the frequency domain.
%
% -- systematic_uncertainty_phase (1-by-N) vector with real numbers of
% type 'double' that represents the systematic uncertainty of the phase in
% the different frequencies of the FFT of the measurement signal with the
% given sensor.
systematic_uncertainty_absolute = sysunc_bias + sysunc_slope + (0.608 * sysunc_linearity);
systematic_uncertainty_phase = atan(0.681 * sysunc_hysteresis ./ abs(FFT_Mean));
end
propagateSystematicUncertainityToFrequencyDomain_package/underlying_functions/cropData.m 0 → 100644
+ 64
− 0
View file @ 3a853f85
function [cropped_value_vector, ...
cropped_time_vector, ...
floored_number_of_cycles] = cropData(...
data_vector,...
time_vector, ...
sampletime, ...
main_frequency, ...
number_of_cycles)
% CROPDATA Crops a harmonic osciallation signal to N (numberofcycles) full
% oscillations beginning at the start of the signal
%
% [cropped_value_vector, cropped_time_vector, floored_numberofcycles] =
% CROPDATA(data_vector,
% time_vector,
% sampletime,
% main_frequency,
% numberofcycles)
%
% -- data_vector needs to be a 1-by-N vector containing the values of
% the signal.
%
% -- time_vector needs to be a 1-by-N vector containing the time values
% of the signal.
%
% -- sampletime needs to be a scalar natural number (1-by-1) of type
% integer that describes the time between steps of the time_vector.
%
% -- main_frequency needs to be a scalar number (1-by-1) of type double
% or int and can be obtained through the getMainFreq(...) function of
% this package this file also belongs to.
%
% -- numberofcycles as a scalar real number of type double and can be
% obtained through the getNumberOfCycles(...) function of this package
% this file also belongs to.
%
%
% Returns:
% -- cropped_value_vector the 1-by-N value vector cropped down to a
% length that only contains a natural number of full harmonic
% oscillations/cycles containing the values of the signal.
%
% -- cropped_time_vector the 1-by-N time vector cropped down to a length
% that only contains a natural number of full harmonic
% oscillations/cycles containing the values of the signal.
%
% -- floored_number_of_cycles as a scalar natural number of type integer.
% It is the rounded down (to a natural number) version of the input
% numberofcycles.
%
%
% See also GETMAINFREQ, GETNUMBEROFCYCLES.
y = data_vector;
t = time_vector;
floored_number_of_cycles = floor(number_of_cycles);
Ndft = round((floored_number_of_cycles * (1/main_frequency)) / sampletime);
% Cut n-periods from the timeseries
cropped_value_vector = y(1:Ndft);
cropped_time_vector = t(1:Ndft);
end
propagateSystematicUncertainityToFrequencyDomain_package/underlying_functions/getMainFreq.m 0 → 100644
+ 30
− 0
View file @ 3a853f85
function main_frequency = getMainFreq(value_vector, time_vector)
% GETMAINFREQ gets the main frequency of a harmonic oscillations signal
%
% main_frequency = GETMAINFREQ(value_vector, time_vector)
% gets the main frequency of a harmonic oscillations signal described by
% a value vector and time vector through regression analysis.
% It uses Mathworks Matlab fit(...) from the 'Mathworks Matlab Curve
% Fitting Toolbox' with the 'fourier1' option.
%
% -- value_vector needs to be a 1-by-N vector containing the values of
% the signal.
%
% -- time_vector needs to be a 1-by-N vector containing the time values
% of the signal.
%
%
% Returns:
% -- main_frequency as a scalar number of type double.
%
%
% See also FIT.
y = value_vector;
t = time_vector;
[c2,~] = fit(t', y', 'fourier1');
main_frequency = c2.w / (2*pi);
end
propagateSystematicUncertainityToFrequencyDomain_package/underlying_functions/getNumberOfCycles.m 0 → 100644
+ 32
− 0
View file @ 3a853f85
function number_of_cycles = getNumberOfCycles(time_vector, main_frequency)
% GETNUMBEROFCYCLES returns the count of full oscillations of a signal
%
% number_of_cycles = GETNUMBEROFCYCLES(time_vector, main_frequency)
% returns the count of full oscillations of the signal given by the
% time_vector and the main frequency of the value vector.
%
% -- time_vector needs to be a 1-by-N vector containing the time values
% of the signal.
%
% -- main_frequency needs to be a scalar number (1-by-1) of type 'double'
% or 'int' and can be obtained through the getMainFreq(...) function of
% this package this file also belongs to.
%
% Returns:
% -- number_of_cycles as a scalar real number of type 'double'.
%
%
% See also GETMAINFREQ.
t = time_vector;
number_of_cycles = (t(end) - t(1)) / (1/main_frequency);
if number_of_cycles < 1
error(append("There is not at least one complete harmonic cycle in", ...
"the signal you have provided. The data cannot be", ...
"analyzed."))
end
end
propagateSystematicUncertainityToFrequencyDomain_package/underlying_functions/getTheStatisticOscillation.m 0 → 100644
+ 55
− 0
View file @ 3a853f85
function [value_mean_vector, value_standartdeviation_vector] = getTheStatisticOscillation(cropped_value_vector, floored_number_of_cycles)
% GETSTATISTICOSCILLATION get the "statistic oscillation" (mean value
% vector and the standart deviation vector) by taking all available cycles
% of the harmonic osciallation, "moving them over each other" and
% calculating the mean and standart deviation for every data points that are
% superimposed upon each other.
% [value_mean_vector,
% value_standartdeviation_vector] =
% getTheStatisticOscillation( ...
% cropped_value_vector
% floored_number_of_cycles)
%
% -- cropped_value_vector needs to be a 1-by-N vector containing the
% values of the signal. It needs to be cropped to a length that
% only a natural number of *full* harmonic oscillations/cycles are
% present. The maximum amount of *full* harmonic oscillations/cycles is
% the optimum you should use.
%
% -- floored_numberofcycles as a scalar natural number of type integer.
% It is the rounded down (to a natural number) version of the input
% number_of_cycles.
%
%
% Returns:
% -- value_mean_vector the 1-by-N value vector of the "statistic
% oscillation" that has the length of
% length(cropped_value_vector)/floored_numberofcycles.
%
% -- value_standartdeviation_vector the 1-by-N value standart deviation
% vector of the "statistic oscillation" that also has the length of
% length(cropped_value_vector)/floored_numberofcycles.
%
%
% See also GETNUMBEROFCYCLES.
data_points_per_cycle = round(length(cropped_value_vector) / floored_number_of_cycles) - 1;
data_matrix = [];
start_index = 1;
for ii = 1:floored_number_of_cycles
end_index = start_index + data_points_per_cycle;
data_matrix(ii,:) = cropped_value_vector(start_index : end_index);
% For every oscillation reset the start_index that only one
% oscillation/cycle gets handled.
start_index = end_index + 1;
end
for jj = 1:(data_points_per_cycle + 1)
value_mean_vector(jj) = mean(data_matrix(:, jj));
value_standartdeviation_vector(jj) = std(data_matrix(:, jj));
end
end
sample_uncertatinty_pagation.m deleted 100644 → 0
+ 0
− 53
View file @ 5c12adcf
clear data
%% Generate sample Data
% sample time series
% getSampleData(amp,freq,nCycles,sampleTime)
data = getSampleData(1,5,10.1,0.001);
% add some noise
data.value= data.value+0.1*randn(1,data.numberpoints);
% add uncertainty
% addSysunc(data,bias,slope,lin,hys)
data = addSysunc(data,'bias',0.01,'slope',0.02,'lin',0.03,'hys',0.04);
% addStatunc(data,sigma)
% data = addStatunc(data,0.001);
%% Test data structure
%testDataStructure(data)
%% Pre Processing
% analyse main frequency
res.excitationfreq=getMainfreq(data);
disp(['Excitation frequency: ',num2str(res.excitationfreq)])
% get Number of cycles
res.numberofcycles = getNumberofcycles(data,res.excitationfreq);
disp(['Number of cycles: ',num2str(res.numberofcycles)])
% crop to Number of cycles
[res.timecrop,res.valuecrop,res.numberofcycles] = cropData(data,res.excitationfreq,res.numberofcycles);
% statistical analysis
[res.valuemean,res.valuestd] = getStatisticoscilation (res.valuecrop,res.numberofcycles);
%% FFT with uncertainty
% fft of mean values
res.FFT.Ndft = length(res.valuemean);
res.FFT.mean = fft(res.valuemean);
res.FFT.uncstat = sqrt(sum(res.valuestd*tinv(1-0.05/2,res.numberofcycles-1)/sqrt(res.numberofcycles))^2/res.FFT.Ndft);
% Scale coeffitients
res.FFT.mean = 2 / res.FFT.Ndft * res.FFT.mean ;
res.FFT.mean(1)=res.FFT.mean(1)/2;
% get frequency vector
res.FFT.frequencies=(0:(res.FFT.Ndft/2))/res.FFT.Ndft/data.time.sampletime;
% propagate systematic uncertainty
[res.FFT.uncsysabs, res.FFT.uncsysphase]= getSysunc(data.unc.sys,res.FFT.mean);
% total uncertainty
% analysis for main frequancy abs / phase and imag / real
%% Plots
figure(1)
plot(res.FFT.frequencies,abs(res.FFT.mean(1:length(res.FFT.frequencies))),'o')
tStatisticOscilation.m deleted 100644 → 0
+ 0
− 0
View file @ 5c12adcf
testDataStructure.m deleted 100644 → 0
+ 0
− 19
View file @ 5c12adcf
function [result] = testDataStructure(data)
%UNTITLED5 Summary of this function goes here
% Detailed explanation goes here
data.value
data.unit
data.date
data.points
data.time.value
data.time.unit
data.time.sampletime
data.unc.sys.bias
data.unc.sys.slope
data.unc.sys.lin
data.unc.sys.hys
data.unc.stat.sigma
%% Results
data.res.excitationfreq
end
\ No newline at end of file
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment