JumpSensor.m

Go to the documentation of this file.

% =========================================================================== %
%> @brief JumpSensor is a simulation of the JUMP sensor of the ZurichMove Project.
%>
%> JumpSensor simulates one JUMP sensor of the ZurichMove Project. The JUMP
%> sensor is a wearable sensor which can be used time synchronized with other
%> JUMP sensors to track the movement of a person.
%>
%> This class is a simulation of the JUMP sensor. It can be configured equally
%> to the real JUMP sensor. It contains the models of the actual sensors of the
%> JUMP sensor according to the datasheets of the sensors. The JUMP sensor
%> simulation can be used to test the algorithms of the ZurichMove Project.
%>
%> The JUMP sensor is equipped with the following sensors:
%> - MPU9250 9-axis IMU containing                       :
%>       - 3-axis accelerometer
%>       - 3-axis gyroscope
%>       - 3-axis magnetometer
%>   - ADXL375 3-axis accelerometer (high-g)
%>   - MS5611 barometer and temperature sensor
%>
%> @code
%> % example code
%> % create a JUMP sensor with default settings
%> jump = JumpSensor();
%>
%> % generate true data (examples, use the \ref trajectoryGenerationApp.m for more realistic data)
%> acc_true = zeros(1000,3); % laying still
%> gyr_true = zeros(1000,3); % not moving
%> ori_true = quaternion.ones(1000,1); % no rotation
%>
%> % generate sensor data
%> data = jump.generate(acc_true, gyr_true, ori_true);
%> @endcode
%>
%> see also the more detailed example in @ref jumpSensorExample.m
%>
%> @copyright see the file @ref LICENSE in the root directory of the repository 
% =========================================================================== %
 
