VocalProductionLearning_P.discolor/scripts/Script2_Analysis.m

clear all;
close all;

% fixed variables
fs = 192000;   % sampling rate in Hz
gain_comp=156;  % mic gain compensation
dur_min=0.005;  % minimum call duration (excluding echolocation calls) in sec
factor_freq=1/1e3;      % converting frequency values into kHz
factor_dur=1000;        % converting duration values into ms

% load matrix with call values
rootpath=['D:\Documents\Manuskripte\PitchShift_artificial_template\scripts\'];
fn=uigetfile([rootpath '*.mat']);
load([rootpath fn]);

% extract batname and start/stop date and day from file name
batnum=char(fn(10));

%% define start and stop dates for the six bats
if batnum=='2',
    batname='DAM';  % Bat ID
    startdate='17/03/08';
    stopdate='17/08/10';
    pre={'17/03/08';'17/03/09';'17/03/10';'17/03/13';'17/03/14'};
    post30={'17/05/09';'17/05/10';'17/05/11';'17/05/12';'17/05/19'};
    post60={'17/07/05';'17/07/13';'17/07/14';'17/07/17';'17/07/18'};
    nofi={'17/07/27';'17/07/28';'17/07/31';'17/08/01';'17/08/02'};
    noshi={'17/08/04';'17/08/07';'17/08/08';'17/08/09';'17/08/10'};
elseif batnum=='4',
    batname='DID';
    startdate='16/12/13';
    stopdate='17/07/28';
    pre={'16/12/13';'16/12/15';'16/12/19';'16/12/20';'16/12/23'};
    post30={'17/02/21';'17/02/22';'17/02/28';'17/03/01';'17/03/06'};
    post60={'17/04/27';'17/04/28';'17/05/08';'17/05/09';'17/05/10'};
    nofi={'17/07/13';'17/07/14';'17/07/17';'17/07/18';'17/07/19'};
    noshi={'17/07/21';'17/07/24';'17/07/25';'17/07/27';'17/07/28'};
elseif batnum=='5',
    batname='MAR';
    startdate='16/12/15';
    stopdate='17/07/28';
    pre={'16/12/15';'16/12/16';'16/12/19';'16/12/20';'16/12/23'};
    post30={'17/02/20';'17/02/21';'17/02/22';'17/02/28';'17/03/01'};
    post60={'17/04/24';'17/04/25';'17/04/26';'17/04/28';'17/05/08'};
    nofi={'17/07/13';'17/07/14';'17/07/17';'17/07/18';'17/07/19'};
    noshi={'17/07/21';'17/07/24';'17/07/25';'17/07/27';'17/07/28'};
elseif batnum=='6',
    batname='PAU';
    startdate='17/02/07';
    stopdate='17/08/10';
    dur_min=0.025;          % excluding all of PAU short calls!!!
    pre={'17/02/07';'17/02/09';'17/02/10';'17/02/14';'17/02/15'};
    post30={'17/04/06';'17/04/07';'17/04/10';'17/04/11';'17/04/12'};
    post60={'17/06/23';'17/06/26';'17/06/27';'17/06/28';'17/06/29'};
    nofi={'17/07/27';'17/07/31';'17/08/01';'17/08/02';'17/08/03'};
    noshi={'17/08/04';'17/08/07';'17/08/08';'17/08/09';'17/08/10'};
elseif batnum=='1',
    batname='ARN';  % Bat ID
    startdate='17/03/10';
    stopdate='17/07/19';
    pre={'17/03/10';'17/03/13';'17/03/14';'17/03/15';'17/03/16';'17/03/17'};
    post30={'17/05/10';'17/05/11';'17/05/12';'17/05/18';'17/05/19';'17/05/22'};
    post60={'17/07/05';'17/07/13';'17/07/14';'17/07/17';'17/07/18';'17/07/19'};
elseif batnum=='3',
    batname='DAV';  % Bat ID
    startdate='17/03/21';
    stopdate='17/08/01';
    pre={'17/03/21';'17/03/22';'17/03/23';'17/03/27';'17/03/28'};
    post30={'17/05/30';'17/05/31';'17/06/01';'17/06/02';'17/06/06'};
    post60={'17/07/25';'17/07/27';'17/07/28';'17/07/31';'17/08/01'};
end


% absolut number of calls found in the recordings
ncalls=size(result.call);
ncalls=ncalls(2);

% calculate the mean fundamental frequency for all calls longer than min.
% call duration
mf0s=zeros(ncalls,1);       % create zeros-vector with length=ncalls
for callnum=1:ncalls,       % go through all calls
    f0m=result.call(callnum).f0;
    
    % extract call day and time info for the calls
    callday(callnum,:)=result.call(callnum).day;
    calltime(callnum,:)=result.call(callnum).time;
    
    % use only calls which are longer than minimum call duration
    if result.call(callnum).dur>=dur_min,
        mf0s(callnum)=nanmean(f0m(:,2));
    else
        mf0s(callnum)=nan;
    end
end

% create pre and post data vectors
pre_level=[];
pre_dur=[];
pre_f0=[];
post30_level=[];
post30_dur=[];
post30_f0=[];
post60_level=[];
post60_dur=[];
post60_f0=[];
nofi_level=[];
nofi_dur=[];
nofi_f0=[];
noshi_level=[];
noshi_dur=[];
noshi_f0=[];

% extract variables from the structure
for cn=1:ncalls,
    current_date=result.call(cn).day;
    if strcmp(current_date,pre{1}) | strcmp(current_date,pre{2}) | strcmp(current_date,pre{3}) | strcmp(current_date,pre{4}) | strcmp(current_date,pre{5}),
        pre_level=[pre_level result.call(cn).level];
        pre_dur=[pre_dur result.call(cn).dur];
        mf0=result.call(cn).f0;
        mf0=nanmean(mf0(:,2));
        pre_f0=[pre_f0 mf0];
    elseif strcmp(current_date,post30{1}) | strcmp(current_date,post30{2}) | strcmp(current_date,post30{3}) | strcmp(current_date,post30{4}) | strcmp(current_date,post30{5}),
        post30_level=[post30_level result.call(cn).level];
        post30_dur=[post30_dur result.call(cn).dur];
        mf0=result.call(cn).f0; % extract matrix
        mf0=nanmean(mf0(:,2));  % generate mean of second column
        post30_f0=[post30_f0 mf0];  % save
    elseif strcmp(current_date,post60{1}) | strcmp(current_date,post60{2}) | strcmp(current_date,post60{3}) | strcmp(current_date,post60{4}) | strcmp(current_date,post60{5}),
        post60_level=[post60_level result.call(cn).level];
        post60_dur=[post60_dur result.call(cn).dur];
        mf0=result.call(cn).f0; % extract matrix
        mf0=nanmean(mf0(:,2));  % generate mean of second column
        post60_f0=[post60_f0 mf0];  % save
    end
end

if batnum == '1' | batnum=='3',
    short=1;
else
    short=2;
    for cn=1:ncalls,
        current_date=result.call(cn).day;
        if strcmp(current_date,nofi{1}) | strcmp(current_date,nofi{2}) | strcmp(current_date,nofi{3}) | strcmp(current_date,nofi{4}) | strcmp(current_date,nofi{5}),
            nofi_level=[nofi_level result.call(cn).level];
            nofi_dur=[nofi_dur result.call(cn).dur];
            mf0=result.call(cn).f0; % extract matrix
            mf0=nanmean(mf0(:,2));  % generate mean of second column
            nofi_f0=[nofi_f0 mf0];  % save
        elseif strcmp(current_date,noshi{1}) | strcmp(current_date,noshi{2}) | strcmp(current_date,noshi{3}) | strcmp(current_date,noshi{4}) | strcmp(current_date,noshi{5}),
            noshi_level=[noshi_level result.call(cn).level];
            noshi_dur=[noshi_dur result.call(cn).dur];
            mf0=result.call(cn).f0; % extract matrix
            mf0=nanmean(mf0(:,2));  % generate mean of second column
            noshi_f0=[noshi_f0 mf0];  % save
        end
    end
end


% exclude calls shorter than call_dur and convert parameters in correct
% units (incl. correcting for different microphone gains)
pre_duration_indices=find(pre_dur>=dur_min);
post30_duration_indices=find(post30_dur>=dur_min);
post60_duration_indices=find(post60_dur>=dur_min);
% pre days
pre_level=pre_level(pre_duration_indices)'+gain_comp;
pre_dur=pre_dur(pre_duration_indices)'.*factor_dur;
pre_f0=pre_f0(pre_duration_indices)'.*factor_freq;
% post 30 days
post30_level=post30_level(post30_duration_indices)'+gain_comp;
post30_dur=post30_dur(post30_duration_indices)'.*factor_dur;
post30_f0=post30_f0(post30_duration_indices)'.*factor_freq;
% post 60 days
post60_level=post60_level(post60_duration_indices)'+gain_comp;
post60_dur=post60_dur(post60_duration_indices)'.*factor_dur;
post60_f0=post60_f0(post60_duration_indices)'.*factor_freq;
if short==2,
    nofi_duration_indices=find(nofi_dur>=dur_min);
    noshi_duration_indices=find(noshi_dur>=dur_min);
    % nofi days
    nofi_level=nofi_level(nofi_duration_indices)'+gain_comp;
    nofi_dur=nofi_dur(nofi_duration_indices)'.*factor_dur;
    nofi_f0=nofi_f0(nofi_duration_indices)'.*factor_freq;
    % noshi days
    noshi_level=noshi_level(noshi_duration_indices)'+gain_comp;
    noshi_dur=noshi_dur(noshi_duration_indices)'.*factor_dur;
    noshi_f0=noshi_f0(noshi_duration_indices)'.*factor_freq;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% This part of the script PLOTS extracted data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% create vectors for boxplots
x_dur=[pre_dur' post30_dur' post60_dur' nofi_dur' noshi_dur'];
y_dur=[ones(size(pre_dur)); 2*ones(size(post30_dur)); 3*ones(size(post60_dur));4*ones(size(nofi_dur));5*ones(size(noshi_dur))];
x_level=[pre_level' post30_level' post60_level' nofi_level' noshi_level'];
y_level=[ones(size(pre_level)); 2*ones(size(post30_level)); 3*ones(size(post60_level));4*ones(size(nofi_level));5*ones(size(noshi_level))];
x_f0=[pre_f0; post30_f0; post60_f0; nofi_f0; noshi_f0];
y_f0=[ones(size(pre_f0)); 2*ones(size(post30_f0)); 3*ones(size(post60_f0));4*ones(size(nofi_f0));5*ones(size(noshi_f0))];

% create Figure with subplots for the different parameters
figure('Name',['Parameter summary boxplots of bat ' batnum],'NumberTitle','off','units','normalized','outerposition',[0 0 1 1]),
% subplot for duration
subplot(222),
boxplot(x_dur, y_dur,'notch','on',...
    'symbol','.','color','k','jitter',0.2,'width',0.75);
set(gca,'xtick',1:length(y_dur),'xticklabel',1:length(y_dur),...
    'XTickLabel',{'pre_filter' 'filter_30' 'filter_60' 'no_filter' 'no_ps'},'fontsize',14,'YGrid','on');
ylabel('duration (ms)','fontsize',16);
ylim([0 95]);
box off
h = findobj(gca,'Tag','Box');
for j=1:length(h)
    patch(get(h(1,1),'XData'),get(h(1,1),'YData'),'r','FaceAlpha',.7);
    patch(get(h(2,1),'XData'),get(h(2,1),'YData'),'g','FaceAlpha',.7);
    patch(get(h(3,1),'XData'),get(h(3,1),'YData'),'c','FaceAlpha',.7);
    if short==2,
        patch(get(h(4,1),'XData'),get(h(4,1),'YData'),'b','FaceAlpha',.7);
        patch(get(h(5,1),'XData'),get(h(5,1),'YData'),'r','FaceAlpha',.7);
    end
end
title(['Bat ' batnum ' - Call duration']);

% subplot for call level
subplot(223),
boxplot(x_level, y_level,'notch','on',...
    'symbol','.','color','k','jitter',0.2,'width',0.75);
set(gca,'xtick',1:length(y_level),'xticklabel',1:length(y_level),...
    'XTickLabel',{'pre_filter' 'filter_30' 'filter_60' 'no_filter' 'no_ps'},'fontsize',14,'YGrid','on');
ylabel('level (dB SPL)','fontsize',16);
ylim([90 130]);
box off
h = findobj(gca,'Tag','Box');
for j=1:length(h)
    patch(get(h(1,1),'XData'),get(h(1,1),'YData'),'r','FaceAlpha',.7);
    patch(get(h(2,1),'XData'),get(h(2,1),'YData'),'g','FaceAlpha',.7);
    patch(get(h(3,1),'XData'),get(h(3,1),'YData'),'c','FaceAlpha',.7);
    if short==2,
        patch(get(h(4,1),'XData'),get(h(4,1),'YData'),'b','FaceAlpha',.7);
        patch(get(h(5,1),'XData'),get(h(5,1),'YData'),'r','FaceAlpha',.7);
    end
end
title(['Bat ' batnum ' - Call level']);

% subplot for mean fundamental frequency
subplot(221),
boxplot(x_f0, y_f0,'notch','on',...
    'symbol','.','color','k','jitter',0.2,'width',0.75);
set(gca,'xtick',1:length(y_f0),'xticklabel',1:length(y_f0),...
    'XTickLabel',{'pre_filter' 'filter_30' 'filter_60' 'no_filter' 'no_ps'},'fontsize',14,'YGrid','on');
ylabel('Mean f0 (kHz)','fontsize',16);
ylim([12 22]);
box off
h = findobj(gca,'Tag','Box');
for j=1:length(h)
    patch(get(h(1,1),'XData'),get(h(1,1),'YData'),'r','FaceAlpha',.7);
    patch(get(h(2,1),'XData'),get(h(2,1),'YData'),'g','FaceAlpha',.7);
    patch(get(h(3,1),'XData'),get(h(3,1),'YData'),'c','FaceAlpha',.7);
    if short==2,
        patch(get(h(4,1),'XData'),get(h(4,1),'YData'),'b','FaceAlpha',.7);
        patch(get(h(5,1),'XData'),get(h(5,1),'YData'),'r','FaceAlpha',.7);
    end
end
title(['Bat ' batnum ' - Mean fundamental frequency (f0)']);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% This part of the script runs the statistics
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% stats for the duration (dur)
PRE_dur=[kstest(pre_dur),size(pre_dur),median(pre_dur),iqr(pre_dur)];
POST30_dur=[kstest(post30_dur),size(post30_dur),median(post30_dur),iqr(post30_dur)];
POST60_dur=[kstest(post60_dur),size(post60_dur),median(post60_dur),iqr(post60_dur)];
p30_dur=ranksum(pre_dur,post30_dur);
p60_dur=ranksum(pre_dur,post60_dur);
p3060_dur=ranksum(post30_dur,post60_dur);
if short ==2,
    NOFI_dur=[kstest(nofi_dur),size(nofi_dur),median(nofi_dur),iqr(nofi_dur)];
    NOSHI_dur=[kstest(noshi_dur),size(noshi_dur),median(noshi_dur),iqr(noshi_dur)];
    p60nofi_dur=ranksum(post60_dur,nofi_dur);
    pnofishi_dur=ranksum(nofi_dur,noshi_dur);
end

% stats for call level
PRE_level=[kstest(pre_level),size(pre_level),median(pre_level),iqr(pre_level)];
POST30_level=[kstest(post30_level),size(post30_level),median(post30_level),iqr(post30_level)];
POST60_level=[kstest(post60_level),size(post60_level),median(post60_level),iqr(post60_level)];
p30_level=ranksum(pre_level,post30_level);
p60_level=ranksum(pre_level,post60_level);
p3060_level=ranksum(post30_level,post60_level);
if short==2,
    NOFI_level=[kstest(nofi_level),size(nofi_level),median(nofi_level),iqr(nofi_level)];
    NOSHI_level=[kstest(noshi_level),size(noshi_level),median(noshi_level),iqr(noshi_level)];
    p60nofi_level=ranksum(post60_level,nofi_level);
    pnofishi_level=ranksum(nofi_level,noshi_level);
end

% stats for mean fundamental frequency (f0)
PRE_f0=[kstest(pre_f0),size(pre_f0),median(pre_f0),iqr(pre_f0)];
POST60_f0=[kstest(post60_f0),size(post60_f0),median(post60_f0),iqr(post60_f0)];
POST30_f0=[kstest(post30_f0),size(post30_f0),median(post30_f0),iqr(post30_f0)];
p30_f0=ranksum(pre_f0,post30_f0);
p60_f0=ranksum(pre_f0,post60_f0);
p3060_f0=ranksum(post30_f0,post60_f0);
if short==2,
    NOFI_f0=[kstest(nofi_f0),size(nofi_f0),median(nofi_f0),iqr(nofi_f0)];
    NOSHI_f0=[kstest(noshi_f0),size(noshi_f0),median(noshi_f0),iqr(noshi_f0)];
    p60nofi_f0=ranksum(post60_f0,nofi_f0);
    pnofishi_f0=ranksum(nofi_f0,noshi_f0);
end

%print tables with stats for comparison pre and post 30 days data
ps1=[p30_dur;p30_dur;p30_level;p30_level;p30_f0;p30_f0];
table1=[PRE_dur;POST30_dur;PRE_level;POST30_level;PRE_f0;POST30_f0];
table1(:,3)=table1(:,4);
table1(:,4)=table1(:,5);
table1(:,5)=ps1;
printmat(table1, 'Summary Stats pre and post30 data', 'pre_dur post30_dur pre_level post30_level pre_f0 post30_f0',...
    'ks-test calls median IQR p-value')

%print tables with stats for comparison pre and post 60 days data
ps2=[p60_dur;p60_dur;p60_level;p60_level;p60_f0;p60_f0];
table2=[PRE_dur;POST60_dur;PRE_level;POST60_level;PRE_f0;POST60_f0];
table2(:,3)=table2(:,4);
table2(:,4)=table2(:,5);
table2(:,5)=ps2;
printmat(table2, 'Summary Stats pre and post60 data', 'pre_dur post60_dur pre_level post60_level pre_f0 post60_f0',...
    'ks-test calls median IQR p-value')

% print tables with stats for comparison post 30 and post 60 days data
ps3=[p3060_dur;p3060_dur;p3060_level;p3060_level;p3060_f0;p3060_f0];
table3=[POST30_dur;POST60_dur;POST30_level;POST60_level;POST30_f0;POST60_f0];
table3(:,3)=table3(:,4);
table3(:,4)=table3(:,5);
table3(:,5)=ps3;
printmat(table3, 'Summary Stats post30 and post60 data', 'post30_dur post60_dur post30_level post60_level post30_f0 post60_f0',...
    'ks-test calls median IQR p-value')

if short==2,
    % print tables with stats for comparison post 60 days and no_filter data
    ps4=[p60nofi_dur;p60nofi_dur;p60nofi_level;p60nofi_level;p60nofi_f0;p60nofi_f0];
    table4=[POST60_dur;NOFI_dur;POST60_level;NOFI_level;POST60_f0;NOFI_f0];
    table4(:,3)=table4(:,4);
    table4(:,4)=table4(:,5);
    table4(:,5)=ps4;
    printmat(table4, 'Summary Stats post60 and no_filter data', 'post60_dur nofi_dur post60_level nofi_level post60_f0 nofi_f0',...
        'ks-test calls median IQR p-value')
    
    % print tables with stats for comparison no_filter and no_shift data
    ps6=[pnofishi_dur;pnofishi_dur;pnofishi_level;pnofishi_level;pnofishi_f0;pnofishi_f0];
    table6=[NOFI_dur;NOSHI_dur;NOFI_level;NOSHI_level;NOFI_f0;NOSHI_f0];
    table6(:,3)=table6(:,4);
    table6(:,4)=table6(:,5);
    table6(:,5)=ps6;
    printmat(table6, 'Summary Stats no_filter and no_shift data', 'nofi_dur noshi_dur nofi_level noshi_level nofi_f0 noshi_f0',...
        'ks-test calls median IQR p-value')
end

%%%%%%%%%%%%%%% mean and standard deviation script

%pre days
pre_level_mean=mean(pre_level);
pre_dur_mean=mean(pre_dur);
pre_f0_mean=mean(pre_f0);
pre_level_std=std(pre_level);
pre_dur_std=std(pre_dur);
pre_f0_std=std(pre_f0);

%post 30 days
post30_level_mean=mean(post30_level);
post30_dur_mean=mean(post30_dur);
post30_f0_mean=mean(post30_f0);
post30_level_std=std(post30_level);
post30_dur_std=std(post30_dur);
post30_f0_std=std(post30_f0);

%post 60 days
post60_level_mean=mean(post60_level);
post60_dur_mean=mean(post60_dur);
post60_f0_mean=mean(post60_f0);
post60_level_std=std(post60_level);
post60_dur_std=std(post60_dur);
post60_f0_std=std(post60_f0);

if short==2,
    %no_filter days
    nofi_level_mean=mean(nofi_level);
    nofi_dur_mean=mean(nofi_dur);
    nofi_f0_mean=mean(nofi_f0);
    nofi_level_std=std(nofi_level);
    nofi_dur_std=std(nofi_dur);
    nofi_f0_std=std(nofi_f0);
    
    %no_pitch-shift days
    noshi_level_mean=mean(noshi_level);
    noshi_dur_mean=mean(noshi_dur);
    noshi_f0_mean=mean(noshi_f0);
    noshi_level_std=std(noshi_level);
    noshi_dur_std=std(noshi_dur);
    noshi_f0_std=std(noshi_f0);
    
    % prints table with mean and std values
    table5=[pre_f0_mean;post30_f0_mean;post60_f0_mean;nofi_f0_mean;noshi_f0_mean;pre_f0_std;post30_f0_std;post60_f0_std;nofi_f0_std;noshi_f0_std];
    table5(1:5,2)=table5(6:10,1);
    table5(6:10,1)=0;
    printmat(table5, 'Summary stats for fundamental freq.', 'pre_f0 post30_f0 post60_f0 no_filter_f0 no_pitch_shift - - - - -',...
        'mean std')
    
    %plots graphic with mean and std values
    y = [pre_f0_mean post30_f0_mean post60_f0_mean nofi_f0_mean noshi_f0_mean];
    x = [1 2 3 4 5];
    figure,
    err = [pre_f0_std post30_f0_std post60_f0_std nofi_f0_std noshi_f0_std];
    errorbar(x,y,err,'-s','LineWidth',3,'MarkerSize',3,...
        'MarkerEdgeColor','k','MarkerFaceColor','k')
    set(gca,'XTickLabel',{'','pre_filter','','filter30','','filter60','','no_filter','','no_pitch-shift',''});
    grid on
    title(['Bat ' num2str(batnum) ' - change of mean fundamental frequency']);
    ylabel('mean f0 (kHz)')
end

if short ==1,
    % prints table with mean and std values
    table5=[pre_f0_mean;post30_f0_mean;post60_f0_mean;pre_f0_std;post30_f0_std;post60_f0_std];
    table5(1:3,2)=table5(4:6,1);
    table5(4:6,1)=0;
    printmat(table5, 'Summary stats for fundamental freq.', 'pre_f0 post30_f0 post60_f0 - - -',...
        'mean std')
    
    %plots graphic with mean and std values
    y = [pre_f0_mean post30_f0_mean post60_f0_mean];
    x = [1 2 3];
    figure,
    err = [pre_f0_std post30_f0_std post60_f0_std];
    errorbar(x,y,err,'-s','LineWidth',3,'MarkerSize',3,...
        'MarkerEdgeColor','k','MarkerFaceColor','k')
    set(gca,'XTickLabel',{'pre_filter','','','','','filter30','','','','','filter60'});
    grid on
    title([ 'Bat ' num2str(batnum) ' - change of mean fundamental frequency']);
    ylabel('mean f0 (kHz)')
end