% Table of Contents
% Stationary conditions
%     !Read Data
%     Chapter 4 - test stand description
%     Chapter 5 - Analysis of measurement system
%         Chapter 5.2 Raw, angle, IAS for PSG
%         Chapter 5.3.1 + 5.3.2 - Measurement accuracy
%         Chapter 5.3.3 - Analysis of Operating Conditions
%         Chapter 5.4 - Tooth-specific signature
%         Chapter 5.5 - Optimisation
%             Influence of fit ellipses
%     Chapter 6 - Damage Detection
%         Generate and save features
%         Chapter 6.2 - at stationary conditions
%         Chapter 6.2.3 - Classification
%         Chapter 6.2.4 - Evaluation: Comparison with other Sensors
%         Chapter 6.5 - Optimisation
%         Chapter 6 - SEW fatigue damage test
%     Chapter 7 - Potentials for Torque, TE, Absolute ... Measurement
%         Absolute position using tooth-specific signature
%         Torque measurement
%         Severity Assessment
% Transient conditions
%     !Read data
%         Chapter 5.3.4 - Investigation of speed and torque ramps
%         CutSpeedRamp and perform equi-angle sampling
%         Chapter 6 - Generate and save feature
%         Chapter 6.6 - Classification at transient conditions
%     Chapter 7.2 - SEW tests at torque ramps
% Functions
%% Stationary conditions
%% !Read Data
% Execute before start

i = 1; %read first test
basefolder = 'C:\Users\sn85zymy\Desktop\Diss'; 
%get file names and extract meta data
dataListing_PSG_GW = getDataListing([basefolder '\Testdata_GSP\stationary\1_GearWheelSensor']);
GW_PSG = AngleSensor_PSG(dataListing_PSG_GW(i).folder, dataListing_PSG_GW(i).name, "Speed"); 

dataListing_PSG_PR = getDataListing([basefolder '\Testdata_GSP\stationary\3_PoleRingSensor']);
PR_PSG = AngleSensor_PSG(dataListing_PSG_PR(i).folder, dataListing_PSG_PR(i).name, "AL 795"); 

dataListing_PSG_Acc = getDataListing([basefolder '\Testdata_GSP\stationary\6_Acceleration']);
% Acc_PSG = AccelerationSensor_PSG(dataListing_PSG_Acc(i).folder, dataListing_PSG_Acc(i).name, "top")

% dataListing_PSG_GW_InFact = getDataListing([basefolder '\Testdata_GSP\InfluencingFactors\1_GearWheelSensor']);

dataListing_SEW_Z3 = getDataListing([basefolder '\Testdata_SEW\MR_Z3']);
% GW_SEW_Z3 = AngleSensor_SEW(dataListing_SEW_Z3(i).folder, dataListing_SEW_Z3(i).name, "Z3_SEW");
dataListing_SEW_Z1 = getDataListing([basefolder '\Testdata_SEW\MR_Z1']);
% GW_SEW_Z1 = AngleSensor_SEW(dataListing_SEW_Z1(i).folder, dataListing_SEW_Z1(i).name, "Z1_SEW");
dataListing_SEW_B = getDataListing([basefolder '\Testdata_SEW\MR_B']);

dataListing_SEW_MR3 = getDataListing([basefolder '\Testdata_SEW\MR_Pos3']);
dataListing_SEW_MR4 = getDataListing([basefolder '\Testdata_SEW\MR_Pos4']);

dataListing_SEW_Acc_HSS = getDataListing([basefolder '\Testdata_SEW\Acc_HSS']);
dataListing_SEW_Acc_LSS = getDataListing([basefolder '\Testdata_SEW\Acc_LSS']);

%
%% Chapter 4 - test stand description

%Plot idealised sensor signal
T = 0.00001; 
t = 0 : T : 0.02; 
sensor1 = sin(1500*t);
sensor2 = cos(1500*t);
angle = atan2(sensor1, sensor2); 
angle_unwrapped = unwrap(angle);
angle_scaled = angle_unwrapped*360/(2*pi()*51);
IAS = (diff(angle_unwrapped)/T)/(2*pi()*51)*60;
%Raw
figure(); title('Idealised: Raw - Angle - Speed'); 
plot(t,sensor1); hold on 
plot(t, sensor2)
xlabel('Time in s')
ylabel('Sensorsignal normalised') 
savePlotWithCustomSettings('Ch4-1_SineCosineExample', 0.48, 0.35, 0)
%Angle & IAS

%

figure(), set(gcf, 'defaultAxesTickLabelInterpreter', 'latex');
plot(t, angle); hold on
plot(t, angle_unwrapped)
ylabel('Angle in rad'), yticklabels(strrep(yticklabels,'-','$-$'));
yyaxis right
plot(t(1:end-1), IAS)
ylabel('IAS in rpm'), ylim([250 300]), xlabel('Time in s');  
legend('$\varphi_{atan2}$', '$\varphi_{Gear wheel}$', 'IAS', 'location', 'northwest', 'Interpreter', 'latex')
savePlotWithCustomSettings('Ch4-1_SignalProcessing', 0.7, 0.4, -0.0003)

%
%% Chapter 5 - Analysis of measurement system
% The Aim of this section is to investigate influences on the measurement system. First, possible influencing factors are identified through modelling of the sensor system, which are then investigated using the test data (mainly GSP tests)
%

ende = 420;
t = t(1:ende); sensor1 =sensor1(1:ende);
figure()
AmplitudeShift = 0.5*sin(1500*t); plot(t, AmplitudeShift); hold on
PhaseShift = 0.9*sin(1500*t*0.75); plot(t, PhaseShift)
plot(t, sensor1, Color='black'), xlim([0 0.0045])
xlabel('Time in ms'), ylabel('Amplitude normalised')
savePlotWithCustomSettings('Ch5-1_AmplitudePhaseShift', 0.48, 0.35, 0)

%
%% Chapter 5.2 Raw, angle, IAS for PSG

start = 1; 
ende = 30000; 
GW_PSG.FilterMethod = 'none';

figure()
plot(GW_PSG.t(start:ende), GW_PSG.sin(start:ende),'DisplayName', 'sine'), hold on  
plot(GW_PSG.t(start:ende), GW_PSG.cos(start:ende), 'DisplayName', 'cosine'), legend, xlabel('Time in s'), ylabel('Voltage in V');
xlim([0 0.05]), ylim([2 3.05])
savePlotWithCustomSettings('Ch5-2_RawSignal005', 0.5, 0.4, 0.002)
figure()
plot(GW_PSG.t(start:ende), GW_PSG.sin(start:ende),'DisplayName', 'sine'), hold on  
plot(GW_PSG.t(start:ende), GW_PSG.cos(start:ende), 'DisplayName', 'cosine'), legend, xlabel('Time in s'), ylabel('Voltage in V');
xlim([0.02 0.028]), ylim([2 3.05])
savePlotWithCustomSettings('Ch5-2_RawSignal002', 0.5, 0.4, 0.002)
figure()
plot(GW_PSG.t(start:ende), GW_PSG.sin(start:ende),'DisplayName', 'sine'), hold on  
plot(GW_PSG.t(start:ende), GW_PSG.cos(start:ende), 'DisplayName', 'cosine')
xlim([0.02295 0.024]), ylim([2.82 3.05])
savePlotWithCustomSettings('Ch5-2_RawSignalZoom', 0.3, 0.2, 0.002)

%

GW_PSG.FilterMethod = 'filter';
GW_PSG = GW_PSG.processSinCos(["raw", "preprocess"]);
start = 1; 
ende = 30000; 
figure()
subplot(3,1,1)
plot(GW_PSG.t(start:ende), GW_PSG.sin(start:ende),'DisplayName', 'sine'), hold on
plot(GW_PSG.t(start:ende), GW_PSG.cos(start:ende), 'DisplayName', 'cosine'), legend, ylabel('Voltage in V');
subplot(3,1,2)
plot(GW_PSG.t(start:ende), GW_PSG.angle_deviation(start:ende)), ylabel('Angle in rad'), ylim([-0.0055 0.007])
subplot(3,1,3)
plot(GW_PSG.t(start:ende), GW_PSG.IAS(start:ende)), xlabel('Time in s'), ylabel('IAS in rpm'), ylim([350 900])
savePlotWithCustomSettings('Ch5-2_Raw_Angle_IAS', 0.5, 0.5, 0.005)
start = 15000; 
ende = 17000; 
figure()
subplot(3,1,1)
plot(GW_PSG.t(start:ende), GW_PSG.sin(start:ende),'DisplayName', 'sine'), hold on
plot(GW_PSG.t(start:ende), GW_PSG.cos(start:ende), 'DisplayName', 'cosine'), legend
subplot(3,1,2)
plot(GW_PSG.t(start:ende), GW_PSG.angle_deviation(start:ende))
subplot(3,1,3)
plot(GW_PSG.t(start:ende), GW_PSG.IAS(start:ende)), xlabel('Time in s'), ylim([400 700])
savePlotWithCustomSettings('Ch5-2_Raw_Angle_IAS_zoom', 0.5, 0.5, 0)

%

figure()
[p,f] = fft_data(GW_PSG.IAS, GW_PSG.fs); f = f/GW_PSG.quant.f1; semilogy(f,p), xlim([0 200]), xlabel('frequency order'), ylabel('Amplitude in V')

% for SEW

disp('SEW: Raw - Angle - Speed')
GW_SEW_Z3.FilterMethod = 'filter';
figure(), hold on 
plot(GW_SEW_Z3.t(start:ende), GW_SEW_Z3.sin(start:ende),'DisplayName', 'sine'), 
plot(GW_SEW_Z3.t(start:ende), GW_SEW_Z3.cos(start:ende), 'DisplayName', 'cosine'), legend, xlabel('Time in s'), ylabel('Voltage in V');

figure()
subplot(3,1,1)
plot(GW_SEW_Z3.t(start:ende), GW_SEW_Z3.sin(start:ende),'DisplayName', 'sine'), hold on
plot(GW_SEW_Z3.t(start:ende), GW_SEW_Z3.cos(start:ende), 'DisplayName', 'cosine'), legend, ylabel('Voltage in V');
subplot(3,1,2)
plot(GW_SEW_Z3.t(start:ende), GW_SEW_Z3.angle_deviation(start:ende)), ylabel('Angle in rad')
subplot(3,1,3)
plot(GW_SEW_Z3.t(start:ende), GW_SEW_Z3.IAS(start:ende)), xlabel('Time in s'), ylabel('IAS in rpm')

figure()
[p, f]  = fft_data(GW_SEW_Z3.IAS, GW_SEW_Z3.fs);  f = f/GW_SEW_Z3.quant.fC; semilogy(f, p), xlim([0 200]), xlabel('frequency order'), ylabel('Amplitude in V')

%
%% Chapter 5.3.1 + 5.3.2 - Measurement accuracy
% Plot Errorbar of Peaks for multiple Revolutions at GSP

GW_PSG = processSinCos(GW_PSG, ["raw", "preprocess"]);
Name = ["1st repetition" "2nd repetition"];
figure() %for repeatability
i = 1;
i_plot = 1;
while i_plot < 3 % step : length(dataListing_PSG_GW)
    %read sensor data
    GW_PSG = AngleSensor_PSG(dataListing_PSG_GW(i).folder, dataListing_PSG_GW(i).name, "Speed"); 
    if GW_PSG.Trq == 20 && GW_PSG.Rev == 1600
        % Name = ['Pinion' num2str(GW_PSG.Pinion) ' Dmg' num2str(GW_PSG.Damage) ' ' num2str(GW_PSG.Trq) 'Nm ' num2str(GW_PSG.Rev) 'rpm, ' num2str(GW_PSG.Temperature_mean) '°C'];
        peaks = get_peaks(GW_PSG, "raw");
        errorbar(peaks.sin.high_mean, 3*peaks.sin.high_std, 'DisplayName', Name(i_plot))
        standard_deviation = 3*mean(peaks.sin.high_std)
        mininum = min(peaks.sin.high_mean)
        maximum = max(peaks.sin.high_mean)
        hold on
        i_plot = i_plot + 1; 
    end
    i = i + 1;
end
xlabel('Teeth'), xlim([0 51]), ylabel('Voltage in V'), legend
savePlotWithCustomSettings('Ch5-3-1_repeatability', 0.7, 0.5, 0.9)

%

GW_PSG_Ref = AngleSensor_PSG(dataListing_PSG_GW(1).folder, dataListing_PSG_GW(1).name, "Speed");
figure() %for different torques
subplot(1,2,1)
Trq2plot = [55, 20, -20, -55]; i = 1; i_plot = 1; peaks_temp = []; peaks_mean = []; 
while i_plot < 5 % step : length(dataListing_PSG_GW)
    %read sensor data
    GW_PSG = AngleSensor_PSG(dataListing_PSG_GW(i).folder, dataListing_PSG_GW(i).name, "Speed"); 
    if GW_PSG.Trq == Trq2plot(i_plot) && GW_PSG.Rev == -1600
        GW_PSG = shiftSignals(GW_PSG, GW_PSG_Ref);
        Name = ['Torque = ' strrep(num2str(GW_PSG.Trq), '-', '$-$'), 'Nm'];
        peaks = get_peaks(GW_PSG, "shifted");
        % errorbar(peaks.sin.high_mean, 3*peaks.sin.high_std, 'DisplayName', Name)
        plot(peaks.cos.high_mean, 'DisplayName', Name)
        if i_plot == 3 
            peaks_temp = []; 
        end
        peaks_temp = [peaks_temp, peaks.cos.high]; 
        peaks_mean = [peaks_mean, peaks.cos.high_mean];
        standardDeviation = mean(3*std(peaks_temp, 0, 2))
        hold on
        i_plot = i_plot + 1; 
    end
    i = i + 1;
end
xlabel('Teeth'), ylabel('Voltage in V'), xlim([0 51]),title(['Speed = ' num2str(GW_PSG.Rev) 'rpm']), legend('location', 'northoutside', 'Interpreter', 'latex')

subplot(1, 2,2) %for different speeds
Rev2plot = [1600, 1300, -1300, -1600]; 
i = 1; i_plot = 1; peaks_temp =[]; peaks_mean = [];
while i_plot < 5 % step : length(dataListing_PSG_GW)
    %read sensor data
    GW_PSG = AngleSensor_PSG(dataListing_PSG_GW(i).folder, dataListing_PSG_GW(i).name, "Speed"); 
    if GW_PSG.Trq == -20 && GW_PSG.Rev == Rev2plot(i_plot)
        GW_PSG = shiftSignals(GW_PSG, GW_PSG_Ref);
        Name = ['Speed = ' strrep(num2str(GW_PSG.Rev), '-', '$-$'), 'rpm'];
        peaks = get_peaks(GW_PSG, "shifted");
        % errorbar(peaks.sin.high_mean, 3*peaks.sin.high_std, 'DisplayName', Name)
        plot(peaks.cos.high_mean, 'DisplayName', Name)
        if i_plot == 3 
            peaks_temp = []; 
        end
        peaks_temp = [peaks_temp, peaks.cos.high]; 
        % peaks_mmean = [peaks_mean, peaks.cos.high_mean];
        standardDeviation = mean(3*std(peaks_temp, 0, 2))
        hold on
        i_plot = i_plot + 1; 
    end
    i = i + 1;
end
xlabel('Teeth'), xlim([0 51]), title(['Torque = ' num2str(GW_PSG.Trq) 'Nm']), legend('location', 'northoutside', 'Interpreter', 'latex')
savePlotWithCustomSettings('Ch5-3-2_Reproducibility', 1, 0.6, 0)

%
% SEW

step = 100; 
GW_SEW_Ref = AngleSensor_SEW(dataListing_SEW_Z3(1).folder, dataListing_SEW_Z3(1).name, "Z3_SEW");
figure()
for i = 1 : 10 % step : length(dataListing_PSG_GW)
    %read sensor data
    GW_SEW_Z3 = AngleSensor_SEW(dataListing_SEW_Z3(i).folder, dataListing_SEW_Z3(i).name, "Z3_SEW"); 
    GW_SEW_Z3 = shiftSignals(GW_SEW_Z3, GW_SEW_Ref); 
    Name = [datetime(GW_SEW_Z3.Date)];
    peaks = get_peaks(GW_SEW_Z3, "shifted");
    errorbar(peaks.sin.high_mean, 3*peaks.sin.high_std, 'DisplayName', Name)
    hold on
end
xlabel('Teeth'), xlim([0 80]), ylabel('Voltage in V'), legend, title('Errorbar of different Peaks Sin High')

%
% Success rate of shift using tooth-specific signature with different references

clear GW_PSG_mat
GW_PSG_Ref_neg = AngleSensor_PSG(dataListing_PSG_GW(1).folder, dataListing_PSG_GW(10).name, "Speed");
GW_PSG_Ref_pos = AngleSensor_PSG(dataListing_PSG_GW(1).folder, dataListing_PSG_GW(77).name, "Speed");

idx_neg = 1; idx_pos = 1; 
for i = 1 : length(dataListing_PSG_GW)
    GW_PSG = AngleSensor_PSG(dataListing_PSG_GW(i).folder, dataListing_PSG_GW(i).name, "Speed"); 
    if GW_PSG.Pinion == 0
        GW_PSG = shiftSignals(GW_PSG, GW_PSG_Ref_neg);
        if GW_PSG.phase_shift_success == 0
            GW_PSG_mat{idx_neg,1} = GW_PSG; 
            idx_neg = idx_neg+1;
        end
        GW_PSG = shiftSignals(GW_PSG, GW_PSG_Ref_pos);
        if GW_PSG.phase_shift_success == 0
            GW_PSG_mat{idx_pos,2} = GW_PSG;
            idx_pos = idx_pos +1;
        end

    end
end
%plot
for i = 1 : size(GW_PSG_mat, 1)
    if ~isempty(GW_PSG_mat{i,1})
        Rev_pos(i) = GW_PSG_mat{i,1}.Rev;
        Trq_pos(i) = GW_PSG_mat{i,1}.Trq;
    end
    if ~isempty(GW_PSG_mat{i,2})
        Rev_neg(i) = GW_PSG_mat{i,2}.Rev;
        Trq_neg(i) = GW_PSG_mat{i,2}.Trq;
    end
end
figure(), plot(Rev_pos, 'DisplayName', '+2000 rpm'), hold on, plot(Rev_neg, 'DisplayName', '-2000 rpm'), ylabel('Speed in rpm'), legend

%
%% Chapter 5.3.3 - Analysis of Operating Conditions
% Create Table with Peak, OC data

t1 = datetime;
idx = 1; 
for i = 1 : length(dataListing_PSG_GW)
    GW_PSG = AngleSensor_PSG(dataListing_PSG_GW(i).folder, dataListing_PSG_GW(i).name, "Speed"); 
%     if GW_PSG.Rev == 1300
        peaks_temp = get_peaks(GW_PSG, "raw"); 
        peaks{idx} = peaks_temp; 
        Torque(idx) = GW_PSG.Trq; 
        Speed(idx) = GW_PSG.Rev;
        Temperature(idx) = GW_PSG.Temperature_mean;
        Pinion(idx) = GW_PSG.Pinion;
        Quadrant(idx) = GW_PSG.Quadrant;
        idx = idx+1; 
%     end
end
varNames = {'Peaks', 'Torque', 'Speed', 'Temperature', 'Pinion', 'Quadrant'};
varTypes = {'struct', 'double', 'double', 'double', 'double', 'double'};
T_data = table('Size', [length(Pinion), numel(varNames)], 'VariableTypes', varTypes, 'VariableNames', varNames);

