小波包分解

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%小波包分解程序%%m=load('300_30.txt');导入文件名为300_30的txt文件N=length(m);fori=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全置0cdd2=zeros(1,length(cd2));%第二层高频系数cd2全置0cdd1=zeros(1,length(cd1));%第一层高频系数cd1全置0c1=[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');%fftN=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('故障信号的功率谱')figuresubplot(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');figuresubplot(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');figureE=s10+s11+s12+s13+s14+s15+s16+s17E0=s10/E;E1=s11/E;E2=s12/E;E3=s13/E;E4=s14/E;E5=s15/E;E6=s16/E;E7=sE;EE=[E0,E1,E2,E3,E4,E5,E6,E7]subplot(211)bar(EE)M=s20+s21+s22+s23+s24+s25+s26+s27M0=s20/M;M1=