p_discolor_auditory_feedback/vocalisations/scripts_for_reference/extract_adult_call_parameters.m

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