Easy way of finding zero crossing of a function

I am trying to find zero-crossings of a function in Matlab and plot the points where zero-crossing occurs. However, i am not able to find an easy way. I tried http://terpconnect.umd.edu/~toh/spectrum/PeakFindingandMeasurement.htm

and Matlab fnzeros, but i can't understand the spmak and x,y used in these function. The function for which i want to find zero crossing is Euclidean distance function. I will really appreciate if someone can tell me an easy way of doing this.

[Walter Roberson]

I do not see the connection to the tags about image analysis that were provided?

[BlueBee77]

Its used in an image processing project. Finding the peak points in an eclosed boundary

Answer 1

[Adam Cook]

Please see the zerocrossrate function from Signal Processing Toolbox. It has a special indices output that you can use to identify the specific indices where zero crossings occur:

https://www.mathworks.com/help/releases/R2024a/signal/ref/zerocrossrate.html?s_tid=doc_srchtitle#mw_62e2c7a7-b2f5-4768-9b4a-9d68e9e7968d

Note that this will just give discrete time indices. Interpolation would be required to find the continuous time indices.

[Star Strider]

The zerocrossrate function was introduced in R2021b.

Answer 2

[Star Strider]

If your function is a vector of values, you can use this little function to approximate them:

zci = @(v) find(v(:).*circshift(v(:), [-1 0]) <= 0);                    % Returns Zero-Crossing Indices Of Argument Vector

It is also helpful if you want to use fzero or interp1 in a loop to get the exact values.

Note that it returns the indices of the zero-crossings, so you have to do the appropriate vector addressing to get the corresponding x and y coordinates of them.

EDIT — (7 Jul 2020 at 2:54)

Another way of defining ‘zci’ is:

zci = @(v) find(diff(sign(v)));

producing the same result.

[BlueBee77]

Thanks alot for your response. I tried to use it, but why is zci returned empty. Can you explain how it works and how can i plot the indices or the zero-crossing peaks.

And yes, my function is a vector and i am using it inside a loop.

[Star Strider]

You take the ‘y’ vector (your dependent variable) and use it as ‘v’:

t = [1:0.01:50;
y = sin(2*pi*t);
zci = @(v) find(v(:).*circshift(v(:), [-1 0]) <= 0);                    % Returns Zero-Crossing Indices Of Argument Vector
zx = zci(y);

[BlueBee77]

Still this function returns zx empty. Can you explain the "t". is it the number of points that we want to check? and what the 0.01?

[Star Strider]

I would have to have your signal and experiment with it to see what the problem is. The first thing you need to do is to plot it to see if it even has zero-crossings.

This works:

t = [1:0.01:5];                                                         % Time Vector
y = sin(2*pi*t);                                                        % Signal
zci = @(v) find(v(:).*circshift(v(:), [-1 0]) <= 0);                    % Returns Zero-Crossing Indices Of Argument Vector
zx = zci(y);                                                            % Approximate Zero-Crossing Indices

figure(1)
plot(t, y, '-r')
hold on
plot(t(zx), y(zx), 'bp')
hold off
grid
legend('Signal', 'Approximate Zero-Crossings')

[BlueBee77]

Yes, for this function it works. But for my function it doesn't. I dont know why because the function has values, it shouldn't be empty. I am finding this zero-crossing so that i can look for how many peak points are there in the signal. Imagine you have an irregular enclosed boundary and you want to look for all the points(peaks). The paper say they found it by looking for zero-crossing. I would really appreciate if you can help me out. This really took alot of time and no results

[Star Strider]

I have no idea what your signal is. It is possible that your signal has no zero crossings, and the paper is looking at the zero-crossings of something else (perhaps the derivative of the signal). If you want to find the peaks instead, and you have the Signal Processing Toolbox, see if the findpeaks function will work for you.

For example, in my illustration, if the original signal were a cosine, the zero-crossings of the sine curve (the negative derivative of the cosine signal) would be the peaks (and valleys) of the cosine signal.

Answer 3

[Isaac Yuki]

I wrote following function based on the zci function of @Star Strider:

function [number_zeros,zero_crossings] = findzeros(array,samplerate)
%FINDZEROS finds zerocrossings

%Finds the zeros or the nearest values to zero in a function and gives back
%as result the number of zerocrossings and an array containing median of the
%array with the positions of the value that are zero or nearst to zero in
%a zero crossing area, so its the middlest value of the zero crossing area

z = find(diff(sign(array)));    
a = 1;
b = 1;
for i=2:1:length(z)
    if z(i) > z(i-1)+round((samplerate/10000))+1
        a = 1;
        if i == 2
            zci(b,a) = z(i-1);
        end
        zci(end+1,a) = z(i); 
        b = b+1;
    else
        zci(b,a) = z(i); 
        a = a+1;
    end
end

number_zeros = b; %output1

zci2 = [];
zb = [];
zc = [];
zero_crossings = [];
for b = 1:1:number_zeros
    zci2 = zci(b,:);
    for j=1:1:length(zci2)
        if zci2(j) == 0 && j ~=1
            break
        end
        zci3(b,j) = array(zci2(j));
    end
    zb = find(abs(zci3(b,:)) == min(abs(zci3(b,:))));
    zb = zci2(zb);
    if length(zb) <= 1
        zero_crossings = [zero_crossings zb]; %output2
    else
        zero_crossings(end+1) = zb(floor(length(zb)/2)); %outpu2
    end
end
end

[Star Strider]

Thank you!

Also, a reference.

.

[María del Rosario]

Crack!

Answer 4

[xszm]

I think that you can interpolate your data. You can find my results as follow. Thanks for Star Strider, I found it for a long time.

NO interp

interp

[Teong Guan Yew]

Hi YL,

Could you share the code that you used in plotting interp?

Thanks.

[xszm]

Iwl2=360:0.001:740;                                                         % interpolating to 0.001-nm resolution
loc_frequ1=0;
for i=1:c
    y = sgf1(:,i);
    y2=interp1(wl,y,Iwl2);                                                  % interpolating to 0.001-nm resolution                  
    zci = @(v) find(v(:).*circshift(v(:), [-1 0]) <= 0);                    % Returns Zero-Crossing Indices Of Argument Vector
    zx = zci(y2);                                                           % Approximate Zero-Crossing Indices
    
    loczeros1=round(Iwl2(zx));
    locfrequ= ismember(Iwl2,loczeros1);                                     % find same data from Iwl2             
    loc_frequ1=loc_frequ1+double(locfrequ);                  
end

[Ruqaiya Ali]

Can you please tell how to perform zero crossing on a Dataset containing multiple EMG signals?

[Mohamed Jamal]

Hi , xszm

what's c? could you please explain and describe your zero-crossings code with interpolation method by some words?

I have signal called y (my signal is like your second photo signal that you attached here - like sinusoidal signal), how can I use your code in order to find my zero-crossings by interpolation method and plot those zero crossings points on my signal?

[Walter Roberson]

c is number of columns of the signal, which would be the number of channels of the signal.

[Mohamed Jamal]

THanks, so what's the signal in the code? what's this row y = sgf1(:,i) ? because in my matlab it shows me this "Undefined function or variable 'sgf1'".

what's w1? it's not defined at all!!!

thanks

[Walter Roberson]

sgf1 is the signal, arranged with one column per channel. (Be careful, that code will not work properly with a row vector signal.)

The particular code posted expects wl to be the input time associated with the signal.

The code interpolates at 1 kHz between times 360 and 740.

[Mohamed Jamal]

Im trying this code below in my matlab: (my signal is y1 -it's implicitly values of samples, my time (x's axis)=1:length(y1); )

but it doesn't work, could you please help me why and how could I correct it? thanks.

Iwl2=360:0.001:740; % interpolating to 0.001-nm resolution
loc_frequ1=0;
w1=1:length(y1);
for i=1:c
    y = y1(:,i);
    y2=interp1(wl,y,Iwl2);                                                  % interpolating to 0.001-nm resolution                  
    zci = @(v) find(v(:).*circshift(v(:), [-1 0]) <= 0);                    % Returns Zero-Crossing Indices Of Argument Vector
    zx = zci(y2);                                                           % Approximate Zero-Crossing Indices
    
    loczeros1=round(Iwl2(zx));
    locfrequ= ismember(Iwl2,loczeros1);                                     % find same data from Iwl2             
    loc_frequ1=loc_frequ1+double(locfrequ);                  
end

[Walter Roberson]

w1 should not be 1:length(y1) . You should be using wl (lower-case-L not digit 1), and it should be the time vector corresponding to your input signal. If you know your sampling rate, Fs, then

if size(y1,1) == 1; y1 = y1.'); end    %ensure columns of signal
c = size(y1,2);

