小波熵图如何获取
本人利用小波变换得到一段信号的小波时频图，直观但是没有确切数字，想通过小波熵图表征信号的变化，求教大侠求小波熵图的程序，多谢~一下是小波时频图，
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); Lorraine 发表于 2014-3-11 09:39
你奏是活雷锋啊~{:3_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}:} 缺少函数：
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:42 编辑
Lorraine 发表于 2014-3-12 10:54
缺少函数：
再次感谢啊！不过我不懂编程，看不懂程序，倒是求出来每段数据的熵值了，只是想获得的是以数据点数为横坐标，熵值为纵坐标的图，比如，一维数据有2000个数字，那横坐标就是0~2000，纵坐标为每个点对应的小波熵值，{:4_74:}，或者是每n个点求一次熵值，刚接触到数据分析这块，不怎么了解，多有麻烦~ 粤语残片 发表于 2014-3-19 08:37
再次感谢啊！不过我不懂编程，看不懂程序，倒是求出来每段数据的熵值了，只是想获得的是以数据点数为横 ...
不明白你的意思 Lorraine 发表于 2014-3-24 14:10
不明白你的意思
就是绘制像这样的图 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:} 额 shannon小波熵是不是指能量熵奥？？shan2-3这个是什么意思奥？？谢谢 韵天之色 发表于 2014-11-1 21:13
额 shannon小波熵是不是指能量熵奥？？shan2-3这个是什么意思奥？？谢谢
这个，我也是菜鸟，熵的类型是依据算法来定的，结合小波特性来使用，2-3表示的是带宽和中心频率，其实说具体了我也说不好 粤语残片 发表于 2014-11-13 17:16
这个，我也是菜鸟，熵的类型是依据算法来定的，结合小波特性来使用，2-3表示的是带宽和中心频率，其实说 ...
谢谢 能不能推荐几本关于小波熵的书 嘿 怎么提取实部奥 粤语残片 发表于 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-12 11:48
for i=1:h for j=1:n;
e(j)=abs(WT(j,i))^2;
P(j)= e(j)/E_total;
恩恩，谢谢你~我也觉得算的有问题 小波熵的介绍在胡昌华老师书籍里有
页:
[1]
2