clear; close all;

%setup
M=3;
type='fft'; %best result - power spectrum, not smoothed
EX='mean';
x=1; %degree of approx

Recordings = uigetdir; %gets directory
Pataka = dir(fullfile(Recordings,'*.wav')); %gets all wav files in struct
Protocols = uigetdir;
Timestamps = dir(fullfile(Protocols,'*.txt'));

data = readtable('data_update.xls','PreserveVariableNames',true);
% data = readtable('data.xls','PreserveVariableNames',true);
data(:,4)=[]; % attribute 'par' deletion
data(19,:)=[]; %missing recording 
data.Gender=categorical(data.Gender);
data=cat(2,data,table(data.Age_m/12,'VariableNames',{'month_decimal'}));

%%
for k = 1 :length(table2array(data(:,1)))
close all;
RecordingName = Pataka(k).name; 
fullFileName = fullfile(Recordings, RecordingName);
fprintf(1, 'Now reading %s\n', fullFileName);
[signal, fs] = audioread(fullFileName);

ProtocolName = Timestamps(k).name; 
fullFileName2 = fullfile(Protocols, ProtocolName);
timestamp=load(fullFileName2);
if size(signal,2)==2
    signal(:,2)=[];
end

%filtr
fm = [100 4000];
[b,a] = butter(5,2*fm/fs,'bandpass');
signal_filt=filtfilt(b,a,signal);
%downsample
% M=2; % M=3 - fs=14700; M=4 - fs=11025; M=5 - fs=8820;
fs=fs/M;
timestamp=uint32(timestamp*fs);
signal_downsampled=downsample(signal_filt,M);

%Timestamps 
segment=signal_downsampled(timestamp(1,1):timestamp(3,3));
siglen=length(segment);

[Pxx,Ff] = pwelch(segment,512,256,siglen-1,fs);

thr=mean(Pxx) +std(Pxx);
[ymax,idx_peak]=findpeaks(Pxx,'MinPeakDistance', 5, 'MinPeakHeight', thr);
plot(Ff,Pxx); hold on; plot(Ff(idx_peak),ymax,'x'); yline(mean(Pxx)+std(Pxx)) %drawnow,
while length(idx_peak)<3
    [ymax,idx_peak]=findpeaks(Pxx,'MinPeakDistance', 5, 'MinPeakHeight', thr);
    thr=thr*0.95;
