%% 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