T_data.Peaks = peaks';
T_data.Torque = Torque';
T_data.Speed = Speed';
T_data.Temperature = Temperature';
T_data.Pinion = Pinion';
T_data.Quadrant = Quadrant';
for i = 1 : height(T_data)
    peaks_temp = T_data.Peaks{i}; 
    sin_high_mean(i) = mean(peaks_temp.sin.high_mean); sin_low_mean(i) = mean(peaks_temp.sin.low_mean); 
    cos_high_mean(i) = mean(peaks_temp.cos.high_mean); cos_low_mean(i) = mean(peaks_temp.cos.low_mean);
end
T_data.("sine high") = sin_high_mean'; T_data.("sine low") = sin_low_mean';
T_data.("cosine high") = cos_high_mean'; T_data.("cosine low") = cos_low_mean';
datetime - t1

%
% Correlation analysis at artificial damage tests

T_data_temp = T_data(:, :);
variables2read = ["Torque", "Speed", "Temperature", "sine high", "cosine high"];
clear idx
idx.Pinion0 = find(T_data_temp.Pinion == 0);
idx.Pinion123 = [find(T_data_temp.Pinion == 1); find(T_data_temp.Pinion == 1); find(T_data_temp.Pinion == 1)];
idx.Pinion56 = [find(T_data_temp.Pinion == 5); find(T_data_temp.Pinion == 6)]; 
idx.Q1 = find(T_data_temp.Quadrant == 1); idx.Q2 = find(T_data_temp.Quadrant == 2); idx.Q3 = find(T_data_temp.Quadrant == 3); idx.Q4 = find(T_data_temp.Quadrant == 4);
idx.Speed300 = T_data_temp.Speed == 300; idx.Speed1300 = T_data_temp.Speed == 1300; idx.Speed2000 = T_data_temp.Speed == 2000;

corr_table = T_data_temp(:, variables2read);
corr_matrix = round(corr(corr_table{:,:}),2);
fig = figure();
subplot(1,2,1)
heatmap(corr_table.Properties.VariableNames, corr_table.Properties.VariableNames, corr_matrix);
title("All Quadrants")

corr_table = T_data_temp(idx.Q4, variables2read);
corr_matrix = round(corr(corr_table{:,:}),2);
subplot(1,2,2)
heatmap(corr_table.Properties.VariableNames, corr_table.Properties.VariableNames, corr_matrix)
title("Quadrant 4")
savePlotWithCustomSettings('Ch5-3-3_CorrArtificialDamages', 1.02, 0.35, 0)

%separate Pinions
fn_names = {"Q1", "Q2", "Q3", "Q4"};
figure()
for i = 1 : length(fn_names)
    corr_table = T_data_temp(idx.(fn_names{i}), variables2read);
    corr_matrix = round(corr(corr_table{:,:}),2);
    subplot(2,2,i)
    heatmap(corr_table.Properties.VariableNames, corr_table.Properties.VariableNames, corr_matrix)
    title(fn_names{i})
end

%correlation with 300 & 1300
corr_table = T_data_temp(idx.Speed300|idx.Speed1300, variables2read);
corr_matrix = round(corr(corr_table{:,:}),2);
figure()
heatmap(corr_table.Properties.VariableNames, corr_table.Properties.VariableNames, corr_matrix)
title('-300 | -1300')
%correlation with 1300 & 2000
corr_table = T_data_temp(idx.Speed1300|idx.Speed2000, variables2read);
corr_matrix = round(corr(corr_table{:,:}),2);
figure()
heatmap(corr_table.Properties.VariableNames, corr_table.Properties.VariableNames, corr_matrix)
title('-1300 | -2000')
%correlation with 1300 & 2000
corr_table = T_data_temp(idx.Speed300|idx.Speed2000, variables2read);
corr_matrix = round(corr(corr_table{:,:}),2);
figure()
heatmap(corr_table.Properties.VariableNames, corr_table.Properties.VariableNames, corr_matrix)
title('-300 | -2000')

%
% Random Test Evaluation
% Evaluation based on the additional tests on the influencing factors. This is aimed to reduce the influence of test design based temperature fluctuations and analyse this influence

GW_PSG_Ref = AngleSensor_PSG(dataListing_PSG_GW_InfFact(26).folder, dataListing_PSG_GW_InfFact(26).name, "Speed");

t1 = datetime
idx = 1; 
for i = 1 : length(dataListing_PSG_GW_InfFact)
    GW_PSG = AngleSensor_PSG(dataListing_PSG_GW_InfFact(i).folder, dataListing_PSG_GW_InfFact(i).name, "Speed"); 
    peaks_temp = get_peaks(GW_PSG, "raw"); 
    peaks{idx} = peaks_temp; 
    Torque(idx) = GW_PSG.Trq; 
    Speed(idx) = GW_PSG.Rev;
    Temperature(idx) = GW_PSG.Temperature_mean;
    Quadrant(idx) = GW_PSG.Quadrant;
    filename = GW_PSG.filename;
    idx_file = find(filename == '_');
    Date{idx} = filename(1:idx_file(1)-1);
    Test{idx} = filename(idx_file(1)+1 : idx_file(2)-1);
    GW_PSG_mat{idx} = GW_PSG;
    idx = idx+1; 
end

varNames = {'Peaks', 'Torque', 'Speed', 'Temperature', 'Quadrant', 'Date', 'Test'};
varTypes = {'struct', 'double', 'double', 'double', 'double', 'string', 'string'};
T_data = table('Size', [length(Torque), numel(varNames)], 'VariableTypes', varTypes, 'VariableNames', varNames);

T_data.Peaks = peaks';
T_data.Torque = Torque'; T_data.Speed = Speed'; T_data.Temperature = Temperature';
T_data.Date = Date'; T_data.Test = Test'; T_data.Quadrant = Quadrant';
for i = 1 : height(T_data)
    peaks_temp = T_data.Peaks{i}; 
    sin_high_mean(i) = mean(peaks_temp.sin.high_mean); sin_low_mean(i) = mean(peaks_temp.sin.low_mean); 
    cos_high_mean(i) = mean(peaks_temp.cos.high_mean); cos_low_mean(i) = mean(peaks_temp.cos.low_mean);
end
T_data.("sine high") = sin_high_mean'; T_data.("sine low") = sin_low_mean';
T_data.("cosine high") = cos_high_mean'; T_data.("cosine low") = cos_low_mean';
datetime - t1

%
% Correlation analysis in four quadrants at artificial tests

variables2read = ["Torque", "Speed", "Temperature", "sine high", "cosine high"];
CorrType = "Pearson";
idx_Heating = cellfun(@(x) isequal(x, 'Heating'), T_data.Test); idx_TempEq  = cellfun(@(x) isequal(x, 'TempEquilibrium'), T_data.Test); idx_Cold = cellfun(@(x) isequal(x, 'Cold'), T_data.Test);
idx_Q(:,1) = T_data.Quadrant == 1; idx_Q(:,2) = T_data.Quadrant == 2; idx_Q(:,3) = T_data.Quadrant == 3; idx_Q(:,4) = T_data.Quadrant == 4; 

corr_table = T_data(idx_TempEq, variables2read);
[corr_matrix, ~] = corr(corr_table{:,:}, 'type', CorrType);
figure()
heatmap(corr_table.Properties.VariableNames, corr_table.Properties.VariableNames, corr_matrix)
title("All rho Temp equilibrium tests")
%get all indices

% Evaluation of different Quadrants without heating

num_Plots = 2;
figure()
for i = 1 : num_Plots*2
    idx = idx_Q(:,i)&idx_TempEq;
    T_temp = T_data(idx, variables2read); 
    corr_matrix = corr(T_temp{:,:}, 'type', CorrType);
    %plot correlation
    subplot(2, num_Plots, i)
    heatmap(T_temp.Properties.VariableNames, T_temp.Properties.VariableNames, corr_matrix)
    title(['Q' num2str(i) ' all'])
%     legend
end

% Random Test Evaluation in one Quadrant

corr_table = T_data(146:277, variables2read);
[corr_matrix, pval] = corr(corr_table{:,:}, 'type', CorrType);
figure()
heatmap(corr_table.Properties.VariableNames, corr_table.Properties.VariableNames, corr_matrix)
% savePlotWithCustomSettings('Ch 5\Heatmap', 0.75, 0.5, 0)
title("Rho in one Quadrant with Temperatre heating")

idx_Trq20 = corr_table.Torque == 20; idx_Trq85 = corr_table.Torque == 85; idx_Rev300 = corr_table.Speed == -300; idx_Rev1300 = corr_table.Speed == -1300;

figure()
scatter(corr_table{idx_Trq20&idx_Rev300, "Temperature"}, corr_table{idx_Trq20&idx_Rev300, "sine high"}, 'DisplayName', 'Torque = 20 Nm, Speed = 300 rpm'), hold on 
scatter(corr_table{idx_Trq85&idx_Rev300, "Temperature"}, corr_table{idx_Trq85&idx_Rev300, "sine high"}, 'DisplayName', 'Torque = 85 Nm, Speed = 300 rpm')
scatter(corr_table{idx_Trq20&idx_Rev1300, "Temperature"}, corr_table{idx_Trq20&idx_Rev1300, "sine high"}, 'DisplayName', 'Torque = 20 Nm, Speed = 1300 rpm')
scatter(corr_table{idx_Trq85&idx_Rev1300, "Temperature"}, corr_table{idx_Trq85&idx_Rev1300, "sine high"}, 'DisplayName', 'Torque = 85 Nm, Speed = 1300 rpm')
xlabel('Temperature in °C'), ylabel('Voltage in V'), title('Sine'), ylim([2.79 2.84])
savePlotWithCustomSettings('Ch5-3-3_HeatingCurveSine', 0.5, 0.5, 0)

figure()
scatter(corr_table{idx_Trq20&idx_Rev300, "Temperature"}, corr_table{idx_Trq20&idx_Rev300, "cosine high"}, 'DisplayName', 'Torque = 20 Nm, Speed = 300 rpm'), hold on 
scatter(corr_table{idx_Trq85&idx_Rev300, "Temperature"}, corr_table{idx_Trq85&idx_Rev300, "cosine high"}, 'DisplayName', 'Torque = 85 Nm, Speed = 300 rpm')
scatter(corr_table{idx_Trq20&idx_Rev1300, "Temperature"}, corr_table{idx_Trq20&idx_Rev1300, "cosine high"}, 'DisplayName', 'Torque = 20 Nm, Speed = 1300 rpm')
scatter(corr_table{idx_Trq85&idx_Rev1300, "Temperature"}, corr_table{idx_Trq85&idx_Rev1300, "cosine high"}, 'DisplayName', 'Torque = 85 Nm, Speed = 1300 rpm')
xlabel('Temperature in °C'), title('Cosine'), ylim([2.97 3.02])
savePlotWithCustomSettings('Ch5-3-3_HeatingCurveCosine', 0.5, 0.5, 0)
legend('NumColumns', 2, 'location', 'northoutside', 'Interpreter', 'latex');
savePlotWithCustomSettings('Ch5-3-3_HeatingCurveLegend', 1, 0.5, 0)

%

%Linear regression 
sineTable = T_data(146:277, ["Temperature" "Torque" "Speed" "sine high"]);
lm_sine = fitlm(sineTable); 
figure(), plot(lm_sine)
disp(lm_sine)
cosineTable = T_data(146:277, ["cosine high" "Temperature"]);
lm_cosine = fitlm(cosineTable.("Temperature"), cosineTable.("cosine high")); 
figure(), scatter(cosineTable.("Temperature"), cosineTable.("cosine high")), plot(lm_cosine)
disp(lm_cosine)

%
%% Chapter 5.4 - Tooth-specific signature
% works for most of the conditions at GSP. Sometimes there are deviations as the tooth-specific signature varies at different measurements. This has to be further investigated, if a method should be used, which works for all data.
% Further analyses showed that the experiments starting with Pinion 1 - 3 need a reference of these experiments. Might be due to the long time between the experiments.

idx = 1;
% ende = 1000; 
figure()
for i = 1 : 400 : length(dataListing_PSG_GW) 
    GW_PSG_vec{idx} = AngleSensor_PSG(dataListing_PSG_GW(i).folder, dataListing_PSG_GW(i).name, "Speed"); 
    GW_PSG_vec{idx}.delta_t = 1/40000; 
    GW_PSG_vec{idx} = shiftSignals(GW_PSG_vec{idx}, GW_PSG_vec{1});
    Peaks = get_peaks(GW_PSG_vec{idx}, "shifted");
    plot(Peaks.cos.high_mean, "DisplayName", ['i=' num2str(i) ' ' dataListing_PSG_GW(i).name]), hold on
    idx = idx +1;
end
legend()

%
% For the SEW results it works very well. Due to the constant speed throughout the tests. At sum measurements a shift of one tooth is present. Signal is way more stable than the one of the PSG test stand

idx = 1;
ende = 1000000; 
figure()
for i = 1 : length(dataListing_SEW_Z3) 
    GW_SEW_vec{idx} = AngleSensor_SEW(dataListing_SEW_Z3(i).folder, dataListing_SEW_Z3(i).name, "Z3_SEW"); 
    GW_SEW_vec{idx} = shiftSignals(GW_SEW_vec{idx}, GW_SEW_vec{1});
    Peaks = get_peaks(GW_SEW_vec{idx}, "shifted");
    plot(Peaks.sin.high_mean-mean(Peaks.sin.high_mean)), hold on    %angle
%     plot(GW_SEW_vec{idx}.angle_shifted(1 : ende), GW_SEW_vec{idx}.angle_deviation_shifted(1: ende), "DisplayName", ['i=' num2str(i) ' ' dataListing_SEW_Z3(i).name]), hold on 
    idx = idx +1;
end
%73/263 wrong identified, but only one tooth wrong to the right or left
%side

%
%% Chapter 5.5 - Optimisation

%Calibration - Ellipsenfit <apply for one revolution> usage of absolute position - Dependent on speed? Angular sampling? 
%IAA implementation 
%angular sampling

%% Influence of fit ellipses
% calculate mean of ellipses fit and save it
% GSP

for i = 1 : length(dataListing_PSG_GW)
    % GW = AngleSensor_PSG(dataListing_PSG_GW(1).folder, dataListing_PSG_GW(i).name, "Speed");
    % [~, Parameters_temp] = fit_ellipse(GW.cos, GW.sin);
    % Parameters_temp.Trq = GW.Trq;
    % Parameters_temp.Rev = GW.Rev; 
    % if i == 1
    %     Parameters = Parameters_temp; 
    % else                        
    %     Parameters(i) = Parameters_temp;
    % end
    GW = AngleSensor_PSG(dataListing_PSG_PR(1).folder, dataListing_PSG_PR(i).name, "AL795");
    [~, Parameters_temp] = fit_ellipse(GW.cos, GW.sin);
    Parameters_temp.Trq = GW.Trq;
    Parameters_temp.Rev = GW.Rev; 
    if i == 1
        Parameters = Parameters_temp; 
    else                        
        Parameters(i) = Parameters_temp;
    end
end
fit_ellipse_ParametersMean_GSP.phase = mean([Parameters.phase]); fit_ellipse_ParametersMean_GSP.X = mean([Parameters.X]); fit_ellipse_ParametersMean_GSP.Y = mean([Parameters.Y]); 
fit_ellipse_ParametersMean_GSP.p = mean([Parameters.p]); fit_ellipse_ParametersMean_GSP.q = mean([Parameters.q]); 
fit_ellipse_ParametersMean_GSP_PR = fit_ellipse_ParametersMean_GSP; 
save('fit_ellipse_ParametersMean_GSP_PR.mat', 'fit_ellipse_ParametersMean_GSP_PR')
% fnNames = fieldnames(Parameters); 
% for i = 1 : length(fnNames)-2
%     figure(), plot([Parameters.(fnNames{i})]), title(fnNames{i})
% end

%

GW = AngleSensor_PSG(dataListing_PSG_GW(1).folder, dataListing_PSG_GW(1).name, "Speed");
GW_processed = GW.processSinCos(["fitEllipses"]);

[Data, fit_ellipse_ParametersMean_GSP] = fit_ellipse(GW.cos, GW.sin);

figure(), hold on
plot(GW.t, GW.angle_deviation), plot(GW.t, GW_processed.angle_deviation), hold off

figure(), hold on,
GW.plotFFT('Hz', 'sin'), GW_processed.plotFFT('Hz', 'sin'), xlim([0 500])

%
% SEW

for i = 1 : length(dataListing_SEW_Z3)
    Z3 = AngleSensor_SEW(dataListing_SEW_Z3(1).folder, dataListing_SEW_Z3(i).name, "Z3_SEW");
    [~, Parameters_temp] = fit_ellipse(Z3.cos, Z3.sin);
    if i == 1
        ParametersZ3 = Parameters_temp; 
    else                        
        ParametersZ3(i) = Parameters_temp;
    end
    Z1 = AngleSensor_SEW(dataListing_SEW_Z1(1).folder, dataListing_SEW_Z1(i).name, "Z1_SEW");
    [~, Parameters_temp] = fit_ellipse(Z1.cos, Z1.sin);
    if i == 1
        ParametersZ1 = Parameters_temp; 
    else                        
        ParametersZ1(i) = Parameters_temp;
    end
end
fit_ellipse_ParametersMean_SEW_Z3.phase = mean([ParametersZ3.phase]); fit_ellipse_ParametersMean_SEW_Z3.X = mean([ParametersZ3.X]); fit_ellipse_ParametersMean_SEW_Z3.Y = mean([ParametersZ3.Y]); 
fit_ellipse_ParametersMean_SEW_Z3.p = mean([ParametersZ3.p]); fit_ellipse_ParametersMean_SEW_Z3.q = mean([ParametersZ3.q]); 
fit_ellipse_ParametersMean_SEW_Z1.phase = mean([ParametersZ1.phase]); fit_ellipse_ParametersMean_SEW_Z1.X = mean([ParametersZ1.X]); fit_ellipse_ParametersMean_SEW_Z1.Y = mean([ParametersZ1.Y]); 
fit_ellipse_ParametersMean_SEW_Z1.p = mean([ParametersZ1.p]); fit_ellipse_ParametersMean_SEW_Z1.q = mean([ParametersZ1.q]); 
save('fit_ellipse_ParametersMean_SEW_Z3.mat', 'fit_ellipse_ParametersMean_SEW_Z3')
save('fit_ellipse_ParametersMean_SEW_Z1.mat', 'fit_ellipse_ParametersMean_SEW_Z1')

%

fnNames = fieldnames(ParametersZ3); 
for i = 1 : length(fnNames)-2
    figure(), plot([ParametersZ3.(fnNames{i})]), title(fnNames{i})
end

%

Z3 = AngleSensor_SEW(dataListing_SEW_Z3(1).folder, dataListing_SEW_Z3(1).name, "Z3_SEW");
Z3_processed = Z3.processSinCos(["fitEllipses"]);

[Data, fit_ellipse_ParametersMean_SEW_Z3] = fit_ellipse(Z3.cos, Z3.sin);

figure(), hold on
plot(Z3.t, Z3.angle_deviation), plot(Z3_processed.t, Z3_processed.angle_deviation), hold off

% Schwingungen von 1 Hz in ersten Versuchen und 0.7 Hz bei späteren Versuchen überlagern, drehzahlunabhängig das Messsignal - was könnte das sein?

