guzhangfenxi.m

matlab中结合小波和偏度对故障信号的处理

%使用小波包变换分析信号
%使用平台：Matlab7.0
%作者：张毅，   北 京 石 油 化 工 学 院
%加axis([0,1600,-2.5,2.5]);

function wavelet
[FileName,PathName] = uigetfile('2(2006_4_21_21_7_52_98).txt','Select the TXT-file');
fid = fopen(FileName, 'r');
[A,COUNT]=fscanf(fid,'%f');
fclose(fid);
% B2=A(176000:2:192000); %正常部分
% B1=A(182000:2:198000); %含故障部分   %(开始：每隔几个取：结束点)

C1=A(24000:2:34000); %含故障部分   %(开始：每隔几个取：结束点)
C2=A(12000:2:22000); %正常部分

D1=C1;
D2=C2;

max1=max(C1);  %发生故障时信号的最大值
max2=max(C2);   %正常信号的最大值
DianShu=10000;  % DianShu 定义为每次一次性读取多少个点
multiple=max1/max2;%  7.916667;   % multiple 定义为放大倍数
threshold =max2+(max1-max2)*0.5; % threshold  定义为放大阈值


%m()=0;
x=1;
for i=1:DianShu/2
    if C1(i)>threshold 
        m(x)=i;
        x=x+1;
        
    end
    if C2(i)<threshold 
        C2(i)=multiple* C2(i);
    end
end

for i=1:m(1)-1000
     C1(i)=multiple* C1(i);
end

B2=C2;
B1=C1;

%维纳滤波器的实现 

c=1;
a=1;
n=1:DianShu/2;
w=0.1*wgn(1,length(B2),2); 
v=0.3*wgn(1,length(B2),2); 
Q=std(w).^2; 
R=std(v).^2; 
P2=0.0358
G=c*P2/(R+P2*c^2); 
f=R*a/(R+P2*c^2); 
[R,P2,K]=residuez(G,[1,-f]); 
h2=R*P2.^n; 
%用维纳滤波器对正常信号进行滤波 
z2=conv(B2,h2); 
s2=1.2*z2(2:DianShu/2);


w1=0.1*wgn(1,length(B1),2); 
v1=0.3*wgn(1,length(B1),2); 
Q1=std(w1).^2; 
R1=std(v1).^2; 
P1=0.0358
G1=c*P1/(R1+P1*c^2); 
f1=R1*a/(R1+P1*c^2); 
[R1,P1,K]=residuez(G1,[1,-f1]); 
h1=R1*P1.^n; 
%用维纳滤波器对含故障信号进行滤波 
z1=conv(B1,h1); 
s1=1.2*z1(2:DianShu/2); 
 

figure(1);
subplot(2,2,1);
plot(D2);
title('原始正常信号');
axis([0,DianShu/2,-0.5,0.5,]);
ylabel('D2');
subplot(2,2,2);
plot(D1);
title('原始的含故障的信号');
axis([0,DianShu/2,-1,1,]);
ylabel('D1');
subplot(2,2,3)
plot(s2);
title('幅值加大后的正常信号');
ylabel('s2');
subplot(2,2,4)
plot(s1);
title('幅值加大后的含故障的信号');
ylabel('s1');

%用db10小波包对含故障信号进行8层分解
wpt=wpdec(s1,8,'db6','shannon');%%wpdec是单尺度二维离散小波分解函数

%对信号第三层各系数进行重构
%s130表示S1的[10，0]节点的重构系数，其余依次类推
s130=wprcoef(wpt,[8,0]);           
s131=wprcoef(wpt,[8,1]);
s132=wprcoef(wpt,[8,2]);
s133=wprcoef(wpt,[8,3]);
s134=wprcoef(wpt,[8,4]);
s135=wprcoef(wpt,[8,5]);
s136=wprcoef(wpt,[8,6]);
s137=wprcoef(wpt,[8,7]);

%计算各重构系数的方差  离均差平方和除以自由度称为方差
%S10表示S130的方差，其余类推
s10=norm(s130);
s11=norm(s131);
s12=norm(s132);
s13=norm(s133);
s14=norm(s134);
s15=norm(s135);
s16=norm(s136);
s17=norm(s137);

%计算各系数的标准方差  %方差开方得标准方差
st10=std(s130);
st11=std(s131);
st12=std(s132);
st13=std(s133);
st14=std(s134);
st15=std(s135);
st16=std(s136);
st17=std(s137);

% 偏度：
% 
% 偏度（Skewness）是描述某变量取值分布对称性的统计量。
% 
% Skewness=0     分布形态与正态分布偏度相同
% Skewness>0     正偏差数值较大，为正偏或右偏。长尾巴拖在右边。
% Skewness<0     负偏差数值较大，为负偏或左偏。长尾巴拖在左边。
% | Skewness| 越大，分布形态偏移程度越大。
% 偏斜程度也可用标准差的三次方除三阶中心矩的方法计算
% 偏度和峭度分别反映信号概率密度函数的不对称程度和陡峭程度

%计算各系数的偏度
sp10=skewness(s130);
sp11=skewness(s131);
sp12=skewness(s132);
sp13=skewness(s133);
sp14=skewness(s134);
sp15=skewness(s135);
sp16=skewness(s136);
sp17=skewness(s137);


% 峰度（Kurtosis）是描述某变量所有取值分布形态陡缓程度的统计量。
% 它是和正态分布相比较的。
% 
% Kurtosis=0       与正态分布的陡缓程度相同。
% Kurtosis>0       比正态分布的高峰更加陡峭——尖顶峰
% Kurtosis<0       比正态分布的高峰来得平台——平顶峰


