123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193 |
- function result = extract_adult_call_parameters(recbuf)
- % fixed variables
- fs=192e3;
- env_flow=500;
- pk_thrs=-45; % dB
- pk_dist=0.005; % in seconds (changed from 0.02 in Meike's extraction script)
- dur_thrs=-20; % dB
- freq_thrs=-20; % dB
- hpf_display=50e3;
- nfft = 128;
- [eb,ea]=butter(2,2*env_flow/fs);
- [hpb,hpa] = butter(8,2*hpf_display/fs,'high'); % 50 (!) kHz high pass
- % parameters for envelope (SAM) analysis
- fftpts=2^16;
- env_f=linspace(0,fs,fftpts);
- amf_min=50;
- amf_max=800;
- amf_ind1=min(find(env_f>amf_min));
- amf_ind2=max(find(env_f<amf_max));
- amf_inds=[amf_ind1 amf_ind2];
- amfs=env_f(amf_inds(1):amf_inds(2));
- am_mindur=10e-3;
- am_pk_thrs=-90; % dB
- am_pk_prom=9; % dB
- recbuf=double(recbuf(:,1:16))/32767;
- % highpass filter
- recbuf=filtfilt(hpb,hpa,recbuf);
- aud_time=(0:size(recbuf,1)-1)'/fs;
- % search for calls
- sumblockbuf=sum(abs(recbuf),2);
- % smooth the envelope
- env_dB=20*log10(max(filtfilt(eb,ea,sumblockbuf),1e-4));
- % find calls
- t=[0:length(env_dB)-1]/fs;
- [pk_val pk_pos]=findpeaks(env_dB,'minpeakheight',pk_thrs,'minpeakdistance',round(pk_dist*fs));
- cn=0;
- if isempty(pk_pos)
- return;
- end
- % call cutting and saving
- call_num=1;
- while call_num<=length(pk_pos),
- % 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)
- cn=cn+1;
-
- call_dur=(stop-start)/fs;
-
- % SET UP RESULTS STRUCTURE
- result.call(cn).call_start_sample=pk_pos(cn);
- result.call(cn).call_dur=call_dur;
-
- call=recbuf(start:stop,1:16);
- call_levels=20*log10(std(call));
- [call_level_max call_level_max_pos]=max(call_levels);
- xcall=call(:,call_level_max_pos);
- xct=(0:length(xcall)-1)'/fs*1e3;
- 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;
-
- % quantify AM in call envelope
- if call_dur>=am_mindur,
- if length(xcall)<fftpts,
- lxcall=[xcall;zeros(fftpts-length(xcall),1)];
- end
-
- [modspec,mfs]=pwelch(abs(lxcall),fftpts,0.5*fftpts,fftpts,fs);
- modspec(1:amf_inds(1))=0;
- modspec(amf_inds(2):end)=0;
- modspecdb=10*log10(modspec);
- [AM_val AM_pos]=findpeaks(modspecdb,'minpeakheight',am_pk_thrs,'MinPeakProminence',am_pk_prom);
- if isempty(AM_pos),
- result.call(cn).amf=0;
- else
- % select the highest peak
- [mist,ind]=max(AM_val);
- maxAMfreq=mfs(AM_pos(ind));
- result.call(cn).amf=maxAMfreq;
- end
-
- else
- result.call(cn).amf=0;
- end
- end
- call_num=call_num+1;
- start_sec = start/fs;
- stop_sec = stop/fs;
- disp(sprintf('DONE: call %u, %f-%f s', cn, start_sec, stop_sec))
- end
|