声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

12
返回列表 发新帖
楼主: chenjc18

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

[复制链接]
 楼主| 发表于 2009-1-5 16:56 | 显示全部楼层

针对含有两(及以上)个零的 离散序列 计算零点算法

针对ChaChing 兄 对上述程序提出的见解
我将问题暂时解决了一下:
请看程序:
此程序只是使用于含有两个零(及以上)的离散序列,对于含有1个或者连续出现两个零点的情形 可以按照我提出的 分段求解 思路 对下面的程序稍加修改 即可使用。

请大家,试试看,还有什么问题。
鉴于我刚刚学习MATLAB,其中出现幼稚的语句也在所难免,还请大家多提出建议。谢谢了
  1. clear;close;

  2. y=[1 2 -3 -1 0 2 8 -1 -1 4 0 -1 1 -4 7 -2];
  3. x=linspace(0,20,length(y));
  4. zero=zeros(length(x),1);
  5. x=x';
  6. y=y';
  7. plot(x,y,'b',x,zero,'r');
  8. ind=find(y==0);
  9. X1=x(ind); %找到自身就有的零点对应的x
  10. %段一 从段首到第一个自身的零点
  11. %以下是对段一求解零点
  12. a=y(1:(ind(1)-1));
  13. b=find(a(1:end-1).*a(2:end)<0);
  14. c=[a(b) a(b+1)];
  15. d=[x(b) x(b+1)];
  16. for k=1:length(b)
  17. X2(k)=interp1([c(k,1), c(k,2)],[d(k,1),d(k,2)],0);
  18. end
  19. %中间的各段
  20. clear a b c d ;
  21. for i=1:(length(ind)-1)
  22. if((ind(i+1)-ind(i)~=1)|(ind(i+1)-ind(i)~=2))
  23. a=y((ind(i)+1):(ind(i+1)-1));
  24. a1=x((ind(i)+1):(ind(i+1)-1));
  25. b=find(a(1:end-1).*a(2:end)<0);
  26. c=[a(b) a(b+1)];
  27. d=[a1(b) a1(b+1)];
  28. for k=1:length(b)
  29. X_3(i,k)=interp1([c(k,1), c(k,2)],[d(k,1),d(k,2)],0);
  30. end
  31. else continue;
  32. end
  33. end
  34. X_30=reshape(X_3,1,[]); %矩阵转向量;
  35. X3=X_30(find(X_30)); %剔除系统自动补的零
  36. %最后一段,最后一个自身的零点 到 y序列末
  37. clear a b c d a1;
  38. a=y((ind(length(ind))+1):length(y));
  39. a1=x((ind(length(ind))+1):length(y));
  40. b=find(a(1:end-1).*a(2:end)<0);
  41. c=[a(b) a(b+1)];
  42. d=[a1(b) a1(b+1)];
  43. X=zeros(size(b));
  44. for k=1:length(b)
  45. X4(k)=interp1([c(k,1), c(k,2)],[d(k,1),d(k,2)],0);
  46. end
  47. X_1=cat(2,X1',X2);
  48. X_2=cat(2,X3,X4);
  49. X=cat(2,X_1,X_2);
  50. X=sort(X) %按升序排列 得到零点序列
复制代码

[ 本帖最后由 chenjc18 于 2009-1-5 21:36 编辑 ]

评分

2

查看全部评分

回复 支持 反对
分享到:

使用道具 举报

发表于 2009-1-5 21:43 | 显示全部楼层

回复 16楼 chenjc18 的帖子

楼主了解涵义了! 人不适, 明天再看看!
 楼主| 发表于 2009-1-5 23:02 | 显示全部楼层

可以处理任意离散序列的程序

我将16楼的做了改到 可以处理任何形式的离散序列:适用于:不含零点 含有一个零点 含有两个(不)连续零点 含有两个以上零点
请大家,看一下程序,有什么问题
谢谢了
我写的这段程序可能有点太乱了,嘿嘿,我已经尽力 欢迎大家批评指正,大家共同进步!

  1. clear;close;
  2. y=[1 2 -3  -1  0 1 -1  0 2 0 8 0 -1 -1  4  1 -4  7 -2];
  3. x=linspace(0,20,length(y));
  4. zero=zeros(length(x),1);
  5. x=x';
  6. y=y';
  7. plot(x,y,'b',x,zero,'r');
  8. ind=find(y==0); %找出自身零点所在下标
  9. if(isempty(ind)) %如果原离散序列不含有零,则进行程序一
  10. %%%%%%%%程序一%%%%%%%%%
  11.     b=find(y(1:end-1).*y(2:end)<0);
  12.     c=[y(b) y(b+1)];
  13.     d=[x(b) x(b+1)];
  14.     for k=1:length(b)
  15.         X(k)=interp1([c(k,1), c(k,2)],[d(k,1),d(k,2)],0);
  16.     end
  17.     X %输出x
  18. else  %否则执行程序二
  19. %%%%%%%%%程序二%%%%%%%%%%%  
  20. X1=x(ind);   %找到自身就有的零点对应的x
  21. %%段一 从段首到第一个自身的零点
  22. a=y(1:(ind(1)-1));
  23. b=find(a(1:end-1).*a(2:end)<0);
  24. c=[a(b) a(b+1)];
  25. d=[x(b) x(b+1)];
  26. X=zeros(size(b));
  27. for k=1:length(b)
  28. X2(k)=interp1([c(k,1), c(k,2)],[d(k,1),d(k,2)],0);
  29. end
  30. %中间的各段
  31. clear a b c d ;
  32. if((length(ind)==2)&(ind(2)-ind(1)==1)|(length(ind)==1))
  33.     X_3=[0];%如果原离散序列含有一个,或者连续两个零点,则执行此句
  34. else        %如果含有两个(不连续)及以上个零点 执行下面的语句
  35. for i=1:(length(ind)-1)
  36.     if((ind(i+1)-ind(i)~=1)|(ind(i+1)-ind(i)~=2))%如果出现相邻的或间隔为1的零点,则跳过
  37.         a=y((ind(i)+1):(ind(i+1)-1));            
  38.         a1=x((ind(i)+1):(ind(i+1)-1));
  39.         b=find(a(1:end-1).*a(2:end)<0);
  40.         c=[a(b) a(b+1)];
  41.         d=[a1(b) a1(b+1)];
  42.         for k=1:length(b)
  43.         X_3(i,k)=interp1([c(k,1), c(k,2)],[d(k,1),d(k,2)],0);
  44.         end
  45.     else continue;
  46.     end
  47. end
  48. end
  49. X_30=reshape(X_3,1,[]); %矩阵转向量;
  50. X3=X_30(find(X_30));    %剔除系统自行添加的零
  51. %%最后一段,最后一个自身的零点 到 y序列末
  52. clear a b c d a1;
  53. a=y((ind(length(ind))+1):length(y));
  54. a1=x((ind(length(ind))+1):length(y));
  55. b=find(a(1:end-1).*a(2:end)<0);
  56. c=[a(b) a(b+1)];
  57. d=[a1(b) a1(b+1)];
  58. for k=1:length(b)
  59. X4(k)=interp1([c(k,1), c(k,2)],[d(k,1),d(k,2)],0);
  60. end
  61. %%将各段计算得出零点放在一起
  62. X_1=cat(2,X1',X2);      
  63. X_2=cat(2,X3,X4);
  64. X=cat(2,X_1,X_2);
  65. X=sort(X)                %按升序排列
  66. end
复制代码

评分

1

查看全部评分

发表于 2009-1-5 23:11 | 显示全部楼层

回复 18楼 chenjc18 的帖子

LZ还真有研究的劲头,佩服!
发表于 2009-1-6 16:02 | 显示全部楼层
真服了楼主的态度! 欣赏!
还记得7F的回覆吗!? 其中说过interp1好像根本不需要!
18F的程序太长了, 逻辑思绪个人感觉太复杂了, 有点懒得看, 不好意思!
欣赏楼主, 列上我的, 不是最佳, 参考参考!
clc; clear;
y=[1 2 -3  -1  0 1 -1  0 2 0 8 0 -1 -1  4  1 -4  7 -2];
x=linspace(0,20,length(y)); x=x'; y=y'; xx=[]; plot(x,y); grid on;
for k=1:length(y)-1
   if y(k)==0, xx=[xx; x(k)];
   elseif y(k)*y(k+1)<0, xx=[xx; x(k)-y(k)*(x(k+1)-x(k))/(y(k+1)-y(k))]; end
end
if y(end)==0, xx=[xx; x(end)]; end
xx

其中x(k)-y(k)*(x(k+1)-x(k))/(y(k+1)-y(k))亦可用interp1([y(k),y(k+1)],[x(k),x(k+1)],0)

[ 本帖最后由 ChaChing 于 2009-1-6 16:09 编辑 ]

评分

1

查看全部评分

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

回复 20楼 ChaChing 的帖子

谢谢 ChaChing
您编写的程序真棒!
又长了很多知识,谢谢了
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

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

GMT+8, 2024-11-25 13:19 , Processed in 0.074121 second(s), 23 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

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