声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 3932|回复: 13

[综合讨论] 有人用过陆振波博士求时间序列关联维数和熵的程序么?

[复制链接]
发表于 2012-12-3 10:29 | 显示全部楼层 |阅读模式

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

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

x
本帖最后由 牛小贱 于 2014-3-7 20:36 编辑

RT,我最近做时间序列的分析,在网上找的陆振波博士的计算时间序列关联维数和熵的程序,对时间序列进行HHT滤波处理之后,导入陆博士的程序,有些得到的是nan和0.00000,这些说明什么情况呢?有人做过吗?求教我一下

下面的是主程序:
  1. clc
  2. clear all
  3. close all
  4. x=importdata('D:\Documents\Desktop\新的要求\10-30%-5.mat'); % 时间序列(列向量)
  5. x=x';
  6. %-----------------------------------  ---------------------------------------
  7. % G-P算法计算关联维
  8. k=length(x);
  9. rr = 0.5;
  10. Log2R = -6:rr:0;        % log2(r)
  11. R = 2.^(Log2R);
  12. fs =100;               % 信号采样频率
  13. t = 10;                 % 时延
  14. dd = 4;                 % 嵌入维间隔
  15. D = 4:dd:72;            % 嵌入维   
  16. p = 50;                 % 限制短暂分离,大于序列平均周期(不考虑该因素时 p = 1)  

  17. tic
  18. Log2Cr = log2(CorrelationIntegral(x,t,D,R,p));   % 输出每一行对应一个嵌入维
  19. toc
  20. %--------------------------------------------------------------------------
  21. % 结果作图
  22. figure
  23. plot(Log2R,Log2Cr','k.-'); axis tight; hold on; grid on;
  24. xlabel('log2(r)');
  25. ylabel('log2(C(r))');
  26. title(['signal, length = ',num2str(k)]);
  27. %--------------------------------------------------------------------------
  28. % 最小二乘拟合
  29. Linear = [6:10];                            % 线性似合区域
  30. [a,B] = LM2(Log2R,Log2Cr,Linear);           % 最小二乘求斜率和截距

  31. disp(sprintf('Correlation Dimension = %.4f',a));
  32. for i = 1:length(D)
  33.     Y = polyval([a,B(i)],Log2R(Linear),1);
  34.     plot(Log2R(Linear),Y,'r');
  35. end
  36. hold off;
  37. %--------------------------------------------------------------------------
  38. % 求梯度

  39. Slope = diff(Log2Cr,1,2)/rr;                % 求梯度
  40. xSlope = Log2R(1:end-1);                    % 梯度所对应的log2(r)

  41. figure;
  42. plot(xSlope,Slope','k.-'); axis tight; grid on;
  43. xlabel('log2(r)');
  44. ylabel('Slope');
  45. title(['signal, length = ',num2str(k)]);

  46. %--------------------------------------------------------------------------
  47. % 差分求K熵

  48. KE = -diff(B)/(dd*t)*fs*log(2);             % 用采样频率 fs 和公式 log(x) = log2(x)*log(2) 将单位转化成 nats/s
  49. D_KE = D(1:end-1);                          % K熵所对应的嵌入维

  50. figure;
  51. plot(D_KE,KE,'k.-'); grid on; hold on;
  52. xlabel('m');
  53. ylabel('Kolmogorov Entropy (nats/s)');
  54. title(['signal, length = ',num2str(k)]);

  55. %--------------------------------------------------------------------------
  56. % 输出显示

  57. disp(sprintf('Kolmogorov Entropy = %.4f',min(KE)));
复制代码
处理一组数据后得到的图形如下:


回复
分享到:

使用道具 举报

发表于 2012-12-3 14:53 | 显示全部楼层
本帖最后由 happy 于 2012-12-3 14:56 编辑

没有数据,没有代码,没有说明输出是什么,这个问题无从回答
建议看看http://forum.vibunion.com/thread-36219-1-1.html
把问题说明清楚
发表于 2012-12-3 18:54 | 显示全部楼层
我试过之后感觉出来的Kolmogorov Entropy图不对,结果并不收敛!!!
发表于 2012-12-4 09:40 | 显示全部楼层
发表于 2012-12-5 20:51 | 显示全部楼层
happy 发表于 2012-12-4 09:40
你的时间序列是仿真结果还是实验结果?

仿真结果,就是Lorenz模型,Kolmogorov Entropy随着m值变化,结果不收敛,图像是先下降然后居然又上升。还有那个程序有使用次数限制吗?
发表于 2012-12-6 10:58 | 显示全部楼层
章紫风 发表于 2012-12-5 20:51
仿真结果,就是Lorenz模型,Kolmogorov Entropy随着m值变化,结果不收敛,图像是先下降然后居然又上升。还 ...

1. 检查一下仿真得到的时间序列是否已经提出了瞬态响应数据
2. 没听说过他的工具包有这方面的限制
发表于 2012-12-6 18:32 | 显示全部楼层
happy 发表于 2012-12-6 10:58
1. 检查一下仿真得到的时间序列是否已经提出了瞬态响应数据
2. 没听说过他的工具包有这方面的限制

非常感谢,我的问题差不多解决了!!!
发表于 2012-12-11 11:22 | 显示全部楼层
章紫风 发表于 2012-12-6 18:32
非常感谢,我的问题差不多解决了!!!

建议解决问题后把相关的情况和大家分享一下
 楼主| 发表于 2012-12-27 15:22 | 显示全部楼层
章紫风 发表于 2012-12-6 18:32
非常感谢,我的问题差不多解决了!!!

同求解决,不知道能否交流一下?
发表于 2013-5-5 22:22 | 显示全部楼层
我想请问相关的问题
头像被屏蔽
发表于 2014-2-21 19:50 | 显示全部楼层
提示: 该帖被管理员或版主屏蔽
发表于 2014-2-21 20:20 | 显示全部楼层
章紫风 发表于 2012-12-6 18:32
非常感谢,我的问题差不多解决了!!!

希望你分享一下你的经验,方便大家学习!谢了
发表于 2014-5-12 15:31 | 显示全部楼层
你好,我是中国矿业大学的研究生,我在网上看到你发的一篇帖子,我对里面的内容很感兴趣。不知道你现在还有没有涉及柯尔莫戈洛夫熵值的领域。我想和你就我们共同感兴趣的地方交流一下。我的QQ是652247469.如果你方便的话请加我QQ好友。谢谢!
发表于 2014-6-3 09:52 | 显示全部楼层
31行程序的LM2是什么?子函数么?matlab没这个命令。如果是子函数,那是求什么的?子函数具体什么样?求解答!
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

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

GMT+8, 2024-12-23 20:15 , Processed in 0.076972 second(s), 18 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

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