小波包变换 matlab

cfs27=wpcoef(wpt1,[n,7]);
wpt2=wpdec(s3,n,ji);
rcfs30=wprcoef(wpt2,[n,0]);
rcfs31=wprcoef(wpt2,[n,1]);
rcfs32=wprcoef(wpt2,[n,2]);
rcfs33=wprcoef(wpt2,[n,3]);
for j=1:2^n
text(x(j),p_node(i),num2str(p_node(j),'%0.2f'),...
'HorizontalAlignment','center',...
'VerticalAlignment','bottom')
end
end
Matlab小波包分解后如何求各频带信号的能量值？
for i=1:2^n
pthi(i)= E(i)/E_total;%求每个节点的概率
end
pthi
%%%%以下提取小波包重构信号的故障诊断
clc
clear
fs=12e3;
loadzhengchang.mat;%加载正常信号
load 98.mat;
s1=X098_FE_time(19001:1:20000,1);% 正常情况下驱动端数据
% s3=X283_FE_time(1:1:10000,1);% 滚动体故障情况下驱动端数据
% load fan_out_12.mat;%加载外圈12点位置故障信号
load out_12.mat;%加载外圈12点位置故障信号
load 158.mat;
s4=X158_DE_time(14001:1:15000,1);% 外圈12点位置情况下驱动端数据
cfs15=wpcoef(wpt,[n,5]);
cfs16=wpcoef(wpt,[n,6]);
cfs17=wpcoef(wpt,[n,7]);
wpt1=wpdec(s2,n,ji);
rcfs20=wprcoef(wpt1,[n,0]);
rcfs21=wprcoef(wpt1,[n,1]);
rcfs22=wprcoef(wpt1,[n,2]);
% E
% disp('小波包分解总能量E_total');
E_total=sum(E); %求总能量
for i=1:2^n
pfir(i)= E(i)/E_total;%求每个节点的概率
end
pfir
% load a1_1-2;
load k302_1-4;
% Data=filter(b,a,Data);
end
% E = wenergy(wpt); only get the last layer
figure;
x=1:2^n;
bar(x,p_node);
title(['第',num2str(n),'层']);
axis([0 2^n 0 100]);
xlabel('结点');
ylabel('能量百分比/%');
Data=Data-mean(Data);
Data=detrend(Data);
wpt2=wpdec(Data,n,wpname); %对数据进行小波包分解
for i=1:2^n %wpcoef(wpt1,[n,i-1])是求第n层第i个节点的系数
E(i)=norm(wpcoef(wpt2,[n,i-1]),2);%求第i个节点的范数平方，其实也就是平方和
rex3=wprcoef(wpt,[n 3]);
;
subplot(5,1,1);
plot(x);
title('原始信号');
subplot(5,1,2);
plot(rex0);
title(['重构结点 ',num2str(n) ' 1',' 系数'])
subplot(5,1,3);
plot(rex1);
subplot(5,1,3);
plot(cfs1);
title(['结点 ',num2str(n) ' 2',' 系数'])
subplot(5,1,4);
plot(cfs2);
title(['结点 ',num2str(n) ' 3',' 系数'])
subplot(5,1,5);
plot(cfs3);
wpt3=wpdec(Data,n,wpname); %对数据进行小波包分解
for i=1:2^n %wpcoef(wpt1,[n,i-1])是求第n层第i个节点的系数
E(i)=norm(wpcoef(wpt3,[n,i-1]),2);%求第i个节点的范数平方，其实也就是平方和
end
E_total=sum(E); %求总能量
% s1=X098_FE_time(1:1:10000,1);% 正常情况下驱动端数据
load in.mat;%加载内圈故障信号
% load fan_in.mat;%加载风扇端内圈故障信号
load 106.mat;
s2=X106_DE_time(19001:1:20000,1);% 内圈故障信号情况下驱动端数据
% Plot wavelet packet tree (binary tree)
plot(wpt)
%% wavelet packet coefficients.default：use the front 4.
cfs0=wpcoef(wpt,[n 0]);
cfs1=wpcoef(wpt,[n 1]);
E(1:2^n )=0;
for i=1:2^n
E(i) = norm(wpcoef(wpt,[n,i-1]),2)^2; %20180604更新 原代码：E(i) = norm(wpcoef(wpt,[n,i-1]),2)
end
%求每个节点的概率
E_total=sum(E);
for i=1:2^n
p_node(i)= 100*E(i)/E_total;
小波包分解与重构
function wpt= wavelet_packetdecomposition_reconstruct( x,n,wpname )
%% 对信号进行小波包分解，得到节点的小波包系数。然后对每个节点系数进行重构。
% Decompose x at depth n with wpname wavelet ing Shannon entropy.
wpt1=wpdec(Data,n,wpname); %对数据进行小波包分解
for i=1:2^n %wpcoef(wpt1,[n,i-1])是求第n层第i个节点的系数
E(i)=norm(wpcoef(wpt1,[n,i-1]),2);%求第i个节点的范数平方，其实也就是平方和
end
% disp('每个节点的能量E(i)');
% s4=X305_FE_time(1:1:10000,1);% 外圈12点位置情况下驱动端数据
n=3;
ji='db5';
wpt=wpdec(s1,n,ji);%使用db5小波包分解正常信号到底三层，使用shannon熵
% plot(wpt);
rcfs10=wprcoef(wpt,[n,0]);%利用相应的频段进行重构
end
能量谱函数
function E = wavelet_energy_spectrum( wpt,n )
%% 计算每一层每一个节点的能量
% wpt-wavelet packet tree
% n-第n层能量
%
% Author hubery_zhang
% Date 20170714
%%
% 求第n层第i个节点的系数
cfs32=wpcoef(wpt2,[n,2]);
cfs33=wpcoef(wpt2,[n,3]);
cfs34=wpcoef(wpt2,[n,4]);
cfs35=wpcoef(wpt2,[n,5]);
cfs36=wpcoef(wpt2,[n,6]);
cfs37=wpcoef(wpt2,[n,7]);
clear;
clc;
n=3;
wpname='db3';
% [b,a]=butter(8,[5/100 99/100]);
% load a1_1-1;
load k301_1-4;
% Data=filter(b,a,Data);
Data=Data-mean(Data);
Data=detrend(Data);
rcfs11=wprcoef(wpt,[n,1]);
rcfs12=wprcoef(wpt,[n,2]);
rcfs13=wprcoef(wpt,[n,3]);
rcfs14=wprcoef(wpt,[n,4]);
rcfs15=wprcoef(wpt,[n,5]);
rcfs16=wprcoef(wpt,[n,6]);
title(['重构结点 ',num2str(n) ' 2',' 系数'])
subplot(5,1,4);
plot(rex2);
title(['重构结点 ',num2str(n) ' 3',' 系数'])
subplot(5,1,5);
plot(rex3);
title(['重构结点 ',num2str(n) ' 4',' 系数'])
rcfs23=wprcoef(wpt1,[n,3]);
rcfs24=wprcoef(wpt1,[n,4]);
rcfs25=wprcoef(wpt1,[n,5]);
rcfs26=wprcoef(wpt1,[n,6]);
rcfs27=wprcoef(wpt1,[n,7]);
cfs20=wpcoef(wpt1,[n,0]);
rcfs17=wprcoef(wpt,[n,7]);
cfs10=wpcoef(wpt,[n,0]);%提取小波包分解系数
cfs11=wpcoef(wpt,[n,1]);
cfs12=wpcoef(wpt,[n,2]);
cfs13=wpcoef(wpt,[n,3]);
cfs14=wpcoef(wpt,[n,4]);
title(['结点 ',num2str(n) ' 4',' 系数'])
%% reconstruct wavelet packet coefficients.
rex0=wprcoef(wpt,[n 0]);
rex1=wprcoef(wpt,[n 1]);
rex2=wprcoef(wpt,[n 2]);
cfs2=wpcoef(wpt,[n 2]);
cfs3=wpcoef(wpt,[n 3]);
figure;
subplot(5,1,1);
plot(x);
title('原始信号');
subplot(5,1,2);
plot(cfs0);
title(['结点 ',num2str(n) ' 1',' 系数'])
cfs21=wpcoef(wpt1,[n,1]);
cfs22=wpcoef(wpt1,[n,2]);
cfs23=wpcoef(wpt1,[n,3]);
cfs24=wpcoef(wpt1,[n,4]);
cfs25=wpcoef(wpt1,[n,5]);
cfs26=wpcoef(wpt1,[n,6]);
%
% x-input signal,列向量。
% n-the number of decomposition layers
% wpname-a particular wavelet.type:string.
%
%Author hubery_zhang
%Date 20170714
%%
wpt=wpdec(x,n,wpname);
end
E_total=sum(E); %求总能量
for i=1:2^n
psec(i)= E(i)/E_total;%求每个节点的概率
end
psec
% load a1_1-4;
load k303_1-4;
% Data=filter(b,a,Data);
Data=Data-mean(Data);
Data=detrend(Data);
% s2=X279_FE_time(1:1:10000,1);% 内圈故障信号情况下驱动端数据
% load fan_ball.mat;%加载滚动体故障信号
load ball.mat;%加载滚动体故障信号
load 119.mat;
s3=X119_DE_time(19001:1:20000,1);% 滚动体故障情况下驱动端数据
wpt3=wpdec(s4,n,ji);
rcfs40=wprcoef(wpt3,[n,0]);
rcfs34=wprcoef(wpt2,[n,4]);
rcfs35=wprcoef(wpt2,[n,5]);
rcfs36=wprcoef(wpt2,[n,6]);
rcfs37=wprcoef(wpt2,[n,7]);
cfs30=wpcoef(wpt2,[n,0]);
cfs31=wpcoef(wpt2,[n,1]);