diff --git a/commonFunctions/sensors/NoiseAnalysis/Functions/NoiseAnalysis_colored.m b/commonFunctions/sensors/NoiseAnalysis/Functions/NoiseAnalysis_colored.m new file mode 100644 index 0000000000000000000000000000000000000000..d525ae57e4b01e1abf43083b7b39d76fa16b018e --- /dev/null +++ b/commonFunctions/sensors/NoiseAnalysis/Functions/NoiseAnalysis_colored.m @@ -0,0 +1,106 @@ +function [sensor_vect] = NoiseAnalysis_colored(sensor_vect, sensor_number, track, plot_val) + +fprintf("Sensor analyzed: %s, track: %d and noise type: %s\n\n", sensor_vect(sensor_number).name, track, sensor_vect(sensor_number).noise_type); + +name = sensor_vect(sensor_number).name; +fs = sensor_vect(sensor_number).fs; + +b_left = sensor_vect(sensor_number).bounds(1); +b_right = sensor_vect(sensor_number).bounds(2); + +white_variance = sensor_vect(sensor_number).colored_data.white_variance; +fcut = sensor_vect(sensor_number).colored_data.fcut; +butterOrder = sensor_vect(sensor_number).colored_data.butterOrder; + + +%% Noise extraction +A = importdata(name).data; +signal = A(:,track); + +cut_left = ceil(length(signal)*b_left); +cut_right = ceil(length(signal)*b_right); +signal = signal(cut_left:cut_right); + +filtered_signal = movmean(signal, fs); + +% Plots +bound = 15; +y_min = min(min(signal), min(filtered_signal)); +y_max = max(max(signal), max(filtered_signal)); +if plot_val + figure(2) + subplot(1,2,1), hold on, grid on, title("Unfiltered signal"), xlabel("Samples [-]"), ylabel("Value"), ylim([y_min y_max]) + plot(signal) + subplot(1,2,2), hold on, grid on, title("Filtered signal"), xlabel("Samples [-]"), ylabel("Value"), ylim([y_min y_max]) + plot(filtered_signal(bound:end-bound)) % removed boundaries +end + +% Noise extraction +noise = signal(bound:end-bound) - filtered_signal(bound:end-bound); + +% fft of noise +tf = length(noise)/fs; +t = 0:1/fs:tf; +[G_n,f] = dft(noise,t); + +if plot_val + figure(3) + subplot(1,2,1), hold on, hold on, grid on, title("Real noise"), xlabel("f [Hz]"), ylabel("|fft(X)|") + plot(f,abs(G_n)) +end + + +%% Colored noise creator +white_noise = sqrt(white_variance)*randn(size(t)); + +[b, a] = butter(butterOrder, fcut, 'low'); +pink_noise = filter(b, a, white_noise); +noise_modeled = pink_noise; + + +%% Final plots +[G_test,f] = dft(noise_modeled,t); +if plot_val + figure(3) + subplot(1,2,2), hold on, hold on, grid on, title("Modeled noise"), xlabel("f [Hz]"), ylabel("|fft(X)|") + plot(f,abs(G_test)) +end + +y_min = min(min(noise), min(noise)); +y_max = max(max(noise_modeled), max(noise_modeled)); + +if plot_val + figure(4) + subplot(1,2,1), hold on, grid on, title("Real noise"), xlabel("Samples [-]"), ylabel("Value"), ylim([y_min y_max]) + plot(noise) + subplot(1,2,2), hold on, grid on, title("Modeled noise"), xlabel("Samples [-]"), ylabel("Value"), ylim([y_min y_max]) + plot(noise_modeled) +end + + +%% Data + +if ~plot_val + tracks = sensor_vect(sensor_number).tracks; + if tracks(1) >= 6 + track = track - (tracks(1) - 2); + end + switch track + case 2 + sensor_vect(sensor_number).track1.white_variance = white_variance*sensor_vect(sensor_number).factor; + sensor_vect(sensor_number).track1.fcut = fcut; + sensor_vect(sensor_number).track1.butterOrder = butterOrder; + case 3 + sensor_vect(sensor_number).track2.white_variance = white_variance*sensor_vect(sensor_number).factor; + sensor_vect(sensor_number).track2.fcut = fcut; + sensor_vect(sensor_number).track2.butterOrder = butterOrder; + case 4 + sensor_vect(sensor_number).track3.white_variance = white_variance*sensor_vect(sensor_number).factor; + sensor_vect(sensor_number).track3.fcut = fcut; + sensor_vect(sensor_number).track3.butterOrder = butterOrder; + end +end + + +end + diff --git a/commonFunctions/sensors/NoiseAnalysis/Functions/NoiseAnalysis_pink.m b/commonFunctions/sensors/NoiseAnalysis/Functions/NoiseAnalysis_pink.m index a1fd16d2606398e88a0f156b9db0a8e8169bf07f..9daf7a7fb24c94bca092a2397dfcb508decb409d 100644 --- a/commonFunctions/sensors/NoiseAnalysis/Functions/NoiseAnalysis_pink.m +++ b/commonFunctions/sensors/NoiseAnalysis/Functions/NoiseAnalysis_pink.m @@ -8,12 +8,13 @@ fs = sensor_vect(sensor_number).fs; b_left = sensor_vect(sensor_number).bounds(1); b_right = sensor_vect(sensor_number).bounds(2); +white_variance = sensor_vect(sensor_number).colored_data.white_variance; +fcut = sensor_vect(sensor_number).colored_data.fcut; +butterOrder = sensor_vect(sensor_number).colored_data.butterOrder; -%% Magic values + +%% Magic value temp_val = 0.1; -butterOrder = 1; -fcut = 0.3; -white_variance = 15; %% Noise extraction diff --git a/commonFunctions/sensors/NoiseAnalysis/main_noiseAnalysis.m b/commonFunctions/sensors/NoiseAnalysis/main_noiseAnalysis.m index 833acbee00b1a18bb7d1d8358dfc3089a42cc307..87c440573e25e3c8e3634f8b474de65accc91048 100644 --- a/commonFunctions/sensors/NoiseAnalysis/main_noiseAnalysis.m +++ b/commonFunctions/sensors/NoiseAnalysis/main_noiseAnalysis.m @@ -16,7 +16,12 @@ clear, clc close all sensor_single.name = "motor_Motor_CCPressureData.csv"; -sensor_single.noise_type = "white"; +sensor_single.noise_type = "colored"; +if strcmp(sensor_single.noise_type, "pink") || strcmp(sensor_single.noise_type, "colored") + sensor_single.colored_data.white_variance = 0.000025; + sensor_single.colored_data.fcut = 0.15; + sensor_single.colored_data.butterOrder = 1; +end sensor_single.track = 2; sensor_single.fs = 100; sensor_single.bound_left = 0.35; @@ -39,6 +44,8 @@ if sensor_single.noise_type == "white" sensor_vect = NoiseAnalysis_white(sensor_single, 1, sensor_single.track, 1); elseif sensor_single.noise_type == "pink" sensor_vect = NoiseAnalysis_pink(sensor_single, 1, sensor_single.track, 1); +elseif sensor_single.noise_type == "colored" + sensor_vect = NoiseAnalysis_colored(sensor_single, 1, sensor_single.track, 1); else error("Noise type is neither 'white' nor 'pink'") end @@ -64,8 +71,10 @@ for sensor_num = 1:length(Lyra_Port_sensor_vect) Lyra_Port_sensor_vect = NoiseAnalysis_white(Lyra_Port_sensor_vect, sensor_num, track, 0); elseif Lyra_Port_sensor_vect(sensor_num).noise_type == "pink" Lyra_Port_sensor_vect = NoiseAnalysis_pink(Lyra_Port_sensor_vect, sensor_num, track, 0); + elseif Lyra_Port_sensor_vect(sensor_num).noise_type == "colored" + Lyra_Port_sensor_vect = NoiseAnalysis_colored(Lyra_Port_sensor_vect, sensor_num, track, 0); else - error("Noise type is neither 'white' nor 'pink'") + error("Noise type is not 'white', 'pink' or 'colored'") end end end @@ -75,7 +84,7 @@ save(save_name, name) %% 3 - fs -name = "payload_Payload_StaticPressureData.csv"; +name = sensor_single.name; track = 1; % timestamp track timestamp = importdata(name).data; diff --git a/commonFunctions/sensors/NoiseAnalysis/mat_creator.m b/commonFunctions/sensors/NoiseAnalysis/mat_creator.m index c24e4ede4d0d2a379b0c8ac5f2ea1bdbfeb3ec2f..502d5e10c9f4a2f90e22e2f569a00355ac5783c4 100644 --- a/commonFunctions/sensors/NoiseAnalysis/mat_creator.m +++ b/commonFunctions/sensors/NoiseAnalysis/mat_creator.m @@ -3,61 +3,133 @@ clear, clc close all -Lyra_Port_sensor_vect(1).name = "main_Main_StaticPressureData1.csv"; -Lyra_Port_sensor_vect(1).info = "barometer1"; -Lyra_Port_sensor_vect(1).fs = 100; -Lyra_Port_sensor_vect(1).noise_type = "pink"; -Lyra_Port_sensor_vect(1).bounds = [0.5 0.6]; -Lyra_Port_sensor_vect(1).tracks = 2; -Lyra_Port_sensor_vect(1).factor = 0.01; - -Lyra_Port_sensor_vect(2).name = "main_Main_StaticPressureData2.csv"; -Lyra_Port_sensor_vect(2).info = "barometer2"; -Lyra_Port_sensor_vect(2).fs = 100; -Lyra_Port_sensor_vect(2).noise_type = "pink"; -Lyra_Port_sensor_vect(2).bounds = [0.6 0.7]; -Lyra_Port_sensor_vect(2).tracks = 2; -Lyra_Port_sensor_vect(2).factor = 0.01; - -Lyra_Port_sensor_vect(3).name = "main_Boardcore_LPS28DFWData.csv"; -Lyra_Port_sensor_vect(3).info = "barometer3"; -Lyra_Port_sensor_vect(3).fs = 50; -Lyra_Port_sensor_vect(3).noise_type = "white"; -Lyra_Port_sensor_vect(3).bounds = [0.13 0.23]; -Lyra_Port_sensor_vect(3).tracks = 2; -Lyra_Port_sensor_vect(3).factor = 0.01; - -Lyra_Port_sensor_vect(4).name = "main_Boardcore_LSM6DSRXData.csv"; -Lyra_Port_sensor_vect(4).info = "accelerometer"; -Lyra_Port_sensor_vect(4).fs = 440; -Lyra_Port_sensor_vect(4).noise_type = "white"; -Lyra_Port_sensor_vect(4).bounds = [0.35 0.5]; -Lyra_Port_sensor_vect(4).tracks = [2 3 4]; -Lyra_Port_sensor_vect(4).factor = 1000/9.81; - -Lyra_Port_sensor_vect(5).name = "main_Boardcore_LSM6DSRXData.csv"; -Lyra_Port_sensor_vect(5).info = "initial gyroscope"; -Lyra_Port_sensor_vect(5).fs = 440; -Lyra_Port_sensor_vect(5).noise_type = "white"; -Lyra_Port_sensor_vect(5).bounds = [0.35 0.5]; -Lyra_Port_sensor_vect(5).tracks = [6 7 8]; -Lyra_Port_sensor_vect(5).factor = 1000; - -Lyra_Port_sensor_vect(6).name = "motor_Motor_CCPressureData.csv"; -Lyra_Port_sensor_vect(6).info = "initial chamber pressure"; -Lyra_Port_sensor_vect(6).fs = 100; -Lyra_Port_sensor_vect(6).noise_type = "white"; -Lyra_Port_sensor_vect(6).bounds = [0.35 0.5]; -Lyra_Port_sensor_vect(6).tracks = 2; -Lyra_Port_sensor_vect(6).factor = 0.01; - -Lyra_Port_sensor_vect(7).name = "payload_Payload_StaticPressureData.csv"; -Lyra_Port_sensor_vect(7).info = "pitot (static + total)"; -Lyra_Port_sensor_vect(7).fs = 100; -Lyra_Port_sensor_vect(7).noise_type = "pink"; -Lyra_Port_sensor_vect(7).bounds = [0.3 0.47]; -Lyra_Port_sensor_vect(7).tracks = 2; -Lyra_Port_sensor_vect(7).factor = 0.01; +% .colored_data.white_variance = 15 +% .colored_data.fcut = 0.3 +% .colored_data.butterOrder = 1 + +ii = 1; + +% Ok, check factor +Lyra_Port_sensor_vect(ii).name = "main_Boardcore_LIS2MDLData.csv"; +Lyra_Port_sensor_vect(ii).info = "Mag"; +Lyra_Port_sensor_vect(ii).fs = 100; +Lyra_Port_sensor_vect(ii).noise_type = "pink"; +Lyra_Port_sensor_vect(ii).colored_data.white_variance = 0.00003; +Lyra_Port_sensor_vect(ii).colored_data.fcut = 0.3; +Lyra_Port_sensor_vect(ii).colored_data.butterOrder = 1; +Lyra_Port_sensor_vect(ii).bounds = [0.5 0.7]; +Lyra_Port_sensor_vect(ii).tracks = [2 3 4]; +Lyra_Port_sensor_vect(ii).factor = 1000; + +ii = ii + 1; + +Lyra_Port_sensor_vect(ii).name = "main_Main_StaticPressureData1.csv"; +Lyra_Port_sensor_vect(ii).info = "barometer1"; +Lyra_Port_sensor_vect(ii).fs = 100; +Lyra_Port_sensor_vect(ii).noise_type = "pink"; +Lyra_Port_sensor_vect(ii).colored_data.white_variance = 15; +Lyra_Port_sensor_vect(ii).colored_data.fcut = 0.3; +Lyra_Port_sensor_vect(ii).colored_data.butterOrder = 1; +Lyra_Port_sensor_vect(ii).bounds = [0.5 0.6]; +Lyra_Port_sensor_vect(ii).tracks = 2; +Lyra_Port_sensor_vect(ii).factor = 0.01; + +ii = ii + 1; + +Lyra_Port_sensor_vect(ii).name = "main_Main_StaticPressureData2.csv"; +Lyra_Port_sensor_vect(ii).info = "barometer2"; +Lyra_Port_sensor_vect(ii).fs = 100; +Lyra_Port_sensor_vect(ii).noise_type = "pink"; +Lyra_Port_sensor_vect(ii).colored_data.white_variance = 15; +Lyra_Port_sensor_vect(ii).colored_data.fcut = 0.3; +Lyra_Port_sensor_vect(ii).colored_data.butterOrder = 1; +Lyra_Port_sensor_vect(ii).bounds = [0.6 0.7]; +Lyra_Port_sensor_vect(ii).tracks = 2; +Lyra_Port_sensor_vect(ii).factor = 0.01; + +ii = ii + 1; + +% Ok, check factor +Lyra_Port_sensor_vect(ii).name = "payload_Payload_StaticPressureData.csv"; +Lyra_Port_sensor_vect(ii).info = "StaticPressure"; +Lyra_Port_sensor_vect(ii).fs = 100; +Lyra_Port_sensor_vect(ii).noise_type = "pink"; +Lyra_Port_sensor_vect(ii).colored_data.white_variance = 30; +Lyra_Port_sensor_vect(ii).colored_data.fcut = 0.3; +Lyra_Port_sensor_vect(ii).colored_data.butterOrder = 1; +Lyra_Port_sensor_vect(ii).bounds = [0.3 0.4]; +Lyra_Port_sensor_vect(ii).tracks = 2; +Lyra_Port_sensor_vect(ii).factor = 0.01; + +ii = ii + 1; + +Lyra_Port_sensor_vect(ii).name = "payload_Payload_DynamicPressureData.csv"; +Lyra_Port_sensor_vect(ii).info = "DynamicPressure"; +Lyra_Port_sensor_vect(ii).fs = 100; +Lyra_Port_sensor_vect(ii).noise_type = "colored"; +Lyra_Port_sensor_vect(ii).colored_data.white_variance = 12; +Lyra_Port_sensor_vect(ii).colored_data.fcut = 0.05; +Lyra_Port_sensor_vect(ii).colored_data.butterOrder = 1; +Lyra_Port_sensor_vect(ii).bounds = [0.06 0.14]; +Lyra_Port_sensor_vect(ii).tracks = 2; +Lyra_Port_sensor_vect(ii).factor = 0.01; + +ii = ii + 1; + +Lyra_Port_sensor_vect(ii).name = "main_Boardcore_LPS28DFWData.csv"; +Lyra_Port_sensor_vect(ii).info = "barometer3"; +Lyra_Port_sensor_vect(ii).fs = 50; +Lyra_Port_sensor_vect(ii).noise_type = "white"; +Lyra_Port_sensor_vect(ii).bounds = [0.13 0.23]; +Lyra_Port_sensor_vect(ii).tracks = 2; +Lyra_Port_sensor_vect(ii).factor = 0.01; + +ii = ii + 1; + +Lyra_Port_sensor_vect(ii).name = "main_Boardcore_LSM6DSRXData.csv"; +Lyra_Port_sensor_vect(ii).info = "accelerometer"; +Lyra_Port_sensor_vect(ii).fs = 440; +Lyra_Port_sensor_vect(ii).noise_type = "white"; +Lyra_Port_sensor_vect(ii).bounds = [0.35 0.5]; +Lyra_Port_sensor_vect(ii).tracks = [2 3 4]; +Lyra_Port_sensor_vect(ii).factor = 1000/9.81; + +ii = ii + 1; + +Lyra_Port_sensor_vect(ii).name = "main_Boardcore_LSM6DSRXData.csv"; +Lyra_Port_sensor_vect(ii).info = "initial gyroscope"; +Lyra_Port_sensor_vect(ii).fs = 440; +Lyra_Port_sensor_vect(ii).noise_type = "white"; +Lyra_Port_sensor_vect(ii).bounds = [0.35 0.5]; +Lyra_Port_sensor_vect(ii).tracks = [6 7 8]; +Lyra_Port_sensor_vect(ii).factor = 1000; + +ii = ii + 1; + +% Ok, check factor +Lyra_Port_sensor_vect(ii).name = "motor_Motor_CCPressureData.csv"; +Lyra_Port_sensor_vect(ii).info = "initial chamber pressure"; +Lyra_Port_sensor_vect(ii).fs = 100; +Lyra_Port_sensor_vect(ii).noise_type = "colored"; +Lyra_Port_sensor_vect(ii).colored_data.white_variance = 0.000025; +Lyra_Port_sensor_vect(ii).colored_data.fcut = 0.15; +Lyra_Port_sensor_vect(ii).colored_data.butterOrder = 1; +Lyra_Port_sensor_vect(ii).bounds = [0.35 0.5]; +Lyra_Port_sensor_vect(ii).tracks = 2; +Lyra_Port_sensor_vect(ii).factor = 0.01; + +ii = ii + 1; + +Lyra_Port_sensor_vect(ii).name = "payload_Payload_StaticPressureData.csv"; +Lyra_Port_sensor_vect(ii).info = "pitot (static + total)"; +Lyra_Port_sensor_vect(ii).fs = 100; +Lyra_Port_sensor_vect(ii).noise_type = "pink"; +Lyra_Port_sensor_vect(ii).colored_data.white_variance = 15; +Lyra_Port_sensor_vect(ii).colored_data.fcut = 0.3; +Lyra_Port_sensor_vect(ii).colored_data.butterOrder = 1; +Lyra_Port_sensor_vect(ii).bounds = [0.3 0.47]; +Lyra_Port_sensor_vect(ii).tracks = 2; +Lyra_Port_sensor_vect(ii).factor = 0.01; save("Lyra_Port_sensor_vect") diff --git a/commonFunctions/sensors/Sensor1D.m b/commonFunctions/sensors/Sensor1D.m index 8470de04db30177b2c4bf318142938a329bb8087..8c7f22615983149811af9f9c854f647cd8698bd4 100644 --- a/commonFunctions/sensors/Sensor1D.m +++ b/commonFunctions/sensors/Sensor1D.m @@ -20,6 +20,7 @@ classdef Sensor1D < handle noiseType; % White (default) or Pink noiseDataTrack1; noiseFactor; + colored_opts; noiseVariance; % Defining gaussian white noise % offset @@ -98,6 +99,13 @@ classdef Sensor1D < handle else error("Sensor not defined") end + + if strcmp(obj.noiseType, "colored") + obj.colored_opts.white_variance = vect(ii).colored_data.white_variance; + obj.colored_opts.fcut = vect(ii).colored_data.fcut; + obj.colored_opts.butterOrder = vect(ii).colored_data.butterOrder; + obj.colored_opts.filterStatus = 0; + end else if strcmp("Sensor1D", class(obj)) || strcmp("SensorFault", class(obj)) obj.noiseDataTrack1 = []; @@ -149,6 +157,10 @@ classdef Sensor1D < handle inputArg = inputArg + obj.noiseDataTrack1.peaks_vect_val(ii)*obj.noiseFactor*sin(2*pi*obj.noiseDataTrack1.peaks_vect_f(ii)*t + randn(1)); end inputArg = inputArg + sqrt(obj.noiseDataTrack1.variance*obj.noiseFactor).*randn(length(inputArg),1); + elseif strcmp(obj.noiseType, "colored") + inputArg = inputArg + sqrt(obj.colored_opts.white_variance*obj.noiseFactor).*randn(length(inputArg),1); + [b, a] = butter(obj.colored_opts.butterOrder, obj.colored_opts.fcut, 'low'); + [inputArg, obj.colored_opts.filterStatus] = filter(b, a, inputArg, obj.colored_opts.filterStatus); else error("This noise is not defined") end diff --git a/commonFunctions/sensors/Sensor3D.m b/commonFunctions/sensors/Sensor3D.m index 14955ff159d208fe8956e981bd03b7b33916d497..bada2453af6fecca427344657d9d96c789eb80d6 100644 --- a/commonFunctions/sensors/Sensor3D.m +++ b/commonFunctions/sensors/Sensor3D.m @@ -142,6 +142,8 @@ classdef Sensor3D < Sensor1D inputArg1 = inputArg1 + sqrt(obj.noiseDataTrack1.variance*obj.noiseFactor).*randn(length(inputArg1),1); inputArg2 = inputArg2 + sqrt(obj.noiseDataTrack2.variance*obj.noiseFactor).*randn(length(inputArg2),1); inputArg3 = inputArg3 + sqrt(obj.noiseDataTrack3.variance*obj.noiseFactor).*randn(length(inputArg3),1); + elseif strcmp(obj.noiseType, "colored") + error("Yet to be implemented") else error("This noise is not defined") end diff --git a/data/2024_Lyra_Portugal_October/Lyra_Port_sensor_vect_res.mat b/data/2024_Lyra_Portugal_October/Lyra_Port_sensor_vect_res.mat index c315505e15c35d9721647dc0d0bbb9264d6f17f5..2c3b98ad9def8becb8ced991dfe7f38b560308d0 100644 Binary files a/data/2024_Lyra_Portugal_October/Lyra_Port_sensor_vect_res.mat and b/data/2024_Lyra_Portugal_October/Lyra_Port_sensor_vect_res.mat differ diff --git a/data/2024_Lyra_Portugal_October/initSensors2024_Lyra_Portugal_October.m b/data/2024_Lyra_Portugal_October/initSensors2024_Lyra_Portugal_October.m index b36ed6b6f36e8067a5659f97f9f9cfb080a6e446..303eb0421c5dd69ba03fbc9f506b4228ecbea70e 100644 --- a/data/2024_Lyra_Portugal_October/initSensors2024_Lyra_Portugal_October.m +++ b/data/2024_Lyra_Portugal_October/initSensors2024_Lyra_Portugal_October.m @@ -155,7 +155,7 @@ sensorSettings.comb_chamber.minMeasurementRange = 0; sensorSettings.comb_chamber.resolution = 1; % random value stolen from baro sensorSettings.comb_chamber.offset = 0; -sensorSettings.comb_chamber.update(Lyra_Port_sensor_vect, "motor_Motor_TopTankPressureData.csv", 1); +sensorSettings.comb_chamber.update(Lyra_Port_sensor_vect, "motor_Motor_CCPressureData.csv", 1); %% pitot @@ -171,7 +171,7 @@ sensorSettings.pitot_total = Sensor1D(); sensorSettings.pitot_total.maxMeasurementRange = 2*1034.21; % mbar (30 psi) sensorSettings.pitot_total.minMeasurementRange = 0; sensorSettings.pitot_total.bit = 12; -sensorSettings.pitot_total.update(Lyra_Port_sensor_vect, "payload_Payload_StaticPressureData.csv", 1); +sensorSettings.pitot_total.update(Lyra_Port_sensor_vect, "payload_Payload_DynamicPressureData.csv", 1); %% total sensor initialization %