设为首页收藏本站联系我们
切换到宽版

声振论坛

 找回密码
 我要加入

QQ登录

只需一步，快速开始

声振论坛»论坛 › 工程技术 › 信号处理方法 › 小波熵图如何获取
12下一页
返回列表 发新帖
查看: 4117|回复: 15

[小波] 小波熵图如何获取

[复制链接]
粤语残片
发表于 2014-3-11 08:50 | 显示全部楼层 |阅读模式

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

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

x
本人利用小波变换得到一段信号的小波时频图，直观但是没有确切数字，想通过小波熵图表征信号的变化，求教大侠求小波熵图的程序，多谢~
一下是小波时频图， a.PNG
回复
分享到：

使用道具 举报

Lorraine
发表于 2014-3-11 09:39 | 显示全部楼层
  1. function [estimate,nbias,sigma,descriptor]=entropy(x,descriptor,approach,base)
  2. %ENTROPY   Estimates the entropy of stationary signals with
  3. %          independent samples using various approaches.
  4. %   [ESTIMATE,NBIAS,SIGMA,DESCRIPTOR] = ENTROPY(X) or
  5. %   [ESTIMATE,NBIAS,SIGMA,DESCRIPTOR] = ENTROPY(X,DESCRIPTOR) or
  6. %   [ESTIMATE,NBIAS,SIGMA,DESCRIPTOR] = ENTROPY(X,DESCRIPTOR,APPROACH) or
  7. %   [ESTIMATE,NBIAS,SIGMA,DESCRIPTOR] = ENTROPY(X,DESCRIPTOR,APPROACH,BASE)
  8. %
  9. %   ESTIMATE    : The entropy estimate
  10. %   NBIAS       : The N-bias of the estimate
  11. %   SIGMA       : The standard error of the estimate
  12. %   DESCRIPTOR  : The descriptor of the histogram, seel alse ENTROPY
  13. %
  14. %   X           : The time series to be analyzed, a row vector
  15. %   DESCRIPTOR  : Where DESCRIPTOR=[LOWERBOUND,UPPERBOUND,NCELL]
  16. %     LOWERBOUND: Lowerbound of the histogram
  17. %     UPPERBOUND: Upperbound of the histogram
  18. %     NCELL     : The number of cells of the histogram      
  19. %   APPROACH    : The method used, one of the following ones:
  20. %     'unbiased': The unbiased estimate (default)
  21. %     'mmse'    : The minimum mean square error estimate
  22. %     'biased'  : The biased estimate
  23. %   BASE        : The base of the logarithm; default e
  24. %
  25. %   See also: http://www.cs.rug.nl/~rudy/matlab/

  27. %   R. Moddemeijer
  28. %   Copyright (c) by R. Moddemeijer
  29. %   $Revision: 1.1 $  $Date: 2001/02/05 08:59:36 $



  33. if nargin <1
  34.    disp('Usage: [ESTIMATE,NBIAS,SIGMA,DESCRIPTOR] = ENTROPY(X)')
  35.    disp('       [ESTIMATE,NBIAS,SIGMA,DESCRIPTOR] = ENTROPY(X,DESCRIPTOR)')
  36.    disp('       [ESTIMATE,NBIAS,SIGMA,DESCRIPTOR] = ENTROPY(X,DESCRIPTOR,APPROACH)')
  37.    disp('       [ESTIMATE,NBIAS,SIGMA,DESCRIPTOR] = ENTROPY(X,DESCRIPTOR,APPROACH,BASE)')
  38.    disp('Where: DESCRIPTOR = [LOWERBOUND,UPPERBOUND,NCELL]')
  39.    return
  40. end

  42. % Some initial tests on the input arguments

  44. [NRowX,NColX]=size(x);

  46. if NRowX~=1
  47.   error('Invalid dimension of X');
  48. end;

  50. if nargin>4
  51.   error('Too many arguments');
  52. end;

  54. if nargin==1
  55.   [h,descriptor]=histogram(x);
  56. end;

  58. if nargin>=2
  59.   [h,descriptor]=histogram(x,descriptor);
  60. end;

  62. if nargin<3
  63.   approach='unbiased';
  64. end;

  66. if nargin<4
  67.   base=exp(1);
  68. end;

  70. lowerbound=descriptor(1);
  71. upperbound=descriptor(2);
  72. ncell=descriptor(3);

  74. estimate=0;
  75. sigma=0;
  76. count=0;

  78. for n=1:ncell
  79.   if h(n)~=0
  80.     logf=log(h(n));
  81.   else
  82.     logf=0;
  83.   end;
  84.   count=count+h(n);
  85.   estimate=estimate-h(n)*logf;
  86.   sigma=sigma+h(n)*logf^2;
  87. end;

  89. % biased estimate

  91. estimate=estimate/count;
  92. sigma   =sqrt( (sigma/count-estimate^2)/(count-1) );
  93. estimate=estimate+log(count)+log((upperbound-lowerbound)/ncell);
  94. nbias   =-(ncell-1)/(2*count);

  96. % conversion to unbiased estimate

  98. if approach(1)=='u'
  99.   estimate=estimate-nbias;
  100.   nbias=0;
  101. end;

  103. % conversion to minimum mse estimate

  105. if approach(1)=='m'
  106.   estimate=estimate-nbias;
  107.   nbias=0;
  108.   lambda=estimate^2/(estimate^2+sigma^2);
  109.   nbias   =(1-lambda)*estimate;
  110.   estimate=lambda*estimate;
  111.   sigma   =lambda*sigma;
  112. end;

  114. % base transformation

  116. estimate=estimate/log(base);
  117. nbias   =nbias   /log(base);
  118. sigma   =sigma   /log(base);
