小波与小波包、小波包分解与信号重构、小波包能量特征提取

本篇为《信号处理》系列博客的第七篇，该系列博客主要记录信号处理相关知识的学习过程和自己的理解，方便以后查阅。

文章原地址：《小波与小波包、小波包分解与信号重构、小波包能量特征提取 暨 小波包分解后实现按频率大小分布重新排列(Matlab 程序详解）》

小波与小波包区别

工程应用中经常需要对一些非平稳信号进行分析，小波分析和小波包分析适合对非平稳信号分析

相比较小波分析，利用小波包分析可以对信号分析更加精细

小波包分析可以将时频平面划分的更为细致，对信号的高频部分的分辨率要好于小波分析

可以根据信号的特征，自适应的选择最佳小波基函数，以便更好的对信号进行分析，所以小波包分析应用更加广泛。

小波分解

小波变换只对信号的低频部分做进一步分解，而对高频部分也即信号的细节部分不再继续分解

所以小波变换能够很好地表征一大类以低频信息为主要成分的信号，不能很好地分解和表示包含大量细节信息（细小边缘或纹理）的信号

如非平稳机械振动信号、遥感图象、地震信号和生物医学信号等。

小波包分解

小波包变换既可以对低频部分信号进行分解，也可以对高频部分进行分解，

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

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

小波包树解读：

以上即是小波包树，其中节点的命名规则是从（1，0）开始，叫1号， （1，1）是2号………依此类推，（3，0）是7号，（3，7）是14号。

每个节点都有对应的小波包系数，这个系数决定了频率的大小，也就是说频率信息已经有了，但是时域信息在哪里呢？

那就是 order。 这个order就是这些节点的顺序，也就是频率的顺序。

Matlab实例

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

fs=1024;  %采样频率
f1=100;   %信号的第一个频率
f2=300;   %信号第二个频率

t=0:1/fs:1;

s=sin(2*pi*f1*t)+sin(2*pi*f2*t);  %生成混合信号

[tt]=wpdec(s,3,'dmey');  %小波包分解，3代表分解3层，'dmey'使用meyr小波

plot(tt)               %画小波包树图

wpviewcf(tt,1);        %画出时间频率图
1
2
3
4
5
6
7
8
9
10
11
12
13

时间频率图主要解释

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是第五段(图片从下往上数第5条)，也就是257-320Hz之间

这样就说通了，正好这个原始信号只有两个频率段，一个100一个300。
如果我们不是分解了3层，而是更多层，那么每个频率段包含的频率也就越窄，图上有颜色的地方也会更细，也就是说更精细了，由于原始信号的频率在整个1秒钟内都没有改变，所以有颜色的地方是一个横线。

下面是分解了5层的时间频率图：

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

在数值分析中，我们学过内积

内积的物理含义：两个图形的相似性，若两个图形完全正交，则内积为0，若两个图形完全一样，则系数为1（相对值）。

小波变换的实质是：原信号与小波基函数的相似性。
小波系数就是小波基函数与原信号相似的系数。

连续小波变换 ：小波函数与原信号对应点相乘，再相加，得到对应点的小波变换系数，平移小波基函数，再计算小波函数与原信号对应点相乘，再相加，这样就得到一系列的小波系数。

对于离散小波变换，由于很多小波函数不是正交函数，因此需要一个尺度函数
所以，原信号函数可以分解成尺度函数和小波函数的线性组合，在这个函数中，尺度函数产生低频部分，小波函数产生高频部分。

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

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

示例1

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

filename = '/home/al007/Matlab/聪姐论文程序/聪数据/肩屈曲总/肩屈曲数据/肩屈曲_Plot_and_Store_Rep_15.csv';
P=csvread(filename,1,0);%从第1行第0列开始读取数据
wave = P(:,2);%获取第一个通道的肌电数据

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

figure(1);%Create figure window
subplot(2,1,1);
plot(wave);
title('原始信号');

subplot(2,1,2);
plot(wave2);
title('小波重构信号');
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

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

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

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

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

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

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

可以预见，因为我们把小波树的节点系数都变为0了，所以信号也就全为0了。所以wave2是一个0向量。

进一步，如果我们只聚合第二层中的某几个节点，比如 3和4，则如下所示：

filename = "E:\研究生论文\基于EEG与sEMG的上肢动作识别研究\梁光金数据存放\sEMG数据存放\1上肢动作\1上肢动作3.csv";
P=csvread(filename,1,0);%从第1行第0列开始读取数据
wave = P(:,2);%获取第一个通道的肌电数据

t=wpdec(wave,3,'dmey');
t2 = wpjoin(t,[3;4;5;6]);
cfs3=wpcoef(t2,3);
cfs4=wpcoef(t2,4);
% cfs5=wpcoef(t2,5);
% cfs6=wpcoef(t2,6);

sNod = read(t,'sizes',[3,4,5,6]);
% cfs3  = zeros(sNod(1,:));
% cfs4  = zeros(sNod(2,:));
cfs5  = zeros(sNod(3,:));
cfs6  = zeros(sNod(4,:));

t3 = write(t2,'cfs',3,cfs3,'cfs',4,cfs4,'cfs',5,cfs5,'cfs',6,cfs6);
wave2=wprec(t3);

plot(wave);
hold on 
plot(wave2);
title('原始信号 与 小波重构信号 对比图');
legend('原始信号','小波重构信号')

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

如果我们只聚合第二层中4和5：

filename = "E:\研究生论文\基于EEG与sEMG的上肢动作识别研究\梁光金数据存放\sEMG数据存放\1上肢动作\1上肢动作3.csv";
P=csvread(filename,1,0);%从第1行第0列开始读取数据
wave = P(:,2);%获取第一个通道的肌电数据

t=wpdec(wave,3,'dmey');
t2 = wpjoin(t,[3;4;5;6]);
% cfs3=wpcoef(t2,3);
cfs4=wpcoef(t2,4);
cfs5=wpcoef(t2,5);
% cfs6=wpcoef(t2,6);

sNod = read(t,'sizes',[3,4,5,6]);
cfs3  = zeros(sNod(1,:));
% cfs4  = zeros(sNod(2,:));
% cfs5  = zeros(sNod(3,:));
cfs6  = zeros(sNod(4,:));

t3 = write(t2,'cfs',3,cfs3,'cfs',4,cfs4,'cfs',5,cfs5,'cfs',6,cfs6);
wave2=wprec(t3);

plot(wave);
hold on 
plot(wave2);
title('原始信号 与 小波重构信号 对比图');
legend('原始信号','小波重构信号')
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

出现上述两个不同的结果的原因，可以结合下面这段代码中的第二层小波包系数图。

示例2

我们首先来看一下wpjoin函数的效果以及第二层分节出来的小波包的系数是什么样的。

filename = "E:\研究生论文\基于EEG与sEMG的上肢动作识别研究\梁光金数据存放\sEMG数据存放\1上肢动作\1上肢动作3.csv";
P=csvread(filename,1,0);%从第1行第0列开始读取数据
wave = P(:,2);%获取第一个通道的肌电数据

figure(1);%Create figure window
t=wpdec(wave,3,'dmey');
plot(t)

figure(2);%Create figure window
t2 = wpjoin(t,[3;4;5;6]);
plot(t2)

figure(3);%Create figure window
cfs3=wpcoef(t2,3);
subplot(4,1,1);
plot(cfs3)
axis([0 1000 -5*10^(-6) 10*10^(-5)])

cfs4=wpcoef(t2,4);
subplot(4,1,2);
plot(cfs4)
axis([0 1000 -5*10^(-6) 10*10^(-5)])

cfs5=wpcoef(t2,5);
subplot(4,1,3);
plot(cfs5)
axis([0 1000 -5*10^(-6) 10*10^(-5)])

cfs6=wpcoef(t2,6);
subplot(4,1,4);
plot(cfs6)
axis([0 1000 -5*10^(-6) 10*10^(-5)])
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32

原始小波分解节点图：

wpjoin作用后的节点图：

第二层各节点小波包系数：

filename = '/home/al007/Matlab/聪姐论文程序/聪数据/肩屈曲总/肩屈曲数据/肩屈曲_Plot_and_Store_Rep_15.csv';
P=csvread(filename,1,0);%从第1行第0列开始读取数据
wave = P(:,2);%获取第一个通道的肌电数据

t=wpdec(wave,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);

figure(1);%Create figure window
subplot(2,1,1);
plot(wave);
title('原始信号');

subplot(2,1,2);
plot(wave2);
title('小波重构信号');
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

解释：

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

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

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

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

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

可以看出，该例子就是对一个小波包展开了，又原封不动的装回去了，所以说 wave2和wave是一样的。
注意：wpjoin命令在这里是必要的，因为write函数只能将最底层的节点写进去。可以回看一下示例2的第一段代码，wpjoin的效果。
也就是说，如果我们将第三层的小波包系数进行修改的话，就不用wpjoin了，具体可以看示例3

示例3

filename = '/home/al007/Matlab/聪姐论文程序/聪数据/肩屈曲总/肩屈曲数据/肩屈曲_Plot_and_Store_Rep_15.csv';
P=csvread(filename,1,0);%从第1行第0列开始读取数据
wave = P(:,2);%获取第一个通道的肌电数据

t=wpdec(wave,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);
wave2=wprec(t3);

figure(1);%Create figure window
subplot(2,1,1);
plot(wave);
title('原始信号');

subplot(2,1,2);
plot(wave2);
title('小波重构信号');
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

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

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

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

接下来进行举例说明：

filename = '/home/al007/Matlab/聪姐论文程序/聪数据/肩屈曲总/肩屈曲数据/肩屈曲_Plot_and_Store_Rep_15.csv';
P=csvread(filename,1,0);%从第1行第0列开始读取数据
wave = P(:,2);%获取第一个通道的肌电数据

plot(wave);
title('输入信号时域图像')   %绘制输入信号时域图像
1
2
3
4
5
6

% 查看频谱范围
x=wave;       
fs=2000;
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
1
2
3
4
5
6
7
8
9
10
11
12
13

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

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

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

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

rex3 = [];
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 50 0 0.03]); 
    title(['小波包第3层',num2str(i),'节点信号频谱']);
