转：https://blog.csdn.net/cqfdcw/article/details/84995904

        小波与小波包、小波包分解与信号重构、小波包能量特征提取   (Matlab 程序详解）

                                                       -----暨 小波包分解后解决频率大小分布重新排列问题

       本人当前对小波理解不是很深入，通过翻阅网络他人博客，进行汇总总结，重新调试Matlab代码，实现对小波与小波包、小波包分解与信号重构、小波包能量特征提取，供大家参考，后续将继续更新！

     本人在分析信号的过程中发现，按照网上所述的小波包分解方法理解，获取每层节点重构后信号频率并不是按照（n,0）、(n,1)...顺序依次由小到大排列的，经过进一步分析研究后发现，需要对节点进行重排序，具体操作见本文分析。

1.小波与小波包区别

        工程应用中经常需要对一些非平稳信号进行，小波分析和小波包分析适合对非平稳信号分析，相比较小波分析，利用小波包分析可以对信号分析更加精细，小波包分析可以将时频平面划分的更为细致，对信号的高频部分的分辨率要好于小波分析，可以根据信号的特征，自适应的选择最佳小波基函数，比便更好的对信号进行分析，所以小波包分析应用更加广泛。

①小波分解

         小波变换只对信号的低频部分做进一步分解，而对高频部分也即信号的细节部分不再继续分解，所以小波变换能够很好地表征一大类以低频信息为主要成分的信号，不能很好地分解和表示包含大量细节信息（细小边缘或纹理）的信号，如非平稳机械振动信号、遥感图象、地震信号和生物医学信号等。

②小波包分解

       小波包变换既可以对低频部分信号进行分解，也可以对高频部分进行分解，而且这种分解既无冗余，也无疏漏，所以对包含大量中、高频信息的信号能够进行更好的时频局部化分析。

2.小波包——小波包树与时频图

小波包树解读：

  以上即是小波包树，其中节点的命名规则是从（1，0）开始，叫1号， （1，1）是2号………依此类推，（3，0）是7号，（3，7）是14号。 每个节点都有对应的小波包系数，这个系数决定了频率的大小，也就是说频率信息已经有了，但是时域信息在哪里呢？ 那就是 order。  这个order就是这些节点的顺序，也就是频率的顺序。

Matlab实例：

采样频率为1024Hz，采样时间是1秒，有一个信号s是由频率100和200Hz的正弦波混合的，我们用小波包来分解。

1.  
   clear all
2.  
   clc
3.  
   fs=1024; %采样频率
4.  
   f1=100; %信号的第一个频率
5.  
   f2=300; %信号第二个频率
6.  
   t=0:1/fs:1;
7.  
   s=sin(2*pi*f1*t)+sin(2*pi*f2*t); %生成混合信号
8.  
   [tt]=wpdec(s,3,\'dmey\'); %小波包分解，3代表分解3层，\'dmey\'使用meyr小波
9.  
   plot(tt) %画小波包树图
10.  
    wpviewcf(tt,1); %画出时间频率图

主要解释：

x轴，就是1024个点，对应1秒，每个点就代表1/1024秒。

y轴，显示的数字对应于小波包树中的节点，从下面开始，顺序是 7号节点，8号，10号，9号，，，，11号节点，这个顺序是这么排列的，这是小波包自动排列的。然后，y轴是频率啊，怎么不是 100Hz和300Hz呢？我们的采样频率是1024Hz，根据采样定理，奈奎斯特采样频率是512Hz，我们分解了3层，最后一层就是 2^3=8个频率段，每个频率段的频率区间是 512/8=64Hz,看图颜色重的地方一个是在8那里，一个在13那里，8是第二段，也就是 65-128Hz之间，13是第五段，也就是257-320Hz之间。这样就说通了，正好这个原始信号只有两个频率段，一个100一个300。如果我们不是分解了3层，而是更多层，那么每个频率段包含的频率也就越窄，图上有颜色的地方也会更细，也就是说更精细了，由于原始信号的频率在整个1秒钟内都没有改变，所以有颜色的地方是一个横线。（引用：http://www.cnblogs.com/welen/articles/5667217.html ）

3.小波包-----小波包分解系数

       在数值分析中，我们学过内积，内积的物理含义：两个图形的相似性，若两个图形完全正交，则内积为0，若两个图形完全一样，则系数为1（相对值）。小波变换的实质是：原信号与小波基函数的相似性。小波系数就是小波基函数与原信号相似的系数。

      连续小波变换：小波函数与原信号对应点相乘，再相加，得到对应点的小波变换系数，平移小波基函数，再计算小波函数与原信号对应点相乘，再相加，这样就得到一系列的小波系数。对于离散小波变换（由于很多小波函数不是正交函数，因此需要一个尺度函数）所以，原信号函数可以分解成尺度函数和小波函数的线性组合，在这个函数中，尺度函数产生低频部分，小波函数产生高频部分。

4.小波包-----信号分解与重构（方法1）

该方法可以实现对任意节点系数选择进行组合重构。

例1

有一个信号，变量名为wave，随便找一个信号load进来就行了。

1.  
   t=wpdec(wave,3,\'dmey\');
2.  
   t2 = wpjoin(t,[3;4;5;6]);
3.  
   sNod = read(t,\'sizes\',[3,4,5,6]);
4.  
   cfs3 = zeros(sNod(1,:));
5.  
   cfs4 = zeros(sNod(2,:));
6.  
   cfs5 = zeros(sNod(3,:));
7.  
   cfs6 = zeros(sNod(4,:));
8.  
   t3 = write(t2,\'cfs\',3,cfs3,\'cfs\',4,cfs4,\'cfs\',5,cfs5,\'cfs\',6,cfs6);
9.  
   wave2=wprec(t3);

第一行：将wave 用 dmey小波进行3层小波包分解，获得一个小波包树 t

第二行：将小波包树的第二行的四个节点收起来，也就是让第二行的节点变为树的最底层节点。因为第一行中小波包树的节点个数是第一层2个，第二层4个，第三层8个。现在t2就是将第三层的节点再聚合回第二层。

第三行：读取第二层四个节点系数的size

第四~七行：将所有四个节点的小波包系数变为0

第八行：将四个节点的系数重组到t3小波树中。

第九行：对t3小波树进行重构，获得信号wave2

       可以预见，因为我们把小波树的节点系数都变为0了，所以信号也就全为0了。所以wave2是一个0向量。读者可以自行plot一下wave和wave2看看。进一步，如果我们只聚合第二层中的某几个节点，比如 4和5，即将第三行到第八行中节点 3 和节点 6的语句删除或修改，那么意思就是将 4  5节点的系数变为0，那么wave2肯定就不是0向量了。

例2

1.  
   t=wpdec(wave,3,\'dmey\');
2.  
   t2 = wpjoin(t,[3;4;5;6]);
3.  
   cfs3=wpcoef(t,3);
4.  
   cfs4=wpcoef(t,4);
5.  
   cfs5=wpcoef(t,5);
6.  
   cfs6=wpcoef(t,6);
7.  
   t3 = write(t2,\'cfs\',3,cfs3,\'cfs\',4,cfs4,\'cfs\',5,cfs5,\'cfs\',6,cfs6);
8.  
   wave2=wprec(t3);

解释：

第一行：将wave 用 dmey小波进行3层小波包分解，获得一个小波包树 t

第二行：将小波包树的第二行的四个节点收起来，也就是让第二行的节点变为树的最底层节点。

第三~六行：获取四个节点的小波包系数 （小波包系数就是一个一维向量）

第七行：将四个节点的系数重组到t3小波树中

第八行：对t3小波树进行重构，获得信号wave2

可以看出，该例子就是对一个小波包展开了，又原封不动的装回去了，所以说 wave2和wave是一样的。

注意，wpjoin命令在这里是必要的，因为write函数只能将最底层的节点写进去。也就是说，如果我们将第三层的小波包系数进行修改的话，就不用wpjoin了，具体可以看例3

例3

1.  
   t=wpdec(wave,3,\'dmey\');
2.  
   cfs7=wpcoef(t,7);
3.  
   cfs8=wpcoef(t,8);
4.  
   cfs9=wpcoef(t,9);
5.  
   cfs10=wpcoef(t,10);
6.  
   cfs11=wpcoef(t,11);
7.  
   cfs12=wpcoef(t,12);
8.  
   cfs13=wpcoef(t,13);
9.  
   cfs14=wpcoef(t,14);
10.  
    t3=write(t,\'cfs\',7,cfs7,\'cfs\',8,cfs8,\'cfs\',9,cfs9,\'cfs\',10,cfs10,\'cfs\',11,cfs11,\'cfs\',...
11.  
    12,cfs12,\'cfs\',13,cfs13,\'cfs\',14,cfs14);
12.  
    y=wprec(t3);

该例子也是对一个小波包展开了，又原封不动的装回去了，只不过这次是直接对第三层节点进行的。

5.小波包-----信号分解与重构（方法2）

该方法只能对某一节点信号系数分别进行重构，不能实现多个节点系数组合进行重构.

接下来进行举例说明：

1.  
   x_input=x_train(:,1,1); %输入数据
2.  
   plot(x_input);title(\'输入信号时域图像\') %绘制输入信号时域图像

1.  
   %% 查看频谱范围
2.  
   x=x_input;
3.  
   fs=128;
4.  
   N=length(x); %采样点个数
5.  
   signalFFT=abs(fft(x,N));%真实的幅值
6.  
   Y=2*signalFFT/N;
7.  
   f=(0:N/2)*(fs/N);
8.  
   figure;plot(f,Y(1:N/2+1));
9.  
   ylabel(\'amp\'); xlabel(\'frequency\');title(\'输入信号的频谱\');grid on

接下来进行4层小波包分解

1.  
   wpt=wpdec(x_input,3,\'dmey\'); %进行3层小波包分解
2.  
   plot(wpt); %绘制小波包树

接下来，我们查看第3层8个节点的频谱分布（我的输入信号采样频率是128，按照采样定理小波包分解根节点（0，0）处的频率应该为0-64Hz，按照这个计算（3，0）节点为0-8Hz，后面依次以8Hz为一个段递增）

首先展示最初的重构方法（频率排布顺序混乱，常见理解错误）

1.  
   for i=0:7
2.  
   rex3(:,i+1)=wprcoef(wpt,[3 i]); %实现对节点小波节点进行重构
3.  
   end
4.  
    
5.  
   figure; %绘制第3层各个节点分别重构后信号的频谱
6.  
   for i=0:7
7.  
   subplot(2,4,i+1);
8.  
   x_sign=rex3(:,i+1);
9.  
   N=length(x_sign); %采样点个数
10.  
    signalFFT=abs(fft(x_sign,N));%真实的幅值
11.  
    Y=2*signalFFT/N;
12.  
    f=(0:N/2)*(fs/N);
13.  
    plot(f,Y(1:N/2+1));
14.  
    ylabel(\'amp\'); xlabel(\'frequency\');grid on
15.  
    axis([0 50 0 0.03]); title([\'小波包第3层\',num2str(i),\'节点信号频谱\']);
16.  
    end

绘制完后你会发现频谱分布并不是按照之前理解的顺序依次递增排列，如下所示

那么问题出在哪里呢？经过仔细深入验证后发现问题出在小波包节点的频谱划分“不是严格按照上述理解的顺序排列的”（可能是一种格雷编码或者其他），要解决这个问题我们需要对节点顺序进行重新编码排序。参考这里https://www.ilovematlab.cn/thread-122226-1-1.html?tdsourcetag=s_pctim_aiomsg

1.  
   nodes=[7;8;9;10;11;12;13;14]; %第3层的节点号
2.  
   ord=wpfrqord(nodes);  %小波包系数重排，ord是重排后小波包系数索引构成的矩阵　如3层分解的[1;2;4;3;7;8;6;5]
3.  
   nodes_ord=nodes(ord); %重排后的小波系数

注意：节点的命名规则是从（1，0）开始，叫1号， （1，1）是2号………依此类推，（3，0）是7号，（3，7）是14号。

那么我们来看看经过上面这段重排序运行后nodes_ord中顺序是什么？

接下来我们再绘制一下第三层8个节点重构信号的频谱如下

1.  
   for i=1:8
2.  
   rex3(:,i)=wprcoef(wpt,nodes_ord(i)); %实现对节点小波节点进行重构
3.  
   end
4.  
    
5.  
   figure; %绘制第3层各个节点分别重构后信号的频谱
6.  
   for i=0:7
7.  
   subplot(2,4,i+1);
8.  
   x_sign= rex3(:,i+1);
9.  
   N=length(x_sign); %采样点个数
10.  
    signalFFT=abs(fft(x_sign,N));%真实的幅值
11.  
    Y=2*signalFFT/N;
12.  
    f=(0:N/2)*(fs/N);
13.  
    plot(f,Y(1:N/2+1));
14.  
    ylabel(\'amp\'); xlabel(\'frequency\');grid on
15.  
    axis([0 50 0 0.03]); title([\'小波包第3层\',num2str(i),\'节点信号频谱\']);
16.  
    end

到这里为止不知道大家有没有明白我要表达的意思，如果没有明白可以再反复看理解，或者自己进行分解后观察每个节点的频谱后或许就理解了。

6.小波包分解------能量特征提取（方法1）

接下来绘制第3层各个频段的能量占比

1.  
   %% wavelet packet coefficients. 求取小波包分解的各个节点的小波包系数
2.  
   cfs3_0=wpcoef(wpt,nodes_ord(1)); %对重排序后第3层0节点的小波包系数0-8Hz
3.  
   cfs3_1=wpcoef(wpt,nodes_ord(2)); %对重排序后第3层0节点的小波包系数8-16Hz
4.  
   cfs3_2=wpcoef(wpt,nodes_ord(3)); %对重排序后第3层0节点的小波包系数16-24Hz
5.  
   cfs3_3=wpcoef(wpt,nodes_ord(4)); %对重排序后第3层0节点的小波包系数24-32Hz
6.  
   cfs3_4=wpcoef(wpt,nodes_ord(5)); %对重排序后第3层0节点的小波包系数32-40Hz
7.  
   cfs3_5=wpcoef(wpt,nodes_ord(6)); %对重排序后第3层0节点的小波包系数40-48Hz
8.  
   cfs3_6=wpcoef(wpt,nodes_ord(7)); %对重排序后第3层0节点的小波包系数48-56Hz
9.  
   cfs3_7=wpcoef(wpt,nodes_ord(8)); %对重排序后第3层0节点的小波包系数56-64Hz
10.  
     
11.  
    E_cfs3_0=norm(cfs3_0,2)^2; %% 1-范数：就是norm(...,1)，即各元素绝对值之和；2-范数：就是norm(...,2)，即各元素平方和开根号；
12.  
    E_cfs3_1=norm(cfs3_1,2)^2;
13.  
    E_cfs3_2=norm(cfs3_2,2)^2;
14.  
    E_cfs3_3=norm(cfs3_3,2)^2;
15.  
    E_cfs3_4=norm(cfs3_4,2)^2;
16.  
    E_cfs3_5=norm(cfs3_5,2)^2;
17.  
    E_cfs3_6=norm(cfs3_6,2)^2;
18.  
    E_cfs3_7=norm(cfs3_7,2)^2;
19.  
    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;
20.  
     
21.  
    p_node(0)= 100*E_cfs3_0/E_total; % 求得每个节点的占比
22.  
    p_node(2)= 100*E_cfs3_1/E_total; % 求得每个节点的占比
23.  
    p_node(3)= 100*E_cfs3_2/E_total; % 求得每个节点的占比
24.  
    p_node(4)= 100*E_cfs3_3/E_total; % 求得每个节点的占比
25.  
    p_node(5)= 100*E_cfs3_4/E_total; % 求得每个节点的占比
26.  
    p_node(6)= 100*E_cfs3_5/E_total; % 求得每个节点的占比
27.  
    p_node(7)= 100*E_cfs3_6/E_total; % 求得每个节点的占比
28.  
    p_node(8)= 100*E_cfs3_7/E_total; % 求得每个节点的占比
29.  
     
30.  
    figure;
31.  
    x=1:8;
32.  
    bar(x,p_node);
33.  
    title(\'各个频段能量所占的比例\');
34.  
    xlabel(\'频率 Hz\');
35.  
    ylabel(\'能量百分比/%\');
36.  
    for j=1:8
37.  
    text(x(j),p_node(j),num2str(p_node(j),\'%0.2f\'),...
38.  
    \'HorizontalAlignment\',\'center\',...
39.  
    \'VerticalAlignment\',\'bottom\')
40.  
    end

7.小波包分解------能量特征提取（方法2）

直接运行matlab自带函数，如下

E = wenergy(wpt);   %该函数只能对最后一层（底层）节点进行能量提取

（引用：https://blog.csdn.net/it_beecoder/article/details/78668273）

这里对比一下，和方法1得到的图形基本上是一致的，不同之处在于排列顺序变了，方法2中的顺序是按照小波包分解函数自动排列顺序，方法1经过重排序后为按照频率段递增顺序依次排列顺序的

1.  
   %% 绘制重构前各个频段小波包系数
2.  
   figure(1);
3.  
   subplot(4,1,1);
4.  
   plot(x_input);
5.  
   title(\'原始信号\');
6.  
   subplot(4,1,2);
7.  
   plot(cfs3_0);
8.  
   title([\'层数 \',num2str(3) \' 节点 0的小波0-8Hz\',\' 系数\'])
9.  
   subplot(4,1,3);
10.  
    plot(cfs3_1);
11.  
    title([\'层数 \',num2str(3) \' 节点 1的小波8-16Hz\',\' 系数\'])
12.  
    subplot(4,1,4);
13.  
    plot(cfs3_2);
14.  
    title([\'层数 \',num2str(3) \' 节点 2的小波16-24Hz\',\' 系数\'])

1.  
   %% 绘制重构后时域各个特征频段的图形
2.  
   figure(3);
3.  
   subplot(3,1,1);
4.  
   plot(rex3(:,1));
5.  
   title(\'重构后0-8Hz频段信号\');
6.  
   subplot(3,1,2);
7.  
   plot(rex3(:,2));
8.  
   title(\'重构后8-16Hz频段信号\')
9.  
   subplot(3,1,3);
10.  
    plot(rex3(:,3));
11.  
    title(\'重构后16-24Hz频段信号\');

以上过程完整代码

1.  
   clear all;
2.  
    
3.  
   x_input=x_train(:,1,1); %输入数据
4.  
   plot(x_input);title(\'输入信号时域图像\') %绘制输入信号时域图像
5.  
    
6.  
   x=x_input; % 查看频谱范围
7.  
   fs=128;
8.  
   N=length(x); %采样点个数
9.  
   signalFFT=abs(fft(x,N));%真实的幅值
10.  
    Y=2*signalFFT/N;
11.  
    f=(0:N/2)*(fs/N);
12.  
    figure;plot(f,Y(1:N/2+1));
13.  
    ylabel(\'amp\'); xlabel(\'frequency\');title(\'输入信号的频谱\');grid on
14.  
     
15.  
    wpt=wpdec(x_input,3,\'dmey\'); %进行4层小波包分解
16.  
    plot(wpt); %绘制小波包树
17.  
     
18.  
    %% 实现对节点顺序按照频率递增进行重排序
19.  
    % nodes=get(wpt,\'tn\'); %小波包分解系数　为什么ｎｏｄｅｓ是[7;8;9;10;11;12;13;14]
20.  
    % N_cfs=length(nodes); %小波包系数个数
21.  
    nodes=[7;8;9;10;11;12;13;14];
22.  
    ord=wpfrqord(nodes); %小波包系数重排，ord是重排后小波包系数索引构成的矩阵　如3层分解的[1;2;4;3;7;8;6;5]
23.  
    nodes_ord=nodes(ord); %重排后的小波系数
24.  
     
25.  
    %% 实现对节点小波节点进行重构
26.  
    for i=1:8
27.  
    rex3(:,i)=wprcoef(wpt,nodes_ord(i));
28.  
    end
29.  
     
30.  
    %% 绘制第3层各个节点分别重构后信号的频谱
31.  
    figure;
32.  
    for i=0:7
33.  
    subplot(2,4,i+1);
34.  
    x_sign= rex3(:,i+1);
35.  
    N=length(x_sign); %采样点个数
36.  
    signalFFT=abs(fft(x_sign,N));%真实的幅值
37.  
    Y=2*signalFFT/N;
38.  
    f=(0:N/2)*(fs/N);
39.  
    plot(f,Y(1:N/2+1));
40.  
    ylabel(\'amp\'); xlabel(\'frequency\');grid on
41.  
    axis([0 50 0 0.03]); title([\'小波包第3层\',num2str(i),\'节点信号频谱\']);
42.  
    end
43.  
     
44.  
    %% wavelet packet coefficients. 求取小波包分解的各个频段的小波包系数
45.  
    cfs3_0=wpcoef(wpt,nodes_ord(1)); %对重排序后第3层0节点的小波包系数0-8Hz
46.  
    cfs3_1=wpcoef(wpt,nodes_ord(2)); %对重排序后第3层0节点的小波包系数8-16Hz
47.  
    cfs3_2=wpcoef(wpt,nodes_ord(3)); %对重排序后第3层0节点的小波包系数16-24Hz
48.  
    cfs3_3=wpcoef(wpt,nodes_ord(4)); %对重排序后第3层0节点的小波包系数24-32Hz
49.  
    cfs3_4=wpcoef(wpt,nodes_ord(5)); %对重排序后第3层0节点的小波包系数32-40Hz
50.  
    cfs3_5=wpcoef(wpt,nodes_ord(6)); %对重排序后第3层0节点的小波包系数40-48Hz
51.  
    cfs3_6=wpcoef(wpt,nodes_ord(7)); %对重排序后第3层0节点的小波包系数48-56Hz
52.  
    cfs3_7=wpcoef(wpt,nodes_ord(8)); %对重排序后第3层0节点的小波包系数56-64Hz
53.  
    %% 求取小波包分解的各个频段的能量
54.  
    E_cfs3_0=norm(cfs3_0,2)^2; %% 1-范数：就是norm(...,1)，即各元素绝对值之和；2-范数：就是norm(...,2)，即各元素平方和开根号；
55.  
    E_cfs3_1=norm(cfs3_1,2)^2;
56.  
    E_cfs3_2=norm(cfs3_2,2)^2;
57.  
    E_cfs3_3=norm(cfs3_3,2)^2;
58.  
    E_cfs3_4=norm(cfs3_4,2)^2;
59.  
    E_cfs3_5=norm(cfs3_5,2)^2;
60.  
    E_cfs3_6=norm(cfs3_6,2)^2;
61.  
    E_cfs3_7=norm(cfs3_7,2)^2;
62.  
    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;
63.  
     
64.  
    p_node(0)= 100*E_cfs3_0/E_total; % 求得每个节点的占比
65.  
    p_node(2)= 100*E_cfs3_1/E_total; % 求得每个节点的占比
66.  
    p_node(3)= 100*E_cfs3_2/E_total; % 求得每个节点的占比
67.  
    p_node(4)= 100*E_cfs3_3/E_total; % 求得每个节点的占比
68.  
    p_node(5)= 100*E_cfs3_4/E_total; % 求得每个节点的占比
69.  
    p_node(6)= 100*E_cfs3_5/E_total; % 求得每个节点的占比
70.  
    p_node(7)= 100*E_cfs3_6/E_total; % 求得每个节点的占比
71.  
    p_node(8)= 100*E_cfs3_7/E_total; % 求得每个节点的占比
72.  
     
73.  
    figure;
74.  
    x=1:8;
75.  
    bar(x,p_node);
76.  
    title(\'各个频段能量所占的比例\');
77.  
    xlabel(\'频率 Hz\');
78.  
    ylabel(\'能量百分比/%\');
79.  
    for j=1:8
80.  
    text(x(j),p_node(j),num2str(p_node(j),\'%0.2f\'),...
81.  <