复制代码
回复 支持 1 反对 0

使用道具 举报

粤语残片
 楼主| 发表于 2014-3-11 10:21 | 显示全部楼层
Lorraine 发表于 2014-3-11 09:39

你奏是活雷锋啊~我自己先跑跑看，非常感谢啊！
回复 支持 反对

使用道具 举报

粤语残片
 楼主| 发表于 2014-3-11 14:53 | 显示全部楼层
Lorraine 发表于 2014-3-11 09:39

喔哦，出现问题解决不好了，我定义了entropy函数，然后输入值load‘X.txt’
entropy(X)
提示错误??? Undefined function or method 'histogram' for input arguments of type
'double'.

Error in ==> entropy at 55
  [h,descriptor]=histogram(x);
请女侠赐教
回复 支持 反对

使用道具 举报

Lorraine
发表于 2014-3-12 10:54 | 显示全部楼层
缺少函数：
  1. function [result,descriptor]=histogram(x,descriptor)
  2. %HISTOGRAM Computes the frequency histogram of the row vector x.
  3. %   [RESULT,DESCRIPTOR] = HISTOGRAM(X) or
  4. %   [RESULT,DESCRIPTOR] = HISTOGRAM(X,DESCRIPTOR) or
  5. %where
  6. %   DESCRIPTOR = [LOWER,UPPER,NCELL]
  7. %
  8. %   RESULT    : A row vector containing the histogram
  9. %   DESCRIPTOR: The used descriptor
  10. %
  11. %   X         : The row vector be analyzed
  12. %   DESCRIPTOR: The descriptor of the histogram
  13. %     LOWER   : The lowerbound of the histogram
  14. %     UPPER   : The upperbound of the histogram
  15. %     NCELL   : The number of cells of the histogram
  16. %
  17. %   See also: http://www.cs.rug.nl/~rudy/matlab/

  19. %   R. Moddemeijer
  20. %   Copyright (c) by R. Moddemeijer
  21. %   $Revision: 1.2 $  $Date: 2001/02/05 09:54:29 $


  24. if nargin <1
  25.    disp('Usage: RESULT = HISTOGRAM(X)')
  26.    disp('       RESULT = HISTOGRAM(X,DESCRIPTOR)')
  27.    disp('Where: DESCRIPTOR = [LOWER,UPPER,NCELL]')
  28.    return
  29. end

  31. % Some initial tests on the input arguments

  33. [NRowX,NColX]=size(x);

  35. if NRowX~=1
  36.   error('Invalid dimension of X');
  37. end;

  39. if nargin>2
  40.   error('Too many arguments');
  41. end;

  43. if nargin==1
  44.   minx=min(x);
  45.   maxx=max(x);
  46.   delta=(maxx-minx)/(length(x)-1);
  47.   ncell=ceil(sqrt(length(x)));
  48.   descriptor=[minx-delta/2,maxx+delta/2,ncell];
  49. end;

  51. lower=descriptor(1);
  52. upper=descriptor(2);
  53. ncell=descriptor(3);

  55. if ncell<1
  56.   error('Invalid number of cells')
  57. end;

  59. if upper<=lower
  60.   error('Invalid bounds')
  61. end;

  63. result(1:ncell)=0;

  65. y=round( (x-lower)/(upper-lower)*ncell + 1/2 );
  66. for n=1:NColX
  67.   index=y(n);
  68.   if index >= 1 & index<=ncell
  69.     result(index)=result(index)+1;
  70.   end;
  71. end;