wl = (0:size(y1,1)-1) / Fs;

Iwl2=360:0.001:740; % interpolating to 0.001-nm resolution
loc_frequ1=0;

for i=1:c
    y = y1(:,i);
    y2=interp1(wl,y,Iwl2);                                                  % interpolating to 0.001-nm resolution                  
    zci = @(v) find(v(:).*circshift(v(:), [-1 0]) <= 0);                    % Returns Zero-Crossing Indices Of Argument Vector
    zx = zci(y2);                                                           % Approximate Zero-Crossing Indices
    
    loczeros1=round(Iwl2(zx));
    locfrequ= ismember(Iwl2,loczeros1);                                     % find same data from Iwl2             
    loc_frequ1=loc_frequ1+double(locfrequ);                  
end

[Star Strider]

Walter — Thank you!

Answer 5

[Nick Hunter]

I have just worked out a quicker bug proof solution, I guess:

clear;
theta = [0:7:360*4,1440]; % Angle Vector (MUST BE a ROW VECTOR)
y = sind(theta);          % Signal Vector (MUST BE a ROW VECTOR)

UpZCi = @(v) find(v(1:end-1) <= 0 & v(2:end) > 0);	% Returns Up Zero-Crossing Indices
DownZCi = @(v) find(v(1:end-1) >= 0 & v(2:end) < 0);    % Returns Down Zero-Crossing Indices
ZeroX = @(x0,y0,x1,y1) x0 - (y0.*(x0 - x1))./(y0 - y1); % Interpolated x value for Zero-Crossing 

ZXi = sort([UpZCi(y),DownZCi(y)]);
ZX = ZeroX(theta(ZXi),y(ZXi),theta(ZXi+1),y(ZXi+1));

% === Checking for zero at the ignored value ===
if y(end)==0
    ZX(end+1) = theta(end);
end
% ==============================================

figure(1)
plot(theta, y, '-b')
hold on;
plot(ZX,zeros(1,length(ZX)),'ro')
grid on;
legend('Signal', 'Interpolated Zero-Crossing')

[Mohamed Jamal]

Hi Nick Hunter I tried your code but unfortunately it doesn't work for me, any help pelase? thanks alot.

it shows me this error "Error using horzcat

Dimensions of arrays being concatenated are not consistent."

the code that I tried is this below, y1 is my sinusoidal signal related to time.

y = y1;          % Signal
theta=1:length(y1); %it's time (x axis)
UpZCi = @(v) find(v(1:end-1) <= 0 & v(2:end) > 0);	% Returns Up Zero-Crossing Indices
DownZCi = @(v) find(v(1:end-1) >= 0 & v(2:end) < 0);    % Returns Down Zero-Crossing Indices
ZeroX = @(x0,y0,x1,y1) x0 - (y0.*(x0 - x1))./(y0 - y1); % Interpolated x value for Zero-Crossing 
ZXi = sort([UpZCi(y),DownZCi(y)]);
ZX = ZeroX(theta(ZXi),y(ZXi),theta(ZXi+1),y(ZXi+1));
% === Checking for zero at the ignored value ===
if y(end)==0
    ZX(end+1) = theta(end);
end
% ==============================================
figure(1)
plot(theta, y, '-b')
hold on;
plot(ZX,zeros(1,length(ZX)),'ro')
grid on;
legend('Signal', 'Interpolated Zero-Crossing')

[Nick Hunter]

The problem is not with the code that I posted. THere is something wrong with the code before or after that. Check your "horzcat" function.

[Walter Roberson]

Nick:

y = [1 0 1 0 1  0 -1].';
UpZCi = @(v) find(v(1:end-1) <= 0 & v(2:end) > 0);	% Returns Up Zero-Crossing Indices
DownZCi = @(v) find(v(1:end-1) >= 0 & v(2:end) < 0);    % Returns Down Zero-Crossing Indices
UpZCi(y),DownZCi(y)

ans =
     2
     4
ans =
     6

However, you cannot [] together [2;4] and [6]

Your code is assuming that the result of UpZCi and DownZCi are row vectors, but that will not be the case if y is a column vector. You did not document an orientation requirement.

[Nick Hunter]

Thank you so much, Walter. This is a good point. I added this condition in comments.

Answer 6

[Mitch Lautigar]

