% AUTHORS: Ana Marta Dias and Pedro Sousa
% Zero Crossing Counter Function

function [NumZC] = func_ZC(emg, low_threshold, high_threshold, Fs)

% 1. EMG Signal Filtering
filt = ZC_HighPass;
filtemg = filt(emg); % High-Pass Filtering 150 Hz

% 2. Add hysteresis to the counter
below = filtemg <= low_threshold;     % Below Threshold
above = filtemg >= high_threshold;      % Above Threshold
Z = ~below & ~above;          % Between Thresholds

for i = 1:size(filtemg, 2)
    z0 = find(Z(:, i) == 0, 1, 'first');    % Find first crossing
    if ~isempty(z0)                         % If any crossing found
        Z(1:z0, i) = false;
        z1 = find(Z(:, i) == 1);
        for j = 1:length(z1)
            below(z1(j), i) = below(z1(j)-1, i);
        end
    end
end

D = [0*below(1, :); abs(diff(below, 1 , 1))];
vec = filtemg;
vec(D == 0) = 0;
vec(vec ~= 0) = 1;

% 3. Count Zero Crossings
% To do so, windows of 1 second, with 75% overlap, will be considered
NumZC = zeros(0);
for n = 1:(0.25*Fs):(length(emg)-Fs)
    NumZC(length(NumZC) + 1) = sum(vec(n:n+Fs));
end

% Then, a threshold will be applied.
% According to the literature, it can be in the range [95, 255].
% Threshold = 145 will be set.

end

