%% Process listening experiment repsonses 
% Demographics of participants
% BvdH: 13-06-2024

clear all; close all; clc

responses = readtable('responses_table.xlsx','VariableNamingRule','preserve');

%% Variables

nParticipants = size(responses, 1);

variables = {
    'age', 
    'gender', 
    'education', 
    'employment', 
    'how_found', 
    'affiliation', 
    'hearing_rating', 
    'accident', 
    'hearing_aid', 
    'wearing_aid', 
    'blast_trauma', 
    'work_protection', 
    'which_protection', 
    'ear_disease', 
    'which_disease', 
    'tinnitus', 
    'cold', 
    'cold_magnitude', 
    'feeling', 
    'tired'};

%% Age

agePivotTable = pivot(responses, Columns = 'age', Rows = 'gender');

figure
boxplot(responses.age, responses.gender)
% violinplot(responses.age,responses.gender)
age.mean = mean(responses.age);
age.std = std(responses.age);


%% Pivot table per variable

for i = 1:length(variables)
    variables{i,2} = pivot(responses, Columns = variables{i});
end



%% Group and average answers from repeat question
% 'bieke-Flyover_Velis/Virus_#flyover_00/01repetition-q1'

Aircraft_Table = table('Size', [nParticipants 6], ...
    'VariableTypes', {'double', 'double', 'double', 'double', 'double', 'double'}, ...
    'VariableNames', {'Velis_1', 'Velis_2', 'Velis_3', 'Virus_1', 'Virus_2', 'Virus_3'});

Difference = table('Size', [nParticipants 6], ...
    'VariableTypes', {'double', 'double', 'double', 'double', 'double', 'double'}, ...
    'VariableNames', {'Velis_1', 'Velis_2', 'Velis_3', 'Virus_1', 'Virus_2', 'Virus_3'});

% Group Velis answers
idxQuestion = contains(responses.Properties.VariableNames, 'Velis_1');
Annoyance_Velis_1 = table2array(responses(:, idxQuestion));
Difference.Velis_1 = diff(Annoyance_Velis_1, [], 2);
Aircraft_Table.Velis_1 = mean(Annoyance_Velis_1, 2);
[h(1,1), p(1,1)] = ttest(Annoyance_Velis_1(:,1), Annoyance_Velis_1(:,2), 'tail', 'both'); % tail=both  --> H0 = means equal, H1 = means not equal

idxQuestion = contains(responses.Properties.VariableNames, 'Velis_2');
Annoyance_Velis_2 = table2array(responses(:, idxQuestion));
Difference.Velis_2 = diff(Annoyance_Velis_2, [], 2);
Aircraft_Table.Velis_2 = mean(Annoyance_Velis_2, 2);
[h(2,1), p(2,1)] = ttest(Annoyance_Velis_2(:,1), Annoyance_Velis_2(:,2), 'tail', 'both'); % tail=both  --> H0 = means equal, H1 = means not equal


idxQuestion = contains(responses.Properties.VariableNames, 'Velis_3');
Annoyance_Velis_3 = table2array(responses(:, idxQuestion));
Difference.Velis_3 = diff(Annoyance_Velis_3, [], 2);
Aircraft_Table.Velis_3 = mean(Annoyance_Velis_3, 2);
[h(3,1), p(3,1)] = ttest(Annoyance_Velis_3(:,1), Annoyance_Velis_3(:,2), 'tail', 'both'); % tail=both  --> H0 = means equal, H1 = means not equal


figure
boxplot([Annoyance_Velis_1, Annoyance_Velis_2, Annoyance_Velis_3])
hold on
plot(1:2, mean(Annoyance_Velis_1, 1), '-k', 'Marker', 'o')
plot(3:4, mean(Annoyance_Velis_2, 1), '-k', 'Marker', 'o')
plot(5:6, mean(Annoyance_Velis_3, 1), '-k', 'Marker', 'o')
set(gca, 'XTickLabel', {'Velis 1 v1', 'Velis 1 v2', 'Velis 2 v1', 'Velis 2 v2', 'Velis 3 v1', 'Velis 3 v2'} )
xline([2.5 4.5])
ylabel('Rated annoyance')