classdef JumpSensor < handle
    
    properties (Access = private)
        % default settings from JUMPconnect
        accRange               = 16 * 9.81; %> default acceleration range (m/s^2)
        gyroRange              = 2000; %> default gyroscope range (°/s)
        sampleRate             = 200; %> default sample rate of the MPU (Hz)
        sampleRatePressure     = 5; %> default sample rate of the pressure and temperature sensor (Hz)
        sampleRateMagnetometer = 10; %> default sample rate of the magnetometer (Hz)
        
        % sensors
        mpu {mustBeA(mpu, 'imuSensor')} = imuSensor; %> model of the MPU9250 9-axis IMU
        hig {mustBeA(hig, 'imuSensor')} = imuSensor; %> model of the ADXL375 3-axis accelerometer
 
        %> @note the barometer is implemented directly in the generate function
        
    end
    
    methods (Access = public)
 
        % ====================================================================== %
        %> @brief JumpSensor is the constructor of the JUMP sensor.
        %>
        %> JumpSensor creates a JUMP sensor with the given settings. The
        %> settings are the same as in the JUMPconnect software. The JUMP sensor
        %> can be used to generate sensor data from true data.
        %>
        %> @param   accRange (optional)           - the acceleration range of the MPU9250 in g
        %> @param   gyroRange (optional)          - the gyroscope range of the MPU9250 in °/s
        %> @param   sampleRate (optional)         - the sample rate of the MPU9250 in Hz
        %> @param   sampleRatePressure (optional) - the sample rate of the pressure and temperature sensor in Hz
        %>
        %> @retval  obj                           - the created JUMP sensor
        % ====================================================================== %
        function obj = JumpSensor(args)
            arguments
                args.accRange           (1,1) double {mustBeMember(args.accRange, [2,4,8,16])}           = 16;
                args.gyroRange          (1,1) double {mustBeMember(args.gyroRange, [250,500,1000,2000])} = 2000;
                args.sampleRate         (1,1) double {mustBeMember(args.sampleRate, [50,200])}           = 200;
                args.sampleRatePressure (1,1) double {mustBeMember(args.sampleRatePressure, [1,5])}      = 5;
            end
            
            obj.accRange           = args.accRange * 9.81; % convert to m/s^2
            obj.gyroRange          = args.gyroRange;
            obj.sampleRate         = args.sampleRate;
            obj.sampleRatePressure = args.sampleRatePressure;
            
            % make checks on the settings
            if obj.gyroRange ~= 2000
                warning("The JUMPconnect software only allows the setting 2000 °/s.");
            end
            if (obj.sampleRate == 50 && obj.sampleRatePressure == 5) || (obj.sampleRate == 200 && obj.sampleRatePressure == 1)
                warning("The JUMPconnect software does not allow the combination of %i Hz sample rate and %i Hz pressure sensor sample rate.", obj.sampleRate, obj.sampleRatePressure);
            end
            
            % create the sensors
            obj.createMpu();
            obj.createHig();
            
            
        end
 
        % ====================================================================== %
        %> @brief generate generates sensor data from true data.
        %>
        %> generate generates sensor data from true data with the defined sensor
        %> models. The data format is identical to the data format of the parsing
        %> function of the function JUMPREADDATA.
        %>
        %> @param acc_true (nx3 double)                     - the true accelerometer data in m/s^2
        %> @param gyr_true (nx3 double)                     - the true gyroscope data in deg/s
        %> @param ori_true (nx1 quaternion)                 - the true orientation data in quaternions
        %> @param options.alt_true (nx1 double, optional)   - the true altitude in m (default: 500 m)
        %> @param options.tmp_true (nx1 double, optional)   - the true temperature in °C (default: 25°C)
        %> @param options.t_start (datetime, optional)      - the start time of the sensor data (default: now)
        %> @param options.simulate_mpu (logical, optional)  - simulate the MPU9250 sensor (default: true)
        %> @param options.simulate_hig (logical, optional)  - simulate the ADXL375 sensor (default: true)
        %> @param options.simulate_prs (logical, optional)  - simulate the MS5611 sensor (default: true)
        %> @param options.simulate_tmp (logical, optional)  - simulate the temperature sensor (default: true)
        %>
        %> @retval data (struct)                           - the sensor data in the format of the parsing function of the function JUMPREADDATA
        %>  - data.tt (timetable)                          - timetable containing the sensor data
        %>  - data.ID (char)                               - sensor ID
        %>  - data.UID (char)                              - sensor UID
        %>  - data.name (char)                             - sensor name (fixed to 'JUMP model')
        %>  - data.sw_version (char)                       - sensor software version
        %>  - data.hw_version (char)                       - sensor hardware version
        %>  - data.acc_range (double)                      - accelerometer range
        %>  - data.gyro_range (double)                     - gyroscope range
        %>  - data.mpu_cal (double)                        - mpu calibration
        %>  - data.act_thrs (double)                       - activity threshold
        %>  - data.inact_dur (double)                      - inactivity duration
        %>  - data.mpu_rate (double)                       - mpu rate
        %>  - data.adxl_rate (double)                      - adxl rate
        %>  - data.pres_rate (double)                      - pressure rate
        %>  - data.magn_en (double)                        - magnetometer enabled
        %>  - data.sync (double)                           - sync mode
        %>  - data.ble_adress (char)                       - ble adress
        %>
        %> @note for whatever reason the MATLAB imuSensor class negates the
        %> accelerometer data. I guess, because they missunderstand how the
        %> gravitational force works. According to their algorithm description,
        %> they negate the input acceleration and then add the gravity vector.
        %> This is not the correct way to do it. The correct way is to negate the
        %> gravity vector. Therefore, the input is negated and the gravity vector
        %> is added two times before the imuSensor class is called. This is
        %> a bug in the MATLAB imuSensor class. In ither words, we have the
        %> actual acceleration \f$ a \f$. to patch the MATLAB bevavior, we do:
        %> \f[ a_{patched} = -a + 2 \cdot [0,0,g], \f] before we call the
        %> imuSensor class.then the imuSensor class does
        %> \f[ totalAcc = -a_{input} + [0,0, g] = -(-a + 2 \cdot [0,0,g]) + 
        %> [0,0, +g] = a - [0,0,g], \f] which brings the correct result.
        %>
        %> @todo check the covariances of the sensors
        %>
        %> @todo check the axis misalignment and sensitivity over time.
        % ====================================================================== %        
        function data = generate(obj, acc_true, gyr_true, ori_true, options)
            
            % check input arguments
            arguments
                obj;
                acc_true (:,3) double;
                gyr_true (:,3) double;
                ori_true (:,1) quaternion;
                options.alt_true     (:,1) double   = repmat(500, length(acc_true), 1);
                options.tmp_true     (:,1) double   = repmat(25, length(acc_true), 1);
                options.t_start      (1,1) datetime = datetime("now");
                options.simulate_mpu (1,1) logical  = true;
                options.simulate_hig (1,1) logical  = true;
                options.simulate_prs (1,1) logical  = true;
                options.simulate_tmp (1,1) logical  = true;
            end
 
            % unpack options
            alt_true = options.alt_true;
            tmp_true = options.tmp_true;
            t_start  = options.t_start;
            
            % check if the input vectors have the same length
            input_length = length(acc_true);
            assert(length(gyr_true) == input_length, "The input vectors must have the same length.");
            assert(length(ori_true) == input_length, "The input vectors must have the same length.");
            assert(length(tmp_true) == input_length, "The input vectors must have the same length.");
 
            % modify the acceleration because of a MATLAB bug
            acc_true = -acc_true + 2*[0,0,9.81]; % negate the acceleration and add 2x the gravity vector
            
            if options.simulate_mpu
                
                % generate the sensor data
                gyr_true = gyr_true * pi / 180; % convert to rad/s
                
                % mpu
                [acc_mpu, gyr_mpu, mag_mpu]        = obj.mpu(acc_true, gyr_true, ori_true);
                gyr_mpu                            = gyr_mpu * 180 / pi; % convert to deg/s
 
                % nownsample and then upsample with previous value to get the same output as the real sensor
                mag_mpu_ds                         = downsample(mag_mpu, obj.sampleRate / obj.sampleRateMagnetometer); % NOTE: the magnetometer is sampled at 10 Hz
                mag_mpu                            = nan(input_length,3);
                mag_mpu(1:obj.sampleRate/obj.sampleRateMagnetometer:end,:) = mag_mpu_ds;
                mag_mpu                           = fillmissing(mag_mpu, 'previous');
 
            else
                acc_mpu = nan(input_length,3);
                gyr_mpu = nan(input_length,3);
                mag_mpu = nan(input_length,3);
            end
            
            if options.simulate_hig
                % adxl
                acc_adxl = obj.hig(acc_true, zeros(input_length,3), ori_true);
            else
                acc_adxl = nan(input_length,3);
            end
            
            if options.simulate_prs
                % pressure
                prs = obj.getPressureFromAltitude(alt_true, tmp_true);
            else
                prs = repmat(1013.25, input_length, 1);
            end
            
 
            if options.simulate_tmp
                % temperature
                tmp = obj.getSensorTemperature(tmp_true);
            else
                tmp = tmp_true;
            end
            
            % create timetable
            t = repmat(t_start, input_length, 1) + ((0:input_length-1) .* seconds(1/obj.sampleRate))';
            data.tt                          = timetable(t, acc_mpu, gyr_mpu, acc_adxl, mag_mpu, prs, tmp);
            data.tt.Properties.VariableNames = {'acc', 'gyr', 'hig', 'mag', 'prs', 'tmp'};
            
            % add sensor information
            data.t_start     = t_start;
            data.ID          = 'S-24-A-SIMUL';
            data.UID         = 'simulation';
            data.name        = 'JUMP simulation';
            data.sw_version  = '132000(s)';
            data.hw_version  = '1627390748(à)';
            data.acc_range   = round(obj.accRange / 9.81); % convert to g
            data.gyro_range  = round(obj.gyroRange * 180 / pi); % convert to deg/s
            data.mpu_cal     = 2;
            data.act_thrs    = 0;
            data.inact_dur   = 0;
            data.mpu_rate    = obj.sampleRate;
            data.adxl_rate   = obj.sampleRate;
            data.pres_rate   = obj.sampleRatePressure;
            data.magn_en     = 255;
            data.sync        = 1;
            data.ble_address = '001122334455';
            
        end
 
        % ====================================================================== %
        %> @brief generates data from a @ref Trajectory object.
        %>
        %> generateFromTrajectory generates sensor data from a @ref Trajectory
        %> object. The sensor data is generated with the defined sensor models.
        %> 
        %> @param obj (JumpSensor) - the JUMP sensor
        %> @param traj (Trajectory) - the trajectory object
        %>
        %> @retval data (struct) - the sensor data in the format defined in @ref generate
        % ====================================================================== %
        function data = generateFromTrajectory(obj, traj)
            
            % check input arguments
            arguments
                obj;
                traj (1,1) Trajectory;
            end
            
            traj = traj.getTrajectory;
 
            % generate the sensor data
            data = obj.generate(traj.acc, traj.gyr, traj.ori);
        end
 
        % ====================================================================== %
        %> @brief returns the parameters of the JUMP sensor.
        %>
        %> getParams returns the parameters of the JUMP sensor. The parameters
        %> are used to generate the erroneous sensor data.
        %>
        %> @param obj (JumpSensor) - the JUMP sensor
        %>
        %> @retval params (struct) - the parameters of the JUMP sensor
        %> - params.mpu.acc (struct) - the parameters of the MPU9250 accelerometer
        %> - params.mpu.gyr (struct) - the parameters of the MPU9250 gyroscope
        %> - params.mpu.mag (struct) - the parameters of the MPU9250 magnetometer
        %> - params.hig.acc (struct) - the parameters of the ADXL375 accelerometer
        % ====================================================================== %
        function params = getParams(obj)
            
            params.mpu.acc = obj.mpu.Accelerometer;
            params.mpu.gyr = obj.mpu.Gyroscope;
            params.mpu.mag = obj.mpu.Magnetometer;
            params.hig.acc = obj.hig.Accelerometer;
            
        end
    end
    
    methods (Access = private)
 
        % ====================================================================== %
        %> @brief createMpu creates the MPU9250 sensor.
        %>
        %> createMpu creates the MPU9250 sensor with the settings specified in
        %> the properties of the JUMP sensor model.
        %>
        %> @param obj (JumpSensor) - the JUMP sensor
        %>
        %> @retval none
        %>
        %> @note in here the actual sensor is modeled with noise and bias. See the
        %> actual code for the parameters and explanations.
        % ====================================================================== %
        function createMpu(obj)
            
            obj.mpu = imuSensor('accel-gyro-mag', ...          % sensor type
                "ReferenceFrame", "NED", ...                   % reference frame for JUMP sensors
                "SampleRate", obj.sampleRate, ...              % sample rate
                "Temperature", 25, ...                         % temperature (tunable)
                "MagneticField", [21.8719 1.2697, 42.8197] ... % magnetic field at 46.8N, 8.23E, 500m (GPS) (NED, according to Magbnetic Field Estimated Values of NOAA)
                );
            
            % accelerometer parameters
            obj.mpu.Accelerometer.MeasurementRange = obj.accRange * 9.81; % m/s^2
            obj.mpu.Accelerometer.Resolution       = obj.accRange / 32768; % (m/s^2) / LSB
            obj.mpu.Accelerometer.ConstantBias     = [obj.rant(0.06), obj.rant(0.06), obj.rant(0.08)] * 9.81; % (m/s^2) (datasheet: 60,60,80 mg)
            % sensitivity is modeled into the misalignment matrix
            obj.mpu.Accelerometer.AxesMisalignment       = obj.getMisalignmentMatrix(2,2,2,3); % (percent cross-axis sensitivity, percent overall sensitivity error)
            obj.mpu.Accelerometer.NoiseDensity           = [obj.rant(0.0025,0.0035,2),obj.rant(0.0025,0.0035,2),obj.rant(0.0055,0.0065,2)]; % (m/s^2) / sqrt(Hz) (datasheet: 300 µg/rtHz, value from allan variance experiments)
            obj.mpu.Accelerometer.BiasInstability        = 6e-3; % (m/s^2) (missing in datasheet, value from allan variance experiments)
            obj.mpu.Accelerometer.RandomWalk             = 3.9e-5; % (m/s^2) * sqrt(Hz) (missing in datasheet, value from allan variance experiments)
            obj.mpu.Accelerometer.TemperatureBias        = [obj.rant(0.0015), obj.rant(0.0015), obj.rant(0.0015)] * 9.81; % (m/s^2) / °C (datasheet: 1.5 mg/°C)
            obj.mpu.Accelerometer.TemperatureScaleFactor = abs(obj.rant(0.026)); % (%/°C) (datasheet: 0.026 %/°C)
            obj.mpu.Accelerometer.NoiseType              = "single-sided";
            
            % gyroscope parameters
            obj.mpu.Gyroscope.MeasurementRange = obj.gyroRange * pi/180; % rad/s
            obj.mpu.Gyroscope.Resolution       = (obj.gyroRange * pi/180) / 32768; % (rad/s) / LSB
            obj.mpu.Gyroscope.ConstantBias     = [obj.rant(5), obj.rant(5), obj.rant(5)] * pi/180; % (rad/s) (datasheet: +-5 dps)
            % sensitivity is modeled into the misalignment matrix
            obj.mpu.Gyroscope.AxesMisalignment       = obj.getMisalignmentMatrix(2,2,2,3); % (percent cross-axis sensitivity)
            obj.mpu.Gyroscope.NoiseDensity           = [obj.rant(4.7e-3,6e-3,2),obj.rant(4.7e-3,6e-3,2),obj.rant(4.7e-3,6e-3,2)] * pi/180; % (rad/s) / sqrt(Hz) (datasheet: 10 mdps/rtHz, value from allan variance experiments)
            obj.mpu.Gyroscope.BiasInstability        = 8.3e-3 * pi/180; % (rad/s) (missing in datasheet, value from allan variance experiments)
            obj.mpu.Gyroscope.RandomWalk             = 155.9e-6 * pi/180; % (rad/s) * sqrt(Hz) (missing in datasheet, value from allan variance experiments)
            obj.mpu.Gyroscope.TemperatureBias        = [obj.rant(0.48), obj.rant(0.48), obj.rant(0.48)] * pi / 180; % (rad/s) / °C (datasheet: (60 °/s) / 125 °C)
            obj.mpu.Gyroscope.TemperatureScaleFactor = abs(obj.rant(3)); % (%/°C) (datasheet: 3 %/°C)
            obj.mpu.Gyroscope.AccelerationBias       = 0.1 * pi/180; % (rad/s) / (m/s^2) (NOTE: missing in datasheet, this is an estimation from https://www.analog.com/en/resources/technical-articles/gyro-mechanical-performance.html)
            obj.mpu.Gyroscope.NoiseType              = "single-sided";
            
            % magnetometer parameters
            obj.mpu.Magnetometer.MeasurementRange = 4800; % µT
            obj.mpu.Magnetometer.Resolution       = 9600 / 2^14; % µT / LSB
            obj.mpu.Magnetometer.ConstantBias     = obj.mpu.Magnetometer.Resolution * randn() * 500; % µT (datasheet: +-500 LSB)
            obj.mpu.Magnetometer.AxesMisalignment = diag([100,100,100]); % (percent cross-axis sensitivity) (NOTE: missing in datasheet)
            obj.mpu.Magnetometer.NoiseDensity     = obj.rant(0.65,0.75,2); % µT / sqrt(Hz) (missing in datasheet, value from allan variance experiments)
            obj.mpu.Magnetometer.BiasInstability  = 0.3673; % µT (missing in datasheet, value from allan variance experiments)
            obj.mpu.Magnetometer.RandomWalk       = 0.01216; % µT * sqrt(Hz) (missing in datasheet, value from allan variance experiments)
            % NOTE: temperature bias and scale factor are missing in the datasheet
            obj.mpu.Magnetometer.NoiseType = "single-sided";
        end
 
        % ====================================================================== %
        %> @brief createHig creates the ADXL375 sensor.
        %>
        %> createHig creates the ADXL375 sensor with the settings specified in
        %> the properties of the JUMP sensor model.
        %>
        %> @param obj (JumpSensor) - the JUMP sensor
        %>
        %> @retval none
        %>
        %> @note in here the actual sensor is modeled with noise and bias. See the
        %> actual code for the parameters and explanations.
        % ====================================================================== %
        function obj = createHig(obj)
            
            obj.hig = imuSensor('accel-gyro', ... % sensor type (NOTE: there is no actual gyro, but the imuSensor class requires it)
                "ReferenceFrame", "NED", ...      % reference frame for JUMP sensors
                "SampleRate", obj.sampleRate, ... % sample rate
                "Temperature", 25 ...             % temperature (tunable)
                );
            
            % accelerometer parameters
            obj.hig.Accelerometer.MeasurementRange       = 200 * 9.81; % m/s^2 (datasheet: min 180g, typ: 200g)
            obj.hig.Accelerometer.Resolution             = 0.049 * 9.81; % (m/s^2) / LSB (datasheet: 49 mg/LSB)
            obj.hig.Accelerometer.ConstantBias           = [obj.rant(-0.4,0.4,1),obj.rant(-0.4,0.4,1),obj.rant(-0.4,0.4,1)] * 9.81; % (m/s^2) (datasheet: +- 0.4g 1 sigma)
            obj.hig.Accelerometer.AxesMisalignment       = obj.getMisalignmentMatrix(2.5); % (percent cross-axis sensitivity) (datasheet: 2.5%)
            obj.hig.Accelerometer.NoiseDensity           = obj.rant(0.0264,0.0381,2); % (m/s^2) / sqrt(Hz) (datasheet: 5 mg/rtHz, value from allan variance experiments)
            obj.hig.Accelerometer.BiasInstability        = 0.0038; % m/s^2 (missing in datasheet, value from allan variance experiments)
            obj.hig.Accelerometer.RandomWalk             = 174.95e-6; % (m/s^2) * sqrt(Hz) (missing in datasheet, value from allan variance experiments)
            obj.hig.Accelerometer.TemperatureBias        = abs(obj.rant(-0.01,0.01,1)) * 9.81; % (m/s^2) / °C (datasheet: 10 mg/°C)
            obj.hig.Accelerometer.TemperatureScaleFactor = abs(obj.rant(-0.02,0.02,1)); % (%/°C) (datasheet: 0.02 %/°C)
            obj.hig.Accelerometer.NoiseType              = "single-sided";
            
        end
        
        % ====================================================================== %
        %> @brief getPressureFromAltitude returns the pressure from the altitude.
        %>
        %> getPressureFromAltitude returns the pressure from the altitude. The
        %> pressure is modeled with the barometric formula and the actual sensor
        %> is modeled with noise and bias.
        %>
        %> @param obj (JumpSensor) - the JUMP sensor
        %> @param alt (nx1 double) - the altitude in m
        %> @param tmp (nx1 double, optional) - the temperature in °C (default: 25°C)
        %> @param prs_sealvl (1x1 double, optional) - the pressure at sea level in mbar (default: 1013.25 mbar)
        %>
        %> @retval mbar (nx1 double) - the pressure in mbar
        %>
        %> @note the experiment 240213_pressure_membrane suggested a mean raise
        %> time of 0.5s. Since this is approx the free fall duration (=0.49 s),
        %> we assume that the pressure sensor is able to measure the pressure with
        %> almost no delay. Therefore, we do not need to model the delay of the
        %> sensor.
        %>
        %> @note formula from: https://www.omnicalculator.com/physics/air-pressure-at-altitude, last accessed: 2024-02
        % ====================================================================== %
        function mbar = getPressureFromAltitude(obj, alt, tmp, prs_sealvl)            
            
            % input validation
            arguments
                obj;
                alt        (:,1) double {mustBeInRange(alt, 0, 2000)}
                tmp        (:,1) double {mustBeInRange(tmp, -20, 40)} = repmat(25, numel(alt),1); % temperature in °C
                prs_sealvl (1,1) double                               = 1013.25; % pressure at sea level in mbar
            end
 
            % downsample the altitude to the sample rate of the pressure sensor
            downsample_ratio = obj.sampleRate / obj.sampleRatePressure;
            alt_ds           = downsample(alt, downsample_ratio);
            tmp              = downsample(tmp, downsample_ratio);
            
            % calculate the pressure from the altitude (source: https://www.omnicalculator.com/physics/air-pressure-at-altitude)
            %
            % Formula: P = P0 * exp(-g * M * (h - h0) / (R * T))
            %
            % h : altitude in m
            % P : pressure at altitude h in Pascal
            % P0: pressure at sea level in Pascal
            % T : temperature at h in Kelvin
            % g : acceleration due to gravity in m/s^2
            % M : molar mass of Earth's air in kg/mol
            % R: universal gas constant in J/(mol*K)
            h  = alt_ds;
            P0 = prs_sealvl * 100; % convert to Pascal
            T  = tmp + 273.15; % convert to Kelvin
            g  = 9.80665; % m/s^2
            M  = 0.0289644; % kg/mol
            R  = 8.31432; % J/(mol*K)
            
            P    = P0 * exp(-g * M .* h ./ (R * T));
            mbar = P / 100; % convert to mbar
            
            % the experiment 240213_pressure_membrane suggested a mean raise time of 0.5s.
            % since this is approx the free fall duration (=0.49 s), we assume that the
            % pressure sensor is able to measure the pressure with almost no delay.
            % Therefore, we do not need to model the delay of the sensor.
            
            % model the actual sensor
            resolution      = 0.027; % mbar / LSB
            noiseDensity    = 0.0109; % mbar / sqrt(Hz) (from allan variance experiments)
            biasInstability = 0.013; % mbar (from allan variance experiments)
            randomWalk      = 0.003; % mbar * sqrt(Hz) (from allan variance experiments)
 
            % generate noise
            noise = obj.generateNoise(length(mbar), noiseDensity, randomWalk, biasInstability, obj.sampleRatePressure);
            
            bias = obj.rant(-1.5,1.5); % mbar
            mbar = mbar + bias + noise; % add noise and bias
            mbar = mbar / resolution; % convert to LSB
            mbar = round(mbar); % round to integer
            mbar = mbar * resolution; % convert back to mbar
 
            % upsample the pressure to the sample rate of the MPU (filling with the last value)
            mbarOut                         = nan(length(alt),1);
            mbarOut(1:downsample_ratio:end) = mbar;
            mbar = mbarOut;
            % mbar                            = fillmissing(mbarOut, 'previous'); % the new jumpReadData function handles this better, therefore not needed anymore.
        end
        
        % ====================================================================== %
        %> @brief getSensorTemperature returns the sensor temperature.
        %>
        %> getSensorTemperature returns the sensor temperature with the actual
        %> sensor noise and bias.
        %>
        %> @param obj (JumpSensor) - the JUMP sensor
        %> @param tmp_true (nx1 double) - the true temperature in °C
        %>
        %> @retval degC (nx1 double) - the sensor temperature in °C
        % ====================================================================== %
        function degC = getSensorTemperature(obj, tmp_true)
            
            % input validation
            arguments
                obj;
                tmp_true (:,1) double {mustBeInRange(tmp_true, -20, 40)} = repmat(25, numel(alt),1); % temperature in °C
            end
 
            % downsample the temperature to the sample rate of the pressure sensor
            downsample_ratio = obj.sampleRate / obj.sampleRatePressure;
            tmp_true_ds      = downsample(tmp_true, downsample_ratio);
            
            % model the actual sensor
            resolution      = 0.01; % °C / LSB
            noiseDensity    = 0.0013; % °C / sqrt(Hz) (from allan variance experiments)
            biasInstability = 0.00158; % °C (from allan variance experiments)
            randomWalk      = 0.0005; % °C * sqrt(Hz) (from allan variance experiments)
 
            % generate noise
            noise = obj.generateNoise(length(tmp_true_ds), noiseDensity, randomWalk, biasInstability, obj.sampleRatePressure);
            
            bias = obj.rant(-0.8,0.8); % °C
            degC = tmp_true_ds + bias + noise; % add noise and bias
 
            % quantize the temperature
            degC = degC / resolution; % convert to LSB
            degC = round(degC); % round to integer
            degC = degC * resolution; % convert back to °C
 
            % upsample the temperature to the sample rate of the MPU (filling with the last value)
            degCOut                         = nan(length(tmp_true),1);
            degCOut(1:downsample_ratio:end) = degC;
            degC = degCOut;
            %degC                            = fillmissing(degCOut, 'previous'); % the new jumpReadData function handles this better, therefore not needed anymore.
        
        end
        
    end
    
 
    methods (Access = private, Static)
 
        % ====================================================================== %
        %> @brief getMisalignmentMatrix returns a misalignment matrix for the
        %> accelerometer or gyroscope.
        %>
        %> getMisalignmentMatrix returns a misalignment matrix for the
        %> accelerometer or gyroscope. The misalignment matrix is a 3x3 matrix
        %> which is used to rotate the sensor axes to the reference frame. The
        %> misalignment matrix is generated randomly.
        %>
        %> @param x (1x1 double) - the typical misalignment of the x-axis in percent
        %> @param y (1x1 double) - the typical misalignment of the y-axis in percent
        %> @param z (1x1 double) - the typical misalignment of the z-axis in percent
        %> @param sensitivity (1x1 double) - the overall sensitivity error in percent
        %>
        %> @retval M (3x3 double) - the misalignment matrix in percent
        %>
        %> @note the misalignment matrix is generated with the following steps:
        %> 1. generate the three axis in perfect alignment, but with sensitivity errors
        %> 2. generate the misalignment angles
        %> 3. generate the misalignment rotation
        %> 4. rotate the axes (attention, the rotx function works with degrees!)
        %> 5. create the misalignment matrix
        % ====================================================================== %
        function M = getMisalignmentMatrix(x, y, z, sensitivity)
            
            % check input arguments
            if nargin      == 1
               y            = x;
               z            = x;
               sensitivity  = 0;
            end
            
            % scale to scalars
            x           = x / 100;
            y           = y / 100;
            z           = z / 100;
            sensitivity = sensitivity / 100;
            
            % define the three axis in perfect alignment, but with sensitivity errors
            ax_x = [1;0;0] * JumpSensor.rant(1-sensitivity, 1+sensitivity);
            ax_y = [0;1;0] * JumpSensor.rant(1-sensitivity, 1+sensitivity);
            ax_z = [0;0;1] * JumpSensor.rant(1-sensitivity, 1+sensitivity);
            
            % generate the misalignment angles
            theta_x = JumpSensor.rant(0,asin(x) * 180/pi);
            theta_y = JumpSensor.rant(0,asin(y) * 180/pi);
            theta_z = JumpSensor.rant(0,asin(z) * 180/pi);
            
            % generate the misalignment rotation
            rot_x = rand() * 360;
            rot_y = rand() * 360;
            rot_z = rand() * 360;
            
            % rotate the axes (attention, the rotx function works with degrees!)
            ax_x = rotx(rot_x) * ( roty(theta_x) * ax_x );
            ax_y = roty(rot_y) * ( rotz(theta_y) * ax_y );
            ax_z = rotz(rot_z) * ( rotx(theta_z) * ax_z );
            
            % create the misalignment matrix
            M = [ax_x'; ax_y'; ax_z'] * 100;
            
        end
        
        % ====================================================================== %
        %> @brief rant returns a typical value in the sense of datasheet typical
        %> values.
        %>
        %> rant(lower) returns a typical random number with a symmetric
        %> distribution around zero, where the range is defined by the lower
        %> bound. The random number is gaussian distributed with lower been the
        %> two standard deviation range.
        %>
        %> rant(lower, upper) returns a typical random number with a gaussian
        %> distribution with the given range.
        %>
        %> rant(lower, upper, nstd) returns a typical random number with a
        %> gaussian distribution with nstd standard deviations.
        %>
        %> @param lower (1x1 double) - the lower bound of the distribution
        %> @param upper (1x1 double) - the upper bound of the distribution
        %> @param nstd (1x1 double) - the number of standard deviations
        %>
        %> @retval number (1x1 double) - the random number
        function number = rant(lower, upper, nstd)
            
            % check input arguments
            arguments
                lower (1,1) double;
                upper (1,1) double = -lower;
                nstd  (1,1) double = 2;
            end
 
            if lower > upper
                [lower, upper] = deal(upper, lower);
            end
            
            % gaussian distribution
            mu     = (lower + upper) / 2;
            sigma  = (upper - lower) / (2 * nstd);
            number = normrnd(mu, sigma);
            
        end
 
        % ====================================================================== %
        %> @brief generateNoise generates noise with the given parameters.
        %>
        %> generateNoise generates noise with the given parameters. The noise is
        %> generated with the following model:
        %> - white noise with the noise density N
        %> - random walk with the coefficient K
        %> - bias instability with the coefficient B
        %> The noise is generated with the sample rate fs.
        %>
        %> @param n (1x1 double) - the number of samples
        %> @param N (1x1 double) - the noise density
        %> @param K (1x1 double) - the random walk coefficient
        %> @param B (1x1 double) - the bias instability coefficient
        %> @param fs (1x1 double) - the sample rate
        %>
        %> @retval noise (nx1 double) - the generated noise
        % ====================================================================== %
        function noise = generateNoise(n,N,K,B,fs)
 
            % check input arguments
            arguments
                n (1,1) double {mustBeInteger, mustBePositive}; % number of samples
                N (1,1) double {mustBePositive}; % noise density
                K (1,1) double {mustBePositive}; % random walk
                B (1,1) double {mustBePositive}; % bias instability
                fs (1,1) double {mustBeInteger, mustBePositive}; % sample rate
                end
 
            % scale the coefficients
            N = N * sqrt(2); % single-sided noise
            K = K * sqrt(2); % single-sided noise
            B = B * sqrt(2); % single-sided noise
 
            % define the PSD of the Noise
            epsilon = 1e-6; % to avoid division by zero
            PSD     = @(f) N^2 + B^2 ./ (2*pi*f + epsilon) + K^2 ./ (2*pi*f + epsilon).^2;
 
            % frequency vector
            f = (0:(n-1)/2) * fs / n;
 
            % generate the noise
            white_noise = randn(n,1);
 
            % one-sided fft of the white noise
            Y = fft(white_noise);
 
            % set the DC component to zero
            Y(1) = 0;
 
            % adjust the amplitude according to the PSD
            Y(1:length(f)) = Y(1:length(f)) .* sqrt(PSD(f))';
 
            % if the signal length is even, the Nyquist frequency component is real
            if mod(n,2) == 0
            Y (n/2+1)    = real(Y(n/2+1));
            end
 
            % inverse fft
            noise = ifft(Y, 'symmetric');
 
        end
        
    end
end