figure(), hold on,
Z3.plotFFT("Hz", "sin"), 
Z3_processed.plotFFT('Hz', 'sin'), xlim([0 500]), legend

%
%% Chapter 6 - Damage Detection
%% Generate and save features
% PSG

MeasuringPeriod = 1; delta_t = 1/40000; FilterMethod = 'Filter'; 
processMethods = ["raw", "preprocess", "fitEllipses"]; SplitSignal = 25; %SplitSignal = 3 at 1 second for all tests "_fitEllipse" overlap 75%; SplitSignal = 25 for Pinion 2, 3, 5 and 6 at 1 second "__fitEllipse_2356" overlap 90% 
Signaltype = ["sin", "cos", "IAS", "IAA"];
%include Trq 20, 55, 85 and Rev 300, 1300, 2000, because they are present
%in all tests
T = struct2table(dataListing_PSG_GW);
idxRev = abs(T{:, "Rev"}) == 300 | abs(T{:, "Rev"}) == 1300 | abs(T{:, "Rev"}) == 2000;
idxTrq = abs(T{:, "Trq"}) == 20 | abs(T{:, "Trq"}) == 55 | abs(T{:, "Trq"}) == 85; 
idxPinion = T{:, "Ritzel"} == 2 | T{:, "Ritzel"} == 3 | T{:, "Ritzel"} == 5 | T{:, "Ritzel"} == 6;
idxRed = idxRev & idxTrq & idxPinion; 
dataListing_PSG_GW = dataListing_PSG_GW(idxRed, :); dataListing_PSG_PR = dataListing_PSG_PR(idxRed, :); dataListing_PSG_Acc = dataListing_PSG_Acc(idxRed, :);  
dataList_temp = struct2table(dataListing_PSG_GW);  

%Gear wheel sensor
t1 = datetime; 
damageLabels = table(); idx = 1; clear statFeatures freqFeatures dataList
for i = 1 : length(dataListing_PSG_GW) 
    for j = 1 : SplitSignal
        Obj = AngleSensor_PSG(dataListing_PSG_GW(i).folder, dataListing_PSG_GW(i).name, "Speed");
        Obj.FilterMethod = FilterMethod; Obj.MeasuringPeriod = MeasuringPeriod; Obj.delta_t = delta_t;
        Obj.PartOfSignal = j;
        Obj = Obj.processSinCos(processMethods);
        statFeatures(idx, :) = calcStatFeatures(Obj, Signaltype); 
        freqFeatures(idx, :) = calcFreqFeatures(Obj, Signaltype); 
        damageLabels.damageLabels(idx) = labelDamages(Obj);
        dataList(idx,:) = dataList_temp(i,:);
        idx = idx +1; 
    end
end
Features_GW_PSG_2356 = [dataList damageLabels statFeatures freqFeatures]; 
save([basefolder '\Testdata_GSP\Features\stationary\Features_GW_PSG_2356.mat'], 'Features_GW_PSG_2356')
datetime - t1 %21 min, 2356 (1s, 10Parts) - 48 min, 2356 (1s, 25 parts 6:46 h)

% %Pole ring sensor
t1 = datetime; 
damageLabels = table(); idx = 1; clear statFeatures freqFeatures dataList 
for i = 1 : length(dataListing_PSG_PR)
     for j = 1 : SplitSignal
        sensorObj = AngleSensor_PSG(dataListing_PSG_PR(i).folder, dataListing_PSG_PR(i).name, "AL795");
        sensorObj.FilterMethod = FilterMethod; sensorObj.MeasuringPeriod = MeasuringPeriod; sensorObj.delta_t = delta_t;
        sensorObj.PartOfSignal = j;
        sensorObj = sensorObj.processSinCos(processMethods);
        statFeatures(idx, :) = calcStatFeatures(sensorObj, Signaltype); 
        freqFeatures(idx, :) = calcFreqFeatures(sensorObj, Signaltype); 
        damageLabels.damageLabels(idx) = labelDamages(sensorObj);
        dataList(idx,:) = dataList_temp(i,:);
        idx = idx +1; 
    end
end
Features_PR_PSG_2356= [dataList damageLabels statFeatures freqFeatures]; 
save([basefolder '\Testdata_GSP\Features\stationary\Features_PR_PSG_2356.mat'], 'Features_PR_PSG_2356')
datetime - t1 %17 min, 2356 (1s, 25 parts 7:51 h)

%

%acceleration sensor
t1 = datetime; 
idx = 1; damageLabels = table(); clear statFeatures freqFeatures dataList 
for i = 1 : length(dataListing_PSG_Acc)
     for j = 1 : SplitSignal
        Obj = AccelerationSensor_PSG(dataListing_PSG_Acc(i).folder, dataListing_PSG_Acc(i).name, "side");
        Obj.FilterMethod = FilterMethod; Obj.MeasuringPeriod = MeasuringPeriod; Obj.delta_t = delta_t;
        Obj.PartOfSignal = j;
        statFeatures(idx, :) = calcStatFeatures(Obj, "Signal"); 
        freqFeatures(idx, :) = calcFreqFeatures(Obj, "Signal"); 
        damageLabels.damageLabels(idx) = labelDamages(Obj);
        dataList(idx,:) = dataList_temp(i,:);
        idx = idx +1; 
    end
end

Features_AccSide_PSG_2356 = [dataList damageLabels statFeatures freqFeatures]; 
save([basefolder '\Testdata_GSP\Features\stationary\Features_AccSide_PSG_2356.mat'], 'Features_AccSide_PSG_2356')
datetime - t1 %16 min, 2356 20 min _2356, 2356 (1s, 25 parts 0:52 h)

%
% SEW

MeasuringPeriod = 50; FilterMethod = 'Filter'; 
processMethods = ["raw", "preprocess", "fitEllipses"]; 
Signaltype = ["sin", "cos", "IAS"];
t1 = datetime; 
for i = 1 : length(dataListing_SEW_Z3)
    Obj = AngleSensor_SEW(dataListing_SEW_Z3(i).folder, dataListing_SEW_Z3(i).name, "Z3_SEW");
    Obj.FilterMethod = FilterMethod; Obj.MeasuringPeriod = MeasuringPeriod;  
    Obj = Obj.processSinCos(processMethods);
    statFeatures_Z3(i, :) = calcStatFeatures(Obj, Signaltype); 
    freqFeatures_Z3(i, :) = calcFreqFeatures(Obj, Signaltype); 
    Obj = AngleSensor_SEW(dataListing_SEW_Z1(i).folder, dataListing_SEW_Z1(i).name, "Z1_SEW");
    Obj.FilterMethod = FilterMethod; Obj.MeasuringPeriod = MeasuringPeriod; 
    statFeatures_Z1(i, :) = calcStatFeatures(Obj, Signaltype); 
    freqFeatures_Z1(i, :) = calcFreqFeatures(Obj, Signaltype);

    Obj = AngleSensor_SEW(dataListing_SEW_MR3(i).folder, dataListing_SEW_MR3(i).name, "Pos3_SEW");
    Obj.FilterMethod = FilterMethod; Obj.MeasuringPeriod = MeasuringPeriod; 
    statFeatures_MR3(i, :) = calcStatFeatures(Obj, Signaltype); 
    freqFeatures_MR3(i, :) = calcFreqFeatures(Obj, Signaltype);

    Obj = AngleSensor_SEW(dataListing_SEW_MR4(i).folder, dataListing_SEW_MR4(i).name, "Pos4_SEW");
    Obj.FilterMethod = FilterMethod; Obj.MeasuringPeriod = MeasuringPeriod; 
    statFeatures_MR4(i, :) = calcStatFeatures(Obj, Signaltype); 
    freqFeatures_MR4(i, :) = calcFreqFeatures(Obj, Signaltype);
end

%

Features_SEW_Z3 = [struct2table(dataListing_SEW_Z3) statFeatures_Z3 freqFeatures_Z3]; 
save([basefolder '\Testdata_SEW\Features\Features_SEW_Z3_50s.mat'], 'Features_SEW_Z3') 
Features_SEW_Z1 = [struct2table(dataListing_SEW_Z1) statFeatures_Z1 freqFeatures_Z1]; 
save([basefolder '\Testdata_SEW\Features\Features_SEW_Z1_50s.mat'], 'Features_SEW_Z1')

Features_SEW_MR3 = [struct2table(dataListing_SEW_Z1) statFeatures_MR3 freqFeatures_MR3]; 
save([basefolder '\Testdata_SEW\Features\Features_SEW_MR3_50s.mat'], 'Features_SEW_MR3')
Features_SEW_MR4 = [struct2table(dataListing_SEW_Z1) statFeatures_MR4 freqFeatures_MR4]; 
save([basefolder '\Testdata_SEW\Features\Features_SEW_MR4_50s.mat'], 'Features_SEW_MR4')
%Acceleration sensor
% for i = 1 : length(dataListing_SEW_Acc_HSS)
%     Obj = AccelerationSensor_SEW(dataListing_SEW_Acc_HSS(i).folder, dataListing_SEW_Acc_HSS(i).name, "Acc_HSS");
%     Obj.FilterMethod = 'Filter'; Obj.MeasuringPeriod = 10;
%     statFeatures_HSS(i, :) = calcStatFeatures(Obj, "Signal"); 
%     freqFeatures_HSS(i, :) = calcFreqFeatures(Obj, "Signal"); 
%     Obj = AccelerationSensor_SEW(dataListing_SEW_Acc_LSS(i).folder, dataListing_SEW_Acc_LSS(i).name, "Acc_LSS");
%     Obj.FilterMethod = 'Filter'; Obj.MeasuringPeriod = 10;
%     statFeatures_LSS(i, :) = calcStatFeatures(Obj, "Signal"); 
%     freqFeatures_LSS(i, :) = calcFreqFeatures(Obj, "Signal"); 
% end
% Features_SEW_Acc_HSS = [struct2table(dataListing_SEW_Acc_HSS) statFeatures_HSS freqFeatures_HSS]; 
% save([basefolder '\Testdata_SEW\Features\Features_SEW_Acc_HSS.mat'], 'Features_SEW_Acc_HSS')
% Features_SEW_Acc_LSS = [struct2table(dataListing_SEW_Acc_HSS) statFeatures_LSS freqFeatures_LSS]; 
% save([basefolder '\Testdata_SEW\Features\Features_SEW_Acc_LSS.mat'], 'Features_SEW_Acc_LSS')
datetime - t1 % 14 min

%
%% Chapter 6.2 - at stationary conditions

%fres of Pspectrum has to be 1.5 Hz at 2000 rpm (divide fres by 111/5) to distinguish the peaks
MeasuringPeriod = 2; delta_t = 1/40000; FilterMethod = 'Filter'; 
processMethods = ["raw", "preprocess", "fitEllipses"];

load([basefolder '\Testdata_GSP\Features\stationary\Features_GW_PSG.mat']);
%visual comparison of sensitec tests between damaged and undamaged so GMF
%and sidebands
Q = 1; Trq = 85; Rev = 2000; Pinion = 0; 
idx_undamaged = [dataListing_PSG_GW.Damage] == 0 & [dataListing_PSG_GW.Quadrant] == Q & abs([dataListing_PSG_GW.Trq]) == Trq & abs([dataListing_PSG_GW.Rev]) == Rev & [dataListing_PSG_GW.Ritzel] == Pinion;
idx_damaged = [dataListing_PSG_GW.Damage] == 6 & [dataListing_PSG_GW.Quadrant] == Q & abs([dataListing_PSG_GW.Trq]) == Trq & abs([dataListing_PSG_GW.Rev]) == Rev & [dataListing_PSG_GW.Ritzel] == Pinion;
GW_PSG_undamaged = AngleSensor_PSG(dataListing_PSG_GW(idx_undamaged).folder, dataListing_PSG_GW(idx_undamaged).name, "Speed"); 
GW_PSG_damaged = AngleSensor_PSG(dataListing_PSG_GW(idx_damaged).folder, dataListing_PSG_GW(idx_damaged).name, "Speed");

GW_PSG_undamaged.FilterMethod = FilterMethod; GW_PSG_undamaged.MeasuringPeriod = MeasuringPeriod; GW_PSG_undamaged.delta_t = delta_t;
GW_PSG_damaged.FilterMethod = FilterMethod; GW_PSG_damaged.MeasuringPeriod = MeasuringPeriod; GW_PSG_damaged.delta_t = delta_t;
GW_PSG_undamaged = GW_PSG_undamaged.processSinCos(processMethods);
GW_PSG_damaged = GW_PSG_damaged.processSinCos(processMethods);

%Raw signals
figure(), hold on
plot(GW_PSG_undamaged.t, GW_PSG_undamaged.sin), plot(GW_PSG_damaged.t, GW_PSG_damaged.sin), title('sine'), hold off

%IAS signals
figure(), hold on
plot(GW_PSG_undamaged.t(1:end-1), GW_PSG_undamaged.IAS), plot(GW_PSG_damaged.t(1:end-1), GW_PSG_damaged.IAS), title('IAS'), hold off
%frequency domain 
limitsX = [16.7 25.3];
figure('Renderer', 'Painters'), subplot(2,1,1)
GW_PSG_undamaged.plotFFT('order', 'sin'), hold on
GW_PSG_damaged.plotFFT('order', 'sin'), xline(limitsX(1)+0.3 : 1 : limitsX(2)-0.3, "LineStyle", "--" ), xline(21), hold off, ylim([-60 0]),  xlim(limitsX); %legend(["undamaged", "damaged"], 'location', 'northoutside')
% subplot(3,1,2)
% GW_PSG_undamaged.plotFFT('order', 'cos'), hold on
% GW_PSG_damaged.plotFFT('order', 'cos'), xline(limitsX(1) : 1 : limitsX(2)), xlim(limitsX); title('Cos - FFT')
subplot(2,1,2)
GW_PSG_undamaged.plotFFT('order', 'IAS'), hold on
GW_PSG_damaged.plotFFT('order', 'IAS'), xline(limitsX(1)+0.3 : 1 : limitsX(2)-0.3, "LineStyle", "--" ), xline(21), hold off, xlim(limitsX); ylim([-25 30])
savePlotWithCustomSettings('Ch6-2-1_FFT_damagedVsundamaged', 1, 0.75, 0)

%

fig = figure(), 
GW_PSG_undamaged.plotFFT('order', 'IAS'), hold on
GW_PSG_damaged.plotFFT('order', 'IAS'), xline(limitsX(1)+0.3 : 1 : limitsX(2)-0.3, "LineStyle", "--" ), xline(21), hold off %legend(["undamaged", "damaged"], 'location', 'northoutside')
fig.Children.YLabel.String =[]; fig.Children.XLabel.String =[];
for i = 1 : 9
    xlim([15.95+i 16.05+i]), savePlotWithCustomSettings(['Ch6-2-1_FFT_damagedVsundamaged_IAS_' num2str(16+i)], 0.05, 0.15, 0)
end

%

%Amplitude quotient value of damaged by undamaged of the GMF and four surrounding sidebands
load('C:\Users\sn85zymy\Desktop\Diss\Testdata_GSP\Features\stationary\Features_GW_PSG_fitEllipse.mat')
load('C:\Users\sn85zymy\Desktop\Diss\Testdata_GSP\Features\stationary\Features_PR_PSG_fitEllipse.mat')
load('C:\Users\sn85zymy\Desktop\Diss\Testdata_GSP\Features\stationary\Features_AccSide_PSG.mat')
Rev = [300 1300 2000]; Trq = [20 55 85]; idx = 1; clear DivisionOfMean SubstractionOfMean
for i = 1 : length(Rev)
    for j = 1 : length(Trq)
        idx_undamaged = [Features_GW_PSG_fitEllipse.Damage] == 0 & [Features_GW_PSG_fitEllipse.Quadrant] == 1 & abs([Features_GW_PSG_fitEllipse.Trq]) == Trq(j) & abs([Features_GW_PSG_fitEllipse.Rev]) == Rev(i) & [Features_GW_PSG_fitEllipse.Ritzel] == 0;
        idx_damaged = [Features_GW_PSG_fitEllipse.Damage] == 6 & [Features_GW_PSG_fitEllipse.Quadrant] == 1 & abs([Features_GW_PSG_fitEllipse.Trq]) == Trq(j) & abs([Features_GW_PSG_fitEllipse.Rev]) == Rev(i) & [Features_GW_PSG_fitEllipse.Ritzel] == 0;
        
        features_GW = Features_GW_PSG_fitEllipse(idx_undamaged | idx_damaged,:);
        features_GW_sin = [features_GW{:,"sin_FFT_GMF12 Order1 Sb_f1"}{:, 4:7} features_GW{:, "sin_FFT_GMF12 Order1"} features_GW{:,"sin_FFT_GMF12 Order1 Sb_f1"}{:, 8:11}];
        features_GW_IAS = [features_GW{:,"IAS_FFT_GMF12 Order1 Sb_f1"}{:, 4:7} features_GW{:, "IAS_FFT_GMF12 Order1"} features_GW{:,"IAS_FFT_GMF12 Order1 Sb_f1"}{:, 8:11}];
        features_PR = Features_PR_PSG_fitEllipse(idx_undamaged | idx_damaged,:);
        features_PR_sin = [features_PR{:,"sin_FFT_GMF12 Order1 Sb_f1"}{:, 4:7} features_PR{:, "sin_FFT_GMF12 Order1"} features_PR{:,"sin_FFT_GMF12 Order1 Sb_f1"}{:, 8:11}];
        features_PR_IAS = [features_PR{:,"IAS_FFT_GMF12 Order1 Sb_f1"}{:, 4:7} features_PR{:, "IAS_FFT_GMF12 Order1"} features_PR{:,"IAS_FFT_GMF12 Order1 Sb_f1"}{:, 8:11}];
        features_Acc = Features_AccSide_PSG(idx_undamaged | idx_damaged,:);
        features_Acc = [features_Acc{:,"Signal_FFT_GMF12 Order1 Sb_f1"}{:, 4:7} features_Acc{:, "Signal_FFT_GMF12 Order1"} features_Acc{:,"Signal_FFT_GMF12 Order1 Sb_f1"}{:, 8:11}];
        
        GW.sin.undamaged.mean = mean(features_GW_sin(1:3,:), 1); 
        GW.sin.damaged.mean = mean(features_GW_sin(4:6,:), 1); 
        GW.sin.divisionOfMean = GW.sin.damaged.mean ./ GW.sin.undamaged.mean; GW.sin.substractionOfMean = GW.sin.damaged.mean - GW.sin.undamaged.mean;
        GW.IAS.undamaged.mean = mean(features_GW_IAS(1:3,:), 1); 
        GW.IAS.damaged.mean = mean(features_GW_IAS(4:6,:), 1); 
        GW.IAS.divisionOfMean = GW.IAS.damaged.mean ./ GW.IAS.undamaged.mean; GW.IAS.substractionOfMean = GW.IAS.damaged.mean - GW.IAS.undamaged.mean;
        PR.sin.undamaged.mean = mean(features_PR_sin(1:3,:), 1); 
        PR.sin.damaged.mean = mean(features_PR_sin(4:6,:), 1); 
        PR.sin.divisionOfMean = PR.sin.damaged.mean ./ PR.sin.undamaged.mean; PR.sin.substractionOfMean = PR.sin.damaged.mean - PR.sin.undamaged.mean;
        PR.IAS.undamaged.mean = mean(features_PR_IAS(1:3,:), 1);
        PR.IAS.damaged.mean = mean(features_PR_IAS(4:6,:), 1); 
        PR.IAS.divisionOfMean = PR.IAS.damaged.mean ./ PR.IAS.undamaged.mean; PR.IAS.substractionOfMean = PR.IAS.damaged.mean - PR.IAS.undamaged.mean;
        Acc.undamaged.mean = mean(features_Acc(1:3,:), 1); Acc.undamaged.std = std(features_Acc(1:3,:), 1);
        Acc.damaged.mean = mean(features_Acc(4:6,:), 1); Acc.damaged.std = std(features_Acc(4:6,:), 1); Acc.stdMean = mean([Acc.damaged.std; Acc.undamaged.std], 1);
        Acc.divisionOfMean = Acc.damaged.mean ./ Acc.undamaged.mean; Acc.substractionOfMean = Acc.damaged.mean - Acc.undamaged.mean;
        
        DivisionOfMean(:,:,idx) = [GW.sin.divisionOfMean; GW.IAS.divisionOfMean; PR.sin.divisionOfMean; PR.IAS.divisionOfMean; Acc.divisionOfMean];
        SubstractionOfMean(:,:,idx) = [GW.sin.substractionOfMean; GW.IAS.substractionOfMean; PR.sin.substractionOfMean; PR.IAS.substractionOfMean; Acc.substractionOfMean]; 
        idx = idx +1;
    end