end
peaks(k,:)=Ff(idx_peak(1:3)');

end
%%
boys_idxs=find(data.Gender=='M');
girls_idxs=find(data.Gender=='F');

boys=zeros(length(boys_idxs),4); % (:,1) - age, 3 peaks
boys(:,1)=table2array(data(boys_idxs, 'Age_y'))+table2array(data(boys_idxs, 'month_decimal'));
boys(:,2:4)=peaks(boys_idxs,:);

[r,p]=corr([boys(:,1),boys(:,2)],'Type','Pearson');
korb(1,1)=r(1,2);
korb(2,1)=p(1,2);
[r,p]=corr([boys(:,1),boys(:,3)],'Type','Pearson');
korb(1,2)=r(1,2);
korb(2,2)=p(1,2);
[r,p]=corr([boys(:,1),boys(:,4)],'Type','Pearson');
korb(1,3)=r(1,2);
korb(2,3)=p(1,2);

girls=zeros(length(girls_idxs),4); % (:,1) - age, 3 peaks
girls(:,1)=table2array(data(girls_idxs, 'Age_y'))+table2array(data(girls_idxs, 'month_decimal'));
girls(:,2:4)=peaks(girls_idxs,:);

[r,p]=corr([girls(:,1),girls(:,2)],'Type','Pearson');
korg(1,1)=r(1,2);
korg(2,1)=p(1,2);
[r,p]=corr([girls(:,1),girls(:,3)],'Type','Pearson');
korg(1,2)=r(1,2);
korg(2,2)=p(1,2);
[r,p]=corr([girls(:,1),girls(:,4)],'Type','Pearson');
korg(1,3)=r(1,2);
korg(2,3)=p(1,2);
%%
figure
subplot 311;
plot(boys(:,1),boys(:,2),'bx'); hold on; plot(girls(:,1),girls(:,2), 'ro'); 
p_pab=polyfit(boys(:,1),boys(:,2),x); y_pab=polyval(p_pab,boys(:,1));
plot(boys(:,1),y_pab,'b');
p_pag=polyfit(girls(:,1),girls(:,2),x); y_pag=polyval(p_pag,girls(:,1));
plot(girls(:,1),y_pag,'r'); hold off;
title('1. vrchol'); xlabel('Věk [rok]'); ylabel('1. spektrální vrchol'); legend('chlapci','děvčata');

subplot 312;
plot(boys(:,1),boys(:,3),'bx'); hold on; plot(girls(:,1),girls(:,3), 'ro'); 
p_pab=polyfit(boys(:,1),boys(:,3),x); y_pab=polyval(p_pab,boys(:,1));
plot(boys(:,1),y_pab,'b');
p_pag=polyfit(girls(:,1),girls(:,3),x); y_pag=polyval(p_pag,girls(:,1));
plot(girls(:,1),y_pag,'r'); hold off;
title('2. vrchol'); xlabel('Věk [rok]'); ylabel('2. spektrální vrchol');

subplot 313;
plot(boys(:,1),boys(:,4),'bx'); hold on; plot(girls(:,1),girls(:,4), 'ro'); 
p_pab=polyfit(boys(:,1),boys(:,4),x); y_pab=polyval(p_pab,boys(:,1));
plot(boys(:,1),y_pab,'b');
p_pag=polyfit(girls(:,1),girls(:,4),x); y_pag=polyval(p_pag,girls(:,1));
plot(girls(:,1),y_pag,'r'); hold off;
title('3. vrchol'); xlabel('Věk [rok]'); ylabel('3. spektrální vrchol');


G1M=boys(boys(:,1)<8,:);
G2M=boys(boys(:,1)>=8 & boys(:,1)<12,:);
G3M=boys(boys(:,1)>=12,:);
G1F=girls(girls(:,1)<8,:);
G2F=girls(girls(:,1)>=8 & girls(:,1)<12,:);
G3F=girls(girls(:,1)>=12,:);

groupGM=cat(2,repmat({'1'},[1 length(G1M)]),repmat({'2'},[1 length(G2M)]),repmat({'3'},[1 length(G3M)]));
groupGF=cat(2,repmat({'1'},[1 size(G1F,1)]),repmat({'2'},[1 size(G2F,1)]),repmat({'3'},[1 size(G3F,1)]));

names = {'Skupina 1'; 'Skupina 2'; 'Skupina 3'};
p1m=[G1M(:,2); G2M(:,2); G3M(:,2)];
p2m=[G1M(:,3); G2M(:,3); G3M(:,3)];
p3m=[G1M(:,4); G2M(:,4); G3M(:,4)];

figure; 
subplot 321; boxplot(p1m,groupGM); title('1. vrchol - chlapci'); ylabel('[Hz]'); set(gca,'xtick',[1:3],'xticklabel',names)
subplot 323; boxplot(p2m,groupGM); title('2. vrchol - chlapci'); ylabel('[Hz]'); set(gca,'xtick',[1:3],'xticklabel',names)
subplot 325; boxplot(p3m,groupGM); title('3. vrchol - chlapci'); ylabel('[Hz]'); set(gca,'xtick',[1:3],'xticklabel',names)
set(gca,'xtick',[1:3],'xticklabel',names)

p1f=[G1F(:,2); G2F(:,2); G3F(:,2)];
p2f=[G1F(:,3); G2F(:,3); G3F(:,3)];
p3f=[G1F(:,4); G2F(:,4); G3F(:,4)];

subplot 322; boxplot(p1f,groupGF); title('1. vrchol - dívky'); ylabel('[Hz]'); set(gca,'xtick',[1:3],'xticklabel',names)
subplot 324; boxplot(p2f,groupGF); title('2. vrchol - dívky'); ylabel('[Hz]'); set(gca,'xtick',[1:3],'xticklabel',names)
subplot 326; boxplot(p3f,groupGF); title('3. vrchol - dívky'); ylabel('[Hz]'); set(gca,'xtick',[1:3],'xticklabel',names)
set(gca,'xtick',[1:3],'xticklabel',names)
%%
%anova
for i=1:3
    M=cat(2,G1M(:,i+1)',G2M(:,i+1)',G3M(:,i+1)');
    F=cat(2,G1F(:,i+1)',G2F(:,i+1)',G3F(:,i+1)');

    p_longan(1,i)=anova1(M,groupGM,'off');
    p_longan(2,i)=anova1(F,groupGF,'off');
    
end
groupan=array2table(p_longan, 'VariableNames', {'1. vrchol','2. vrchol','3. vrchol'}, 'RowNames', {'chlapci','dívky'})

%%
figure; subplot 211;
p1m=([mean(G1M(:,2)) mean(G2M(:,2)) mean(G3M(:,2))]);
p2m=([mean(G1M(:,3)) mean(G2M(:,3)) mean(G3M(:,3))]);
p3m=([mean(G1M(:,4)) mean(G2M(:,4)) mean(G3M(:,4))]);
plot(p1m,'--.b','MarkerSize',15); hold on; plot(p2m,'--.r','MarkerSize',15); 
plot(p3m,'--.k','MarkerSize',15); title('Spektrální vrcholy - chlapci'); legend('1. vrchol','2. vrchol','3. vrchol'); %legend('P','T','K');
names = {'Skupina 1'; 'Skupina 2'; 'Skupina 3'};
set(gca,'xtick',[1:3],'xticklabel',names)

subplot 212;
p1f=([mean(G1F(:,2)) mean(G2F(:,2)) mean(G3F(:,2))]);
p2f=([mean(G1F(:,3)) mean(G2F(:,3)) mean(G3F(:,3))]);
p3f=([mean(G1F(:,4)) mean(G2F(:,4)) mean(G3F(:,4))]);
plot(p1f,'--.b','MarkerSize',15); hold on; plot(p2f,'--.r','MarkerSize',15); 
plot(p3f,'--.k','MarkerSize',15); title('Spektrální vrcholy - děvčata'); %legend('P','T','K');
names = {'Skupina 1'; 'Skupina 2'; 'Skupina 3'}; 
set(gca,'xtick',[1:3],'xticklabel',names)
load('class.mat')

%%
for i=1:3
    [h(1,i),pp(1,i)]=ttest2(G1M(:,i+1),G1F(:,i+1));
    [h(2,i),pp(2,i)]=ttest2(G2M(:,i+1),G2F(:,i+1));
    [h(3,i),pp(3,i)]=ttest2(G3M(:,i+1),G3F(:,i+1));
end

statB1(1,:)=mean(G1M); statB1(2,:)=var(G1M); statB1(3,:)=std(G1M);
statB2(1,:)=mean(G2M); statB2(2,:)=var(G2M); statB2(3,:)=std(G2M);
statB3(1,:)=mean(G3M); statB3(2,:)=var(G3M); statB3(3,:)=std(G3M);
statG1(1,:)=mean(G1F); statG1(2,:)=var(G1F); statG1(3,:)=std(G1F);
statG2(1,:)=mean(G2F); statG2(2,:)=var(G2F); statG2(3,:)=std(G2F);
statG3(1,:)=mean(G3F); statG3(2,:)=var(G3F); statG3(3,:)=std(G3F);

all=[boys; girls];
% statistical analysis
statboys(1,:)=mean(boys);
statgirls(1,:)=mean(girls);
statall(1,:)=mean(all);

statboys(2,:)=var(boys);
statgirls(2,:)=var(girls);
statall(2,:)=var(all);

statboys(3,:)=std(boys);
statgirls(3,:)=std(girls);
statall(3,:)=std(all);