end
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

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

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

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

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

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

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

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

%% 
rex3 = [];
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(:,ord(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 50 0 0.03]); 
    title(['小波包第3层',num2str(i),'节点信号频谱']);
end
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

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

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

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

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

wpt=wpdec(wave,3,'dmey');        %进行3层小波包分解

cfs3_0=wpcoef(wpt,nodes_ord(1));  %对重排序后第3层0节点的小波包系数0-8Hz
cfs3_1=wpcoef(wpt,nodes_ord(2));  %对重排序后第3层0节点的小波包系数8-16Hz
cfs3_2=wpcoef(wpt,nodes_ord(3));  %对重排序后第3层0节点的小波包系数16-24Hz
cfs3_3=wpcoef(wpt,nodes_ord(4));  %对重排序后第3层0节点的小波包系数24-32Hz
cfs3_4=wpcoef(wpt,nodes_ord(5));  %对重排序后第3层0节点的小波包系数32-40Hz
cfs3_5=wpcoef(wpt,nodes_ord(6));  %对重排序后第3层0节点的小波包系数40-48Hz
cfs3_6=wpcoef(wpt,nodes_ord(7));  %对重排序后第3层0节点的小波包系数48-56Hz
cfs3_7=wpcoef(wpt,nodes_ord(8));  %对重排序后第3层0节点的小波包系数56-64Hz
 
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;           % 求得每个节点的占比
 
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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45

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

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