end

% visual assessment of whole frequency domain is hard to differentiate --> sidebands carry damage information. Extract sidebands and plot them for different damage states at 2000 rpm and 85 Nm
%
% Analysis of features - carried out with FFT

load('C:\Users\sn85zymy\Desktop\Diss\Testdata_GSP\Features\stationary\Features_GW_PSG_fitEllipse.mat')
idxRev2000Trq85Q1 = [Features_GW_PSG_fitEllipse.Quadrant] == 1 & abs([Features_GW_PSG_fitEllipse.Trq]) == 85 & abs([Features_GW_PSG_fitEllipse.Rev]) == 2000 & [Features_GW_PSG_fitEllipse.Ritzel] == 0 & [Features_GW_PSG_fitEllipse.WDH] == 1;
idxRev2000Trq85Q2 = [Features_GW_PSG_fitEllipse.Quadrant] == 2 & abs([Features_GW_PSG_fitEllipse.Trq]) == 85 & abs([Features_GW_PSG_fitEllipse.Rev]) == 2000 & [Features_GW_PSG_fitEllipse.Ritzel] == 0 & [Features_GW_PSG_fitEllipse.WDH] == 1; 
FeaturesRev2000Trq85Q1 = Features_GW_PSG_fitEllipse(idxRev2000Trq85Q1,:);% FeaturesRev2000Trq85Q1 = FeaturesRev2000Trq85Q1(1:3:end,:);
FeaturesRev2000Trq85Q2 = Features_GW_PSG_fitEllipse(idxRev2000Trq85Q2,:); %FeaturesRev2000Trq85Q2 = FeaturesRev2000Trq85Q2(1:3:end,:);
lLabels = length(FeaturesRev2000Trq85Q1.damageLabels)/3; 
%statistical features 14 : 23
figure('Renderer', 'Painters'), 
idx2plot = 34:42;
for i = 1 : size(FeaturesRev2000Trq85Q1,1) / 3
    idx = i*3 - 2; 
    Stat2plot(i,:) = mean(FeaturesRev2000Trq85Q1{:, idx2plot}(idx : idx+2, :),1);
end
lgdnames = cellfun(@(x) x(5:end), FeaturesRev2000Trq85Q1.Properties.VariableNames(idx2plot), 'UniformOutput', false);
subplot(1,2,1), plot(Stat2plot ./ Stat2plot(1,:)), legend(lgdnames, "Location", "southwest")
ax = gca; ax.XTick = 1 : lLabels; ax.XTickLabel = FeaturesRev2000Trq85Q1.damageLabels(1:3:end); ax.XLim = [1, lLabels]; xlabel("Damage condition"); ylabel("Normalised feature value"), %title('Statistical features at 2000 rpm and 85 Nm in Quadrant 4')
for i = 1 : size(FeaturesRev2000Trq85Q2,1) / 3
    idx = i*3 - 2; 
    Stat2plot(i,:) = mean(FeaturesRev2000Trq85Q2{:, idx2plot}(idx : idx+2, :),1);
end
subplot(1,2,2), plot(Stat2plot ./ Stat2plot(1,:))
ax = gca; ax.XTick = 1 : lLabels; ax.XTickLabel = FeaturesRev2000Trq85Q1.damageLabels(1:3:end); ax.XLim = [1 lLabels]; xlabel("Damage condition") %title('Statistical features at 1300 rpm and 85 Nm in Quadrant 3')
savePlotWithCustomSettings('Ch6-2-1_statisticalFeatures', 1, 0.6, 0)

% %sidebands and GMF
SBs2plot = "sin_FFT_GMF12 Order1 Sb_f1"; GMF2plot = "sin_FFT_GMF12 Order1"; clear Feature2plot
legendnames = ["SB -4" "SB -3" "SB -2" "SB -1" "GMF" "SB +1" "SB +2" "SB +3" "SB +4"]; colornames = ["#77AC30" "#7E2F8E" "#EDB120" "#D95319" "#0072BD" "#D95319" "#EDB120" "#7E2F8E" "#77AC30"]; linestyles = ["--" "--" "--" "--" "-" "-." "-." "-." "-."];
SBsAroundGMF = [FeaturesRev2000Trq85Q1{:,SBs2plot}{:, 4:7} FeaturesRev2000Trq85Q1{:, GMF2plot} FeaturesRev2000Trq85Q1{:,SBs2plot}{:, 8:11}]; 
%for each damage label three features sets are available, through split of
%the signal to calculate the features -> meaning of this three feature sets
for i = 1 : size(SBsAroundGMF,1) / 3
    idx = i*3 - 2; 
    Feature2plot(i,:) = mean(SBsAroundGMF(idx : idx+2, :),1);
end
fig = figure('Renderer', 'Painters'), subplot(1,2,1), hold on
for i = 1 : size(SBsAroundGMF, 2)
    plot(Feature2plot(:,i) ./ Feature2plot(1,i), "Color", colornames(i), "lineStyle", linestyles(i)) 
end, hold off
legend(legendnames, "Location", "northwest"), ylim([0 18.5])
ax = gca; ax.XTick = 1 : lLabels; ax.XTickLabel = FeaturesRev2000Trq85Q1.damageLabels(1:3:end); ax.XLim = [1 lLabels]; xlabel("Damage condition"), ylabel("Relative increase in amplitude")%, title('Sidebands around GMF at 2000 rpm and 85 Nm in Quadrant 4')

SBsAroundGMF = [FeaturesRev2000Trq85Q2{:, SBs2plot}{:, 4:7} FeaturesRev2000Trq85Q2{:, GMF2plot} FeaturesRev2000Trq85Q2{:, SBs2plot}{:, 8:11}]; 
for i = 1 : size(SBsAroundGMF,1) / 3
    idx = i*3 - 2; 
    Feature2plot(i,:) = mean(SBsAroundGMF(idx : idx+2, :),1);
end
subplot(1,2,2), hold on
for i = 1 : size(SBsAroundGMF, 2)
    plot(Feature2plot(:,i) ./ Feature2plot(1,i), "Color", colornames(i), "lineStyle", linestyles(i)) 
end, hold off
% legend(legendnames, "Location", "northwest")
ax = gca; ax.XTick = 1 : lLabels; ax.XTickLabel = FeaturesRev2000Trq85Q2.damageLabels(1:3:end); ax.XLim = [1 lLabels]; xlabel("Damage condition")%, title('Sidebands around GMF at 2000 rpm and 85 Nm in Quadrant 2')
savePlotWithCustomSettings('Ch6-2-1_frequencyFeatures', 1, 0.6, 0)

SBs2plot = "sin_FFT_GMF12 Order1 Sb_f1"; GMF2plot = "sin_FFT_GMF12 Order1"; clear Feature2plot
legendnames = ["SB -4" "SB -3" "SB -2" "SB -1" "SB +1" "SB +2" "SB +3" "SB +4"]; colornames = ["#77AC30" "#7E2F8E" "#EDB120" "#D95319" "#D95319" "#EDB120" "#7E2F8E" "#77AC30"]; linestyles = ["--" "--" "--" "--" "-." "-." "-." "-."];
SBsAroundGMF = [FeaturesRev2000Trq85Q1{:,SBs2plot}{:, 4:11}]; 
%for each damage label three features sets are available, through split of
%the signal to calculate the features -> meaning of this three feature sets
for i = 1 : size(SBsAroundGMF,1) / 3
    idx = i*3 - 2; 
    Feature2plot(i,:) = mean(SBsAroundGMF(idx : idx+2, :),1);
end
fig = figure('Renderer', 'Painters'), subplot(1,2,1), hold on
for i = 1 : size(SBsAroundGMF, 2)
    plot(Feature2plot(:,i) - Feature2plot(1,i), "Color", colornames(i), "lineStyle", linestyles(i)) 
end, hold off
legend(legendnames, "Location", "northwest")
ax = gca; ax.XTick = 1 : lLabels; ax.XTickLabel = FeaturesRev2000Trq85Q1.damageLabels(1:3:end); ax.XLim = [1 lLabels]; 
xlabel("Damage condition"), ylabel("Sideband amplitudes in V $(x10^{-4})$", 'Interpreter', 'latex')%, title('Sidebands around GMF at 2000 rpm and 85 Nm in Quadrant 4')

SBsAroundGMF = [FeaturesRev2000Trq85Q2{:, SBs2plot}{:, 4:11}]; 
for i = 1 : size(SBsAroundGMF,1) / 3
    idx = i*3 - 2; 
    Feature2plot(i,:) = mean(SBsAroundGMF(idx : idx+2, :),1);
end
subplot(1,2,2), hold on
for i = 1 : size(SBsAroundGMF, 2)
    plot(Feature2plot(:,i) - Feature2plot(1,i), "Color", colornames(i), "lineStyle", linestyles(i)) 
end, hold off
% legend(legendnames, "Location", "northwest")
ax = gca; ax.XTick = 1 : lLabels; ax.XTickLabel = FeaturesRev2000Trq85Q2.damageLabels(1:3:end); ax.XLim = [1 lLabels]; xlabel("Damage condition")%, title('Sidebands around GMF at 2000 rpm and 85 Nm in Quadrant 2')
savePlotWithCustomSettings('Ch6-2-1_frequencyFeaturesSubstraction', 1, 0.6, 0)

%
% Acc

Acc_undamaged = AccelerationSensor_PSG(dataListing_PSG_Acc(idx_undamaged).folder, dataListing_PSG_Acc(idx_undamaged).name, "side");
Acc_undamaged.FilterMethod = FilterMethod; Acc_undamaged.MeasuringPeriod = MeasuringPeriod; Acc_undamaged.delta_t = delta_t;
Acc_damaged = AccelerationSensor_PSG(dataListing_PSG_Acc(idx_damaged).folder, dataListing_PSG_Acc(idx_damaged).name, "side");
Acc_damaged.FilterMethod = FilterMethod; Acc_damaged.MeasuringPeriod = MeasuringPeriod; Acc_damaged.delta_t = delta_t;

figure(), 
Acc_undamaged.plotFFT('order'), hold on
Acc_damaged.plotFFT('order'), xline(15 : 1 : 27), xlim([14 28]); title('Acc - FFT')

%
%% Chapter 6.2.3 - Classification

Features = load([basefolder '\Testdata_GSP\Features\stationary\Features_GW_PSG_fitEllipse.mat']);
Features = Features.(cell2mat(fieldnames(Features))); 
idxRed = (Features.Quadrant == 1 | Features.Quadrant == 4) & Features.damageLabels ~= 4.1 & Features.WDH ~= 3 & Features.Ritzel ~= 0; 
FeatureTable = Features(idxRed, 14:end); %end); %9-Quadrant, 11-Rev
allFeatures = table2array(bringAllTableDataInFirstLevel(FeatureTable));
allLabels = Features(idxRed, 1:13); 

labels = allLabels{:, ["Ritzel", "Damage"]};
[labels, labelNames, ~] = combineLabels(labels); 
accuracy = classifyRandomForest(allFeatures, labels, labelNames, 'KFold', []); title(['GW - all test data - 5-fold Accuracy: ', num2str(round(accuracy,2) * 100), '%']);

idxDmg0 = allLabels.damageLabels == 0; 
labels = allLabels{idxDmg0, ["Ritzel", "damageLabels"]};
[labels, labelNames, ~] = combineLabels(labels); 
accuracy = classifyRandomForest(allFeatures(idxDmg0,:), labels, labelNames, 'KFold', []); title(['All Pinions at Dmg 0 - 5-fold - Accuracy: ', num2str(round(accuracy,2) * 100), '%']);

Labels_temp = allLabels(idxDmg0, :);
trainIdx = Labels_temp.WDH == 1; testIdx = Labels_temp.WDH == 2; 
accuracy = classifyRandomForest(allFeatures(idxDmg0, :), labels, labelNames, trainIdx, testIdx); title(['All Pinions at Dmg 0 - Wdh1/Wdh2 - split Accuracy: ', num2str(round(accuracy,2) * 100), '%']);

labels = allLabels{:, "damageLabels"};
[labels, labelNames, ~] = combineLabels(labels); 
accuracy = classifyRandomForest(allFeatures, labels, labelNames, 'KFold', []); title(['all test data - 5-fold - damagedLabels Accuracy: ', num2str(round(accuracy,2) * 100), '%']);

%
% Pinions 2356

FeatureTableGW = load('C:\Users\sn85zymy\Desktop\Diss\Testdata_GSP\Features\stationary\Features_GW_PSG_fitEllipse_2356.mat');  
FeatureTableGW = FeatureTableGW.(cell2mat(fieldnames(FeatureTableGW))); 
cv = cvpartition(FeatureTableGW.damageLabels == 0, 'Holdout', 0.5, 'Stratify', true); 
idxRed = FeatureTableGW.WDH ~= 3 & (FeatureTableGW.Quadrant == 1 | FeatureTableGW.Quadrant == 4) & FeatureTableGW.damageLabels ~= 4.1; idxDmg345 = FeatureTableGW.damageLabels == 3.1 | FeatureTableGW.damageLabels == 3.2 | FeatureTableGW.damageLabels == 4.2 |FeatureTableGW.damageLabels == 5;
idxDmg012 = FeatureTableGW.damageLabels == 0 | FeatureTableGW.damageLabels == 1.4 | FeatureTableGW.damageLabels == 2;
idxRed = ((test(cv) & idxDmg012 | idxDmg345) & idxRed);

FeatureTableGW = FeatureTableGW(idxRed, :);
labels2356 = FeatureTableGW(:, 1:13); features2356 = table2array(bringAllTableDataInFirstLevel(FeatureTableGW(:, 14:end)));

idxDmg0 = labels2356.damageLabels == 0; 
labels_temp = labels2356{idxDmg0, "Ritzel"}; features = features2356(idxDmg0,:); 
[labels_temp, ~, ~] = combineLabels(labels_temp); labelNames = {'Pinion 2', 'Pinion 3', 'Pinion 5', 'Pinion 6'};
accuracy = classifyRandomForest(features, labels_temp, labelNames, 'KFold', []); title(['5-fold, Accuracy: ', num2str(round(accuracy,2) * 100), '%']);
savePlotWithCustomSettings('Ch6-2-3_GW_Dmg0_KFold', 0.4, 0.3, 0)

Labels_temp = labels2356(idxDmg0, :);
trainIdx = [Labels_temp.WDH == 1, Labels_temp.WDH == 2]; testIdx = ~trainIdx; 
accuracy = classifyRandomForest(features, labels_temp, labelNames, trainIdx, testIdx); title(['Repetition split, Accuracy: ', num2str(round(accuracy,2) * 100), '%']);
savePlotWithCustomSettings('Ch6-2-3_GW_Dmg0_RepSplit', 0.4, 0.3, 0)

createColorBarForConfChart(size(Labels_temp,1) / 4, [0 0.4470 0.7410])
savePlotWithCustomSettings('Ch6-2-3_GW_Dmg0_ColorBarCorrect', 0.4, 0.248, 0)
createColorBarForConfChart(size(Labels_temp,1) / 4, [0.8500 0.3250 0.0980])
savePlotWithCustomSettings('Ch6-2-3_GW_Dmg0_ColorBarFalse', 0.4, 0.248, 0)

%


[labels, labelNames, ~] = combineLabels(labels2356(:, "damageLabels"));
accuracy = classifyRandomForest(features2356, labels, labelNames, 'KFold', []); title(['5-fold, Accuracy: ', num2str(round(accuracy,2) * 100), '%']);
savePlotWithCustomSettings('Ch6-2-3_GW_Dmg_KFold', 0.6, 0.566, 0)

trainIdx = [labels2356.WDH == 1, labels2356.WDH == 2]; testIdx = ~trainIdx; 
accuracy = classifyRandomForest(features2356, labels, labelNames, trainIdx, testIdx); title(['Repetition split, Accuracy: ', num2str(round(accuracy,2) * 100), '%']);
savePlotWithCustomSettings('Ch6-2-3_GW_Dmg_RepSplit', 0.53, 0.5, 0)

[labels, labelNames, ~] = combineLabels(labels2356(:, "damageLabels"));
trainIdx = [labels2356.Ritzel == 2 | labels2356.Ritzel == 5, labels2356.Ritzel == 2 | labels2356.Ritzel == 6, labels2356.Ritzel == 3 | labels2356.Ritzel == 5, labels2356.Ritzel == 3 | labels2356.Ritzel == 6]; 
testIdx = ~trainIdx; 
accuracy = classifyRandomForest(features2356, labels, labelNames, trainIdx, testIdx); title(['Pinion split, Accuracy: ', num2str(round(accuracy,2) * 100), '%']);
savePlotWithCustomSettings('Ch6-2-3_GW_Dmg_PinionSplit', 0.53, 0.5, 0)

%

createColorBarForConfChart(size(labels,1) / 7, [0 0.4470 0.7410])
savePlotWithCustomSettings('Ch6-2-3_GW_Dmg_ColorBarCorrect', 0.4, 0.46, 0)
createColorBarForConfChart(size(labels,1) / 7, [0.8500 0.3250 0.0980])
savePlotWithCustomSettings('Ch6-2-3_GW_Dmg_ColorBarFalse', 0.4, 0.46, 0)

% idx56 = labels2356.Ritzel == 5 | labels2356.Ritzel == 6; 
% labels56 = labels2356(idx56, :); features56 = features2356(idx56, :); 
% [labels, labelNames, ~] = combineLabels(labels2356(idx56, "damageLabels"));
% trainIdx = [labels56.Ritzel == 5, labels56.Ritzel == 6]; 
% testIdx = [labels56.Ritzel == 6, labels56.Ritzel == 5];
% accuracy = classifyRandomForest(features56, labels, labelNames, trainIdx, testIdx); title(['Pinion split, Accuracy: ', num2str(round(accuracy,2) * 100), '%']);
% savePlotWithCustomSettings('Ch6-2-3_GW_Dmg_Pinion56Split', 0.53, 0.5, 0)

%
%% Chapter 6.2.4 - Evaluation: Comparison with other Sensors

% %Pole ring sensor
FeatureTablePR = load([basefolder '\Testdata_GSP\Features\stationary\Features_PR_PSG_fitEllipse_2356.mat']);
FeatureTablePR = FeatureTablePR.(cell2mat(fieldnames(FeatureTablePR))); 
FeatureTablePR = FeatureTablePR(idxRed, :); %9-Quadrant, 11-Rev
features2356 = table2array(bringAllTableDataInFirstLevel(FeatureTablePR(:, 14:end)));
labels2356 = FeatureTablePR(:, 1:13);