% Group Virus answers
idxQuestion = contains(responses.Properties.VariableNames, 'Virus_1');
Annoyance_Virus_1 = table2array(responses(:, idxQuestion));
Difference.Virus_1 = diff(Annoyance_Virus_1, [], 2);
Aircraft_Table.Virus_1 = mean(Annoyance_Virus_1, 2);
[h(1,2), p(1,2)] = ttest(Annoyance_Virus_1(:,1), Annoyance_Virus_1(:,2), 'tail', 'both'); % tail=both  --> H0 = means equal, H1 = means not equal

idxQuestion = contains(responses.Properties.VariableNames, 'Virus_2');
Annoyance_Virus_2 = table2array(responses(:, idxQuestion));
Difference.Virus_2 = diff(Annoyance_Virus_2, [], 2);
Aircraft_Table.Virus_2 = mean(Annoyance_Virus_2, 2);
[h(2,2), p(2,2)] = ttest(Annoyance_Virus_2(:,1), Annoyance_Virus_2(:,2), 'tail', 'both'); % tail=both  --> H0 = means equal, H1 = means not equal

idxQuestion = contains(responses.Properties.VariableNames, 'Virus_3');
Annoyance_Virus_3 = table2array(responses(:, idxQuestion));
Difference.Virus_3 = diff(Annoyance_Virus_3, [], 2);
Aircraft_Table.Virus_3 = mean(Annoyance_Virus_3, 2);
[h(3,2), p(3,2)] = ttest(Annoyance_Virus_3(:,1), Annoyance_Virus_3(:,2), 'tail', 'both'); % tail=both  --> H0 = means equal, H1 = means not equal

hConsistency = array2table(h, 'VariableNames', {'Velis', 'Virus'}, 'RowNames', {'f-o 1','f-o 2','f-o 3'});
pConsistency = array2table(p, 'VariableNames', {'Velis', 'Virus'}, 'RowNames', {'f-o 1','f-o 2','f-o 3'});

figure
boxplot([Annoyance_Virus_1, Annoyance_Virus_2, Annoyance_Virus_3])
hold on
plot(1:2, mean(Annoyance_Virus_1, 1), '-k', 'Marker', 'o')
plot(3:4, mean(Annoyance_Virus_2, 1), '-k', 'Marker', 'o')
plot(5:6, mean(Annoyance_Virus_3, 1), '-k', 'Marker', 'o')
set(gca, 'XTickLabel', {'Virus 1 v1', 'Virus 1 v2', 'Virus 2 v1', 'Virus 2 v2', 'Virus 3 v1', 'Virus 3 v2'} )
xline([2.5 4.5])
ylabel('Rated annoyance')


%% Responses for questionnaire Bieke

% The next part of this survey inquires about the annoyance of aircraft, both electric and non-electric.
% This block contains 12 questions.
% While listening, imagine that this is the sound situation in your garden.
% What number from 0 to 10 best shows how much you would be bothered, disturbed or annoyed by this sound?

%% Check consistency of answers per participant
Difference = abs(Difference);
maxDifference = max(table2array(max(Difference)));
minDifference = min(table2array(min(Difference)));
range = minDifference : maxDifference;
idxGoodRange = find(range<2.05);

colorBad = [1 0 0];
colorGood = [0 1 0];

colormapBieke = ones(length(range), 3) .* colorBad;
colormapBieke(idxGoodRange, :) = ones(length(idxGoodRange), 3) .* colorGood ;
% cannot say anything about the positive/negative difference, because you don't know which one came first per person.  Might be that repetition 2 in responses was actually the first
% So maybe I need to do abs() on everythinig in general

figure
imagesc(1:6, 1:nParticipants, table2array(Difference)) % pcolor not an option "Because of this relationship between the vertex colors and face colors, none of the values in the last row and column of C are represented in the plot"
colormap( colormapBieke )
colorbar
axis square
xlabel('Flyover measurement')
ylabel('Participant number')
set(gca, 'XTickLabel', {'Velis 1', 'Velis 2', 'Velis 3', 'Virus 1', 'Virus 2', 'Virus 3'} )

% How many percent has offset of larger than 2?
Difference = Difference{:,:};
largeOffset = sum(sum(Difference>2 | Difference<-2));
percentageOffset = largeOffset/numel(Difference)*100;


