在 由离散点绘制的曲线上 找出零点对应的横坐标

chenjc18

有以下代码生成了 图一
y=[-0.5572   -3.3084   -1.6488   -1.3941   -0.0379   -0.2872   -0.7380......];%大小为256 的向量(在附件中)
x=linspace(0,20,256);
zero_1=(zeros(1,length(x)));  %生成零线点
plot(x,y,'b',x,zero_1,'r');

我想要得到，图中曲线与零线交点的横坐标。

我的思路是：如果能找到，零点y（i）=0的下标值 i ，即可根据此下标找到横坐标值x（i）（因为，画图时，每一个纵坐标点 y，均一一对应着一个横坐标点 x）。
按照此思路，有：
使用find（y==0）但得到的是：Empty matrix: 1-by-0，
我个人认为可能是，曲线是由离散点用直线直接相连画出的，而与零线的交点完全是在 离散点之间的连线上。
所以我曾经尝试插值的方式，找零点 代码如下
   
   xi=0:0.0001:20;
   yi=interp1(x,y,xi);
yi即是经线性插值获得的纵坐标点，但是使用find 函数 还是返回Empty matrix: 1-by-0 ；

我的问题是：如何获得 图中曲线与零线交点的横坐标？
                     我的思路可不可行？
先谢谢了

ch_j1985

按照您的思路，定义一个比较小的数，比如t=0.00001;然后find(y<t&y>t);试试

chenjc18

我尝试了a=0.0001;
               find(abs(yi)<a)  但只返回三个下标，而从图上可以看出零点将近三十个
我想过：
       虽然能通过设定a，如增大a  可以返回更多的点，但这些点不一定保证就是靠近某一零点的唯一的点（会出现这种情况  ÷ ，靠近同一零点的两个离散点均在[-a,a]里)。
      
我曾经尝试继续加大插值密度
  xi=0:0.000001:20;
   yi=interp1(x,y,xi);
但是计算机，一算就死机了。

请问，我的这个问题可不可以从别的角度来解决呢？
      我觉得我的这个思路实现起来，太困难了

wwbeyondww

插值了又何必要那么小步长,
换个思路哈
xi=interp1(y,x,0)

chenjc18

首先 谢谢您
尝试了 您提供的思路方法
但是 我只能得到一个值 xi = 13.5833
我看到help中xi=interp1(y,x,a) 要求a是向量或是数组，才能得到一系列的零点对应的横坐标。我尝试将a定义成全零向量，但得到的结果只是重复上面的结果：  13.5833  13.5833  13.5833 ...

如何得到曲线与零线交点的横坐标？还请wwbeyondww您多多指教

friendchj

可以用数学的知识，因为y是连续函数(假设如此)，对相邻的y值作比较，如果有正有负，则说明必与x轴相交，即存在零点，此时可用直线拟合方法求出y=0的x值，好像有个零点定理吧。

ChaChing

总感觉楼主的appoach方式不太好, 不管取多密, 最後还不是得找到穿越零点前後两点, 然後进行线性插分! 即然如此而不一开始即如此? 进行interp1好像根本不需要!

wwbeyondww

想了一下,我是找不到好办法了,搞了个循环,
y=[-0.5572   -3.3084   -1.6488   -1.3941   -0.0379   -0.2872   -0.7380......];%大小为256 的向量(在附件中)
x=linspace(0,20,256);
x=x';
[a,b]=find(y(1:255).*y(2:256)>0);
c=[y(b);y(b)+1];
d=[x(b);x(b)+1];
for k=1:length(b),
X(k)=interp1([c(1,k), c(2,k)],[d(1,k),d(2,k)],0);
end

X =
  Columns 1 through 4
   0.60005859334129   0.64380698910120   0.87037318032210   1.01835006625789
  Columns 5 through 8
   1.17264125060563   1.98487178693200   2.18191884898283   2.28237073503462
  Columns 9 through 12
   4.00352835325873   6.36007130124777   7.17001600086784   8.06632939853850
  Columns 13 through 16
   9.06452415112386   9.36163861870207   9.51294812848588  10.43708783962015
  Columns 17 through 20
  10.96006869901245  11.13835762776112  11.33807569539444  11.43431372549020
  Columns 21 through 24
  11.45713423831071  11.80171611702678  13.76608873755757  13.88224490898288
  Columns 25 through 28
  14.26936338839998  14.30467571644042  15.44961107624427  15.63123635562499
  Columns 29 through 32
  15.76971510070616  17.79951854823776  17.80694904224316  18.23896592005688
  Columns 33 through 34
  18.37779743706101  18.89690920161748

