


%% matlab script zf_contours.m

%This is the single channel soundfile - keep this file small  - unless you have huge
%system memory.

signal=wavread('finch_original.wav');
SAMPLING=44100; %sampling rate of the signal
FREQUENCYLIMIT=12; %Upper Range of plot in kHz



%% Parameters


Nfft = 2048;%FFT size 
Nshift = 5;% overlap in the windowed FFT. 

%these parameter settings are more practical, for speed.
Nfft = 1024;
Nshift = 20;

angles = pi/8:pi/8:pi; % compute contours in all these angles - recommended not to change.
TScale = 1.5:.3:2.0; %time scales in milliseconds for this analysis
sonogramindex=2; %which of the time-scales above should be used for the sonogram?
ARThreshold=95; %keep only contours longer than this length percentile 98 or 99 recommended.

NtScale = length(TScale);
Nangle = length(angles);
 %Tested only for files shorter than one second. Memory issues may arise for longer files.

tic
            %to make optimum use of the next line, type 'matlabpool' before
            %running your matlab session, to use all processors in the
            %computation.
parfor sigmacount=1:NtScale %use the matlab parallel computing toolbox if multiple processors are available.
   
            tScale=TScale(sigmacount);
            t=-Nfft/2+0.5:Nfft/2-0.5;
            tScale=(tScale/1000)*SAMPLING;
            w=exp(-(t/tScale).^2);
            dw=(w).*((t)/(tScale^2))*-2;

            q=specgram(signal,Nfft,[],w,Nfft-Nshift)+eps; %gaussian windowed spectrogram
            q2=specgram(signal,Nfft,[],dw,Nfft-Nshift)+eps; %deriv gaussian windowed spectrogram
            dx=(q2./q)/(2*pi); %displacement according to the remapping algorithm
            sonogramrecord{sigmacount}=q;
            
            for angle_variable=1:Nangle
                [sigmacount angle_variable]
                
                theta = angles(angle_variable)+pi/8;
                
                %a quick approximation to the contour detection
                s=-1*(imag(dx*exp(1j*theta))<0)+(imag(dx*exp(1j*theta))>0);
                [gx gy]=gradient(s);
                BW=((-gx*cos(theta+pi/2)+gy*sin(theta+pi/2))>.1);
             
                
                CC = bwconncomp(BW); %build a strucure containing a separate entry for each contour (connected component of BW).
                cc_pix=regionprops(CC,'Area'); %this is the length of each contour
                weightv=[];
               
                for i = 1:length(cc_pix),
                    weightv(i)=cc_pix(i).Area;
                end
   
                a=find(weightv>=prctile(weightv,ARThreshold));
               
                tempv=zeros(size(dx));
                
                for(ik=1:length(a)),
                    ind=CC.PixelIdxList(a(ik));
                    tempv(ind{1})=1;
                end          
                BWallk{sigmacount}{angle_variable}=sparse(tempv);
                
            end
   
    end
                
 % Superimpose all contours from all angles and time scales to produce a single image
 % This final image contains no amplitude information
 
consensus=zeros(size(BWallk{1}{1}));
consensuspower=consensus;

  for sigmacount=1:NtScale-1, 
        
            temp=zeros(size(consensus));
            
           
            
            for angle_variable=1:Nangle,
                
                 if(angle_variable==1),
                
                cv= BWallk{sigmacount}{1}+BWallk{sigmacount+1}{angle_variable}+BWallk{sigmacount}{8};
                consensus=consensus+(cv>1);
                temp=temp+consensus;
                 else
                
             
                cv= BWallk{sigmacount}{angle_variable}+BWallk{sigmacount+1}{angle_variable}+BWallk{sigmacount}{angle_variable-1};
                consensus=consensus+(cv>1); %keep only contour points that show some agreement across neighboring angles or time-scales.
                temp=temp+consensus;
                 end
            end
            consensuspower=temp.*abs(sonogramrecord{sigmacount}); %multiply consenus score by sonogram power for final image
                
  end

  

PIXELRANGE=floor((Nfft/2)*1000*FREQUENCYLIMIT/(SAMPLING/2));
FONTSIZE=15;

[a,b]=size(consensus);
timescale=1000*(b*Nshift)/SAMPLING;
figure(1);
subplot(2,1,1)
imagesc([0 timescale],[0 FREQUENCYLIMIT],log(consensuspower(1:PIXELRANGE,:)+30))
ylabel('frequency [kHz]','FontSize', FONTSIZE);
set(gca,'YDir','normal');
xlabel('time [ms]','FontSize', FONTSIZE);
title('Consensus contour representation','FontSize', FONTSIZE);
colormap(hot)
subplot(2,1,2);


imagesc([0 timescale],[0 FREQUENCYLIMIT],log(abs(sonogramrecord{sonogramindex}(1:PIXELRANGE,:))+30))

titlev=sprintf('Sonogram computed with timescale %6.3g ms',TScale(sonogramindex));
title(titlev,'FontSize', FONTSIZE);
ylabel('frequency [kHz]','FontSize', FONTSIZE);
set(gca,'YDir','normal');
xlabel('time [ms]','FontSize', FONTSIZE);
colormap(hot)

toc

print('-dtiff','-r100','zf_no_pulses.tiff')




