小波包能量谱分析

%% 假设对信号f(x)进行三层小波变换，得到小波系数d1,d2,d3和逼近系数a3，我想求出各个频率段的能量（能量=小波系数的平方和），请问如何编程呢？ 谢谢！

figure(1)
fs=2000;
t=0:1/fs:1;
for i=1:length(t)
if t(i)<=0.2
x(i)=8*sin(60*pi*t(i)+pi/12)+6*sin(240*pi*t(i)+pi/2)+sin(300*pi*t(i)+pi/6);
else
x(i)=8*sin(60*pi*t(i)+pi/12)+6*sin(240*pi*t(i)+pi/2)+sin(300*pi*t(i)+pi/6)+20*sin(500*pi*t(i)+pi/3).*exp(-20*(t(i)-0.2));
end
end
plot(t,x)
grid on
figure(2)
wpt=wpdec(x,4,'db3');
plot(wpt);
for i=0:15
figure(3);
rfs=wprcoef(wpt,[4 i]);
subplot(4,4,i+1);
plot(rfs);
grid on;
title(['[4 ',num2str(i),']']);
axis([0 2000 min(rfs) max(rfs)]);
end
for i=0:3
figure(4);
rfs=wprcoef(wpt,[4 i]);
subplot(4,1,i+1);
plot(rfs);
grid on;
title(['[4 ',num2str(i),']']);
axis([0 2000 min(rfs) max(rfs)]);
end
for i=1:4
wpt=wpdec(x,i,'db3');
e=wenergy(wpt);
E=zeros(1,length(e));
for j=1:2^i
E(j)=sum(abs(wprcoef(wpt,[i,j-1])).^2);
end
figure(5)
subplot(4,1,i);
bar(e);
axis([0 length(e) 0 130]);
title(['第 ',num2str(i), ' 层']);
for j=1:length(e)
text(j-0.2,e(j)+20,num2str(e(j),'%2.2f'));
end
end