matlab小波分频,我对一个小波包分频之后进行样本熵计算matlab出现点问题

load D:\105.mat;%内圈故障数据

x=X105_DE_time;%驱动计数段的内圈故障

fs=12000;%采样率

N=2048;%采样点数(100倍)

M=0;%采样数据段的起始位置

n=M:N-1;

t=n/fs;%信号时间序列

X=X105_DE_time(1:N);%装载 驱动计数端的内圈故障数据

%X=X107_FE_time(1:N);%装载 风扇计数端的内圈故障数据

%X=X107_BA_time(1:N);%装载 基础计数端的内圈故障数据

%X=X107_DE_time(1:N)-X107_BA_time(1:N);

y=X';%信号幅值序列

k_in=kurtosis(y);%峭度系数,正常轴承为3左右

figure(1)

plot(t,y);

ylabel('振幅 m/v');

xlabel('时间 t/s');

t=wpdec(y,3,'db10');

figure(2)

plot(t);

s130=wprcoef(t,[3,0]);

s131=wprcoef(t,[3,1]);

s132=wprcoef(t,[3,2]);

s133=wprcoef(t,[3,3]);

s134=wprcoef(t,[3,4]);

s135=wprcoef(t,[3,5]);

s136=wprcoef(t,[3,6]);

s137=wprcoef(t,[3,7]);

figure(3)

subplot(4,1,1);plot(s130);

Ylabel('s130');

subplot(4,1,2);plot(s131);

Ylabel('s131');

subplot(4,1,3);plot(s132);

Ylabel('s132');

subplot(4,1,4);plot(s133);

Ylabel('s133');

figure(4)

subplot(4,1,1);plot(s134);

Ylabel('s134');

subplot(4,1,2);plot(s135);

Ylabel('s135');

subplot(4,1,3);plot(s136);

Ylabel('s136');

subplot(4,1,4);plot(s137);

Ylabel('s137');

s=hilbert(s136);

s1=abs(s);

s2=fft(s1,N);

s3=abs(s2);

f=n*fs/N;

figure(5)

plot(f,s3);

xlim([0 500]);

ylim([0 500]);

ylabel('功率谱 p/w');

ss=hilbert(s134);

ss1=abs(ss);

ss2=fft(ss1,N);

ss3=abs(ss2);

f=n*fs/N;

figure(6)

plot(f,ss3);

xlim([0 500]);

ylim([0 30]);

ylabel('功率谱 p/w');

for i=1:8

E(i)=norm(wpcoef(t,[3,i-1]),2);

end

E_max=max(E);

E_total=sum(E);

for i=1:8

x(i)= E(i)/E_total;

end

y=[x(1);x(2);x(3);x(4);x(5);x(6);x(7);x(8)]

figure(7)

bar(y,0.8);

ssss=SampEn(x(1), 2, 0.2*std(x(1)));

这是主程序,下面是样本熵程序

function SampEnVal = SampEn(data, m, r)

%SAMPEN  计算时间序列data的样本熵

%        data为输入数据序列

%        m为初始分段,每段的数据长度

%        r为阈值

% $Author: lskyp

% $Date:   2010.6.20

% Orig Version: V1.0--------分开计算长度为m的序列和长度为m+1的序列

%                           这一版的计算有些问题,需要注意两个序列总数都要为N-m

% Modi Version: V1.1--------综合计算,计算距离时通过矩阵减法完成,避免重循环

% V1.1 Modified date: 2010.6.23

data = data(:)';

N = length(data);

Nkx1 = 0;

Nkx2 = 0;

% 分段计算距离,x1为长度为m的序列,x2为长度为m+1的序列

for k = N - m:-1:1

x1(k, :) = data(k:k + m - 1);

x2(k, :) = data(k:k + m);

end

for k = N - m:-1:1

% x1序列计算

% 统计距离,由于每行都要与其他行做减法,因此可以先将该行复制为N-m的矩阵,然后

% 与原始x1矩阵做减法,可以避免两重循环,增加效率

x1temprow = x1(k, :);

x1temp    = ones(N - m, 1)*x1temprow;

% 可以使用repmat函数完成上面的语句,即

% x1temp = repmat(x1temprow, N - m, 1);

% 但是效率不如上面的矩阵乘法

% 计算距离,每一行元素相减的最大值为距离

dx1(k, :) = max(abs(x1temp - x1), [], 2)';

% 模板匹配数

Nkx1 = Nkx1 + (sum(dx1(k, :) < r) - 1)/(N - m - 1);

% x2序列计算,和x1同样方法

x2temprow = x2(k, :);

x2temp    = ones(N - m, 1)*x2temprow;

dx2(k, :) = max(abs(x2temp - x2), [], 2)';

Nkx2      = Nkx2 + (sum(dx2(k, :) < r) - 1)/(N - m - 1);

end

% 平均值

Bmx1 = Nkx1/(N - m);

Bmx2 = Nkx2/(N - m);

% 样本熵

SampEnVal = -log(Bmx2/Bmx1);运行结果为ssss=NaN

y =

0.0917

0.1425

0.2247

0.1089

0.0113

0.0335

0.2802

0.1074

求大神指点有什么毛病怎样改正

  • 0
    点赞
  • 3
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值