小波包分解与重构能量熵特征提取MATLAB代码

%% 该方法可以实现对任意节点系数选择进行组合重构。 x=faultsignal1(:,12); t=wpdec(x,3,'db3','shannon');%将wave 用 shannondb3小波进行3层小波包分解，获得一个小波包树 t t2 = wpjoin(t,[3;4;5;6]);%将小波包树的第二行的四个节点收起来，也就是让第二行的节点变为树的最底层节点。因为第一行中小波包树的节点个数是第一层2个，第二层4个，第三层8个。现在t2就是将第三层的节点再聚合回第二层。 sNod = read(t,'sizes',[3,4,5,6]);%读取第二层四个节点系数的size cfs3 = zeros(sNod(1,:));%将所有四个节点的小波包系数变为0 cfs4 = zeros(sNod(2,:));%将四个节点的系数重组到t3小波树中。 cfs5 = zeros(sNod(3,:)); cfs6 = zeros(sNod(4,:)); t3 = write(t2,'cfs',3,cfs3,'cfs',4,cfs4,'cfs',5,cfs5,'cfs',6,cfs6);%对t3小波树进行重构，获得信号wave2 wave2=wprec(t3); %% t=wpdec(x,3,'dmey'); t2 = wpjoin(t,[3;4;5;6]);%将小波包树的第二行的四个节点收起来，也就是让第二行的节点变为树的最底层节点。 cfs3=wpcoef(t,3);%获取节点小波包系数，小波包系数就是一个一维向量 cfs4=wpcoef(t,4); cfs5=wpcoef(t,5); cfs6=wpcoef(t,6); t3 = write(t2,'cfs',3,cfs3,'cfs',4,cfs4,'cfs',5,cfs5,'cfs',6,cfs6); wave2=wprec(t3); %注意，wpjoin命令在这里是必要的， %因为write函数只能将最底层的节点写进去。 %也就是说，如果我们将第三层的小波包系数进行修改的话，就不用wpjoin了，具体可以看下面 %% t=wpdec(x,3,'dmey'); cfs7=wpcoef(t,7); cfs8=wpcoef(t,8); cfs9=wpcoef(t,9); cfs10=wpcoef(t,10); cfs11=wpcoef(t,11); cfs12=wpcoef(t,12); cfs13=wpcoef(t,13); cfs14=wpcoef(t,14); t3=write(t,'cfs',7,cfs7,'cfs',8,cfs8,'cfs',9,cfs9,'cfs',10,cfs10,'cfs',11,cfs11,'cfs',... 12,cfs12,'cfs',13,cfs13,'cfs',14,cfs14); y=wprec(t3); %% 该方法只能对某一节点信号系数分别进行重构，不能实现多个节点系数组合进行重构. % x_input=x ; %输入数据 plot(x_input);title('输入信号时域图像') %绘制输入信号时域图像 %% 查看频谱范围 x=x_input; fs=20000;%采样频率要大于等于原信号中最高频率的二倍 N=length(x); %采样点个数 signalFFT=abs(fft(x,N)); %显示真实的幅值 Y=2*signalFFT/N; f=(0:N/2)*(fs/N); figure; plot(f,Y(1:N/2+1)); ylabel('amp'); xlabel('frequency');title('输入信号的频谱');grid on %% wpt=wpdec(x_input,3,'db1','shannon'); %进行3层小波包分解 plot(wpt); %绘制小波包树 %% %接下来，我们查看第3层8个节点的频谱分布（假设我的输入信号采样频率是128，按照采样定理小波包分解根节点（0，0）处的频率应该为0-64Hz， %按照这个计算（3，0）节点为0-9Hz，后面依次以8Hz为一个段递增） for i=0:7 rex3(:,i+1)=wprcoef(wpt,[3 i]); %实现对节点小波节点进行重构 end figure; %绘制第3层各个节点分别重构后信号的频谱 for i=0:7 subplot(2,4,i+1); x_sign=rex3(:,i+1); N=length(x_sign); %采样点个数 signalFFT=abs(fft(x_sign,N));%真实的幅值 Y=2*signalFFT/N; f=(0:N/2)*(fs/N); plot(f,Y(1:N/2+1)); ylabel('amp'); xlabel('frequency');grid on axis([0 10000 0 100]); title(['小波包第3层',num2str(i),'节点信号频谱']); end %% %需要对节点顺序进行重新编码排序，按照频率递增 nodes=[7:14]'; %第3层的节点号 ord=wpfrqord(nodes); %小波包系数重排，ord是重排后小波包系数索引构成的矩阵　如3层分解的[1;2;4;3;7;8;6;5] nodes_ord=nodes(ord); %重排后的小波系数 %% for i=1:8 rex3(:,i)=wprcoef(wpt,nodes_ord(i)); %实现对节点小波节点进行重构 end figure; %绘制第3层各个节点分别重构后信号的频谱 for i=0:7 subplot(2,4,i+1); x_sign= rex3(:,i+1); N=length(x_sign); %采样点个数 signalFFT=abs(fft(x_sign,N));%真实的幅值 Y=2*signalFFT/N; f=(0:N/2)*(fs/N); plot(f,Y(1:N/2+1)); ylabel('幅值'); xlabel('频率');grid on axis([0 50 0 100]); title(['小波包第3层',num2str(i),'节点信号频谱']); end %% %小波包分解------能量特征提取（方法1） %% wavelet packet coefficients. 求取小波包分解的各个节点的小波包系数 cfs3_0=wpcoef(wpt,nodes_ord(1)); %对重排序后第3层1节点的小波包系数0-8Hz cfs3_1=wpcoef(wpt,nodes_ord(2)); %对重排序后第3层2节点的小波包系数8-16Hz cfs3_2=wpcoef(wpt,nodes_ord(3)); %对重排序后第3层3节点的小波包系数16-24Hz cfs3_3=wpcoef(wpt,nodes_ord(4)); %对重排序后第3层4节点的小波包系数24-32Hz cfs3_4=wpcoef(wpt,nodes_ord(5)); %对重排序后第3层5节点的小波包系数32-40Hz cfs3_5=wpcoef(wpt,nodes_ord(6)); %对重排序后第3层6节点的小波包系数40-48Hz cfs3_6=wpcoef(wpt,nodes_ord(7)); %对重排序后第3层7节点的小波包系数48-56Hz cfs3_7=wpcoef(wpt,nodes_ord(8)); %对重排序后第3层8节点的小波包系数56-64Hz fori=0:7 cfs3_i=wpcoef(wpt,nodes_ord(1)); E_cfs3_0=norm(cfs3_0,2)^2; %% 1-范数：就是norm(...,1)，即各元素绝对值之和；2-范数：就是norm(...,2)，即各元素平方和开根号； E_cfs3_1=norm(cfs3_1,2)^2; E_cfs3_2=norm(cfs3_2,2)^2; E_cfs3_3=norm(cfs3_3,2)^2; E_cfs3_4=norm(cfs3_4,2)^2; E_cfs3_5=norm(cfs3_5,2)^2; E_cfs3_6=norm(cfs3_6,2)^2; E_cfs3_7=norm(cfs3_7,2)^2; E_total=E_cfs3_0+E_cfs3_1+E_cfs3_2+E_cfs3_3+E_cfs3_4+E_cfs3_5+E_cfs3_6+E_cfs3_7; p_node(1)= 100*E_cfs3_0/E_total; % 求得每个节点的占比 p_node(2)= 100*E_cfs3_1/E_total; % 求得每个节点的占比 p_node(3)= 100*E_cfs3_2/E_total; % 求得每个节点的占比 p_node(4)= 100*E_cfs3_3/E_total; % 求得每个节点的占比 p_node(5)= 100*E_cfs3_4/E_total; % 求得每个节点的占比 p_node(6)= 100*E_cfs3_5/E_total; % 求得每个节点的占比 p_node(7)= 100*E_cfs3_6/E_total; % 求得每个节点的占比 p_node(8)= 100*E_cfs3_7/E_total; % 求得每个节点的占比 %求小波包熵 disp('小波熵的值S_wt'); % dim=length( E ); % for i=1:dim for i=1:8 m(i)=p_node(i)*log(p_node(i)) end S_wt=sum(m)*(-1) ; figure; x=1:8; bar(x,p_node); title('各个频段能量所占的比例'); xlabel('频率 Hz'); ylabel('能量百分比/%'); for j=1:8 text(x(j),p_node(j),num2str(p_node(j),'%0.2f'),... 'HorizontalAlignment','center',... 'VerticalAlignment','bottom') end %% %小波包分解------能量特征提取（方法2） %直接运行matlab自带函数，如下 E = wenergy(wpt); %该函数只能对最后一层（底层）节点进行能量提取 %这里对比一下，和方法1得到的图形基本上是一致的，不同之处在于排列顺序变了， %方法2中的顺序是按照小波包分解函数自动排列顺序，方法1经过重排序后为按照频率段递增顺序依次排列顺序的 %% %% 绘制重构前各个特征频段小波包系数的图形 figure(1); subplot(4,1,1); plot(x_input); title('原始信号'); subplot(4,1,2); plot(cfs3_0); title(['层数 ',num2str(3) ' 节点 0的小波0-8Hz',' 系数']) subplot(4,1,3); plot(cfs3_1); title(['层数 ',num2str(3) ' 节点 1的小波8-16Hz',' 系数']) subplot(4,1,4); plot(cfs3_2); title(['层数 ',num2str(3) ' 节点 2的小波16-24Hz',' 系数']) %% 绘制重构后时域各个特征频段的图形 figure(3); subplot(3,1,1); plot(rex3(:,1)); title('重构后0-8Hz频段信号'); subplot(3,1,2); plot(rex3(:,2)); title('重构后8-16Hz频段信号') subplot(3,1,3); plot(rex3(:,3)); title('重构后16-24Hz频段信号'); %% %Wpcoef 是求解某个节点的小波包系数，数据长度是1/(2^n)（分解n层的话），其实就是将原始信号化成2^N段，每段的长度是相等的且比原信号短 %wprcoef