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
求大神指点有什么毛病怎样改正