[labels, labelNames, ~] = combineLabels(labels2356(:, "damageLabels"));
accuracy = classifyRandomForest(features2356, labels, labelNames, 'KFold', []); title(['5-fold, Accuracy: ', num2str(round(accuracy,2) * 100), '%']);
savePlotWithCustomSettings('Additional\Ch6-2-3_PR_Dmg_5-fold', 0.53, 0.5, 0)

trainIdx = [labels2356.WDH == 1, labels2356.WDH == 2]; testIdx = ~trainIdx; 
accuracy = classifyRandomForest(features2356, labels, labelNames, trainIdx, testIdx); title(['Repetition split, Accuracy: ', num2str(round(accuracy,2) * 100), '%']);
savePlotWithCustomSettings('Ch6-2-3_PR_Dmg_RepSplit', 0.57, 0.57, 0) %0.57 for big, 0.47 for small

[labels, labelNames, ~] = combineLabels(labels2356(:, "damageLabels"));
trainIdx = [labels2356.Ritzel == 2 | labels2356.Ritzel == 5, labels2356.Ritzel == 2 | labels2356.Ritzel == 6, labels2356.Ritzel == 3 | labels2356.Ritzel == 5, labels2356.Ritzel == 3 | labels2356.Ritzel == 6]; 
testIdx = ~trainIdx; 
accuracy = classifyRandomForest(features2356, labels, labelNames, trainIdx, testIdx); title(['Pinion split, Accuracy: ', num2str(round(accuracy,2) * 100), '%']);
savePlotWithCustomSettings('Ch6-2-3_PR_Dmg_PinionSplit', 0.57, 0.57, 0)

createColorBarForConfChart(size(labels,1) / 7, [0 0.4470 0.7410])
savePlotWithCustomSettings('Ch6-2-3_RepSplit_ColorBarCorrect', 0.4, 0.25, 0)
createColorBarForConfChart(size(labels,1) / 7, [0.8500 0.3250 0.0980])
savePlotWithCustomSettings('Ch6-2-3_RepSplit_ColorBarFalse', 0.4, 0.25, 0)

%Acceleration sensor
FeatureTableAcc = load([basefolder '\Testdata_GSP\Features\stationary\Features_AccSide_PSG_2356.mat']);
FeatureTableAcc = FeatureTableAcc.(cell2mat(fieldnames(FeatureTableAcc)));
FeatureTableAcc = FeatureTableAcc(idxRed, :);
labels2356 = FeatureTableAcc(:, 1:13); features2356 = table2array(bringAllTableDataInFirstLevel(FeatureTableAcc(:, 14:end)));

[labels, labelNames, ~] = combineLabels(labels2356(:, "damageLabels"));
accuracy = classifyRandomForest(features2356, labels, labelNames, 'KFold', []); title(['5-fold, Accuracy: ', num2str(round(accuracy,2) * 100), '%']);
savePlotWithCustomSettings('Additional\Ch6-2-3_Acc_Dmg_5-fold', 0.53, 0.5, 0)

trainIdx = [labels2356.WDH == 1, labels2356.WDH == 2]; testIdx = ~trainIdx; 
accuracy = classifyRandomForest(features2356, labels, labelNames, trainIdx, testIdx); title(['Repetition split, Accuracy: ', num2str(round(accuracy,2) * 100), '%']);
savePlotWithCustomSettings('Ch6-2-3_Acc_Dmg_RepSplit', 0.57, 0.57, 0)

[labels, labelNames, ~] = combineLabels(labels2356(:, "damageLabels"));
trainIdx = [labels2356.Ritzel == 2 | labels2356.Ritzel == 5, labels2356.Ritzel == 2 | labels2356.Ritzel == 6, labels2356.Ritzel == 3 | labels2356.Ritzel == 5, labels2356.Ritzel == 3 | labels2356.Ritzel == 6]; 
testIdx = ~trainIdx; 
accuracy = classifyRandomForest(features2356, labels, labelNames, trainIdx, testIdx); title(['Pinion split, Accuracy: ', num2str(round(accuracy,2) * 100), '%']);
savePlotWithCustomSettings('Ch6-2-3_Acc_Dmg_PinionSplit', 0.57, 0.57, 0)

createColorBarForConfChart(3600, [0 0.4470 0.7410])
savePlotWithCustomSettings('Ch6-2-3_PinionSplit_ColorBarCorrect', 0.4, 0.25, 0)
createColorBarForConfChart(3600, [0.8500 0.3250 0.0980])
savePlotWithCustomSettings('Ch6-2-3_PinionSplit_ColorBarFalse', 0.4, 0.25, 0)

%
%% Chapter 6.5 - Optimisation

Features = load([basefolder '\Testdata_GSP\Features\stationary\Features_GW_PSG_fitEllipse_2356.mat']);
Features = Features.(cell2mat(fieldnames(Features))); 
%Definition of Pinions to investigate and norm parameters
idxRed = (Features.Quadrant == 1 | Features.Quadrant == 4) & (Features.Ritzel == 2 | Features.Ritzel == 3 | Features.Ritzel == 5 | Features.Ritzel == 6) & Features.damageLabels ~= 4.1 & Features.WDH ~= 3 & abs(Features.Rev) ~= 300; 
Revs = [300 1300 2000 -300 -1300 -2000]; Pinion = [2 3 5 6]; Dmg = 0; 
%delete 50% samples of Dmg 0, 1.4 and 2 to get equal amount of samples for
%all classes
idxDmg0 = Features.damageLabels == 0; idxDmg1 = Features.damageLabels == 1.4; idxDmg2 = Features.damageLabels == 2;
cv = cvpartition(idxDmg0, 'Holdout', 0.5, 'Stratify', true);
idxRedDmg = ~((idxDmg0 & test(cv)) | (idxDmg1 & test(cv)) | (idxDmg2 & test(cv))); 
idxRed = idxRed & idxRedDmg; 
%norm features by damage 0 to have a higher repeatability when training
%with one pinion and testing with an other
FeatureTable = Features(idxRed, 14:end); %9-Quadrant, 11-Rev
all_Features = table2array(bringAllTableDataInFirstLevel(FeatureTable));
allLabels = Features(idxRed, 1:13); 

% [all_Features, allLabels] = normFeatures(all_Features, allLabels, Revs, Pinion, Dmg);

[labels, labelNames, ~] = combineLabels(allLabels(:, "damageLabels"));
trainIdx = [allLabels.Ritzel == 2 | allLabels.Ritzel == 5, allLabels.Ritzel == 2 | allLabels.Ritzel == 6, allLabels.Ritzel == 3 | allLabels.Ritzel == 5, allLabels.Ritzel == 3 | allLabels.Ritzel == 6]; 
testIdx = ~trainIdx;
accuracy = classifyRandomForest(all_Features, labels, labelNames, trainIdx, testIdx); title(['GW - MT+BT Accuracy: ', num2str(round(accuracy,2) * 100), '% - normed - first/second']);

Features = load([basefolder '\Testdata_GSP\Features\stationary\Features_AccSide_PSG_2356.mat']);
Features = Features.(cell2mat(fieldnames(Features))); 
FeatureTable = Features(idxRed, 14:end); %9-Quadrant, 11-Rev
all_Features = table2array(bringAllTableDataInFirstLevel(FeatureTable));
allLabels = Features(idxRed, 1:13); 

% [all_Features, allLabels] = normFeatures(all_Features, allLabels, Revs, Pinion, Dmg);

[labels, labelNames, ~] = combineLabels(allLabels(:, "damageLabels"));
trainIdx = [allLabels.Ritzel == 2 | allLabels.Ritzel == 5, allLabels.Ritzel == 2 | allLabels.Ritzel == 6, allLabels.Ritzel == 3 | allLabels.Ritzel == 5, allLabels.Ritzel == 3 | allLabels.Ritzel == 6]; 
testIdx = ~trainIdx;
accuracy = classifyRandomForest(all_Features, labels, labelNames, trainIdx, testIdx); title(['Acc - MT+BT Accuracy: ', num2str(round(accuracy,2) * 100), '% - normed Features - first/second split']);

%

idxPinion2 = allLabels.Ritzel == 2; idxPinion3 = allLabels.Ritzel == 3; idxPinion5 = allLabels.Ritzel == 5; idxPinion6 = allLabels.Ritzel == 6; 
labels2356 = allLabels((idxPinion2 | idxPinion3 | idxPinion5 | idxPinion6), :); 
features2356 = all_Features((idxPinion2 | idxPinion3 | idxPinion5 | idxPinion6), :);
[labels, labelNames, ~] = combineLabels(labels2356(:, "damageLabels"));
trainIdx = labels2356.Ritzel == 2 | labels2356.Ritzel == 5; testIdx = labels2356.Ritzel == 3 | labels2356.Ritzel == 6; 
accuracy = classifyRandomForest(features2356, labels, labelNames, trainIdx, testIdx); title(['MT+BT - first/second split Accuracy: ', num2str(round(accuracy,2) * 100), '%']);

%
%% Chapter 6 - SEW fatigue damage test
% investigate raw signal, angle, IAS and influence of ellipses fit

MeasuringPeriod = 4; 
SEW_Z3 = AngleSensor_SEW(dataListing_SEW_Z3(1).folder, dataListing_SEW_Z3(1).name, "Z3_SEW");
SEW_Z1 = AngleSensor_SEW(dataListing_SEW_Z1(1).folder, dataListing_SEW_Z1(1).name, "Z1_SEW");
SEW_Z3.MeasuringPeriod = MeasuringPeriod; SEW_Z1.MeasuringPeriod = MeasuringPeriod;
SEW_Z3 = SEW_Z3.processSinCos("preprocess"); SEW_Z1 = SEW_Z1.processSinCos("preprocess"); 
[~, EllipsesfitParameters] = fit_ellipse(SEW_Z3.cos, SEW_Z3.sin);
SEW_Z3_Elips = SEW_Z3.fitEllipses(EllipsesfitParameters); 
%raw signal
figure(), hold on
plot(SEW_Z3.t, SEW_Z3.sin), plot(SEW_Z3.t, SEW_Z3.cos); 
plot(SEW_Z1.t, SEW_Z1.sin), plot(SEW_Z1.t, SEW_Z1.cos); 
%angle
figure(), hold on
plot(SEW_Z3.t, SEW_Z3.angle_deviation)
plot(SEW_Z1.t, SEW_Z1.angle_deviation) 
%IAS
figure(), hold on
plot(SEW_Z3.t(1:end-1), SEW_Z3.IAS)
plot(SEW_Z1.t(1:end-1), SEW_Z1.IAS) 
%ellipses fit
figure(), hold on 
scale = mean(abs(SEW_Z3_Elips.sin)) / mean(abs(SEW_Z3.sin));
plot(SEW_Z3.t, SEW_Z3.sin*scale), plot(SEW_Z3_Elips.t, SEW_Z3_Elips.sin); legend(["raw", "Ellipsesfit"]), title('Ellipsesfit')
figure(), hold on 
plot(SEW_Z3.t, SEW_Z3.angle_deviation), plot(SEW_Z3_Elips.t, SEW_Z3_Elips.angle_deviation); legend(["raw", "Ellipsesfit"])
figure(), hold on 
plot(SEW_Z3.t(1:end-1), SEW_Z3.IAS), plot(SEW_Z3_Elips.t(1:end-1), SEW_Z3_Elips.IAS); legend(["raw", "Ellipsesfit"])

% Ellipses fit has no to little influence on the sensor-specific deviations --> more specific approach. Revolution specific compensation
%
% compare undamaged vs damaged

FilterMethod = "filter"; processMethods = ["preprocess", "fitEllipses"];
MeasuringPeriod = 10;
SEW_Z3_undamaged = AngleSensor_SEW(dataListing_SEW_Z3(1).folder, dataListing_SEW_Z3(1).name, "Z3_SEW");
SEW_Z3_undamaged.FilterMethod = FilterMethod; SEW_Z3_undamaged.MeasuringPeriod = MeasuringPeriod; SEW_Z3_undamaged.PartOfSignal = 2;
SEW_Z3_undamaged = SEW_Z3_undamaged.processSinCos(processMethods);
SEW_Z3_damaged = AngleSensor_SEW(dataListing_SEW_Z3(1).folder, dataListing_SEW_Z3(260).name, "Z3_SEW");
SEW_Z3_damaged.FilterMethod = FilterMethod; SEW_Z3_damaged.MeasuringPeriod = MeasuringPeriod; SEW_Z3_damaged.PartOfSignal = 2; 
SEW_Z3_damaged = SEW_Z3_damaged.processSinCos(processMethods);
SEW_Acc_undamaged = AccelerationSensor_SEW(dataListing_SEW_Acc_HSS(1).folder, dataListing_SEW_Acc_HSS(1).name, "Acc_HSS")
SEW_Acc_undamaged.MeasuringPeriod = MeasuringPeriod;
SEW_Acc_damaged = AccelerationSensor_SEW(dataListing_SEW_Acc_HSS(1).folder, dataListing_SEW_Acc_HSS(85).name, "Acc_HSS")
SEW_Acc_damaged.MeasuringPeriod = MeasuringPeriod;
%frequency domain 
limitsX = [10.7 19.3];
figure(), subplot(2,1,1)
SEW_Z3_undamaged.plotFFT('order', 'sin'), hold on
SEW_Z3_damaged.plotFFT('order', 'sin'), xline(limitsX(1)+0.3 : 1 : limitsX(2), "LineStyle", "--" ), xline(15), xlim(limitsX); %ylim([-70 0])%legend(["undamaged", "damaged"], 'location', 'northoutside') 
subplot(2,1,2)
SEW_Z3_undamaged.plotFFT('order', 'IAS'), hold on
SEW_Z3_damaged.plotFFT('order', 'IAS'), xline(limitsX(1)+0.3 : 1 : limitsX(2), "LineStyle", "--" ), xline(15), xlim(limitsX); ylim([-45 20])
savePlotWithCustomSettings('Ch6-3_FFT_damagedVsundamaged', 1, 0.75, 0)

%

fig = figure(), 
SEW_Z3_undamaged.plotFFT('order', 'sin'), hold on
SEW_Z3_damaged.plotFFT('order', 'sin'), xline(limitsX(1)+0.3 : 1 : limitsX(2)-0.3, "LineStyle", "--" ), xline(15), hold off 
fig.Children.YLabel.String =[]; fig.Children.XLabel.String =[];
for i = 1 : 9
    xlim([9.95+i 10.05+i]), savePlotWithCustomSettings(['Ch6-3_FFT_damagedVsundamaged_sin_' num2str(10+i)], 0.03, 0.15, 0)
end

%

%Bearing frequencies
figure(), hold on
SEW_Z3_undamaged.plotFFT('Hz', 'IAS'), SEW_Z3_damaged.plotFFT('Hz', 'IAS')
SEW_Acc_undamaged.plotFFT('Hz'), SEW_Acc_damaged.plotFFT('Hz')
xline([SEW_Z3_undamaged.quant.Bearing.A1{:,:}, SEW_Z3_undamaged.quant.Bearing.A2{:,:}, SEW_Z3_undamaged.quant.Bearing.B{:,:}]), xlim([0 max([SEW_Z3_undamaged.quant.Bearing.A1{:,:}, SEW_Z3_undamaged.quant.Bearing.A2{:,:}, SEW_Z3_undamaged.quant.Bearing.B{:,:}])])

%
% compare features

load([basefolder '\Testdata_SEW\Features\Features_SEW_Z3_50s.mat']); load([basefolder '\Testdata_SEW\Features\Features_SEW_Z1_50s.mat']);
load([basefolder '\Testdata_SEW\Features\Features_SEW_MR3_50s.mat']); load([basefolder '\Testdata_SEW\Features\Features_SEW_MR4_50s.mat']);
load([basefolder '\Testdata_SEW\Features\Features_SEW_Acc_HSS_10s.mat']); load([basefolder '\Testdata_SEW\Features\Features_SEW_Acc_LSS_10s.mat']);
testdaysMR = linspace(0, 88, 263); %85 test days at 263 data points acquired every 8 hours

%

%fatigue pitting test - show rise of features throughout the tests
idxEMV = [1:199 207:263];
idxStat = find(strcmp(Features_SEW_Z3.Properties.VariableNames, 'IAS_RMS')) : find(strcmp(Features_SEW_Z3.Properties.VariableNames, 'IAS_NA4')); 
figure(), 
% subplot(3,1,1), plot(Features_SEW_Z3{:, idxStat} - Features_SEW_Z3{1, idxStat}), ylim([0.4 1.3]), title('Z3 statistical features')
% subplot(3,1,2), plot(Features_SEW_Z1{:, idxStat} - Features_SEW_Z1{1, idxStat}), ylim([0 2.5]), title('Z1 statistical features')
% subplot(3,1,3), plot(Features_SEW_Acc_HSS{:, 13:22} - Features_SEW_Acc_HSS{1, 13:22}), title('Acc HSS statistical features')
ColMin = min(Features_SEW_Z3{idxEMV, idxStat}); ColMax = max(Features_SEW_Z3{idxEMV, idxStat});
subplot(3,1,1), plot(rescale(Features_SEW_Z3{idxEMV, idxStat}, "InputMin", ColMin, "InputMax", ColMax)), title('Z3 statistical features')
ColMin = min(Features_SEW_Z1{idxEMV, idxStat}); ColMax = max(Features_SEW_Z1{idxEMV, idxStat});
subplot(3,1,2), plot(rescale(Features_SEW_Z1{idxEMV, idxStat}, "InputMin", ColMin, "InputMax", ColMax)), title('Z1 statistical features')
ColMin = min(Features_SEW_Acc_HSS{:, 13:21}); ColMax = max(Features_SEW_Acc_HSS{:, 13:21});
subplot(3,1,3), plot(rescale(Features_SEW_Acc_HSS{:, 13:21}, "InputMin", ColMin, "InputMax", ColMax)), title('Acc HSS statistical features')
legend(Features_SEW_Z3.Properties.VariableNames{idxStat})

%

testdaysMR = linspace(0, 88, 263); testdaysAcc = linspace(0, 88, 85);
legendnames = ["SB -4" "SB -3" "SB -2" "SB -1" "SB +1" "SB +2" "SB +3" "SB +4"]; colornames = ["#77AC30" "#7E2F8E" "#EDB120" "#D95319" "#D95319" "#EDB120" "#7E2F8E" "#77AC30"]; linestyles = ["--" "--" "--" "--" "-." "-." "-." "-."];
%GMF of Z3, Z1, Acc
fig = figure('Renderer', 'Painters'), hold on, 
plot(testdaysMR, Features_SEW_Z3{:, "IAS_FFT_GMF34 Order1"} ./ Features_SEW_Z3{1, "IAS_FFT_GMF34 Order1"}, 'DisplayName', 'GWS Z3')
plot(testdaysMR, Features_SEW_Z1{:, "IAS_FFT_GMF34 Order1"} ./ Features_SEW_Z1{1, "IAS_FFT_GMF34 Order1"}, 'DisplayName', 'GWS Z1')
plot(testdaysMR, Features_SEW_MR3{:, "IAS_FFT_GMF34 Order1"} ./ Features_SEW_MR3{1, "IAS_FFT_GMF34 Order1"}, 'DisplayName', 'GE 3')
plot(testdaysMR, Features_SEW_MR4{:, "IAS_FFT_GMF34 Order1"} ./ Features_SEW_MR4{1, "IAS_FFT_GMF34 Order1"}, 'DisplayName', 'GE 4')
plot(testdaysAcc, Features_SEW_Acc_HSS{:, "Signal_FFT_GMF34 Order1"} ./ Features_SEW_Acc_HSS{1, "Signal_FFT_GMF34 Order1"}, 'DisplayName', 'Acc HSS'), legend('location', 'northwest'), xlim([0 88]), hold off
xlabel('Test days'), ylabel('Relative increase in amplitude')
savePlotWithCustomSettings('Ch6-3_GMF_Z3_Z1_Acc', 1, 0.4, 0)
legend off, fig.Children.XLabel.String = []; fig.Children.YLabel.String = [];
xlim([0 55]); ylim([0.8 1.8]);
savePlotWithCustomSettings('Ch6-3_GMF_Z3_Z1_Acc_zoom', 0.35, 0.2, 0)

