% featurext_wk.m - updated 14/3/05 
%
% updated 31/3/05 - added feature(12)
%
% Function to perform festure extraction on calls
% For inputted spectrogram data (spec) program returns 
% 10 features of the signal which can be selected to be calculated by array (sel):
% 1.  maximum frequency     2.  minimum frequency
% 3.  bandwidth             4.  entropy of call
% 5.  length of call        6.  peak frequency
% 7.  peak amplitude        8.  time to max freq.
% 9.  time to min freq      10. curve fitting (accurate)
% 11. curve fitting (fast)  12. fraction of time with call content above
%                               threshold (50kHz) 

function [feature] = featurext_wk(wave, sel, timeint,prefs)

feature = zeros(1,13);
    
q=sum(wave.spec/length(wave.spec),2);

rescale = find(q > 0);

% max and min frequency and bandwidth
f_max = ((wave.fs*0.5)/size(wave.spec,1))*max(rescale);
f_min = ((wave.fs*0.5)/size(wave.spec,1))*(min(rescale)-1);
q1=sum(wave.spec/length(wave.spec),1);
if find(sel == 1) > 0
    feature(1) = f_max;
end

if find(sel == 2) > 0
    feature(2) = f_min;
end

if find(sel == 3) > 0
    feature(3) = f_max - f_min;
end

% entropy of call
if find(sel == 4) > 0
	counts = hist(wave.spec,size(wave.spec,1));
	counts = counts/sum(counts(:));
	k = find(counts);
	feature(4) = -sum(counts(k).*log2(counts(k)));
end

% call length calculation
if find(sel == 5) > 0
    q1=sum(wave.spec/length(wave.spec),1);
    feature(5) = timeint*(max(find(q1 > 0)) - min(find(q1 > 0)));
end

% peak frequency and amplitude calculation
if find(sel == 6) > 0
    feature(6) = ((wave.fs*0.5)/size(wave.spec,1))*mean(find(q==max(q)));
end

if find(sel == 7) > 0
    feature(7) = q(round(mean(find(q==max(q)))));
end

% times to max and min frequencies
if find(sel == 8) > 0
    feature(8) = timeint*min(find(wave.spec(max(rescale),:) > 0));
end

if find(sel == 9) > 0
    feature(9) = timeint*min(find(wave.spec(min(rescale),:) > 0));
end

% curve fitting data  
% slower more accurate method
if find(sel == 10) > 0
	spec_hold = wave.spec(:,min(find(q1 > 0)):max(find(q1 > 0)));
	reord = sort(reshape(spec_hold,size(spec_hold,1)*size(spec_hold,2),1));
	minval = reord(min(find(reord)));
	spec_hold_r = round(spec_hold/minval);
	bw_spec = im2bw(spec_hold_r,0.1);
	no_pixs = sum(sum(bw_spec));
	[ny,nx] = find(spec_hold_r);
	
	x = [];
	y = [];
	m = 0;
	
	for l = 1:no_pixs
        for k = 1:spec_hold_r(ny(l),nx(l))
            y(l+m) = ny(l);
            x(l+m) = nx(l);
            m = m + 1;
        end
        m = m - 1;
	end
	
	sx = (x - mean(x))./std(x);  
	[p,s] = polyfit(sx,y,5);  
	[fit,delta] = polyval(p,sx,s);
    if mean(delta)>50
        feature(10)=50;
    else
        feature(10) = mean(delta);
    end
end

% quicker less accurate method
if find(sel == 11) > 0
    if length(find(wave.spec))<2
        feature(11)=0;
    else
        spec_hold = wave.spec(:,min(find(q1 > 0)):max(find(q1 > 0)));

        reord = sort(reshape(spec_hold,size(spec_hold,1)*size(spec_hold,2),1));
        minval = reord(min(find(reord)));
        spec_hold1 = round(spec_hold/minval);

        [y,x] = find(spec_hold1);

        sx = (x - mean(x))./std(x);
        [p,s] = polyfit(sx,y,3);
        [fit,delta] = polyval(p,sx,s);
        if mean(delta)>50
            feature(11)=50;
        else
            feature(11) = mean(delta);
        end
    end
end

if find(sel == 12) > 0
    %prefs
    threshold=round( (prefs.f12thresh/(0.5*wave.fs)) * (size(wave.spec,1)) );
    spec_chop=wave.spec(threshold:129,:);
    len=max(find(sum(wave.spec,1)))-min(find(sum(wave.spec,1)));
    len_thresh=length(find(sum(spec_chop,1)));
    if isempty(len_thresh) || len==0 || isempty(len) || len_thresh==0
        feature(12)=0;
    else
        feature(12)=100*len_thresh/len;
    end    
end    

if find(sel==13) > 0
    feature(13)=sum(wave.spec(:))/size(wave.spec,2);
end