复制代码
回复 支持 反对

使用道具 举报

粤语残片
 楼主| 发表于 2014-3-19 08:37 | 显示全部楼层
本帖最后由 粤语残片 于 2014-3-19 08:42 编辑
Lorraine 发表于 2014-3-12 10:54
缺少函数：


再次感谢啊！不过我不懂编程，看不懂程序，倒是求出来每段数据的熵值了，只是想获得的是以数据点数为横坐标，熵值为纵坐标的图，比如，一维数据有2000个数字，那横坐标就是0~2000，纵坐标为每个点对应的小波熵值，，或者是每n个点求一次熵值，刚接触到数据分析这块，不怎么了解，多有麻烦~
回复 支持 反对

使用道具 举报

Lorraine
发表于 2014-3-24 14:10 | 显示全部楼层
粤语残片 发表于 2014-3-19 08:37
再次感谢啊！不过我不懂编程，看不懂程序，倒是求出来每段数据的熵值了，只是想获得的是以数据点数为横 ...

不明白你的意思
回复 支持 反对

使用道具 举报

粤语残片
 楼主| 发表于 2014-3-24 16:05 | 显示全部楼层
Lorraine 发表于 2014-3-24 14:10
不明白你的意思

就是绘制像这样的图 求教.PNG
回复 支持 反对

使用道具 举报

粤语残片
 楼主| 发表于 2014-8-29 11:12 | 显示全部楼层
Lorraine 发表于 2014-3-24 14:10
不明白你的意思

您好，想请教您一个基本的小程序，就是小波熵尺度熵的算法，算法是这样的，求出一段数据多尺度小波系数矩阵，行表示数据个数，列表示尺度序列长度，然后求矩阵一列各小波系数的平方和，而该列各小波系数的平方除以平方和就是该小波系数的概率p,代入公式s=-p*log（p），然后将这一列的s求和即为该列小波熵，依次求出每列小波熵，绘制这段数据的熵图，希望麻烦您给编一下，我自己编的为clc;
clear all;
load 'f320010110yeweibodon.txt';
%t1=f320000929(:,1);
y1=f320010110yeweibodon(:,1);



t1=2:2:954;

%figure(1);
subplot(211);
plot(t1,y1);
%axis([0 1000,0 80]);
%set(gca, 'Fontname', 'Times newman', 'Fontsize', 12)
%xlabel('Time(s)');ylabel('Friction force(kN)');
wavename='shan2-3';
      totalscal=512;                    %尺度序列的长度，即scal的长度
      wcf=centfrq(wavename);            %小波的中心频率
      cparam=2*wcf*totalscal;           %为得到合适的尺度所求出的参数
      a=totalscal:-1:256;  
      Scales=cparam./a;                   %得到各个尺度，以使转换得到频率序列为等差序列
