123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222 |
- clear all
- close all
- rootpath='vpl_behav';
- bats = {'vpl_b2b3', 'vpl_b2b5', 'vpl_b2b8', 'vpl_b2mh', 'vpl_b3mt', ...
- 'vpl_b4b3', 'vpl_b4b5', 'vpl_b4b8', 'vpl_b4mb', 'vpl_b5mt', ...
- 'vpl_b7b3', 'vpl_b7b5', 'vpl_b7b8', 'vpl_b7my', 'vpl_b8mp'};
- % fixed variables
- fs=192e3;
- env_mad=0.005; % duration of the boxcar kernel(for moving average) to smooth envelope (s)
- env_flow=500;
- pk_thrs=-45; % dB
- pk_dist=0.02; % in seconds
- dur_thrs=-20; % dB
- freq_thrs=-20; % dB
- [hpb,hpa] = butter(8,2*50e3/fs,'high'); % 10 kHz high pass
- blockdur=60;
- mph=0.002;
- mpd=0.01*fs;
- nfft=128;
- for batnum=1:length(bats),
- filepath=[rootpath 'recordings\' bats{batnum} '\'];
- d=dir(filepath);
- wavnames=strvcat(d.name);
- for wavnum=3:size(wavnames,1);
- wavname=deblank(wavnames(wavnum,:));
- if 1-isempty(findstr(wavname,'m1.wav'));
- wavname
- cn=0;
- result=[];
- savname=[rootpath 'results\' bats{batnum} '\' wavname(1:end-4) '.mat'];
- fid=fopen(savname);
- if fid<0, % this file has not been analysed
- filesize=wavread([filepath wavname],'size');
- filename1=wavname(1:end-5);
- nblocks=ceil(filesize(1)/(blockdur*fs));
-
- for bn=1:nblocks,
- disp(['working on block ' num2str(bn) ' of ' num2str(nblocks)]);
- blockbuf=zeros(round(blockdur*fs),1);
- block_start=(bn-1)*round(blockdur*fs)+1;
-
- if bn < nblocks
- block_stop=bn*blockdur*fs;
- else
- block_stop=filesize(1);
- end
-
- blockbuf=zeros(block_stop-block_start+1,6);
-
- for micnum=1:6,
- micfile=[filename1 num2str(micnum) '.wav'];
- blockbuf(:,micnum)=wavread([filepath micfile],[block_start block_stop]);
- end
- compblockbuf=filtfilt(hpb,hpa,blockbuf);
-
- % call finder
- sumblockbuf=sum(abs(compblockbuf),2);
-
- % smooth the envelope
- [eb,ea]=butter(2,2*env_flow/fs);
- env_dB=20*log10(max(filtfilt(eb,ea,sumblockbuf),1e-4));
-
- %% find calls
- t=[0:length(env_dB)-1]/fs;
-
-
- % figure(2),
- % plot(env_dB)
- % hold on
- % xs=get(gca,'xlim');
- % ys=[pk_thrs pk_thrs];
- % lh=line(xs,ys);
- % set(lh,'color','r')
-
- [pk_val pk_pos]=findpeaks(env_dB,'minpeakheight',pk_thrs,'minpeakdistance',round(pk_dist*fs));
-
- if 1-isempty(pk_pos)
- % draw found peaks positions
- % hold on
- % plot(pk_pos,pk_val,'r*')
- % hold off
-
- % call cutting and saving
- for call_num=1:size(pk_pos)
- cn=cn+1;
-
- % 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
- start=max(find(env_dB(1:pk_pos(call_num))<dur_thrs));
- stop=min(find(env_dB(pk_pos(call_num):end)<dur_thrs))+pk_pos(call_num);
-
- if 1-isempty(stop) & 1-isempty(start)
-
- call_starts=(bn-1)*blockdur*fs+start;
- call_dur=(stop-start)/fs;
-
- result.call(cn).call_start_sample=call_starts;
- result.call(cn).call_dur=call_dur;
-
- call=blockbuf(start:stop,1:6);
- call_levels=20*log10(std(call));
- [call_level_max call_level_max_pos]=max(call_levels);
- xcall=call(:,call_level_max_pos);
- % figure(3)
- % subplot(2,1,1)
- % plot(xcall)
- % subplot(2,1,2)
- % spectrogram(xcall,nfft,120,nfft,fs,'yaxis');
-
- result.call(cn).call_levels=call_levels;
- result.call(cn).loudest_call=xcall;
- result_call(cn).loudest_mic=call_level_max_pos;
-
- %% --- call frequency parameter ----
- %% spectrogram method
- spects = pwelch(xcall,256,250,nfft,fs);
- spectsdB = 10*log10(spects);
- [ampl_pf pos_pf]= max(spectsdB);
- r = find(spectsdB >= ampl_pf+freq_thrs);
- faxis = (0:1/nfft:1/2).*fs;
- PF = faxis(pos_pf); % peak frequency
- Fmin = faxis(r(1)); % min frequency
- Fmax= faxis(r(end)); % max frequency
- BW= Fmax - Fmin; % bandwidth
-
- %% --- spectral centroid frequency
- SCF = sum(spects.*faxis')/sum(spects);
- SCF_dB = sum(spectsdB.*faxis')/sum(spectsdB);
-
- result.call(cn).PF=PF;
- result.call(cn).Fmin=Fmin;
- result.call(cn).Fmax=Fmax;
- result.call(cn).BW=BW;
- result.call(cn).SCF=SCF;
- result.call(cn).SCF_dB=SCF_dB;
-
-
- % use yin to calculate f0 as a function of time
- p.sr = fs;
-
- R=yin(xcall,p,0);
- f0s=R.f0;
- nframes=length(f0s);
- f0ts=(1:nframes)/(fs/R.hop);
- aperiodicity=R.ap0;
-
- % correct for f0 jumps
- fcin=[f0ts' f0s' aperiodicity'];
- fdiff=abs(diff(fcin(:,2)));
- finds=find(fdiff>9000);
- if 1-isempty(finds) % corrections are necessary
- f0_corrected=1;
- 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
- % figure(10), plot(fcin(:,1),fcin(:,2),'b',fcin(:,1),f0_corr,'r');
- %corr_ok=input('correction (1) ok or (0) not ok');
- corr_ok=1;
- if corr_ok==1,
- fcin(:,2)=f0_corr;
- end
- else
- f0_corrected=0;
- end
-
- result.call(cn).f0=fcin;
- result.call(cn).f0_corrected=f0_corrected;
-
- % fit f0(t)
- f0inds=find(1-isnan(fcin(:,2)));
- f0data=fcin(f0inds,1:2);
- [f0slope,sfm_amp,sfm_frq,sfm_phase]=fit_f0(f0data,0);
-
- result.call(cn).f0_slope=f0slope;
- result.call(cn).sfm_amp=sfm_amp;
- result.call(cn).sfm_frq=sfm_frq;
- result.call(cn).sfm_phase=sfm_phase;
- end
- end
- % pause;
- end
- drawnow;
- end
- save(savname,'result');
- else
- fclose(fid);
- end
- end
- end
- end
|