粤语残片 发表于 2014-3-11 08:50

小波熵图如何获取

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

Lorraine 发表于 2014-3-11 09:39

function =entropy(x,descriptor,approach,base)
%ENTROPY   Estimates the entropy of stationary signals with
%          independent samples using various approaches.
%    = ENTROPY(X) or
%    = ENTROPY(X,DESCRIPTOR) or
%    = ENTROPY(X,DESCRIPTOR,APPROACH) or
%    = ENTROPY(X,DESCRIPTOR,APPROACH,BASE)
%
%   ESTIMATE    : The entropy estimate
%   NBIAS       : The N-bias of the estimate
%   SIGMA       : The standard error of the estimate
%   DESCRIPTOR: The descriptor of the histogram, seel alse ENTROPY
%
%   X         : The time series to be analyzed, a row vector
%   DESCRIPTOR: Where DESCRIPTOR=
%   LOWERBOUND: Lowerbound of the histogram
%   UPPERBOUND: Upperbound of the histogram
%   NCELL   : The number of cells of the histogram      
%   APPROACH    : The method used, one of the following ones:
%   'unbiased': The unbiased estimate (default)
%   'mmse'    : The minimum mean square error estimate
%   'biased': The biased estimate
%   BASE      : The base of the logarithm; default e
%
%   See also: http://www.cs.rug.nl/~rudy/matlab/

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



if nargin <1
   disp('Usage: = ENTROPY(X)')
   disp('       = ENTROPY(X,DESCRIPTOR)')
   disp('       = ENTROPY(X,DESCRIPTOR,APPROACH)')
   disp('       = ENTROPY(X,DESCRIPTOR,APPROACH,BASE)')
   disp('Where: DESCRIPTOR = ')
   return
end

% Some initial tests on the input arguments

=size(x);

if NRowX~=1
error('Invalid dimension of X');
end;

if nargin>4
error('Too many arguments');
end;

if nargin==1
=histogram(x);
end;

if nargin>=2
=histogram(x,descriptor);
end;

if nargin<3
approach='unbiased';
end;

if nargin<4
base=exp(1);
end;

lowerbound=descriptor(1);
upperbound=descriptor(2);
ncell=descriptor(3);

estimate=0;
sigma=0;
count=0;

for n=1:ncell
if h(n)~=0
    logf=log(h(n));
else
    logf=0;
end;
count=count+h(n);
estimate=estimate-h(n)*logf;
sigma=sigma+h(n)*logf^2;
end;

% biased estimate

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

% conversion to unbiased estimate

if approach(1)=='u'
estimate=estimate-nbias;
nbias=0;
end;

% conversion to minimum mse estimate

if approach(1)=='m'
estimate=estimate-nbias;
nbias=0;
lambda=estimate^2/(estimate^2+sigma^2);
nbias   =(1-lambda)*estimate;
estimate=lambda*estimate;
sigma   =lambda*sigma;
end;

% base transformation

estimate=estimate/log(base);
nbias   =nbias   /log(base);
sigma   =sigma   /log(base);

粤语残片 发表于 2014-3-11 10:21

Lorraine 发表于 2014-3-11 09:39


你奏是活雷锋啊~{:3_53:}我自己先跑跑看,非常感谢啊!

粤语残片 发表于 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
=histogram(x);
请女侠赐教{:{05}:}

Lorraine 发表于 2014-3-12 10:54

缺少函数:
function =histogram(x,descriptor)
%HISTOGRAM Computes the frequency histogram of the row vector x.
%    = HISTOGRAM(X) or
%    = HISTOGRAM(X,DESCRIPTOR) or
%where
%   DESCRIPTOR =
%
%   RESULT    : A row vector containing the histogram
%   DESCRIPTOR: The used descriptor
%
%   X         : The row vector be analyzed
%   DESCRIPTOR: The descriptor of the histogram
%   LOWER   : The lowerbound of the histogram
%   UPPER   : The upperbound of the histogram
%   NCELL   : The number of cells of the histogram
%
%   See also: http://www.cs.rug.nl/~rudy/matlab/

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


if nargin <1
   disp('Usage: RESULT = HISTOGRAM(X)')
   disp('       RESULT = HISTOGRAM(X,DESCRIPTOR)')
   disp('Where: DESCRIPTOR = ')
   return
end

% Some initial tests on the input arguments

=size(x);

if NRowX~=1
error('Invalid dimension of X');
end;

if nargin>2
error('Too many arguments');
end;

if nargin==1
minx=min(x);
maxx=max(x);
delta=(maxx-minx)/(length(x)-1);
ncell=ceil(sqrt(length(x)));
descriptor=;
end;

lower=descriptor(1);
upper=descriptor(2);
ncell=descriptor(3);

if ncell<1
error('Invalid number of cells')
end;

if upper<=lower
error('Invalid bounds')
end;

result(1:ncell)=0;

y=round( (x-lower)/(upper-lower)*ncell + 1/2 );
for n=1:NColX
index=y(n);
if index >= 1 & index<=ncell
    result(index)=result(index)+1;
end;
end;

粤语残片 发表于 2014-3-19 08:37

本帖最后由 粤语残片 于 2014-3-19 08:42 编辑

Lorraine 发表于 2014-3-12 10:54
缺少函数:

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

Lorraine 发表于 2014-3-24 14:10

粤语残片 发表于 2014-3-19 08:37
再次感谢啊!不过我不懂编程,看不懂程序,倒是求出来每段数据的熵值了,只是想获得的是以数据点数为横 ...

不明白你的意思

粤语残片 发表于 2014-3-24 16:05

Lorraine 发表于 2014-3-24 14:10
不明白你的意思
就是绘制像这样的图

粤语残片 发表于 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();
%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);
感觉不对,谢谢了啊~{:3_49:}

韵天之色 发表于 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表示的是带宽和中心频率,其实说 ...

谢谢 能不能推荐几本关于小波熵的书 嘿 怎么提取实部奥

周文静 发表于 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

小波熵的介绍在胡昌华老师书籍里有
页: [1] 2
查看完整版本: 小波熵图如何获取