%


figure('Renderer', 'Painters'), hold on, 
plot(testdaysMR, Features_SEW_Z3{:, "IAS_FFT_GMF12 Order1"} ./ Features_SEW_Z3{1, "IAS_FFT_GMF12 Order1"}, 'DisplayName', 'GWS Z3')
plot(testdaysMR, Features_SEW_Z1{:, "IAS_FFT_GMF12 Order1"} ./ Features_SEW_Z1{1, "IAS_FFT_GMF12 Order1"}, 'DisplayName', 'GWS Z1')
plot(testdaysMR, Features_SEW_MR3{:, "IAS_FFT_GMF12 Order1"} ./ Features_SEW_MR3{1, "IAS_FFT_GMF12 Order1"}, 'DisplayName', 'GE 3')
plot(testdaysMR, Features_SEW_MR4{:, "IAS_FFT_GMF12 Order1"} ./ Features_SEW_MR4{1, "IAS_FFT_GMF12 Order1"}, 'DisplayName', 'GE 4')
plot(testdaysAcc, Features_SEW_Acc_HSS{:, "Signal_FFT_GMF12 Order1"} ./ Features_SEW_Acc_HSS{1, "Signal_FFT_GMF12 Order1"}, 'DisplayName', 'Acc HSS'),
plot(testdaysAcc, Features_SEW_Acc_LSS{:, "Signal_FFT_GMF12 Order1"} ./ Features_SEW_Acc_LSS{1, "Signal_FFT_GMF12 Order1"}, 'DisplayName', 'Acc HSS')
legend('location', 'northwest'), xlim([0 88]), hold off
xlabel('Test days'), ylabel('Relative increase in amplitude')

%Sidebands of Z3
fig = figure('Renderer', 'Painters'), hold on,
Features2plot = Features_SEW_Z3{:, "IAS_FFT_GMF34 Order1 Sb_fC"}{:, 4:11};
for i = 1 : size(Features2plot, 2)
    plot(testdaysMR, Features2plot(:,i) - Features2plot(1, i), 'DisplayName', legendnames(i), 'Color', colornames(i), 'LineStyle', linestyles(i))
end
xlabel('Test days'), ylabel('Amplitude in rpm'), legend('location', 'northwest'), xlim([0 88])
savePlotWithCustomSettings('Ch6-3_Sidebands_Z3_absolute', 0.53, 0.45, 0)
% legend off, fig.Children.XLabel.String = []; fig.Children.YLabel.String = []; ylim([0 5])
% savePlotWithCustomSettings('Ch6-3_Sidebands_Z3_absolute_zoom', 0.3, 0.15, 0)

%Sidebands of MR3
figure('Renderer', 'Painters'), hold on,
Features2plot = Features_SEW_MR3{:, "IAS_FFT_GMF34 Order1 Sb_fC"}{:, 4:11};
for i = 1 : size(Features2plot, 2)
    plot(testdaysMR, Features2plot(:,i) ./ Features2plot(1, i), 'DisplayName', legendnames(i), 'Color', colornames(i), 'LineStyle', linestyles(i))
end
xlabel('Testd ays'), ylabel('Relative increase in amplitude'), legend('location', 'northwest'), xlim([0 88])
%Sidebands of Acc
figure('Renderer', 'Painters'), hold on,
Features2plot = Features_SEW_Acc_HSS{:, "Signal_FFT_GMF34 Order1 Sb_fC"}{:, 4:11};
for i = 1 : size(Features2plot, 2)
    plot(testdaysAcc, Features2plot(:,i) ./ Features2plot(1, i), 'DisplayName', legendnames(i), 'Color', colornames(i), 'LineStyle', linestyles(i))
end
xlabel('Test days'), ylabel('Relative increase in amplitude'), legend('location', 'northwest'), xlim([0 88])
%Sum of Sidebands of Z3, Z1, Acc
figure('Renderer', 'Painters'), hold on
plot(testdaysMR, Features_SEW_Z3{:, "IAS_FFT_GMF34 Order1 Sb_fC Sum"} ./ Features_SEW_Z3{1, "IAS_FFT_GMF34 Order1 Sb_fC Sum"})
plot(testdaysMR, Features_SEW_Z1{:, "IAS_FFT_GMF34 Order1 Sb_fC Sum"} ./ Features_SEW_Z1{1, "IAS_FFT_GMF34 Order1 Sb_fC Sum"})
plot(testdaysMR, Features_SEW_MR3{:, "IAS_FFT_GMF34 Order1 Sb_fC Sum"} ./ Features_SEW_MR3{1, "IAS_FFT_GMF34 Order1 Sb_fC Sum"})
plot(testdaysMR, Features_SEW_MR4{:, "IAS_FFT_GMF34 Order1 Sb_fC Sum"} ./ Features_SEW_MR4{1, "IAS_FFT_GMF34 Order1 Sb_fC Sum"})
plot(testdaysAcc, Features_SEW_Acc_HSS{:, "Signal_FFT_GMF34 Order1 Sb_fC Sum"} ./ Features_SEW_Acc_HSS{1, "Signal_FFT_GMF34 Order1 Sb_fC Sum"})
xlabel('Testdays'), ylabel('Amplitude in rpm and m/s^2'), legend(["GW Z3", "GW Z1", "GE 3", "GE 4", "Acc HSS", "Acc LSS"], 'location', 'northwest'), %xlim([0 88])
title('Sum of Sidebands GMF34 sb fC - Substraction')
%Sum of Sidebands
fig = figure('Renderer', 'Painters'), hold on
plot(testdaysMR, sum(Features_SEW_Z3{:, "IAS_FFT_GMF34 Order1 Sb_fC"}{:,:} ./ Features_SEW_Z3{:, "IAS_FFT_GMF34 Order1 Sb_fC"}{1,:}, 2)/14)
plot(testdaysMR, sum(Features_SEW_Z1{:, "IAS_FFT_GMF34 Order1 Sb_fC"}{:,:} ./ Features_SEW_Z1{:, "IAS_FFT_GMF34 Order1 Sb_fC"}{1,:}, 2)/14)
plot(testdaysMR, sum(Features_SEW_MR3{:, "IAS_FFT_GMF34 Order1 Sb_fC"}{:,:} ./ Features_SEW_MR3{:, "IAS_FFT_GMF34 Order1 Sb_fC"}{1,:}, 2)/14)
plot(testdaysMR, sum(Features_SEW_MR4{:, "IAS_FFT_GMF34 Order1 Sb_fC"}{:,:} ./ Features_SEW_MR4{:, "IAS_FFT_GMF34 Order1 Sb_fC"}{1,:}, 2)/14)
plot(testdaysAcc, sum(Features_SEW_Acc_HSS{:, "Signal_FFT_GMF34 Order1 Sb_fC"}{:,:} ./ Features_SEW_Acc_HSS{1, "Signal_FFT_GMF34 Order1 Sb_fC"}{:,:}, 2)/14)
xlabel('Test days'), ylabel('Relative increase in amplitude'),  xlim([0 88]), legend(["GW Z3", "GW Z1", "GE 3", "GE 4", "Acc HSS"], 'location', 'northwest')
savePlotWithCustomSettings('Ch6-3_SumOfSidebands_Z3_Z1_Acc_relative', 0.53, 0.45, 0)
legend off, fig.Children.XLabel.String = []; fig.Children.YLabel.String = []; fig.Children.XTickLabel = {'10', '15', '20', '25'};
xlim([8 29]); ylim([0.8 2.5]); 
savePlotWithCustomSettings('Ch6-3_SumOfSidebands_Z3_Z1_Acc_relative_zoom', 0.27, 0.13, 0)
fig = figure('Renderer', 'Painters'), hold on
plot(testdaysMR, sum(Features_SEW_Z3{:, "IAS_FFT_GMF34 Order1 Sb_fC"}{:,:} - Features_SEW_Z3{:, "IAS_FFT_GMF34 Order1 Sb_fC"}{1,:}, 2)/14)
plot(testdaysMR, sum(Features_SEW_Z1{:, "IAS_FFT_GMF34 Order1 Sb_fC"}{:,:} - Features_SEW_Z1{:, "IAS_FFT_GMF34 Order1 Sb_fC"}{1,:}, 2)/14)
plot(testdaysMR, sum(Features_SEW_MR3{:, "IAS_FFT_GMF34 Order1 Sb_fC"}{:,:} - Features_SEW_MR3{:, "IAS_FFT_GMF34 Order1 Sb_fC"}{1,:}, 2)/14)
plot(testdaysMR, sum(Features_SEW_MR4{:, "IAS_FFT_GMF34 Order1 Sb_fC"}{:,:} - Features_SEW_MR4{:, "IAS_FFT_GMF34 Order1 Sb_fC"}{1,:}, 2)/14)
ylabel('Amplitude in rpm'),ylim([-0.0016 0.016]), yyaxis right
plot(testdaysAcc, sum(Features_SEW_Acc_HSS{:, "Signal_FFT_GMF34 Order1 Sb_fC"}{:,:} - Features_SEW_Acc_HSS{1, "Signal_FFT_GMF34 Order1 Sb_fC"}{:,:}, 2)/14)
xlabel('Test days'), ylabel('Amplitude in {m/s^2}'),  xlim([0 88]), ylim([-0.0006 0.006]), legend(["GW Z3", "GW Z1", "GE 3", "GE 4", "Acc HSS"], 'location', 'northwest')
savePlotWithCustomSettings('Ch6-3_SumOfSidebands_Z3_Z1_Acc_absolute', 0.53, 0.45, 0)
legend off, fig.Children.XLabel.String = []; fig.Children.YLabel.String = []; fig.Children.XTickLabel = {'10', '15', '20', '25'};
xlim([8 29]), yyaxis left, ylim([-0.00008 0.0004]), fig.Children.YLabel.String = []; yyaxis right, ylim([-0.00008 0.0004])
savePlotWithCustomSettings('Ch6-3_SumOfSidebands_Z3_Z1_Acc_absolute_zoom', 0.29, 0.13, 0)

%

figure('Renderer', 'Painters'), hold on
plot(testdaysMR, sum(Features_SEW_Z3{:, "IAS_FFT_GMF34 Order1 Sb_fB"}{:,:} ./ Features_SEW_Z3{:, "IAS_FFT_GMF34 Order1 Sb_fB"}{1,:}, 2)/14)
plot(testdaysMR, sum(Features_SEW_Z1{:, "IAS_FFT_GMF34 Order1 Sb_fB"}{:,:} ./ Features_SEW_Z1{:, "IAS_FFT_GMF34 Order1 Sb_fB"}{1,:}, 2)/14)
plot(testdaysMR, sum(Features_SEW_MR3{:, "IAS_FFT_GMF34 Order1 Sb_fB"}{:,:} ./ Features_SEW_MR3{:, "IAS_FFT_GMF34 Order1 Sb_fB"}{1,:}, 2)/14)
plot(testdaysMR, sum(Features_SEW_MR4{:, "IAS_FFT_GMF34 Order1 Sb_fB"}{:,:} ./ Features_SEW_MR4{:, "IAS_FFT_GMF34 Order1 Sb_fB"}{1,:}, 2)/14)
plot(testdaysAcc, sum(Features_SEW_Acc_HSS{:, "Signal_FFT_GMF34 Order1 Sb_fB"}{:,:} ./ Features_SEW_Acc_HSS{1, "Signal_FFT_GMF34 Order1 Sb_fB"}{:,:}, 2)/14)
% plot(testdaysAcc, sum(Features_SEW_Acc_LSS{:, "Signal_FFT_GMF34 Order1 Sb_fB"}{:,:} ./ Features_SEW_Acc_LSS{1, "Signal_FFT_GMF34 Order1 Sb_fB"}{:,:}, 2)/14)
xlabel('Testdays'), ylabel('Relative increase in amplitude'), legend(["GW Z3", "GW Z1", "GE 3", "GE 4", "Acc HSS", "Acc LSS"], 'location', 'northwest'), xlim([0 88])
title('Sum of Sidebands GMF34 sb fB - Division')

figure('Renderer', 'Painters'), hold on
plot(testdaysMR, sum(Features_SEW_Z3{:, "IAS_FFT_GMF12 Order1 Sb_fB"}{:,:} ./ Features_SEW_Z3{:, "IAS_FFT_GMF12 Order1 Sb_fB"}{1,:}, 2)/14)
plot(testdaysMR, sum(Features_SEW_Z1{:, "IAS_FFT_GMF12 Order1 Sb_fB"}{:,:} ./ Features_SEW_Z1{:, "IAS_FFT_GMF12 Order1 Sb_fB"}{1,:}, 2)/14)
plot(testdaysMR, sum(Features_SEW_MR3{:, "IAS_FFT_GMF12 Order1 Sb_fB"}{:,:} ./ Features_SEW_MR3{:, "IAS_FFT_GMF12 Order1 Sb_fB"}{1,:}, 2)/14)
plot(testdaysMR, sum(Features_SEW_MR4{:, "IAS_FFT_GMF12 Order1 Sb_fB"}{:,:} ./ Features_SEW_MR4{:, "IAS_FFT_GMF12 Order1 Sb_fB"}{1,:}, 2)/14)
plot(testdaysAcc, sum(Features_SEW_Acc_HSS{:, "Signal_FFT_GMF12 Order1 Sb_fB"}{:,:} ./ Features_SEW_Acc_HSS{1, "Signal_FFT_GMF12 Order1 Sb_fB"}{:,:}, 2)/14)
% plot(testdaysAcc, sum(Features_SEW_Acc_LSS{:, "Signal_FFT_GMF12 Order1 Sb_fB"}{:,:} ./ Features_SEW_Acc_LSS{1, "Signal_FFT_GMF12 Order1 Sb_fB"}{:,:}, 2)/14)
xlabel('Testdays'), ylabel('Relative increase in amplitude'), legend(["GW Z3", "GW Z1", "GE 3", "GE 4", "Acc HSS", "Acc LSS"], 'location', 'northwest'), xlim([0 88])
title('Sum of Sidebands GMF12 sb fB - Division')

%
%% Chapter 7 - Potentials for Torque, TE, Absolute ... Measurement
%% Absolute position using tooth-specific signature

idx = 1;
ende = 1000000; 
PreProcessMethod = ["preprocess"];
figure()
t1 = datetime; 
for i = 1 : 1 : length(dataListing_SEW_Z3) 
    GW_SEW_vec{idx} = AngleSensor_SEW(dataListing_SEW_Z3(i).folder, dataListing_SEW_Z3(i).name, "Z3_SEW"); 
    GW_SEW_vec{idx} = processSinCos(GW_SEW_vec{idx}, PreProcessMethod);
    GW_SEW_vec{idx} = shiftSignals(GW_SEW_vec{idx}, GW_SEW_vec{1});
    Peaks{idx} = get_peaks(GW_SEW_vec{idx}, "shifted");
    plot(Peaks{idx}.sin.high_mean-mean(Peaks{idx}.sin.high_mean), 'DisplayName', GW_SEW_vec{idx}.filename), hold on 
    idx = idx +1;
end
phase_shift_Z3 = cellfun(@(x) x.phase_shift, GW_SEW_vec);
Peaks_Z3_sin_high_mean = cellfun(@(x) x.sin.high_mean, Peaks, 'UniformOutput', false);
Peaks_Z3_sin_high_mean = [Peaks_Z3_sin_high_mean{:,:}];

idx = 1;
figure()
for i = 1 : 1 :  length(dataListing_SEW_Z1) 
    GW_SEW_vec{idx} = AngleSensor_SEW(dataListing_SEW_Z1(i).folder, dataListing_SEW_Z1(i).name, "Z1_SEW"); 
    GW_SEW_vec{idx} = shiftSignals(GW_SEW_vec{idx}, GW_SEW_vec{1});
    Peaks{idx} = get_peaks(GW_SEW_vec{idx}, "shifted");
    plot(Peaks{idx}.sin.high_mean-mean(Peaks{idx}.sin.high_mean), 'DisplayName', GW_SEW_vec{idx}.filename), hold on    %angle
    idx = idx +1;
end
phase_shift_Z1 = cellfun(@(x) x.phase_shift, GW_SEW_vec);
Peaks_Z1_sin_high_mean = cellfun(@(x) x.sin.high_mean, Peaks, 'UniformOutput', false);
Peaks_Z1_sin_high_mean = [Peaks_Z1_sin_high_mean{:,:}];

idx = 1;
figure()
for i = 1 : 1 :  length(dataListing_SEW_B) 
    GW_SEW_B_vec{idx} = AngleSensor_SEW(dataListing_SEW_B(i).folder, dataListing_SEW_B(i).name, "EncoderMagnet_B_SEW"); 
    GW_SEW_B_vec{idx} = shiftSignals(GW_SEW_B_vec{idx}, GW_SEW_B_vec{1});
    plot(GW_SEW_B_vec{idx}.sin_shifted, 'DisplayName', GW_SEW_B_vec{idx}.filename), hold on    %angle
    idx = idx +1;
end
phase_shift_B = cellfun(@(x) x.phase_shift, GW_SEW_B_vec);

%

Peaks2plot = Peaks_Z3_sin_high_mean;
x = 1 : size(Peaks2plot, 1); 
testdaysMR = linspace(0, 88, 262);
Peaks_mean = mean(Peaks2plot - mean(Peaks2plot, 1), 2)*1000; 
Peaks_std = std(Peaks2plot - mean(Peaks2plot, 1), 0, 2)*1000; 
mean_Peaks_std = mean(Peaks_std)
Peaks = (Peaks2plot - mean(Peaks2plot, 1)) * 1000;

fig = figure(), hold on
for i = 1:size(Peaks, 2)
    y = Peaks(:, i)';     
    patch([x, fliplr(x)], [y, fliplr(y)], 'r', ...
        'FaceColor', 'none', 'EdgeColor', "#4DBEEE", ...
        'FaceAlpha', 0.1, 'EdgeAlpha', 0.1);
end
errorbar(Peaks_mean, 3*Peaks_std, 'Color', 'k'), hold off
xlabel('Teeth'), ylabel('Voltage in mV'), xlim([1 size(Peaks2plot, 1)])
savePlotWithCustomSettings('Ch7-1_PhaseShift_ToothSpecificSignatureDuringLongTermTest', 0.5, 0.5, 0)
fig.Children.XLabel.String = []; fig.Children.YLabel.String = []; 
xlim([32 45]); ylim([-0.56 0.15]);
savePlotWithCustomSettings('Ch7-1_PhaseShift_ToothSpecificSignatureDuringLongTermTest_zoom', 0.25, 0.25, 0)

