小波包分解
m=load('300_30.txt');
N=length(m);
for i=1:N-1 ;
q(i,1)=m(i,1);
end;
d=q';
s1=d;
change=1000;
[c,l] = wavedec(d,3,'db4');
%提取小波分解后的低频系数
ca3=appcoef(c,l,'db4',3);
%提取各层小波分解后的高频系数
cd3=detcoef(c,l,3);
cd2=detcoef(c,l,2);
cd1=detcoef(c,l,1);
%对信号强制消噪
cdd3=zeros(1,length(cd3));%第三层高频系数cd3全置0
cdd2=zeros(1,length(cd2));%第二层高频系数cd2全置0
cdd1=zeros(1,length(cd1));%第一层高频系数cd1全置0
c1=[ca3,cdd3,cdd2,cdd1];%建立新的系数矩阵
s2=waverec(c1,l,'db4')%为新的分解结构
%[thr,sorh,keepapp]=ddencmp('den','wv',d);
%s2=wdencmp('gbl',c,l,'db4',4,thr,sorh,keepapp);
%subplot(413)
%plot(1:change,s2(1:change));
%title('默认软阈值消噪后信号')
figure(1)
subplot(9,2,1)
plot(1:change,s1(1:change))
title('原始信号')
ylabel('S1')
subplot(9,2,2)
plot(1:change,s2(1:change))
title('强制消噪后信号')
ylabel('S2')
wpt=wpdec(s1,3,'db1','shannon');
%plot(wpt);
%重构第三层个节点小波系数
s130=wprcoef(wpt,[3,0]);
s131=wprcoef(wpt,[3,1]);
s132=wprcoef(wpt,[3,2]);
s133=wprcoef(wpt,[3,3]);
s134=wprcoef(wpt,[3,4]);
s135=wprcoef(wpt,[3,5]);
s136=wprcoef(wpt,[3,6]);
s137=wprcoef(wpt,[3,7]);
%计算第三层个节点小波能量
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);
disp('正常信号的特征向量');
snorm1=[s10,s11,s12,s13,s14,s15,s16,s17];
std1=[st10,st11,st12,st13,st14,st15,st16,st17];
%显示三层个节点小波系数
subplot(9,2,3);plot(1:change,s130(1:change));
ylabel('S130');
subplot(9,2,5);plot(1:change,s131(1:change));
ylabel('S131');
subplot(9,2,7);plot(1:change,s132(1:change));
ylabel('S132');
subplot(9,2,9);plot(1:change,s133(1:change));
ylabel('S133');
subplot(9,2,11);plot(1:change,s134(1:change));
ylabel('S134');
subplot(9,2,13);plot(1:change,s135(1:change));
ylabel('S135');
subplot(9,2,15);plot(1:change,s136(1:change));
ylabel('S136');
subplot(9,2,17);plot(1:change,s137(1:change));
ylabel('S137');
wpt=wpdec(s2,3,'db1','shannon');
%plot(wpt);
s230=wprcoef(wpt,[3,0]);
s231=wprcoef(wpt,[3,1]);
s232=wprcoef(wpt,[3,2]);
s233=wprcoef(wpt,[3,3]);
s234=wprcoef(wpt,[3,4]);
s235=wprcoef(wpt,[3,5]);
s236=wprcoef(wpt,[3,6]);
s237=wprcoef(wpt,[3,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);
disp('故障信号的特征向量');
snorm2=[s20,s21,s22,s23,s24,s25,s26,s27];
std2=[st20,st21,st22,st23,st24,st25,st26,st27];
subplot(9,2,4);plot(1:change,s230(1:change));
ylabel('S230');
subplot(9,2,6);plot(1:change,s231(1:change));
ylabel('S231');
subplot(9,2,8);plot(1:change,s232(1:change));
ylabel('S232');
subplot(9,2,10);plot(1:change,s233(1:change));
ylabel('S233');
subplot(9,2,12);plot(1:change,s234(1:change));
ylabel('S234');
subplot(9,2,14);plot(1:change,s235(1:change));
ylabel('S235');
subplot(9,2,16);plot(1:change,s236(1:change));
ylabel('S236');
subplot(9,2,18);plot(1:change,s237(1:change));
ylabel('S237');
%fft
N=1024;
figure(2)
y1=fft(s1,N);
py1=abs(y1)
%py1=y1.*conj(y1)/N;
y2=fft(s2,N);
py2=abs(y2)
%py2=y2.*conj(y2)/N;
y130=fft(s130,N);
py130=y130.*conj(y130)/N;
y131=fft(s131,N);
py131=y131.*conj(y131)/N;
y132=fft(s132,N);
py132=y132.*conj(y132)/N;
y133=fft(s133,N);
py133=y133.*conj(y133)/N;
y134=fft(s134,N);
py134=y134.*conj(y134)/N;
y135=fft(s135,N);
py135=y135.*conj(y135)/N;
y136=fft(s136,N);
py136=y136.*conj(y136)/N;
y137=fft(s137,N);
py137=y137.*conj(y137)/N;
y230=fft(s230,N);
py230=y230.*conj(y230)/N;
y231=fft(s231,N);
py231=y231.*conj(y231)/N;
y232=fft(s232,N);
py232=y232.*conj(y232)/N;
y233=fft(s233,N);
py233=y233.*conj(y233)/N;
y234=fft(s234,N);
py234=y234.*conj(y234)/N;
y235=fft(s235,N);
py235=y235.*conj(y235)/N;
y236=fft(s236,N);
py236=y236.*conj(y236)/N;
y237=fft(s237,N);
py237=y237.*conj(y237)/N;
para=512
f=1000*(0:(para-1))/N;
subplot(1,2,1);
plot(f,py1(1:para));
ylabel('P1');
title('原始信号的功率谱')
subplot(1,2,2);
plot(f,py2(1:para));
ylabel('P2');
title('故障信号的功率谱')
figure
subplot(4,2,1);
plot(f,py130(1:para));
ylabel('P130');
title('S130的功率谱')
subplot(4,2,2);
plot(f,py131(1:para));
ylabel('P131');
title('S131的功率谱')
subplot(4,2,3);
plot(f,py132(1:para));
ylabel('P132');
subplot(4,2,4);
plot(f,py133(1:para));
ylabel('P133');
subplot(4,2,5);
plot(f,py134(1:para));
ylabel('P134');
subplot(4,2,6);
plot(f,py135(1:para));
ylabel('P135');
subplot(4,2,7);
plot(f,py136(1:para));
ylabel('P136');
subplot(4,2,8);
plot(f,py137(1:para));
ylabel('P137');
figure
subplot(4,2,1);
plot(f,py230(1:para));
ylabel('P230');
title('S230的功率谱')
subplot(4,2,2);
plot(f,py231(1:para));
ylabel('P231');
title('S231的功率谱')
subplot(4,2,3);
plot(f,py232(1:para));
ylabel('P232');
subplot(4,2,4);
plot(f,py233(1:para));
ylabel('P233');
subplot(4,2,5);
plot(f,py234(1:para));
ylabel('P234');
subplot(4,2,6);
plot(f,py235(1:para));
ylabel('P235');
subplot(4,2,7);
plot(f,py236(1:para));
ylabel('P236');
subplot(4,2,8);
plot(f,py237(1:para));
ylabel('P237');
figure
E=s10+s11+s12+s13+s14+s15+s16+s17
E0=s10/E;E1=s11/E;E2=s12/E;E3=s13/E;E4=s14/E;E5=s15/E;E6=s16/E;E7=s17/E;
EE=[E0,E1,E2,E3,E4,E5,E6,E7]
subplot(211)
bar(EE)
M=s20+s21+s22+s23+s24+s25+s26+s27
M0=s20/M;M1=s21/M;M2=s22/M;M3=s23/M;M4=s24/M;M5=s25/M;M6=s26/M;M7=s27/M;
MM=[M0,M1,M2,M3,M4,M5,M6,M7]
subplot(212)
bar(MM)
评论0