%计算各系数的峭度
sq10=kurtosis(s130);
sq11=kurtosis(s131);
sq12=kurtosis(s132);
sq13=kurtosis(s133);
sq14=kurtosis(s134);
sq15=kurtosis(s135);
sq16=kurtosis(s136);
sq17=kurtosis(s137);
spiandu1=[sp10,sp11,sp12,sp13,sp14,sp15,sp16,sp17];
sqiaodu1=[sq10,sq11,sq12,sq13,sq14,sq15,sq16,sq17,];


% figure(3);
% subplot(1,2,1);
% bar(spiandu1);
% subplot(1,2,2);
% bar(sqiaodu1);
%提取信号的特征向量
disp('含故障的信号的特征向量');

snorm1=[s10,s11,s12,s13,s14,s15,s16,s17];
std1=[st10,st11,st12,st13,st14,st15,st16,st17];

%将信号的特征向量保存到文件中
fid=fopen('含故障的信号的特征向量.txt', 'w');
fprintf(fid,'%f\r\n',snorm1);
fclose(fid);


%求正常信号的特征向量
wpt=wpdec(s2,8,'db10','shannon');
%plot(wpt);
s230=wprcoef(wpt,[8,0]);
s231=wprcoef(wpt,[8,1]);
s232=wprcoef(wpt,[8,2]);
s233=wprcoef(wpt,[8,3]);
s234=wprcoef(wpt,[8,4]);
s235=wprcoef(wpt,[8,5]);
s236=wprcoef(wpt,[8,6]);
s237=wprcoef(wpt,[8,7]);

s20=norm(s230);
s21=norm(s231);
s22=norm(s232);
s23=norm(s233);
s24=norm(s234);
s25=norm(s235);
s26=norm(s236);
s27=norm(s237);

st20=std(s230);
st21=std(s231);
st22=std(s232);
st23=std(s233);
st24=std(s234);
st25=std(s235);
st26=std(s236);
st27=std(s237);


%计算各系数的偏度
sp20=skewness(s230);
sp21=skewness(s231);
sp22=skewness(s232);
sp23=skewness(s233);
sp24=skewness(s234);
sp25=skewness(s235);
sp26=skewness(s236);
sp27=skewness(s237);

%计算各系数的峭度
sq20=kurtosis(s230);
sq21=kurtosis(s231);
sq22=kurtosis(s232);
sq23=kurtosis(s233);
sq24=kurtosis(s234);
sq25=kurtosis(s235);
sq26=kurtosis(s236);
sq27=kurtosis(s237);
spiandu2=[sp20,sp21,sp22,sp23,sp24,sp25,sp26,sp27];
sqiaodu2=[sq20,sq21,sq22,sq23,sq24,sq25,sq26,sq27,];

figure(2);
subplot(2,2,1);
bar(spiandu1);
title('含故障信号的偏度图');
subplot(2,2,2);
bar(sqiaodu1);
title('含故障信号的峭度图');
subplot(2,2,3);
bar(spiandu2);
title('正常信号的偏度图');
subplot(2,2,4);
bar(sqiaodu2);
title('正常信号的峭度图');


disp('正常信号的特征向量');
snorm2=[s20,s21,s22,s23,s24,s25,s26,s27];
std2=[st20,st21,st22,st23,st24,st25,st26,st27];

fid=fopen('正常信号的特征向量.txt', 'w');
fprintf(fid,'%f\r\n',snorm2);
fclose(fid);

% figure(3);
% subplot(1,2,1);
% bar(snorm2+std2);
% %axis([0,8,0,2]);
% title('正常信号的特征向量图');
% subplot(122);
% bar(snorm1+std1);
% %axis([0,8,0,20]);
% title('含故障的信号的特征向量图');



%求含故障信号的功率普
y1=fft(s1,1024);
py1=y1.*conj(y1)/1024;
y130=fft(s130,1024);
py130=y130.*conj(y130)/1024;
y131=fft(s131,1024);
py131=y131.*conj(y131)/1024;
y132=fft(s132,1024);
py132=y132.*conj(y132)/1024;
y133=fft(s133,1024);
py133=y133.*conj(y133)/1024;
y134=fft(s134,1024);
py134=y134.*conj(y134)/1024;
y135=fft(s135,1024);
py135=y135.*conj(y135)/1024;
y136=fft(s136,1024);
py136=y136.*conj(y136)/1024;
y137=fft(s137,1024);
py137=y137.*conj(y137)/1024;



%求正常信号的功率普
y2=fft(s2,1024);
py2=y2.*conj(y2)/1024;
y230=fft(s230,1024);
py230=y230.*conj(y230)/1024;
y231=fft(s231,1024);
py231=y231.*conj(y231)/1024;
y232=fft(s232,1024);
py232=y232.*conj(y232)/1024;
y233=fft(s233,1024);
py233=y233.*conj(y233)/1024;
y234=fft(s234,1024);
py234=y234.*conj(y234)/1024;
y235=fft(s235,1024);
py235=y235.*conj(y235)/1024;
y236=fft(s236,1024);
py236=y236.*conj(y236)/1024;
y237=fft(s237,1024);
py237=y237.*conj(y237)/1024;

% figure(4);
% f=1000*(0:511)/1024;
% subplot(1,2,1);
% plot(f,py2(1:512));
% %axis([0,0.35,0,200]);
% %axis([0,200,0,0.35,]);
% ylabel('P2');
% title('正常信号的功率谱');
% 
% subplot(1,2,2);
% plot(f,py1(1:512));
% %axis([0,9,0,200]);
% %axis([0,100,0,3,]);
% ylabel('P1');
% title('含故障的信号的功率谱');
%