x = [1:2 4:263]; %delte third test from results, because there is something wrong with sin and cos raw data of sensor B
figure(), hold on
plot(testdaysMR, phase_shift_Z3(x) - phase_shift_B(x)), xlim([0 88])
xlabel('Testdays'), ylabel('Difference in phase shift')
savePlotWithCustomSettings('Ch7-1_PhaseShiftDifference_Z3-B', 0.5, 0.5, 0)

%
% Peaks edgealpha according to test day

Peaks2plot = Peaks_Z3_sin_high_mean;
Peaks = (Peaks2plot) * 1000;
x = 1 : size(Peaks, 1); 
fig = figure(), hold on
for i = 1:size(Peaks, 2)
    y = Peaks(:, i)';
    alpha = 0.2; %i / size(Peaks,2); 
    patch([x, fliplr(x)], [y, fliplr(y)], 'r', ...
        'FaceColor', 'none', 'EdgeColor', "#4DBEEE", ...
        'FaceAlpha', 0.1, 'EdgeAlpha', alpha);
    % plot(x, y)
end
xlabel('Teeth'), ylabel('Voltage in mV')
savePlotWithCustomSettings('Appendix\ToothSpecificSignature_LongTermTest', 0.8, 0.6, 0)

%
%% Torque measurement

%Measurement of torque using GWS Z1 and GE 4 
processMethods = ["preprocess"];
l = 300000; idx = 1; 
Z1_ref = AngleSensor_SEW(dataListing_SEW_Z1(1).folder, dataListing_SEW_Z1(1).name, "Z1_SEW");
GE4_ref = AngleSensor_SEW(dataListing_SEW_MR4(1).folder, dataListing_SEW_MR4(1).name, "Pos4_SEW");
GE3_ref = AngleSensor_SEW(dataListing_SEW_MR3(1).folder, dataListing_SEW_MR3(1).name, "Pos3_SEW");
figure(), hold on
for i = 1 : 1 : length(dataListing_SEW_Z3) 
    Z1 = AngleSensor_SEW(dataListing_SEW_Z1(i).folder, dataListing_SEW_Z1(i).name, "Z1_SEW");
    GE4 = AngleSensor_SEW(dataListing_SEW_MR4(i).folder, dataListing_SEW_MR4(i).name, "Pos4_SEW");
    GE3 = AngleSensor_SEW(dataListing_SEW_MR3(i).folder, dataListing_SEW_MR3(i).name, "Pos3_SEW");
    Z1 = Z1.processSinCos(processMethods);
    GE4 = GE4.processSinCos(processMethods);
    GE3 = GE3.processSinCos(processMethods);
    %Identification of phase shiftwith GWS Z1
    Z1 = Z1.shiftSignals(Z1_ref); 
    GE4 = GE4.shiftSignals(GE4_ref);
    GE3 = GE3.shiftSignals(GE3_ref);
    %Phase shift GWS Z1 and GE 4
    if GE3.phase_shift > 0 
        Z1_angle = Z1.angle_shifted(1:l);
        GE4_angle = GE4.angle_shifted(1:l);
        GE3_angle = GE3.angle_shifted(1:l);
        % plot(Z1.t(1:l), Z1_angle - GE4_angle)
        meanAngleDiffZ1GE4(idx) = mean(Z1_angle) - mean(GE4_angle);
        meanAngleDiffGE3GE4(idx) = mean(GE3_angle) - mean(GE4_angle);
        idx = idx + 1;  
    end
end
figure(), plot(meanAngleDiffZ1GE4), hold on, plot(meanAngleDiffGE3GE4)

% Identification of phase shift does not work properly for GE3 - Why?
%
%% Severity Assessment

Z1 = AngleSensor_SEW(dataListing_SEW_Z1(1).folder, dataListing_SEW_Z1(1).name, "Z1_SEW");
Z3 = AngleSensor_SEW(dataListing_SEW_Z3(1).folder, dataListing_SEW_Z3(1).name, "Z3_SEW");
HTP = 1 / Z3.quant.HTF12; 
idx_HTP = find(Z3.t > HTP); idx_HTP = idx_HTP(1); 

processMethods = ["preprocess", "fitEllipses"];
start = 1; 
ende = start + idx_HTP-1; 
step = 100;
window = 9001; 
t = Z3.t(1:idx_HTP); 

figure(), hold on
idx = 1; clear TE mean_TE std_TE
for i = 1 : length(dataListing_SEW_Z3)
    Z1 = AngleSensor_SEW(dataListing_SEW_Z1(i).folder, dataListing_SEW_Z1(i).name, "Z1_SEW");
    Z3 = AngleSensor_SEW(dataListing_SEW_Z3(i).folder, dataListing_SEW_Z3(i).name, "Z3_SEW");
    Z1 = Z1.processSinCos(processMethods);
    Z3 = Z3.processSinCos(processMethods); 
    Z3_angle = Z3.angle(start:ende) - mean(Z3.angle(start:start+step));
    Z1_angle = Z1.angle(start:ende) - mean(Z1.angle(start:start+step));
    TE(idx,:) = -Z3_angle/Z3.quant.i12 - Z1_angle;
    mean_TE(idx) = mean(TE(idx,:));
    std_TE(idx) = std(TE(idx,:));
    plot(t, TE(idx,:), 'Displayname', num2str(idx))
    idx = idx+1;
end, hold off; 
figure(), plot(mean_TE), title('Mean TE')
figure(), plot(std_TE), title('Std TE')

%
% HTP check

idx_HTP = find(Z1_ref.t > (1/Z1_ref.quant.HTF12)); idx_HTP = idx_HTP(1); 

figure(), hold on 
plot(Z1_ref.sin), plot(circshift(Z1_ref.sin, idx_HTP))

%

i = 2; processMethods = ["preprocess"]; 
Z1_ref = AngleSensor_SEW(dataListing_SEW_Z1(1).folder, dataListing_SEW_Z1(1).name, "Z1_SEW");
Z3_ref = AngleSensor_SEW(dataListing_SEW_Z3(1).folder, dataListing_SEW_Z3(1).name, "Z3_SEW");
Z1_i = AngleSensor_SEW(dataListing_SEW_Z1(i).folder, dataListing_SEW_Z1(i).name, "Z1_SEW");
Z3_i = AngleSensor_SEW(dataListing_SEW_Z3(i).folder, dataListing_SEW_Z3(i).name, "Z3_SEW");
MeasPer = 1 / Z1_ref.quant.HTF12; 
Z1_ref.MeasuringPeriod = MeasPer; Z1_i.MeasuringPeriod = MeasPer; 
Z3_ref.MeasuringPeriod = MeasPer; Z3_i.MeasuringPeriod = MeasPer;
Z1_ref = Z1_ref.processSinCos(processMethods);
Z3_ref = Z3_ref.processSinCos(processMethods);
Z1_i = Z1_i.processSinCos(processMethods);
Z3_i = Z3_i.processSinCos(processMethods);
%shift signals
Z1_i = Z1_i.shiftSignals(Z1_ref);
Z3_i = Z3_i.shiftSignals(Z3_ref); 

Z1_i_angle_phase_shifted = Z1_i.angle - Z1_i.angle(1) + Z1_i.phase_shift_rad;
Z3_i_angle_phase_shifted = -Z3_i.angle - Z3_i.angle(1) + Z3_i.phase_shift_rad; 
figure(),hold on
plot(Z1_ref.angle-Z1_ref.angle(1)), plot(Z1_i_angle_phase_shifted)
plot(-Z3_ref.angle-Z3_ref.angle(1)), plot(Z3_i_angle_phase_shifted), hold off

%wrap Z1 annd Z3
Z1_i_angle_wrapped = wrapToPi(Z1_i_angle_phase_shifted);
Z3_i_angle_wrapped = wrapToPi(Z3_i_angle_phase_shifted); 

%add wrapped Z1 and Z3 and movmean
Z1_Z3_angle_wrapped = Z1_i_angle_wrapped + Z3_i_angle_wrapped;

figure(), hold on
plot(Z1_i_angle_wrapped), plot(Z3_i_angle_wrapped), plot(Z1_Z3_angle_wrapped), hold off
%min of wrapped signals is point where both are 0 thus equal to the
%reference signal
delta_HTP_ref = find(Z1_Z3_angle_wrapped == min(Z1_Z3_angle_wrapped));

Z1_i_sin_HTP_ref = circshift(Z1_i.sin, delta_HTP_ref);
Z1_i_cos_HTP_ref = circshift(Z1_i.cos, delta_HTP_ref);
figure(), hold on
plot(Z1_ref.t, Z1_ref.sin), plot(Z1_ref.t, [Z1_i.sin(delta_HTP_ref:end) Z1_i.sin(1:delta_HTP_ref-1)])
plot(Z3_ref.t, Z3_ref.sin), plot(Z3_ref.t, [Z3_i.sin(delta_HTP_ref:end) Z3_i.sin(1:delta_HTP_ref-1)])

%wrong detection of HTP reference? How to check it? Shifted sin and cos
%signal and correct start? 
%jump in circshifted sine and cosine signal - HTP reference is not correct? to less to much?

%



%
%% Transient conditions
%% !Read data

basefolder = 'C:\Users\sn85zymy\Desktop\Diss'; 
folder2read = [basefolder '\Testdata_GSP\transient\1_GearWheelSensor'];
dataListing_PSG_GW_trans = getDataListing(folder2read);
folder2read = [basefolder '\Testdata_GSP\transient\3_PoleRingSensor'];
dataListing_PSG_PR_trans = getDataListing(folder2read);
folder2read = [basefolder '\Testdata_GSP\transient\6_Acceleration'];
dataListing_PSG_Acc_trans = getDataListing(folder2read);

%
%% Chapter 5.3.4 - Investigation of speed and torque ramps

idx2plot = [7 9];
for i = 1 : length(idx2plot)
    data2plot = AngleSensor_PSG(dataListing_PSG_GW_trans(1).folder, dataListing_PSG_GW_trans(idx2plot(i)).name, "Speed"); 
    data2plot = processSinCos(data2plot, "preprocess");
    figure()
    subplot(3,1,1)
    plot(data2plot.t, data2plot.cos), hold on, plot(data2plot.t, data2plot.sin), title(data2plot.filename)
    subplot(3,1,2)
    plot(data2plot.t,data2plot.angle)
    subplot(3,1,3)
    plot(data2plot.t(1:end-1), data2plot.IAS)
    %phase shift
    figure()
    plot(data2plot.t, data2plot.angle_deviation), ylabel('angle in rad'), xlabel('Time in s'), title('angle - linear regression of angle')
    %Peaks
    peaks = get_peaks(data2plot, "raw");
    figure(), plot(mean(peaks.sin.high, 1)), ylabel('Voltage in V'), xlabel('Number of revolutions')
    figure(), hold on
    for j = 1 : size(peaks.sin.high, 2)
        alpha = j / size(peaks.sin.high, 2);
        color = [0, 0, 1, alpha]; 
        plot(peaks.sin.high(:,j), 'Color', color)
    end
    ylabel('Voltage in V'), xlabel('Number of teeth'), hold off
end

%
%% CutSpeedRamp and perform equi-angle sampling

idx2plot = [9 10 517 519 520];
% idx2plot = 9;
for i = 172: 172 %length(dataListing_PSG_GW_trans)
    if strcmp(dataListing_PSG_GW_trans(i).Ramp, 'Speed')
        data2plot = AngleSensor_PSG(dataListing_PSG_GW_trans(1).folder, dataListing_PSG_GW_trans(i).name, "Speed"); 
        data2plot = data2plot.processSinCos("preprocess");
        Signal = movmean(data2plot.IAS, 1000); 
        t = data2plot.t(1:end-1); 
        SamplesPerSineWave = data2plot.SamplesPerSineWave;
    
        [idx_Start, idx_End, idx_0] = getSpeedRamp(Signal, t, data2plot.SamplesPerSineWave); 
    
        figure(), hold on
        plot(t(1:idx_Start), Signal(1:idx_Start), 'Color', "#D95319"), plot(t(idx_Start:idx_End), Signal(idx_Start:idx_End), 'Color', "#0072BD")
        plot(t(idx_End:end), Signal(idx_End:end), 'Color', "#D95319"), %xline(t(idx_0))
        xlabel('time in s'), ylabel('Speed in rpm')
    end
    savePlotWithCustomSettings('Ch6-2-3_SpeedRamp', 0.5, 0.3, 0)
    
% %Angular Sampling 
    phi = data2plot.angle(idx_Start : idx_End)-data2plot.angle(idx_Start);
    IAS = data2plot.IAS(idx_Start : idx_End); 
    window_length = 1001;
    delta_phi = 0.005; 
    [p, f] = fft_data(IAS, data2plot.fs);
    phi_equi = getPhi_equi(phi, delta_phi); 
    for j = 1 : size(phi_equi, 1)
        IAS_equi(j) = mean(IAS(phi_equi(j, 2) : phi_equi(j,3))); 
    end
    [p_equi, f_equi] = fft_data(IAS_equi, 1/delta_phi);
    f_equi = f_equi * pi() * data2plot.quant.fShaftRef / data2plot.quant.fShaftDmg;
    phi_plot = 0 : delta_phi : max(abs(phi_equi(:,1))) *2.5;
    figure(), 
    subplot(1,2,1), plot(t(1:length(IAS)), IAS), xlabel('Time in s')
    subplot(1,2,2), plot(abs(phi_equi(1:length(IAS_equi))), IAS_equi), ylabel('Speed in rpm'), xlabel('Angle in rad'), hold off
    savePlotWithCustomSettings('Ch6-2-3_SpeedRamp_IAS_equi', 0.7, 0.4, 790)
    % figure(), hold on, plot(f, pow2db(p)), xline(data2plot.quant.GMF), hold off
    % figure(), hold on, plot(f_equi, pow2db(p_equi)), xline(21), hold off
end

%

idx2plot = [172 172];  
figure(), hold on, 
for i = 1 : length(idx2plot)
    Obj = AngleSensor_PSG(dataListing_PSG_GW_trans(1).folder, dataListing_PSG_GW_trans(idx2plot(i)).name, "Speed"); 
    Obj = Obj.processSinCos(["preprocess", "fitEllipses"]);
    [idx_Start, idx_End, idx_0] = getSpeedRamp(Obj.IAS, Obj.t, Obj.SamplesPerSineWave); 
    Obj.PartOfSignal = [idx_Start, idx_0 - 5000]; Obj.FilterMethod = "none"; 
    Obj = Obj.processSinCos("preprocess");
    
    if i == 1 
        Signaltype = "IAS";
    else
        Signaltype = "IAS_equi"; 
    end
    [p, ~, f_ord] = calcFFT(Obj, Signaltype);
    plot(f_ord, pow2db(p))
end, hold off,  xlim([14 28]), xline([21 42 63 84]), xline([17 : 1 : 27], 'LineStyle', '--')
xlabel('Frequency order of shaft with damaged pinion'), ylabel(['Amplitude in dB']), legend(["time sampling", "equidistant angle sampling"]), xlim([16.8 25.2])
savePlotWithCustomSettings('Ch6-2-3_TimeVsEquiAnlge_FFT', 1, 0.6, 0)

%

idx2plot = [172 222]; %[19 163]; 11 155 
figure(), hold on
for i = 1 : length(idx2plot)
    Obj = AngleSensor_PSG(dataListing_PSG_GW_trans(1).folder, dataListing_PSG_GW_trans(idx2plot(i)).name, "Speed"); 
    Obj = Obj.processSinCos(["preprocess", "fitEllipses"]);
    [idx_Start, idx_End, idx_0] = getSpeedRamp(Obj.IAS, Obj.t, Obj.SamplesPerSineWave); 
    Obj.PartOfSignal = [idx_Start, idx_0 - 5000]; Obj.FilterMethod = ["none"]; 
    Obj = Obj.processSinCos("preprocess");
    Signaltype = "IAS_equi"; 
    [p, ~, f_ord] = calcFFT(Obj, Signaltype);
    plot(f_ord, pow2db(p))
    Amplitudes(i, :) = getAmplitudes(Obj, p, f_ord);
end, hold off, 
Amp_undamaged = [Amplitudes{1, 2}{1, 4:7} Amplitudes{1,1} Amplitudes{1, 2}{1, 8:11}];
Amp_damaged = [Amplitudes{2, 2}{1, 4:7} Amplitudes{2,1} Amplitudes{2, 2}{1, 8:11}];
SubstractionOfAmp = Amp_damaged - Amp_undamaged 
DivisionOfAmp = Amp_damaged ./ Amp_undamaged

%

xline([17 : 1 : 27], 'LineStyle', '--'), xlim([16.8 25.2]), xline(21), ylim([-35 25])
xlabel('Frequency order of shaft with damaged pinion'), ylabel(['Amplitude in dB'])
legend(["undamaged", "damaged"], 'location', 'northoutside')
savePlotWithCustomSettings('Ch6-2-3_Trans_17-25', 1, 0.6, 0)
% legend('off')
% for i = 1 : 9
%     xlim([15.95+i 16.05+i]), savePlotWithCustomSettings(['Ch6-2-3_Trans_' num2str(16+i)], 0.075, 0.3, 0)
% end

%
%% Chapter 6 - Generate and save feature

Overlap = 0.5; %overlap is 1- overlap 
delta_t = 1/40000; Sample2split = 40000; %1 second 
FilterMethod = 'Filter'; processMethods = ["raw", "preprocess"];
T = struct2table(dataListing_PSG_GW_trans); 
idxRed = T{:, "Ramp"} == "Speed" & (T{:, "Ritzel"} == 2 | T{:, "Ritzel"} == 3 | T{:, "Ritzel"} == 5 | T{:, "Ritzel"} == 6); 
T = T(idxRed, :); 
dataList_PSG_GW_trans = dataListing_PSG_GW_trans(idxRed, :);
dataList_PSG_Acc_trans = dataListing_PSG_Acc_trans(idxRed, :);  
dataList_PSG_PR_trans = dataListing_PSG_PR_trans(idxRed, :); 
%Gear wheel sensor
t1 = datetime; 
damageLabels = table(); clear statFeatures_GW freqFeatures_GW statFeatures_Acc freqFeatures_Acc dataList 
idx2save = 1; 
for i = 1 : length(dataList_PSG_GW_trans)
%gear wheel sensor
    ObjGW = AngleSensor_PSG(dataList_PSG_GW_trans(1).folder, dataList_PSG_GW_trans(i).name, "Speed");
    ObjGW.FilterMethod = FilterMethod; ObjGW.delta_t = delta_t;
    ObjGW = ObjGW.processSinCos(processMethods);

    ObjPR = AngleSensor_PSG(dataList_PSG_PR_trans(i).folder, dataList_PSG_PR_trans(i).name, "AL795");
    ObjPR.FilterMethod = FilterMethod; ObjPR.delta_t = delta_t;

    ObjAcc = AccelerationSensor_PSG(dataList_PSG_Acc_trans(i).folder, dataList_PSG_Acc_trans(i).name, "side");
    ObjAcc.FilterMethod = FilterMethod; ObjAcc.delta_t = delta_t;

    [idx_Start, idx_End, ~] = getSpeedRamp(ObjGW.IAS, ObjGW.t, ObjGW.SamplesPerSineWave);
    deltaRamp = idx_End - idx_Start;
    n_split = floor(deltaRamp/Sample2split); 
    for j = 1 : n_split/Overlap - 1 %split signal to length of 1 s with overlap of signal
        idx_Start_temp = idx_Start+(j-1)*Sample2split*Overlap;
        idx_End_temp = idx_Start_temp + Sample2split; 
        ObjGW.PartOfSignal = [idx_Start_temp, idx_End_temp];
        ObjGW = ObjGW.processSinCos(["raw", "preprocess"]); 
        
        ObjPR.PartOfSignal = [idx_Start_temp, idx_End_temp];
        ObjPR = ObjPR.processSinCos(["raw", "preprocess"]);
            
        if (ObjGW.quant.SamplesPerRevDamagedPinion / ObjGW.fs) < 1 %exclude signal around 0 rpm, with less then one rev in one second         
            statFeatures_GW(idx2save, :) = calcStatFeatures(ObjGW, ["sin", "cos", "IAS"]); 
            freqFeatures_GW(idx2save, :) = calcFreqFeatures(ObjGW, ["sin_equi", "cos_equi", "IAS_equi"]); 
            damageLabels.damageLabels(idx2save) = labelDamages(ObjGW);
            dataList(idx2save, :) = T(i, :); 

            statFeatures_PR(idx2save, :) = calcStatFeatures(ObjPR, ["sin", "cos", "IAS"]); 
            freqFeatures_PR(idx2save, :) = calcFreqFeatures(ObjPR, ["sin_equi", "cos_equi", "IAS_equi"]); 

            % ObjAcc.PartOfSignal = [idx_Start_temp, idx_End_temp];
            % statFeatures_Acc(idx2save, :) = calcStatFeatures(ObjAcc, "Signal"); 
            % freqFeatures_Acc(idx2save, :) = calcFreqFeatures(ObjAcc, "Signal"); 

            idx2save = idx2save + 1;
        end
    end