There are many ways to try and skin this problem, many people have tried and find varying different levels of success. Here's a function you can use that requires a sinusoidal waveform to approximate 0 intersections.

[out_array] = signal_manip(s_in)
sign_array = []; %predeclaration
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Following for loop checks every spot to create the following numbers:
%For positive numbers, a "1" is placed in an array.
%For negative numbers, a "-1" is placed in an array.
%For a zero, a "0" is placed in an array.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i = 1:length(s_in)
        if s_in(1,i) > 0
            curr_sign = 1;
        elseif  s_in(1,i) < 0
            curr_sign = -1;
        else
            curr_sign = 0;
        end  %end "if s_in > 0"
sign_array = [sign_array,curr_sign]; %gives an output array that shows you all negative and positive numbers
end %end for i = 1:length(s_in)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%following for loop looks for the change points and handles the spots where a 0 occurs.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
curr = sign_array(1,1); %Starting comparison Point
out_array = []               %Predeclaration of array
for i = 2:length(s_in)
    if ( sign_array(1,i) ~= curr) %if number change detected, mark it.
        out_array = [out_array,i];
    elseif  (sign_array(1,i) == 0) %if we have a zero on the graph, mark it and skip i forward by 1 so as to not have a false hit.
        out_array = [out_array,i];
        i = i + 1;
    end
curr = sign_array(1,i);
end

[Mohamed Jamal]

Hi Mitch Lautigar