%% Correlation between participants
% high correlation means they score similarly
responses_BvdH = responses(:,44:55);
correlation_BvdH = corrcoef(responses_BvdH{:,:}.');

correlation_participants_average = mean(mean(correlation_BvdH));

figure
heatmap(correlation_BvdH)
colormap(flipud(hot))
clim([-1 1])

%% Check mean, std, boxplot of responses Velis and Virus
meanAnnoyance = mean(Aircraft_Table, 1);
stdAnnoyance = std(Aircraft_Table, 1);

figure
boxplot(table2array(Aircraft_Table))
hold on
plot(1:3, table2array(meanAnnoyance(1,1:3)), '-k', 'Marker', 'o')
xline(3.5,'-k')
plot(4:6, table2array(meanAnnoyance(1,4:6)), '-k', 'Marker', 'o')
set(gca, 'XTickLabel', {'Velis 1', 'Velis 2', 'Velis 3', 'Virus 1', 'Virus 2', 'Virus 3'} )
ylabel('Rated annoyance')

%% Normalised annoyance for flyovers of Velis and Virus
participant_mean_RA = mean(responses_BvdH{:,:},2);
Aircraft_Table_Normalised = Aircraft_Table./participant_mean_RA;

normalised_meanAnnoyance = mean(Aircraft_Table_Normalised, 1);

figure
boxplot(table2array(Aircraft_Table_Normalised))
hold on
plot(1:3, table2array(normalised_meanAnnoyance(1,1:3)), '-k', 'Marker', 'o')
xline(3.5,'-k')
plot(4:6, table2array(normalised_meanAnnoyance(1,4:6)), '-k', 'Marker', 'o')
set(gca, 'XTickLabel', {'Velis 1', 'Velis 2', 'Velis 3', 'Virus 1', 'Virus 2', 'Virus 3'} )
ylabel('Normalised rated annoyance (-)')
h1 = get(gcf,'Children');
set(h1,'FontSize',16)


%% Statistical testing
% Is the Velis scored significantly less annoying by the participants?
% One-sided (only 'less' annoying, not 'more annoying' as well), t-test (between two groups of results)
Annoyance_Velis_all = mean(Aircraft_Table{:,1:3},2);
Annoyance_Virus_all = mean(Aircraft_Table{:,4:6},2);

figure
boxplot([Annoyance_Velis_all, Annoyance_Virus_all])
set(gca, 'XTickLabel', {'Velis', 'Virus'} )
ylabel('Rated annoyance')
h1 = get(gcf,'Children');
set(h1,'FontSize',14)

[h_Velis_Virus, p_Velis_Virus] = ttest(Annoyance_Velis_all, Annoyance_Virus_all, 'tail', 'left'); % tail = left --> H0 = means equal, H1 = means mu_x < mu_y


% Check also all individual combinations
[h_VV(1,1), p_VV(1,1)] = ttest(Aircraft_Table{:,1}, Aircraft_Table{:,4}, 'tail', 'left'); % tail = left --> H0 = means equal, H1 = means mu_x < mu_y
[h_VV(1,2), p_VV(1,2)] = ttest(Aircraft_Table{:,1}, Aircraft_Table{:,5}, 'tail', 'left'); % tail = left --> H0 = means equal, H1 = means mu_x < mu_y
[h_VV(1,3), p_VV(1,3)] = ttest(Aircraft_Table{:,1}, Aircraft_Table{:,6}, 'tail', 'left'); % tail = left --> H0 = means equal, H1 = means mu_x < mu_y

[h_VV(2,1), p_VV(2,1)] = ttest(Aircraft_Table{:,2}, Aircraft_Table{:,4}, 'tail', 'left'); % tail = left --> H0 = means equal, H1 = means mu_x < mu_y
[h_VV(2,2), p_VV(2,2)] = ttest(Aircraft_Table{:,2}, Aircraft_Table{:,5}, 'tail', 'left'); % tail = left --> H0 = means equal, H1 = means mu_x < mu_y
[h_VV(2,3), p_VV(2,3)] = ttest(Aircraft_Table{:,2}, Aircraft_Table{:,6}, 'tail', 'left'); % tail = left --> H0 = means equal, H1 = means mu_x < mu_y

[h_VV(3,1), p_VV(3,1)] = ttest(Aircraft_Table{:,3}, Aircraft_Table{:,4}, 'tail', 'left'); % tail = left --> H0 = means equal, H1 = means mu_x < mu_y
[h_VV(3,2), p_VV(3,2)] = ttest(Aircraft_Table{:,3}, Aircraft_Table{:,5}, 'tail', 'left'); % tail = left --> H0 = means equal, H1 = means mu_x < mu_y
[h_VV(3,3), p_VV(3,3)] = ttest(Aircraft_Table{:,3}, Aircraft_Table{:,6}, 'tail', 'left'); % tail = left --> H0 = means equal, H1 = means mu_x < mu_y


hDifferenceVelisVirus = array2table(h_VV, 'VariableNames', {'Virus 1', 'Virus 2', 'Virus 3'}, 'RowNames', {'Velis 1','Velis 2','Velis 3'});
pDifferenceVelisVirus = array2table(p_VV, 'VariableNames', {'Virus 1', 'Virus 2', 'Virus 3'}, 'RowNames', {'Velis 1','Velis 2','Velis 3'});

% figure
% imagesc(pDifferenceVelisVirus{:,:})
% colormap(summer)

figure
heatmap(pDifferenceVelisVirus{:,:})
colormap(summer)
clim([0 0.05])
axs = struct(gca); %ignore warning that this should be avoided
cb = axs.Colorbar;
cb.Ticks = {0, 0.1, 0.2, 0.3, 0.4, 0.5};
cb.Labels = {'p-value'};
axs.XDisplayLabels('Virus 1', 'Virus 2', 'Virus 3')
axs.YDisplayLabels('Velis 1', 'Velis 2', 'Velis 3')


%% Check trends - PA models vs. responses - best match?
load('Pipistrel_SQM.mat','table')
% table order: N5, K5, S5, FS5, R5, I, PADi, PAWillemsen, PAMore, PAZwicker

entriesTable = {'N_5 (sone)', 'K_5 (tu)', 'S_5 (acum)', 'FS_5 (vacil)', 'R_5 (asper)', 'I (sone)', 'PA_{Di,2016} (-)', 'PA_{Willemsen,2010} (-)', 'PA_{More,2010} (-)', 'PA_{Zwicker,1999} (-)'};

for i = 1:length(entriesTable)
    figure
    scatter(table(1:3,i), meanAnnoyance{1,1:3},'filled','Marker','o')
    hold on
    scatter(table(4:6,i), meanAnnoyance{1,4:6},'filled','Marker','d')
    xlabel(entriesTable{i},'FontSize',14)
    ylabel('Mean noise annoyance (-)','FontSize',14)
    grid on
    h1 = get(gcf,'Children');
    set(h1,'FontSize',14)
    
    p = polyfit(table(:,i), meanAnnoyance{:,:},1); % order 1 indicates a linear fit
    new_line = polyval(p, sort(table(:,i)));
    plot(sort(table(:,i)), new_line ,'--k','LineWidth',1.5)

    legend('Velis','Virus','Least-squares fit','FontSize',14,'location','best')

    saveas(gcf,['Annoyance - ' entriesTable{i}],'fig')
    saveas(gcf,['Annoyance - ' entriesTable{i}],'png')

    [rho_intermediate, p_intermediate] = corrcoef(table(:,i), meanAnnoyance{:,:});

    rho(i) = rho_intermediate(2,1);
    p_value(i) = p_intermediate(2,1);

    close gcf
end

%% Check trends - which SQM predicts best the annoyance? With correction for loudness?
% ANOVA with all 6 SQMs and all 4 PA models with response data being "annoyance" --> only have 6 samples unfortunately, not sufficient?

%% Pivot table for questionnaire Bieke
Q_Bieke = {'flyover', % Do you hear aircraft passing over your house on a daily basis due to close proximity to an airport?
          'airport_type', % What best describes the size of the airport?
          'electric', % Do you think it is likely that you can board an electric aircraft for your holiday trip in 2035?
          'aviation_attitude'}; % What is your general attitude towards aviation?

% check for biases
for i = 1:length(Q_Bieke)
    Q_Bieke{i,2} = pivot(responses, Columns = Q_Bieke{i});
end

delta_annoyance_participant = Annoyance_Virus_all - Annoyance_Velis_all;
delta_normalised_participant = delta_annoyance_participant./(mean(responses_BvdH,2));

%% Check bias - based on questionnaire
% if people answered that they live close to an airport, do they score it differently? 
% if people think electric aviation is unlikely, do they score it differently?

figure
boxplot(delta_normalised_participant{:,:},responses.flyover,'GroupOrder',["Yes","No"])
title('Flyover everyday')
ylabel('Normalised difference in annoyance (RA_{Virus} - RA_{Velis})')
hold on
plot(1,table2array(mean(delta_normalised_participant((responses.flyover=="Yes"),1))) ...
    ,'Marker','d','MarkerEdgeColor','k','MarkerFaceColor','k')
plot(2,table2array(mean(delta_normalised_participant((responses.flyover=="No"),1))) ...
    ,'Marker','d','MarkerEdgeColor','k','MarkerFaceColor','k')

figure
boxplot(delta_normalised_participant{:,:},responses.airport_type,"GroupOrder",["Major (Rotterdam, Schiphol, Eindhoven, etc.)","Small regional / leisure","n/a"])
title('Airport type')
ylabel('Normalised difference in annoyance (RA_{Virus} - RA_{Velis})')
hold on
plot(1,table2array(mean(delta_normalised_participant((responses.airport_type=="Major (Rotterdam, Schiphol, Eindhoven, etc.)"),1))) ...
    ,'Marker','d','MarkerEdgeColor','k','MarkerFaceColor','k')
plot(2,table2array(mean(delta_normalised_participant((responses.airport_type=="Small regional / leisure"),1))) ...
    ,'Marker','d','MarkerEdgeColor','k','MarkerFaceColor','k')
plot(3,table2array(mean(delta_normalised_participant((responses.airport_type=="n/a"),1))) ...
    ,'Marker','d','MarkerEdgeColor','k','MarkerFaceColor','k')

figure
boxplot(delta_normalised_participant{:,:},responses.electric,"GroupOrder",["Very unlikely","Not likely","Not sure", "Likely, but later", "Very likely"])
title('Electric likely')
ylabel('Normalised difference in annoyance (RA_{Virus} - RA_{Velis})')
hold on
plot(1,table2array(mean(delta_normalised_participant((responses.electric=="Very unlikely"),1))) ...
    ,'Marker','d','MarkerEdgeColor','k','MarkerFaceColor','k')
plot(2,table2array(mean(delta_normalised_participant((responses.electric=="Not likely"),1))) ...
    ,'Marker','d','MarkerEdgeColor','k','MarkerFaceColor','k')
plot(3,table2array(mean(delta_normalised_participant((responses.electric=="Not sure"),1))) ...
    ,'Marker','d','MarkerEdgeColor','k','MarkerFaceColor','k')
plot(4,table2array(mean(delta_normalised_participant((responses.electric=="Likely, but later"),1))) ...
    ,'Marker','d','MarkerEdgeColor','k','MarkerFaceColor','k')
plot(5,table2array(mean(delta_normalised_participant((responses.electric=="Very likely"),1))) ...
    ,'Marker','d','MarkerEdgeColor','k','MarkerFaceColor','k')

figure
boxplot(delta_normalised_participant{:,:},responses.aviation_attitude,"GroupOrder",["Negative", "Neutral", "Positive","Very positive"])
title('Attitude to aviation ')
ylabel('Normalised difference in annoyance (RA_{Virus} - RA_{Velis})')
hold on
plot(1,table2array(mean(delta_normalised_participant((responses.aviation_attitude=="Negative"),1))) ...
    ,'Marker','d','MarkerEdgeColor','k','MarkerFaceColor','k')
plot(2,table2array(mean(delta_normalised_participant((responses.aviation_attitude=="Neutral"),1))) ...
    ,'Marker','d','MarkerEdgeColor','k','MarkerFaceColor','k')
plot(3,table2array(mean(delta_normalised_participant((responses.aviation_attitude=="Positive"),1))) ...
    ,'Marker','d','MarkerEdgeColor','k','MarkerFaceColor','k')
plot(4,table2array(mean(delta_normalised_participant((responses.aviation_attitude=="Very positive"),1))) ...
    ,'Marker','d','MarkerEdgeColor','k','MarkerFaceColor','k')



%% How much do participants weigh it differently?
idxVirus = contains(responses.Properties.VariableNames, 'Virus');
idxVelis = contains(responses.Properties.VariableNames, 'Velis');

for i = 1:57
    difference_participant = responses{i,idxVirus}- responses{i,idxVelis}.'; % rows: Velis, columns: Virus
    max_difference(i) = max(max(difference_participant));
    min_difference(i) = min(min(difference_participant));
end