p_discolor_auditory_feedback/vocalisations/scripts_for_reference/extract_pup_call_parameters.m

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