end
Features_GW_PSG_trans = [dataList damageLabels statFeatures_GW freqFeatures_GW]; 
save([basefolder '\Testdata_GSP\Features\transient\Features_GW_PSG_trans_2356_deltaPhi005.mat'], 'Features_GW_PSG_trans')

Features_PR_PSG_trans = [dataList damageLabels statFeatures_PR freqFeatures_PR]; 
save([basefolder '\Testdata_GSP\Features\transient\Features_PR_PSG_trans_2356_deltaPhi005.mat'], 'Features_PR_PSG_trans')

% Features_AccSide_PSG_trans = [dataList damageLabels statFeatures_Acc freqFeatures_Acc]; 
% save([basefolder '\Testdata_GSP\Features\transient\Features_AccSide_PSG_trans_2356.mat'], 'Features_AccSide_PSG_trans')
datetime - t1 %1:07 h

%
%% Chapter 6.6 - Classification at transient conditions

%Gear wheel sensor
Features_GW = load([basefolder '\Testdata_GSP\Features\transient\Features_GW_PSG_trans_2356_deltaPhi005.mat']);
Features_GW = Features_GW.(cell2mat(fieldnames(Features_GW))); 

splitFeaturesLabels = 15; %equi angle frequency features = 90
idx = Features_GW.WDH ~= 3 & Features_GW.damageLabels ~= 4.1 & Features_GW.Rev ~= 300; 
cv = cvpartition(Features_GW.damageLabels == 0, 'Holdout', 0.5, 'Stratify', true); 
idxDmg012 = Features_GW.damageLabels == 0 | Features_GW.damageLabels == 1.4 | Features_GW.damageLabels == 2; idxDmg345 = ~idxDmg012; 
idx = ((test(cv) & idxDmg012 | idxDmg345) & idx);

FeatureTable = Features_GW(idx, splitFeaturesLabels:end); 
allFeatures = table2array(bringAllTableDataInFirstLevel(FeatureTable));
allLabels = Features_GW(idx, 1:14);

% [allFeatures, allLabels] = normFeatures(allFeatures, allLabels, [1300, 2000], [2, 3, 5, 6], 0);

labels = allLabels{:, "damageLabels"};
[labels, labelNames, ~] = combineLabels(labels); 
accuracy = classifyRandomForest(allFeatures, labels, labelNames, 'KFold', []); 
title({'Gear wheel sensor', ['Accuracy: ', num2str(round(accuracy,2) * 100), '%']});
savePlotWithCustomSettings('Ch6-2-4_Trans_GW_5-fold', 0.53, 0.5, 0)

% trainIdx = [allLabels.Ritzel == 2 | allLabels.Ritzel == 5, allLabels.Ritzel == 2 | allLabels.Ritzel == 6, allLabels.Ritzel == 3 | allLabels.Ritzel == 5, allLabels.Ritzel == 3 | allLabels.Ritzel == 6]; 
% testIdx = ~trainIdx;
% accuracy = classifyRandomForest(allFeatures, labels, labelNames, trainIdx, testIdx);
% title(['Pinion split, Accuracy: ', num2str(round(accuracy,2) * 100), '%']);

%Pole ring sensor
Features_PR = load([basefolder '\Testdata_GSP\Features\transient\Features_PR_PSG_trans_2356_deltaPhi005.mat']);
Features_PR = Features_PR.(cell2mat(fieldnames(Features_PR))); 

FeatureTable = Features_PR(idx, splitFeaturesLabels:end); %9-Quadrant, 11-Rev
allFeatures = table2array(bringAllTableDataInFirstLevel(FeatureTable));
allLabels = Features_PR(idx, 1:14);

% [allFeatures, allLabels] = normFeatures(allFeatures, allLabels, [1300, 2000], [2, 3, 5, 6], 0);

labels = allLabels{:, "damageLabels"};
[labels, labelNames, ~] = combineLabels(labels); 
accuracy = classifyRandomForest(allFeatures, labels, labelNames, 'KFold', []); 
title({'Pole ring sensor', ['Accuracy: ', num2str(round(accuracy,2) * 100), '%']});
savePlotWithCustomSettings('Ch6-2-4_Trans_PR_5-fold', 0.53, 0.5, 0)

% trainIdx = [allLabels.Ritzel == 2 | allLabels.Ritzel == 5, allLabels.Ritzel == 2 | allLabels.Ritzel == 6, allLabels.Ritzel == 3 | allLabels.Ritzel == 5, allLabels.Ritzel == 3 | allLabels.Ritzel == 6]; 
% testIdx = ~trainIdx;
% accuracy = classifyRandomForest(allFeatures, labels, labelNames, trainIdx, testIdx);
% title(['Pinion split, Accuracy: ', num2str(round(accuracy,2) * 100), '%']);

%Acceleration sensor
Features_Acc = load([basefolder '\Testdata_GSP\Features\transient\Features_AccSide_PSG_trans_2356.mat']);
Features_Acc = Features_Acc.(cell2mat(fieldnames(Features_Acc))); 

FeatureTable = Features_Acc(idx, 15:end); %9-Quadrant, 11-Rev
allFeatures = table2array(bringAllTableDataInFirstLevel(FeatureTable));
allLabels = Features_Acc(idx, 1:14);

labels = allLabels{:, "damageLabels"};
[labels, labelNames, ~] = combineLabels(labels); 
accuracy = classifyRandomForest(allFeatures, labels, labelNames, 'KFold', []); 
title({'Acceleration sensor', ['Accuracy: ', num2str(round(accuracy,2) * 100), '%']});
savePlotWithCustomSettings('Ch6-2-4_Trans_Acc_5-fold', 0.53, 0.5, 0)

% trainIdx = [allLabels.Ritzel == 2 | allLabels.Ritzel == 5, allLabels.Ritzel == 2 | allLabels.Ritzel == 6, allLabels.Ritzel == 3 | allLabels.Ritzel == 5, allLabels.Ritzel == 3 | allLabels.Ritzel == 6]; 
% testIdx = ~trainIdx;
% accuracy = classifyRandomForest(allFeatures, labels, labelNames, trainIdx, testIdx);
% title(['Pinion split, Accuracy: ', num2str(round(accuracy,2) * 100), '%']);

createColorBarForConfChart(360, [0 0.4470 0.7410])
savePlotWithCustomSettings('Ch6-2-4_Trans_ColorBarCorrect', 0.4, 0.46, 0)
createColorBarForConfChart(360, [0.8500 0.3250 0.0980])
savePlotWithCustomSettings('Ch6-2-4_Trans_ColorBarFalse', 0.4, 0.46, 0)

%
%% Chapter 7.2 - SEW tests at torque ramps

directory = 'C:\Users\sn85zymy\Desktop\Diss\Testdata_SEW\03_Vorversuche_Drehzahl_Drehmomentrampen\';
name = {'2021-10-12 14-52-47 (6)_Drehmomentrampe_4_kalt', '2021-10-12 17-14-05 (4)_Drehmomentrampe_4_warm'};
f_cutoff = 100; 
torqueNorm = 7000; 
cutidx = 80; %Signal of MayrRoba and Speed signal is shifted by approximately 0.08 s -> 80 points
for i = 1 : length(name)
    folder = fullfile(directory, name{i}); 
    MayrRoba = load(fullfile(folder, 'Mayr_Roba_DSM_abtrieb.mat')); MayrRoba = MayrRoba.(cell2mat(fieldnames(MayrRoba)));
    Speed = load(fullfile(folder, 'Istdrehzahl_FU.mat')); Speed = Speed.(cell2mat(fieldnames(Speed))); 
    figure()
    subplot(2,1,1)
    plot(Speed.t(cutidx:end) - Speed.t(cutidx), Speed.Data(cutidx:end)), ylabel('Speed in rpm')
    yyaxis right
    plot(MayrRoba.t(cutidx:end) - MayrRoba.t(cutidx), MayrRoba.Data(cutidx:end)/torqueNorm), ylabel('Torque normalised'), xlim([0 30]),% title(name{i})
    grid on
    Z1 = AngleSensor_SEW(folder, 'MR_Z1.mat', "Z1_SEW"); Z1.MeasuringPeriod = 30; %Z1.f_cutoff = f_cutoff;
    Z1 = processSinCos(Z1, "preprocess"); 
    Z3 = AngleSensor_SEW(folder, 'MR_Z3.mat', "Z3_SEW"); Z3.MeasuringPeriod = 30; %Z3.f_cutoff = f_cutoff;
    Z3 = processSinCos(Z3, "preprocess"); 
    MR3 = AngleSensor_SEW(folder, 'MR_Pos3.mat', "Pos3_SEW"); MR3.MeasuringPeriod = 30; %MR3.f_cutoff = f_cutoff; 
    MR3 = processSinCos(MR3, "preprocess"); 
    MR4 = AngleSensor_SEW(folder, 'MR_Pos4.mat', "Pos4_SEW"); MR4.MeasuringPeriod = 30; %MR4.f_cutoff = f_cutoff; 
    MR4 = processSinCos(MR4, "preprocess");
    
    norm = 1/3.81; 
    Z1_angle = sgolayfilt(Z1.angle, 3, 9001)*1000*norm;
    Z3_angle = sgolayfilt(Z3.angle, 3, 9001)*1000*norm;
    MR3_angle = MR3.angle*1000*norm;
    MR4_angle = MR4.angle*1000*norm;
    subplot(2,1,2), hold on
    % plot(Z1.t, Z1_angle - MR3_angle - mean(Z1_angle - MR3_angle), 'Displayname', 'GWS Z1 - GE 3')
    plot(Z1.t, Z1_angle - MR4_angle - mean(Z1_angle(1:10000) - MR4_angle(1:10000)), 'Displayname', 'GE 4 - GWS Z1')
    % plot(Z1.t, MR3_angle - MR4_angle - mean(MR3_angle - MR4_angle), 'Displayname', 'GE 3 - GE 4')
    plot(Z1.t, -Z3_angle/Z3.quant.i12 - MR4_angle - mean(-Z3_angle(1:10000)/Z3.quant.i12 - MR4_angle(1:10000)), 'Displayname', 'GE 4 - GWS Z3')
    plot(Z1.t, -Z3_angle/Z3.quant.i12 - Z1_angle - mean(-Z3_angle(1:10000)/Z3.quant.i12 - Z1_angle(1:10000)), 'Displayname', 'GWS Z1 - GWS Z3')
    % plot(Z3.t, -Z3_angle/Z3.quant.i12 - MR3_angle - mean(-Z3_angle/Z3.quant.i12 - MR3_angle), 'Displayname', 'GE 3 - GWS Z3')
    ylabel({'Normalsied rotational', 'angle difference'}), xlabel('Time in s'), %legend('location', 'northeastoutside')
    grid on
end
savePlotWithCustomSettings('Ch7-2_TorqueRamp_AngleDifference', 1, 0.75, 0)

%
% Cut signal and use absolute reference for the angle calculation for GWS Z3 and GE 4

directory = fullfile(basefolder, '\Testdata_SEW\03_Vorversuche_Drehzahl_Drehmomentrampen\');
name = {'2021-10-12 14-52-47 (6)_Drehmomentrampe_4_kalt'};
folder = fullfile(directory, name{1}); 
MayrRoba = load(fullfile(folder, 'Mayr_Roba_DSM_abtrieb.mat')); MayrRoba = MayrRoba.(cell2mat(fieldnames(MayrRoba)));
Z1 = AngleSensor_SEW(folder, 'MR_Z1.mat', "Z1_SEW"); 
MR4 = AngleSensor_SEW(folder, 'MR_Pos4.mat', "Pos4_SEW"); 
MR3 = AngleSensor_SEW(folder, 'MR_Pos3.mat', "Pos3_SEW"); 


t_start = [0; 3; 7; 11; 15; 19; 23]; %in s ; 27
t_end = t_start + 9.1; %+ 4.1s to have at least one full revolution
t = Z1.t; 
figure(), hold on
for i = 1 : length(t_start)
    Z1_i = Z1; 
    MR4_i = MR4; 
    MR3_i = MR3;
    %Cut signal for each stage with at least one revolution --> more than 4.1 s
    idx_start = find(t > t_start(i)); 
    idx_end = find(t > t_end(i)); 
    Z1_i.PartOfSignal = [idx_start(1) idx_end(1)];
    MR4_i.PartOfSignal = [idx_start(1) idx_end(1)];
    MR3_i.PartOfSignal = [idx_start(1) idx_end(1)];
    Z1_i = Z1_i.processSinCos("preprocess"); 
    MR4_i = MR4_i.processSinCos("preprocess");
    MR3_i = MR3_i.processSinCos("preprocess");
    %set first as reference
    if i == 1
        Z1_ref = Z1_i;
        MR4_ref = MR4_i; 
        MR3_ref = MR3_i;
    end
    
    %calcualte phase shift for all data, shift signals and calc phase
    %difference
    Z1_i = Z1_i.shiftSignals(Z1_ref); 
    MR4_i = MR4_i.shiftSignals(MR4_ref); 
    MR3_i = MR3_i.shiftSignals(MR3_ref); 

    %plot results
    subplot(3,1,1)
    plot((Z1_i.angle -Z1_i.angle(1) - Z1_i.phase_shift_rad) - (MR4_i.angle - MR4_i.angle(1) - MR4_i.phase_shift_rad)), hold on
    
    subplot(3,1,2)
    plot((Z1_i.angle -Z1_i.angle(1) - Z1_i.phase_shift_rad) - (MR3_i.angle - MR3_i.angle(1) - MR3_i.phase_shift_rad)), hold on

    subplot(3,1,3)
    plot((MR4_i.angle - MR4_i.angle(1) - MR4_i.phase_shift_rad) - (MR3_i.angle - MR3_i.angle(1) - MR3_i.phase_shift_rad)), hold on

    Z1_phase_shift(i) = Z1_i.phase_shift;
    MR4_phase_shift(i) = MR4_i.phase_shift; 
    MR3_phase_shift(i) = MR3_i.phase_shift; 
end 
figure(), hold on
plot(Z1_phase_shift), plot(MR4_phase_shift), plot(MR3_phase_shift)

%
%% Functions

function dataListing = getDataListing(folder2read)
    dataListing = dir(fullfile(folder2read));
    dataListing = dataListing(~ismember({dataListing.name}, {'.', '..'}));
    dataListing = extractInformationOfFilename(dataListing); %contains all names in specified folder
    dataListing = dataListing([dataListing.Trq] ~= 0); 
end

function dataListing = extractInformationOfFilename(dataListing)
    variableNames = {'Ritzel', 'Damage', 'Quadrant', 'Trq', 'Rev', 'WDH', 'Ramp'};
    for i = 1 : length(variableNames)
        for j = 1 : length(dataListing)
            filename = dataListing(j).name;
            searchname = variableNames{i}; 

            idxstart = strfind(filename, searchname)+length(searchname);
            if searchname == "Ritzel" || searchname == "Damage" || searchname == "Quadrant"
                idxend = strfind(filename(idxstart:end), '_'); %for Ritzel, Damage, Quadrant
            end
            if searchname == "Trq"
                idxend = strfind(filename(idxstart:end), 'Nm'); %for Trq
            end
            if searchname == "Rev"
                idxend = strfind(filename(idxstart:end), 'rpm'); %for Rev
            end
            if searchname == "WDH"
                % idxend = strfind(filename(idxstart:end), '.');
                idxend = 2; %for Wdh (one digit)
            end

            if searchname == "Ramp"
                idxend = strfind(filename, searchname)-1;
                if ~isempty(idxend)
                    idxTemp = strfind(filename(1:idxend), '_');
                    idxstart = idxTemp(end)+1;
                    dataListing(j).(searchname) = filename(idxstart:idxend); 
                end
            else
                idxend = idxstart + idxend -2;
                dataListing(j).(searchname) = str2double(filename(idxstart:idxend)); 
            end            
        end
        if searchname == "WDH" %change NaN to 1 for first Ramp tests
            for k = 1 : length(dataListing)
                if isnan(dataListing(k).WDH)
                    dataListing(k).WDH = 1; 
                end
            end
        end
    end
end

function [labels, labelNames, labelClasses] = combineLabels(labelsBasic)
    %combine  labels and output labels with all labels, labelNames, with the combined names and
    %labelClasses with the classes from 1 : n. The rows match between these two
    %classes
    if isa(labelsBasic, 'table')
        labelsBasic = table2array(labelsBasic);
    end
    idx = 0; labelClasses = [0; 0]; labels = [0;0]; labelNames = {0};
    for i = 1 : size(labelsBasic, 1)
        if size(labelsBasic, 2) == 1
            name = ['Dmg ' num2str(labelsBasic(i, 1))];
        elseif size(labelsBasic, 2) == 2
            name = ['Pinion ' num2str(labelsBasic(i, 1)) ', Damage ' num2str(labelsBasic(i, 2))];
        else
            error(['Size ' num2str(size(labelsBasic,2)) ' of labels to combine not implemented'])
        end
        if i == 1 %initialise
            idx = idx+1;
            labelNames{idx} = name; 
            labelClasses(idx) = idx; 
            labels(i) = labelClasses(idx); 
            idx = idx+1; 
        end
        if any(strcmp(labelNames, name)) %assign label class
            labels(i) = labelClasses(strcmp(labelNames, name));
        else %add new label class and assign it
            labelNames{idx} = name; 
            labelClasses(idx) = idx; 
            labels(i) = labelClasses(idx); 
            idx = idx+1; 
        end
    end
end



% function data = openDataForAngleSensor_PSG(folder, filename)
%     %opens data and assigns the variables to data so it can be used as
%     %input for AngleSensor_PSG
%     path = fullfile(folder, name);
%     data = open(path); data = data.data;  
%     (dataListing_GW(i).name, data.oben, data.unten, ...
%         data.t, data.wf_increment, data.wf_start_time, "Speed"); 
% 
% end