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