1## Copyright (C) 2008 Carlo de Falco <carlo.defalco@gmail.com>
2##
3## This program is free software: you can redistribute it and/or modify
4## it under the terms of the GNU General Public License as published by
5## the Free Software Foundation, either version 3 of the License, or
6## (at your option) any later version.
7##
8## This program is distributed in the hope that it will be useful,
9## but WITHOUT ANY WARRANTY; without even the implied warranty of
10## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11## GNU General Public License for more details.
12##
13## You should have received a copy of the GNU General Public License
14## along with this program; see the file COPYING. If not, see
15## <https://www.gnu.org/licenses/>.
16
17## -*- texinfo -*-
18## @deftypefn {Function File} {@var{x0} =} zerocrossing (@var{x}, @var{y})
19## Estimates the points at which a given waveform y=y(x) crosses the
20## x-axis using linear interpolation.
21## @seealso{fzero, roots}
22## @end deftypefn
23
24function  retval  = zerocrossing (x,y)
25
26  x = x(:);y = y(:);
27  crossing_intervals = (y(1:end-1).*y(2:end) <= 0);
28  left_ends = (x(1:end-1))(crossing_intervals);
29  right_ends = (x(2:end))(crossing_intervals);
30  left_vals = (y(1:end-1))(crossing_intervals);
31  right_vals = (y(2:end))(crossing_intervals);
32  mid_points = (left_ends+right_ends)/2;
33  zero_intervals = find(left_vals==right_vals);
34  retval1 = mid_points(zero_intervals);
35  left_ends(zero_intervals) = [];
36  right_ends(zero_intervals) = [];
37  left_vals(zero_intervals) = [];
38  right_vals(zero_intervals) = [];
39  retval2=left_ends-(right_ends-left_ends).*left_vals./(right_vals-left_vals);
40  retval = union(retval1,retval2);
41
42endfunction
43
44%!test
45%! x = linspace(0,1,100);
46%! y = rand(1,100)-0.5;
47%! x0= zerocrossing(x,y);
48%! y0 = interp1(x,y,x0);
49%! assert(norm(y0,inf), 0, 100*eps)
50
51%!test
52%! x = linspace(0,1,100);
53%! y = rand(1,100)-0.5;
54%! y(10:20) = 0;
55%! x0= zerocrossing(x,y);
56%! y0 = interp1(x,y,x0);
57%! assert(norm(y0,inf), 0, 100*eps)
58
59%!demo
60%! x = linspace(0,1,100);
61%! y = rand(1,100)-0.5;
62%! x0= zerocrossing(x,y);
63%! y0 = interp1(x,y,x0);
64%! plot(x,y,x0,y0,'x')
65
66%!demo
67%! x = linspace(0,1,100);
68%! y = rand(1,100)-0.5;
69%! y(10:20) = 0;
70%! x0= zerocrossing(x,y);
71%! y0 = interp1(x,y,x0);
72%! plot(x,y,x0,y0,'x')
73