E = wenergy(wpt);%该函数只能对最后一层（底层）节点进行能量提取
bar(E);
title('各个频段能量所占的比例');
xlabel('频率 Hz');
ylabel('能量百分比/%');
1
2
3
4
5

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

重构前和重构后小波包系数对比

filename = '/home/al007/Matlab/聪姐论文程序/聪数据/肩屈曲总/肩屈曲数据/肩屈曲_Plot_and_Store_Rep_15.csv';
P=csvread(filename,1,0);%从第1行第0列开始读取数据
wave = P(:,2);%获取第一个通道的肌电数据

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

wpt=wpdec(wave,3,'dmey');        %进行3层小波包分解

cfs3_0=wpcoef(wpt,nodes_ord(1));  %对重排序后第3层0节点的小波包系数0-8Hz
cfs3_1=wpcoef(wpt,nodes_ord(2));  %对重排序后第3层0节点的小波包系数8-16Hz
cfs3_2=wpcoef(wpt,nodes_ord(3));  %对重排序后第3层0节点的小波包系数16-24Hz
cfs3_3=wpcoef(wpt,nodes_ord(4));  %对重排序后第3层0节点的小波包系数24-32Hz
cfs3_4=wpcoef(wpt,nodes_ord(5));  %对重排序后第3层0节点的小波包系数32-40Hz
cfs3_5=wpcoef(wpt,nodes_ord(6));  %对重排序后第3层0节点的小波包系数40-48Hz
cfs3_6=wpcoef(wpt,nodes_ord(7));  %对重排序后第3层0节点的小波包系数48-56Hz
cfs3_7=wpcoef(wpt,nodes_ord(8));  %对重排序后第3层0节点的小波包系数56-64Hz

% 绘制重构前各个频段小波包系数
figure(1);
subplot(4,1,1);
plot(wave);
title('原始信号');
subplot(4,1,2);
plot(cfs3_0);
title(['层数 ',num2str(3) '  节点 0的小波0-125Hz',' 系数'])
subplot(4,1,3);
plot(cfs3_1);
title(['层数 ',num2str(3) '  节点 1的小波125-250Hz',' 系数'])
subplot(4,1,4);
plot(cfs3_2);
title(['层数 ',num2str(3) '  节点 2的小波375-500Hz',' 系数'])

rex3 = [];
for i=1:8
    rex3(:,i)=wprcoef(wpt,nodes_ord(i));  %实现对节点小波节点进行重构        
