% directory where the WAV-files of the birds are stored:
dpath = fullfile(mildatapath,'birds');

fname = 'PC13_20090531_050000_0.wav';
%fname = 'PC8_20090531_050000_99.wav';
[signal,fs] = wavread(fullfile(dpath,'wavs',fname));

[S,f,t] = spectrogram(signal,512,256,[],fs);
fmax = 256;%Defined in paper

figure(2); clf; showspec(S,f,t);
[S,f,t] = spectrogram(signal,512,256,fmax,fs);

figure(2); imagesc(t,f,abs(S));
axis xy; xlabel('time (s)'); ylabel('freq.');

figure(3); clf; showspec(imag(S), f,t);

%%
% smooth and threshold the spectrogram?
G = fspecial('gaussian',[5 5],2); % window 5x5, sigma=2
I = imfilter(abs(S),G,'same');
t = dd_threshold(I(:),0.8);
mask = (I>t);
mask(f<2000,:) = 0;

figure(4); clf; imagesc(mask); axis xy;
% measure something:
x = regionprops(bwlabel(mask),abs(S));
y1 = x([x.Area]>10)



%%
s2 = specwhiten(S,fmax); 
%figure(5); imagesc(t,f,s2);
%axis xy; xlabel('time (s)'); ylabel('freq.');

% do the same smoothing but change the threshold
G = fspecial('gaussian',[5 5],2); 
I = imfilter(abs(s2),G,'same');
t = dd_threshold(I(:),0.95);
mask = (I>t);
mask(f<2000,:) = 0;

figure(6); clf; imagesc(mask); axis xy;    %Looks pretty similar, but less small blobs 

labeledx = bwlabel(mask);
Nseg = max(unique(labeledx)); %Segment labels, segment 0 is background

absim = abs(S);
realim = real(S);
imagim = imag(S);


maskfeats = nan(Nseg,3);
absfeats = nan(Nseg, 7);


for i=1:Nseg
   
    ix = (labeledx == i);
    
    pixtotal = sum(sum(ix));
    pixheight = max(sum(ix,1));
    pixwidth = max(sum(ix,2));

    maskfeats(i,:) = [pixtotal pixheight pixwidth];
    
    
    seg = absim(ix);
    absfeats(i,1) = mean(seg); 
    absfeats(i,2) = std(seg); 
    
    absfeats(i,3) = quantile(seg(:),0);
    absfeats(i,4) = quantile(seg(:),0.25);
    absfeats(i,5) = quantile(seg(:),0.5);
    absfeats(i,6) = quantile(seg(:),0.75);
    absfeats(i,7) = quantile(seg(:),1);
    
end

mildata = [maskfeats absfeats];