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