end
% 绘制重构后时域各个特征频段的图形
figure(2);
subplot(4,1,1);
plot(wave);
title('原始信号');
subplot(4,1,2);
plot(rex3(:,1));
title('重构后0-125Hz频段信号');
subplot(4,1,3);
plot(rex3(:,2));
title('重构后125-250Hz频段信号')
subplot(4,1,4);
plot(rex3(:,3));
title('重构后375-500Hz频段信号');
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52

小波----常见基函数

与标准的傅里叶变换相比，小波分析中使用到的小波函数具有不唯一性，即小波函数具有多样性。
小波分析在工程应用中，一个十分重要的问题就是 最优小波基 的选择问题，因为用不同的小波基分析同一个问题会产生不同的结果。
目前我们主要是通过用小波分析方法处理信号的结果与理论结果的误差来判定小波基的好坏，由此决定小波基。

常用小波基有Haar小波、Daubechies(dbN)小波、Mexican Hat(mexh)小波、Morlet小波、Meyer小波等。

Haar小波

Haar函数是小波分析中最早用到的一个具有紧支撑的正交小波函数，也是最简单的一个小波函数，它是支撑域在范围内的单个矩形波。

Haar函数的定义如下：Haar小波在时域上是不连续的。
所以作为基本小波性能不是特别好。但它也有自己的优点:计算简单。

不但正交，而且与自己的整数位移正交，因此，在多分辨率系统中，Haar小波构成一组最简单的正交归一的小波族。

[phi,g1,xval] = wavefun('haar',20);

subplot(2,1,1);
plot(xval,g1,'LineWidth',2);
xlabel('t');
title('haar 时域');

g2=fft(g1);
g3=abs(g2); 

subplot(2,1,2);
plot(g3,'LineWidth',2);
xlabel('f');
title('haar 频域');
1
2
3
4
5
6
7
8
9
10
11
12
13
14

Daubechies(dbN)小波

Daubechies小波是世界著名的小波分析学者Inrid·Daubechies构造的小波函数，简写为dbN，N是小波的阶数。

% db4的时域和频域波形：
[phi,g1,xval] = wavefun('db4',10);

subplot(2,1,1);
plot(xval,g1,'LineWidth',2);
xlabel('t');
title('db4 时域');

g2=fft(g1);
g3=abs(g2);

subplot(2,1,2);
plot(g3,'LineWidth',2);
xlabel('f');
title('db4 频域');
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

[Lo_D,Hi_D,Lo_R,Hi_R] = wfilters('db4'); %计算该小波的4个滤波器

subplot(2,2,1); 
stem(Lo_D,'LineWidth',2);
title('分解低通滤波器');

subplot(2,2,2); 
stem(Hi_D,'LineWidth',2);
title('分解高通滤波器');

subplot(2,2,3); 
stem(Lo_R,'LineWidth',2);
title('重构低通滤波器');

subplot(2,2,4); 
stem(Hi_R,'LineWidth',2);
title('重构高通滤波器');
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

Mexican

Hat(mexh)小波Mexican Hat函数为Gauss函数的二阶导数：因为它的形状像墨西哥帽的截面，所以也称为墨西哥帽函数。Mexihat小波的时域和频域波形：

% Mexihat小波的时域和频域波形：
d=-6; h=6; n=100;   
[g1,x]=mexihat(d,h,n);

subplot(2,1,1);
plot(x,g1,'LineWidth',2);
xlabel('t');
title('Mexihat 时域');

g2=fft(g1);
g3=(abs(g2));

subplot(2,1,2);
plot(g3,'LineWidth',2);
xlabel('f');
title('mexihat 频域');
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

Mexihat小波的特点：在时间域与频率域都有很好的局部化，并且满足。不存在尺度函数，所以Mexihat小波函数不具有正交性。

Morlet小波

它是高斯包络下的单频率副正弦函数：其中C是重构时的归一化常数。Morlet小波没有尺度函数，而且是非正交分解。Morlet小波的时域和频域波形图：d=-6; h=6; n=100;

[g1,x]=morlet(d,h,n);

subplot(2,1,1);
plot(x,g1,'LineWidth',2);
xlabel('t');
title('morlet 时域');  

g2=fft(g1);
g3=(abs(g2));

subplot(2,1,2);
plot(g3,'LineWidth',2);
xlabel('f');
title('mexihat 频域');
1
2
3
4
5
6
7
8
9
10
11
12
13
14

注:

获取小波系数的两个函数理解：

wpcoef 是求解某个节点的小波包系数，数据长度是$1/(2^n$ （分解n层的话），其实就是将原始信号化成2^N段，每段的长度是相等的且比原信号短

wprcoef是把某个节点的小波包系数重构，得到的是和原信号一样长度的信号。