I tried your code but Im not getting the zero crossings points on my signal (it's vector of values called y)

the problem when I try your code in matlab is this : "Undefined function or variable 'signal_manip'."

could you please attach/share me your complete code for zero crossings by interpolation? thanks alot

[Walter Roberson]

The first line should be

function [out_array] = signal_manip(s_in)

That is, the code is really the definition for a function named signal_manip that takes an input signal that is expected to be a row vector.

[Mohamed Jamal]

OOOPS didn't notice that there's keyword function, my bad, but whatever it still gives me bugs:

signal_manip(y1) %y1 is my signal (it's vector of values)

Index in position 2 exceeds array bounds (must not exceed 1).

Error in signal_manip (line 10)

if s_in(1,i) > 0

[Walter Roberson]

The person who posted it left out the function keyword accidentally, so it is not surprising you missed it!

In order to get the error message you are seeing when you are passing in a vector of values, then you must have passed in a column vector of values -- but as I was careful to point out before, it expects a row vector of values.

[Mohamed Jamal]

Yes it works fine now but this function is returnning array of zero crossings points right? because Im not getting my signal plot with zero crossings points marked in .. not getting a plot at all , just array!

Answer 7

[UMAIR RASOOL]

close all

clear all

clc

n=-10:10;

a=.9;

y=zeros(size(n));

for k=-5:5

temp=find(n-(2*k)==0); %find xaxis where it is zero(zero crossing)

y(temp)=1; %palce the value at the zero crossing

x=n.*(sum(a^k)*y);

end

plot(n,x)

xlabel('n')

ylabel('x(n)')

Maximum=max(x);

Answer 8

[Mitch Lautigar]

So I took the code previously used and modified it to what I needed (finding zero point crossing for a sin/cos wave. The code is below

function [freq_val] = peakfind(dm_allow)

%Step 1: Load file and find zero crossings

load('s_m.mat','s_m2','t')

zci = @(v) find(v(:) .* circshift(v(:), [-1 0]) <= 0);

zx = zci(s_m2)';

rp_zx = [];

%Step 2: manipulate above values so as to find unique values.

zx2 = round(zx ./ max(zx),2);

for i = 1:length(zx2)-1

if round(zx2(1,i),3) ~= round(zx2(1,i+1),3)

rp_zx = [rp_zx,zx2(1,i+1)];

end

end

end_array = [1];

%Eliminate redundant numbers

for i = 1:length(rp_zx)

zeta = find(zx2 == rp_zx(1,i));

end_array = [end_array,zeta(1,1)];

end

actual_spots = zx(end_array);

t_val = [];

%convert above numbers to to time values.

for i = 1:length(end_array)-1

low_spot = actual_spots(1,i);

high_spot = actual_spots(1,i+1);

max_spot = max(s_m2(low_spot:high_spot));

min_spot = abs(min(s_m2(low_spot:high_spot)));

if max_spot > min_spot

beta = find(s_m2(low_spot:high_spot) == max_spot);

t_val = [t_val,t(beta+low_spot)];

elseif max_spot < min_spot

%do nothing

end

end

%convert time values to frequency values.

if length(t_val) > 1

freq_val = [];

for i = 1:length(t_val)-1

freq_val = [freq_val,1/(t_val(i+1)-t_val(i))];

end

else

error('Cycle time not long enough')

end

end

If you use my code, simply change the following:

zx = zci(s_m2)' %replace s_m2 with whatever signal you are wanting to find zero crossings for.

[Mathieu NOE]

dear all

FYI - this is the function I use since long time - have always worked perfect for all type of signals

function [ind,t0,s0] = crossing_V6(S,t,level,imeth)

% [ind,t0,s0,t0close,s0close] = crossing_V6(S,t,level,imeth)

% CROSSING find the crossings of a given level of a signal
%   ind = CROSSING(S) returns an index vector ind, the signal
%   S crosses zero at ind or at between ind and ind+1
%   [ind,t0] = CROSSING(S,t) additionally returns a time
%   vector t0 of the zero crossings of the signal S. The crossing
%   times are linearly interpolated between the given times t
%   [ind,t0] = CROSSING(S,t,level) returns the crossings of the
%   given level instead of the zero crossings
%   ind = CROSSING(S,[],level) as above but without time interpolation
%   [ind,t0] = CROSSING(S,t,level,par) allows additional parameters
%   par = {'none'|'linear'}.
%	With interpolation turned off (par = 'none') this function always
%	returns the value left of the zero (the data point thats nearest
%   to the zero AND smaller than the zero crossing).
%
%	[ind,t0,s0] = ... also returns the data vector corresponding to 
%	the t0 values.
%
%	[ind,t0,s0,t0close,s0close] additionally returns the data points
%	closest to a zero crossing in the arrays t0close and s0close.
%
%	This version has been revised incorporating the good and valuable
%	bugfixes given by users on Matlabcentral. Special thanks to
%	Howard Fishman, Christian Rothleitner, Jonathan Kellogg, and
%	Zach Lewis for their input. 

% Steffen Brueckner, 2002-09-25
% Steffen Brueckner, 2007-08-27		revised version

% Copyright (c) Steffen Brueckner, 2002-2007
% [e-mail]

% check the number of input arguments
error(nargchk(1,4,nargin));

% check the time vector input for consistency
if nargin < 2 | isempty(t)
	% if no time vector is given, use the index vector as time
    t = 1:length(S);
elseif length(t) ~= length(S)
	% if S and t are not of the same length, throw an error
    error('t and S must be of identical length!');    
end

% check the level input
if nargin < 3
	% set standard value 0, if level is not given
    level = 0;
end

% check interpolation method input
if nargin < 4
    imeth = 'linear';
end

% make row vectors
t = t(:)';
S = S(:)';

% always search for zeros. So if we want the crossing of 
% any other threshold value "level", we subtract it from
% the values and search for zeros.
S   = S - level;

% first look for exact zeros
ind0 = find( S == 0 ); 

% then look for zero crossings between data points
S1 = S(1:end-1) .* S(2:end);
ind1 = find( S1 < 0 );

% bring exact zeros and "in-between" zeros together 
ind = sort([ind0 ind1]);

% and pick the associated time values
t0 = t(ind); 
s0 = S(ind);

if strcmp(imeth,'linear')
    % linear interpolation of crossing
    for ii=1:length(t0)
        %if abs(S(ind(ii))) > eps(S(ind(ii)))    % MATLAB V7 et +
        if abs(S(ind(ii))) > eps*abs(S(ind(ii)))    % MATLAB V6 et -    EPS * ABS(X)
            
            % interpolate only when data point is not already zero
            NUM = (t(ind(ii)+1) - t(ind(ii)));
            DEN = (S(ind(ii)+1) - S(ind(ii)));
            slope =  NUM / DEN;
            terme = S(ind(ii)) * slope;
            t0(ii) = t0(ii) - terme;
            % I'm a bad person, so I simply set the value to zero
            % instead of calculating the perfect number ;)
            s0(ii) = 0;
        end
    end
end

% % Addition:
% % Some people like to get the data points closest to the zero crossing,
% % so we return these as well
% [CC,II] = min(abs([S(ind-1) ; S(ind) ; S(ind+1)]),[],1); 
% ind2 = ind + (II-2); %update indices 
% 
% t0close = t(ind2);
% s0close = S(ind2);