Dateien hochladen nach 'scripts'

scripts/Script1_Extraction.m

@@ -0,0 +1,298 @@
+%% This script extracts calls from recordings for a specific bat in a
+%% specific timeframe. Parameters of all calls are extracted and saved in a
+%% structure.
+
+%set up/clear the workspace
+clear all;
+close all;
+
+plot_on=0;  %select if results should be plotted
+save_on=1;  %select if results should be saved
+
+%set rootpath of the raw recordings
+rootpath='D:\Documents\Manuskripte\PitchShift_artificial_template\recordings';
+
+%% fixed variables
+fs=192000;      %sampling rate (Hz)
+fhigh=1000;     %highpass filter frequency (Hz)
+nfft=256;       %FFT samples
+
+env_mad=0.005;   %duration of the boxcar kernel(for moving average) to smooth envelope (s)
+env_flow=500;
+
+pk_thrs=-60;    %peak threshold (dB)
+pk_dist=0.06;   %peak-to-peak distance (seconds)
+dur_thrs=-20;   %call duration threshold (-20 dB from peak)
+freq_thrs=-20;  %dB
+
+%select bat and timeframe to be analysed
+batnum = input('Enter bat number to be analysed! (i.e. 1-6) ','s');
+if batnum=='2',
+    batname='DAM';  % Bat ID
+    startdate='17/03/08';
+    stopdate='17/08/10';
+elseif batnum=='4',
+    batname='DID';
+    startdate='16/12/13';
+    stopdate='17/07/28';
+elseif batnum=='5',
+    batname='MAR';
+    startdate='16/12/15';
+    stopdate='17/07/28';
+elseif batnum=='6',
+    batname='PAU';
+    startdate='17/02/07';
+    stopdate='17/08/10';
+elseif batnum=='1',
+    batname='ARN';  % Bat ID
+    startdate='17/03/10';
+    stopdate='17/07/19';
+elseif batnum=='3',
+    batname='DAV';  % Bat ID
+    startdate='17/03/21';
+    stopdate='17/08/01';
+end
+
+startday=datenum(startdate);    %reformat start date
+stopday=datenum(stopdate);      %reformat stop date
+current_call=0;                 %call number counter
+
+for day=startday:stopday, %repeat within selected time frame
+    folder=[rootpath '\' strrep(datestr(day,25),'/','')]; %access folder with raw recordings
+    if isdir(folder),   
+        d=dir(folder);  %get filenames in folder
+        filenames=strvcat(d.name); 
+        for nn=3:size(filenames,1), %for all files in the selected folder
+            filename=deblank(filenames(nn,:));  
+            if 1-isempty(findstr(filename,batname)),    %if the folder contains files from the selected bat
+                clear('call_samples');
+                abort_analysis=0;                   %do not abort analysis
+                load([folder '\' filename]);        %load file from folder
+                prev_ana=0;                         %set indicator: not previously analysed
+                clear('call_samples');
+                disp(filename)                      %show filename
+                disp(prev_ana)                      %show if previously analysed
+                
+                %% high pass filtering of full recording
+                [b,a]=butter(2,2*fhigh/fs,'high'); %generate filter
+                if isinteger(xrec),         %convert recording if not integer
+                fxrec=double(xrec)/32767; 
+                else
+                    fxrec=xrec;             %do not overwrite variable with recording
+                end
+                fxrec=filtfilt(b,a,fxrec);   %filter the recording
+                
+                %if plotting is selected the next paragraph is active
+                if plot_on,
+                    figure(1),
+                    spectrogram(fxrec,64,60,nfft,fs,'yaxis');              
+                    set(gcf,'Position',[7 758 1317 220]);
+                    title('Spectrogram');
+                end
+                
+                if prev_ana==0,     %if recording has not previously been analysed
+                    %smooth the envelope of the recording
+                    [eb,ea]=butter(2,2*env_flow/fs);
+                    env_dB=20*log10(max(filtfilt(eb,ea,abs(fxrec)),1e-5));    
+                    
+                    %find calls in the recording
+                    t=[0:length(env_dB)-1]/fs;
+                    
+                    if plot_on,
+                        figure(2),
+                        plot(t,env_dB)
+                        set(gcf,'Position',[7   463   1317   220]);
+                        xlabel('Rec time (s)');
+                        ylabel('Amplitude (arb. units)');
+                        hold on
+                        ylim([min(env_dB) max(env_dB)]);
+                        xs=get(gca,'xlim');
+                        ys=[pk_thrs pk_thrs];
+                        lh=line(xs,ys);
+                        set(lh,'color','r')
+                        hold off
+                    end
+                    
+                    %identify peaks and save value and location in recording
+                    [pk_val pk_pos]=findpeaks(env_dB,'minpeakh',pk_thrs,'minpeakdist',round(pk_dist*fs));
+                    
+                    %check if last detected call is still fully recorded, otherwise ignore it
+                    if env_dB(end)>pk_thrs,
+                        if length(pk_pos)>1,
+                            pk_pos=pk_pos(1:end-1);
+                            pk_val=pk_val(1:end-1);
+                        else
+                            pk_pos=[];
+                            pk_val=[];
+                        end
+                    end
+                    
+                    % check if first detected call is fully recorded, otherwise ignore it
+                    if env_dB(1)>pk_thrs,
+                        if length(pk_pos)>1,
+                            pk_pos=pk_pos(2:end);
+                            pk_val=pk_val(2:end);
+                        else
+                            pk_pos=[];
+                            pk_val=[];
+                        end
+                    end
+                end
+                
+                if prev_ana==1,  %if recording has previously been analysed
+                    if isempty(call_samples)==0,
+                        pk_pos=call_samples(:,1);
+                        abort_analysis=0;
+                    else %file was analysed but no calls found
+                        abort_analysis=1;
+                    end
+                end
+                
+                %if you analysis was not marked as 'to be aborted'
+                if abort_analysis==0, 
+                    for call_num=1:length(pk_pos),
+                        if prev_ana==0,
+                            % determine start and stop of each call as -x dB duration
+                            env_dB=env_dB-env_dB(pk_pos(call_num));
+                            % find start of call
+                            startsample=max(find(env_dB(1:pk_pos(call_num))<dur_thrs+(pk_thrs-pk_val(call_num))));
+                            % calculate IPI to previous call
+                            if call_num>1,
+                                ipi=(startsample-stopsample)/fs;
+                            else
+                                ipi=nan;
+                            end
+                            stopsample=min(find(env_dB(pk_pos(call_num):end)<dur_thrs+(pk_thrs-pk_val(call_num))))+pk_pos(call_num);
+                        else
+                            startsample=call_samples(call_num,1);
+                            stopsample=call_samples(call_num,2);
+                        end
+                        if startsample>1000 & stopsample>startsample & stopsample < length(fxrec),
+                            call=fxrec(startsample:stopsample); 
+                            t=0:length(call)-1;
+                            t=t/fs;
+                           
+                            if plot_on,
+                                figure(3)
+                                spectrogram(call,64,60,nfft,fs,'yaxis');
+                                set(gcf,'Position',[415   243   400   220]);
+                                title(['Record Start time = '  sprintf('%4.2f',(startsample+recpos)/fs) ' s']);
+                                xlabel('Time (s)');
+                                ylabel('Frequency (Hz)');
+                            end
+
+                            call_ok=1;
+                            if call_ok,
+                                if prev_ana==0,
+                                    call_samples(call_num,1)=startsample;
+                                    call_samples(call_num,2)=stopsample;
+                                end
+                                
+                                current_call=current_call+1; %add 1 to call counter
+                                
+                                % call levels for the call on all 45 channels
+                                call_level=20*log10(std(call));
+                                
+                                % max dB value and microphone number with the max. level
+                                call_dur=max(t);
+                                
+                                %% --- extract call frequency parameter ----
+                                %% spectrogram method
+                                spects=pwelch(call,256,250,nfft,fs);  %PSD estimate of current call
+                                spectsdB=10*log10(spects);            %PSD estimate in dB
+                                [ampl_pf pos_pf]=max(spectsdB);
+                                r=find(spectsdB >= ampl_pf+freq_thrs);
+                                faxis=(0:1/nfft:1/2).*fs;
+                                
+                                %% --- extract fundamental frequency
+                                % use yin to calculate f0 as a function of time
+                                p.sr=fs;
+                                
+                                R=yin(call,p,0);
+                                f0s=R.f0;
+                                nframes=length(f0s);
+                                f0ts=(1:nframes)/(fs/R.hop);
+                                aperiodicity=R.ap0;
+                                
+                                if plot_on,
+                                    figure(4)
+                                    set(gcf,'position',[ 822    36   438   429])
+                                    subplot(3,1,1), plot(f0ts,f0s);
+                                    title('Fundamental frequency analysis');
+                                    xlabel('Time (s)');
+                                    ylabel('Fundamental frequency (Hz)');
+                                    subplot(3,1,2), plot(f0ts,aperiodicity);
+                                    title('Aperiodicity analysis');
+                                    xlabel('Time (s)');
+                                    ylabel('Gross aperiodicity');
+                                end
+                                
+                                % check if f0 corrections are necessary
+                                % ('harmonic jumps' of yin algorithm)
+                                fcin=[f0ts' f0s' aperiodicity'];
+                                fdiff=abs(diff(fcin(:,2)));
+                                finds=find(fdiff>9000);
+                                if 1-isempty(finds)         % corrections are necessary
+                                    f0_corr=fcin(:,2);
+                                    % check start frequency of the irregularity
+                                    if f0_corr(finds(1))>30e3;
+                                        % correct first jump
+                                        f0_corr(1:finds(1))=f0_corr(1:finds(1))/2;
+                                        if length(finds)>1,
+                                            finds=finds(2:end);
+                                        end
+                                    end
+                                    while length(finds)>1
+                                        xfinds=finds(1:2);
+                                        f0_corr(xfinds(1)+1:xfinds(2))=f0_corr(xfinds(1)+1:xfinds(2))/2;
+                                        if length(finds)>2,
+                                            finds=finds(3:end);
+                                        elseif length(finds)==2,
+                                            finds=[];
+                                            break
+                                        else
+                                            break
+                                        end
+                                    end
+                                    if length(finds)>0,
+                                        f0_corr(finds(1)+1:end)=f0_corr(finds(1)+1:end)/2;
+                                    end
+                                    if plot_on,
+                                        figure(8),
+                                        subplot(3,1,3),
+                                        plot(fcin(:,1),fcin(:,2),'b',fcin(:,1),f0_corr,'r');
+                                    end
+                                    corr_ok=1;
+                                    if corr_ok==1,
+                                        fcin(:,2)=f0_corr;
+                                    end
+                                end
+                                drawnow;
+                                
+                                result.call(current_call).level=call_level;
+                                result.call(current_call).dur=call_dur;
+                                result.call(current_call).f0=fcin;
+                                result.call(current_call).day=datestr(day,25);
+                                result.call(current_call).time=filename(end-10:end-4);
+                            end
+                        end
+                    end
+                    if prev_ana==0 & isempty(pk_pos)==0,
+                        % overwrite raw data file with version including
+                        % sample information of calls
+                        if exist('call_samples','var')==0,
+                            call_samples=[];
+                        end
+                        if save_on==1,
+                            save([folder '\' filename],'xrec','recpos','call_samples');
+                        end
+                    end
+                end
+            end
+        end
+    end
+end
+result_filename=['D:\Documents\Manuskripte\PitchShift_artificial_template\scripts\Calls_bat' num2str(batnum) '_' strrep(startdate,'/','') '_to_' strrep(stopdate,'/','')];
+if save_on==1,
+    save(result_filename,'result');
+end

scripts/Script2_Analysis.m

@@ -0,0 +1,442 @@
+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