%%%%%%%%%%%%%计算小波熵%%%%%%%%%%%%%%%%%
WT=cwt(y1,Scales,wavename); %对数据进行连续小波变换
n=length(Scales);
h=length(y1);

for i=1:h
      
      for j=1:n;
          E(j)=0;
      E(j)=E(j)+abs(WT(j,i))^2;
      end
      %求第i个节点的范数平方，其实也就是平方和
end
E_total=sum(E);  %求总能量


for i=1:h
      for j=1:n;
      e(j)=abs(WT(j,i))^2;
      P(j)= e(j)/E_total  ;
      m(j)=-(P(j)*log(P(j)));
      end
      s(i)=sum(m)
      %求第i个节点的范数平方，其实也就是平方和
end


%以下计算小波熵，即－sum（pj*lnpj），
%disp('小波熵的值s(j)');
%for i=1:h
  %m(i)=-(P(i)*log(P(i)));  
%end
%S_wt=sum(m);
%S_wt
subplot(212);
plot(t1,s);
感觉不对，谢谢了啊~
回复 支持 反对

使用道具 举报

韵天之色
发表于 2014-11-1 21:13 | 显示全部楼层
额   shannon小波熵是不是指能量熵奥？？shan2-3这个是什么意思奥？？谢谢
回复 支持 反对

使用道具 举报

粤语残片
 楼主| 发表于 2014-11-13 17:16 | 显示全部楼层
韵天之色 发表于 2014-11-1 21:13
额   shannon小波熵是不是指能量熵奥？？shan2-3这个是什么意思奥？？谢谢

这个，我也是菜鸟，熵的类型是依据算法来定的，结合小波特性来使用，2-3表示的是带宽和中心频率，其实说具体了我也说不好
回复 支持 反对

使用道具 举报

韵天之色
发表于 2014-11-14 20:56 | 显示全部楼层
粤语残片 发表于 2014-11-13 17:16
这个，我也是菜鸟，熵的类型是依据算法来定的，结合小波特性来使用，2-3表示的是带宽和中心频率，其实说 ...

谢谢 能不能推荐几本关于小波熵的书 嘿 怎么提取实部奥
回复 支持 0 反对 1

使用道具 举报

周文静
发表于 2014-12-12 11:48 | 显示全部楼层
粤语残片 发表于 2014-8-29 11:12
您好，想请教您一个基本的小程序，就是小波熵尺度熵的算法，算法是这样的，求出一段数据多尺度小波系数矩 ...

for i=1:h    for j=1:n;
      e(j)=abs(WT(j,i))^2;
      P(j)= e(j)/E_total  ;
      m(j)=-(P(j)*log(P(j)));
      end
      s(i)=sum(m)
      %求第i个节点的范数平方，其实也就是平方和
end

楼主，你的尺度熵中循环中  m(j)=-(P(j)*log(P(j)))得到的一个点的在所有尺度上的能量熵，最后得到熵m(j)是各点在所有尺度上的熵和，不是一个尺度上所有点的熵和。我觉得是应该改成 for j=1:n;   for i=1:h  
。
不知道我说的对不对
回复 支持 反对

使用道具 举报

粤语残片
 楼主| 发表于 2014-12-15 13:53 | 显示全部楼层
周文静 发表于 2014-12-12 11:48
for i=1:h    for j=1:n;
      e(j)=abs(WT(j,i))^2;
      P(j)= e(j)/E_total  ;

恩恩，谢谢你~我也觉得算的有问题
回复 支持 反对

使用道具 举报

HHT
发表于 2014-12-27 21:22 | 显示全部楼层
小波熵的介绍在胡昌华老师书籍里有
回复 支持 反对

使用道具 举报

下一页 »
12下一页
返回列表 发新帖
高级模式
B Color Image Link Quote Code Smilies
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

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

GMT+8, 2024-11-5 20:41 , Processed in 0.075710 second(s), 21 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

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