function zcr = get_zcr(SIGNAL, FRAME_LEN, TIME_OVERLAP, THRESH)
% ENERGY_ZCR computes the energy and zero-crossings in SIGNAL
%   SIGNAL = input signal
%   FRAME_LEN = length of window over which the energy is summed and zero-crossings are computed
%   TIME_OVERLAP = number of samples for frame windows to overlap
%   THRESH = threshhold percentage of signal amplitude
%   energy and zrc are arrays of length length(signal)/(frame_len - time_overlap)
% Harriet Fell 2002, using ideas from Childers' ti_e_zcr.m

% Normalize the signal amplitude
%signal = 14000/max(abs(signal))*signal;

% Remove the components around d.c.
signal = filtfilt([1 -1],[1 -.99], SIGNAL);

% Calculate energy and zeros
window_len = FRAME_LEN - TIME_OVERLAP;
signal_len = length(SIGNAL);

if (FRAME_LEN > signal_len)
    disp(['Please use a frame_len < ',int2str(signal_len)]);
else
    num_windows = floor(signal_len/window_len);
    energy = zeros(1, num_windows);
    zcr = zeros(1, num_windows);

    % Find start and end for next sample
    for z = 1:num_windows
        sample_start = (z-1)*window_len + 1;
        if (sample_start + FRAME_LEN -1 > signal_len)
            sample_end = signal_len;
        else
            sample_end = sample_start + FRAME_LEN - 1;
        end
        
        frame = SIGNAL(sample_start:sample_end);
        zcr(z) = zcr_cnt(frame, THRESH);
    end
end
return;
 
 %---------------------------------------------------------
 %  Local function
 %---------------------------------------------------------
function count = zcr_cnt(DATA, THRESH);
% DATA is a signal array
% THRESH is a percentage, 0 <= THRESH <= 100.
% DATA values are set to 0 if their magnitude is less than
% (THRESH/100)*max(abs(DATA))
% count is the number of zero crossings
% Harriet Fell 2004
nmax = THRESH * max(abs(DATA))/100;

data_above = (DATA > nmax);     % 1 if data > nmax, 0 otherwise
data_below = (DATA < -nmax);     % 1 if data < -nmax, 0 otherwise
data_above = DATA.*data_above;
data_below = DATA.*data_below;
data = data_above + data_below; % values between nmin and nmax are now 0
%data = data';
% for i=2:length(data) 
%     if (data(i)==0)data(i)=data(i-1);
%     end
% end
data=data(data~=0);
data_new = (data.*circshift(data,1) < 0); % 1 if data(i) & data(i+1) have opposite sign
npoint = length(data_new);
count = sum(data_new(2:npoint-1));

return