clear all
clc
X= load('D:\管道数据\重复堵塞（空管）\dry_55mmBlockageBefore_40mmBlockageAfter_150mm-1');%正常信号
SampEn=zeros(7,4);
for i=1:4
zh=X.gdAv(1:4410,i);
data(i,:) =zh;
end
x=data(1,:);
Fs=44100;
wlen=1000; inc=500; % 给出帧长和帧移
win=hanning(wlen); % 给出海宁窗
N=length(x); % 信号长度
X=enframe(x,win,inc)'; % 分帧
fn=size(X,2); % 求出帧数
time=(0:N-1)/Fs; % 计算出信号的时间刻度
for i = 1:fn
u = X(:,i);% 取出一帧
zhen(:,i)=u;
end
%****************样本熵**********************************
for iii=1:7
x=zhen(:,iii);
N = length(x);
p= 0.15; % p 可以选择 在 0.1~0.25
r=p*std(x);
%% 求解m=2
m=2;
B2=0;
num=zeros(N-m+1);
c=zeros(N-m+1);
for i=1:N-m+1
for j=1:N-m+1
if j~=i
d(i,j)=max([abs(x(i+0)-x(j+0)) abs(x(i+1)-x(j+1))]);
end
end
end
for i=1:N-m+1
for j=1:N-m+1
if d(i,j)<r
num(i)=num(i)+1;
end
end
c(i)=num(i)/(N-m+1);
B2=B2+c(i);
end
B2=B2/(N-m+1);
%% 求解m=3
m=3;
B3=0;
num=zeros(N-m+1);
c=zeros(N-m+1);
for i=1:N-m+1
for j=1:N-m+1
if j~=i
d(i,j)=max([abs(x(i+0)-x(j+0)) abs(x(i+1)-x(j+1)) abs(x(i+2)-x(j+2))]);
end
end
end
for i=1:N-m+1
for j=1:N-m+1
if d(i,j)<r
num(i)=num(i)+1;
end
end
c(i)=num(i)/(N-m+1);
B3=B3+c(i);
end
B3=B3/(N-m+1);
%% 求得样本熵
SampEn=log(B2/B3);
YBS(iii,:)=SampEn;
end
x2=data(2,:);
Fs=44100;
wlen=1000; inc=500; % 给出帧长和帧移
win=hanning(wlen); % 给出海宁窗
N=length(x2); % 信号长度
X2=enframe(x2,win,inc)'; % 分帧
fn=size(X2,2); % 求出帧数
time=(0:N-1)/Fs; % 计算出信号的时间刻度
for i = 1:fn
u = X2(:,i);% 取出一帧
zhen2(:,i)=u;
end
%****************样本熵**********************************
for iii=1:7
x2=zhen2(:,iii);
N = length(x2);
p= 0.15; % p 可以选择 在 0.1~0.25
r=p*std(x2);
%% 求解m=2
m=2;
B2=0;
num=zeros(N-m+1);
c=zeros(N-m+1);
for i=1:N-m+1
for j=1:N-m+1
if j~=i
d(i,j)=max([abs(x(i+0)-x(j+0)) abs(x(i+1)-x(j+1))]);
end
end
end
for i=1:N-m+1
for j=1:N-m+1
if d(i,j)<r
num(i)=num(i)+1;
end
end
c(i)=num(i)/(N-m+1);
B2=B2+c(i);
end
B2=B2/(N-m+1);
%% 求解m=3
m=3;
B3=0;
num=zeros(N-m+1);
c=zeros(N-m+1);
for i=1:N-m+1
for j=1:N-m+1
if j~=i
d(i,j)=max([abs(x(i+0)-x(j+0)) abs(x(i+1)-x(j+1)) abs(x(i+2)-x(j+2))]);
end
end
end
for i=1:N-m+1
for j=1:N-m+1
if d(i,j)<r
num(i)=num(i)+1;
end
end
c(i)=num(i)/(N-m+1);
B3=B3+c(i);
end
B3=B3/(N-m+1);
%% 求得样本熵
SampEn=log(B2/B3);
YBS(iii,2)=SampEn;
end
x3=data(3,:);
Fs=44100;
wlen=1000; inc=500; % 给出帧长和帧移
win=hanning(wlen); % 给出海宁窗
N=length(x3); % 信号长度
X3=enframe(x3,win,inc)'; % 分帧
fn=size(X3,2); % 求出帧数
time=(0:N-1)/Fs; % 计算出信号的时间刻度
for i = 1:fn
u = X3(:,i);% 取出一帧
zhen3(:,i)=u;
end
%****************样本熵**********************************
for iii=1:7
x3=zhen3(:,iii);
N = length(x3);
p= 0.15; % p 可以选择 在 0.1~0.25
r=p*std(x3);
%% 求解m=2
m=2;
B2=0;
num=zeros(N-m+1);
c=zeros(N-m+1);
for i=1:N-m+1
for j=1:N-m+1
if j~=i
d(i,j)=max([abs(x(i+0)-x(j+0)) abs(x(i+1)-x(j+1))]);
end
end
end
for i=1:N-m+1
for j=1:N-m+1
if d(i,j)<r
num(i)=num(i)+1;
end
end
c(i)=num(i)/(N-m+1);
B2=B2+c(i);
end
B2=B2/(N-m+1);
%% 求解m=3
m=3;
B3=0;
num=zeros(N-m+1);
c=zeros(N-m+1);
for i=1:N-m+1
for j=1:N-m+1
if j~=i
d(i,j)=max([abs(x(i+0)-x(j+0)) abs(x(i+1)-x(j+1)) abs(x(i+2)-x(j+2))]);
end
end
end
for i=1:N-m+1
for j=1:N-m+1
if d(i,j)<r
num(i)=num(i)+1;
end
end
c(i)=num(i)/(N-m+1);
B3=B3+c(i);
end
B3=B3/(N-m+1);
%% 求得样本熵
SampEn=log(B2/B3);
YBS(iii,3)=SampEn;
end
x4=data(4,:);
Fs=44100;
wlen=1000; inc=500; % 给出帧长和帧移
win=hanning(wlen); % 给出海宁窗
N=length(x4); % 信号长度
X4=enframe(x4,win,inc)'; % 分帧
fn=size(X4,2); % 求出帧数
time=(0:N-1)/Fs; % 计算出信号的时间刻度
for i = 1:fn
u = X4(:,i);% 取出一帧
zhen4(:,i)=u;
end
%****************样本熵**********************************
for iii=1:7
x4=zhen4(:,iii);
N = length(x4);
p= 0.15; % p 可以选择 在 0.1~0.25
r=p*std(x4);
%% 求解m=2
m=2;
B2=0;
num=zeros(N-m+1);
c=zeros(N-m+1);
for i=1:N-m+1
for j=1:N-m+1
if j~=i
d(i,j)=max([abs(x(i+0)-x(j+0)) abs(x(i+1)-x(j+1))]);
end
end
end
for i=1:N-m+1
for j=1:N-m+1
if d(i,j)<r
num(i)=num(i)+1;
end
end
c(i)=num(i)/(N-m+1);
B2=B2+c(i);
end
B2=B2/(N-m+1);
%% 求解m=3
m=3;
B3=0;
num=zeros(N-m+1);
c=zeros(N-m+1);
for i=1:N-m+1
for j=1:N-m+1
if j~=i
d(i,j)=max([abs(x(i+0)-x(j+0)) abs(x(i+1)-x(j+1)) abs(x(i+2)-x(j+2))]);
end
end
end
for i=1:N-m+1
for j=1:N-m+1
if d(i,j)<r
num(i)=num(i)+1;
end
end
c(i)=num(i)/(N-m+1);
B3=B3+c(i);
end
B3=B3/(N-m+1);
%% 求得样本熵
SampEn=log(B2/B3);
YBS(iii,4)=SampEn;
end
