声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 5145|回复: 20

[综合讨论] 在 由离散点绘制的曲线上 找出零点对应的横坐标

[复制链接]
发表于 2009-1-3 21:53 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
有以下代码生成了 图一
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 ;

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

y.txt

3 KB, 下载次数: 11

回复
分享到:

使用道具 举报

发表于 2009-1-3 22:16 | 显示全部楼层

回复 楼主 chenjc18 的帖子

按照您的思路,定义一个比较小的数,比如t=0.00001;然后find(y<t&y>t);试试
 楼主| 发表于 2009-1-4 09:03 | 显示全部楼层

回复 沙发 ch_j1985 的帖子

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

请问,我的这个问题可不可以从别的角度来解决呢?
      我觉得我的这个思路实现起来,太困难了
发表于 2009-1-4 09:30 | 显示全部楼层
插值了又何必要那么小步长,
换个思路哈
xi=interp1(y,x,0)

评分

1

查看全部评分

 楼主| 发表于 2009-1-4 10:14 | 显示全部楼层

回复 地板 wwbeyondww 的帖子

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

如何得到曲线与零线交点的横坐标?还请wwbeyondww您多多指教
发表于 2009-1-4 11:25 | 显示全部楼层
可以用数学的知识,因为y是连续函数(假设如此),对相邻的y值作比较,如果有正有负,则说明必与x轴相交,即存在零点,此时可用直线拟合方法求出y=0的x值,好像有个零点定理吧。
发表于 2009-1-4 11:26 | 显示全部楼层
总感觉楼主的appoach方式不太好, 不管取多密, 最後还不是得找到穿越零点前後两点, 然後进行线性插分! 即然如此而不一开始即如此? 进行interp1好像根本不需要!
发表于 2009-1-4 13:40 | 显示全部楼层

回复 5楼 chenjc18 的帖子

想了一下,我是找不到好办法了,搞了个循环,
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 编辑 ]

评分

1

查看全部评分

发表于 2009-1-4 21:04 | 显示全部楼层
原帖由 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 编辑 ]
 楼主| 发表于 2009-1-4 21:11 | 显示全部楼层

回复 8楼 wwbeyondww 的帖子

谢谢您
我根据您的思想 把您给的程序改正了一下
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
谢谢您

评分

1

查看全部评分

 楼主| 发表于 2009-1-4 21:14 | 显示全部楼层

回复 9楼 ChaChing 的帖子

谢谢您的关注
8楼给的方法,我从7点就开始研究,到现在才改正过来。
呵呵 谢谢您给我们这些菜鸟的帮助
 楼主| 发表于 2009-1-4 21:19 | 显示全部楼层

回复 6楼 friendchj 的帖子

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

我的问题已经解决了 呵呵:@)
发表于 2009-1-4 21:58 | 显示全部楼层

回复 10楼 chenjc18 的帖子

不好意思,哈哈哈
我开始写程序的时候好像把y取了转置
写上来的时候搞成对x取转置了,
所以后面都不对了,哈哈。
贴出来之后没有重新验算一遍
发表于 2009-1-5 08:56 | 显示全部楼层

回复 10楼 chenjc18 的帖子

本想算了, 反正楼主的例子没遇到, 但还是提醒下!
若原有资料就刚好有零点, 楼主的判断式就会少掉!! 但改为<=同样会有些问题(多解了)!!
 楼主| 发表于 2009-1-5 14:49 | 显示全部楼层

回复 14楼 ChaChing 的帖子

我验证了您的想法,没错确实出现了重解,甚至会报错,请看下面的例子
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);

我正在尝试,将原离散序列 以各个零点为界进行拆分 然后对每段进行求取零点
呵呵 如果可行 我会发到论坛上的
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-25 11:37 , Processed in 0.083944 second(s), 27 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表