[ 本帖最后由 wwbeyondww 于 2009-1-4 13:42 编辑 ]

ChaChing

原帖由 wwbeyondww 于 2009-1-4 13:40 发表
想了一下,我是找不到好办法了,搞了个循环,
y=[-0.5572   -3.3084   -1.6488   -1.3941   -0.0379   -0.2872   -0.7380......];%大小为256 的向量(在附件中)
x=linspace(0,20,256); ...

怎麽没人发现异样! 楼上提供的答案是对的, 但所给的程序是错的!
我想可能复制错了吧!
不好意思! 帮忙改了下!

load y.txt; x=linspace(0,20,256); x=x'; zero_1=(zeros(1,length(x))); plot(x,y,'b',x,zero_1,'r');
a=y(1:end-1).*y(2:end); b=find(a<=0); xx=zeros(size(b));
c=[y(b),y(b+1)]; d=[x(b),x(b+1)];
for k=1:length(b), xx(k)=interp1([c(k,1), c(k,2)],[d(k,1),d(k,2)],0); end

[ 本帖最后由 ChaChing 于 2009-1-5 08:48 编辑 ]

chenjc18

谢谢您
我根据您的思想 把您给的程序改正了一下
load y.txt;
x=linspace(0,20,256);
x=x';
b=find(y(1:255).*y(2:256)<0);
c=[y(b) y(b+1)];
d=[x(b) x(b+1)];
for k=1:length(b)
X(k)=interp1([c(k,1), c(k,2)],[d(k,1),d(k,2)],0);
end
谢谢您

chenjc18

谢谢您的关注
8楼给的方法，我从7点就开始研究，到现在才改正过来。
呵呵 谢谢您给我们这些菜鸟的帮助

chenjc18

您好，您的思路很正确 8楼的wwbeyondww给出了实现此算法思想的代码（但是，代码有点错误，正确的在9楼 10楼 已经给出）
谢谢您

我的问题已经解决了 呵呵:@)

wwbeyondww

不好意思，哈哈哈
我开始写程序的时候好像把y取了转置
写上来的时候搞成对x取转置了，
所以后面都不对了，哈哈。
贴出来之后没有重新验算一遍

ChaChing

本想算了, 反正楼主的例子没遇到, 但还是提醒下!
若原有资料就刚好有零点, 楼主的判断式就会少掉!! 但改为<=同样会有些问题(多解了)!!

chenjc18

我验证了您的想法，没错确实出现了重解，甚至会报错，请看下面的例子
y=[1 2 -3 0 -1 2 0 -1 4 -1];
x=linspace(0,20,10);
x=x';
y=y';
plot(x,y);
b=find(y(1:end-1).*y(2:end)<=0);
c=[y(b) y(b+1)];
d=[x(b) x(b+1)];
X=zeros(size(b));
for k=1:length(b)
X(k)=interp1([c(k,1), c(k,2)],[d(k,1),d(k,2)],0);
end

X=
    3.1111
        6.6667
        6.6667
    9.6296
        13.3333
        13.3333
   16.0000
   19.5556
在本身就存在零点的坐标点，获得重解。
当y=[1 2 -3 0 -1 2 0 0 -1 4 -1]; %两个零点相邻
此时，报错
程序运行结果：
插值区间：
c =  2    -3
    -3     0
     0    -1
    -1     2
     2     0
     0     0
     0    -1
    -1     4
     4    -1
??? Error using ==> interp1
The values of X should be distinct.

Error in ==> x_y_zero at 15
X(k)=interp1([c(k,1), c(k,2)],[d(k,1),d(k,2)],0);

我正在尝试，将原离散序列 以各个零点为界进行拆分 然后对每段进行求取零点
呵呵 如果可行 我会发到论坛上的