小波包分解、重构、去噪与matlab函数使用

1 小波包分解、重构原理

1.1 小波变换

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

1.2 小波包变换

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

2 小波包使用到的函数及示例代码

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

2.1 小波包分解

2.1.1 T = wpdec(X,N,‘wname’)

T = wpdec(X,N,‘wname’) 小波分解函数
T为小波包分解树，X为待分解信号，N分解层数，'wname’小波基选择

% Load signal. 载入原始信号
load noisdopp; 
x = noisdopp;

figure(1); subplot(211); 
plot(x); title('Original signal');

% Decompose x at depth 3 with db1 wavelet packets 
% using Shannon entropy. 
wpt = wpdec(x,3,'db1');
plot(wpt) % 打印小波包分解树

2.1.2 X = wpcoef(T,N)

X = wpcoef(T,N) 小波包分解系数函数
X为返回的小波包系数，T为小波包分解树，N为小波树节点

% Read packet (2,1) coefficients. 读取小波包（2,1）节点系数
cfs = wpcoef(wpt,[2 1]);

figure(1); subplot(212); 
plot(cfs); title('Packet (2,1) coefficients');

需要注意的是：小波包节点系数的长度等于250，为原始信号长度1000/N = 250。其中N为2^2,为第N层小波包树的节点个数

2.2 小波包重构

2.2.1 X = wprcoef(T,N)

X = wprcoef(T,N) 小波包重构系数函数
X为重构后的小波包节点系数，T为小波包分解树，N为小波分解节点数

% Reconstruct packet (2,1).  重构小波包节点系数(2,1)
rcfs = wprcoef(wpt,[2 1]);

figure(1); subplot(212); 
plot(rcfs); title('Reconstructed packet (2,1)');

需要注意的是：小波重构系数的信号长度与原始信号长度一致

2.2.2 [T,X] = wpjoin(T,N)

[T,X] = wpjoin(T,N) 用于合并小波包树的节点
返回T：为修改后的小波包树T
返回X：合并后的节点系数
T：输入的小波包数T
N：为合并的节点

% Recompose packet (1,1)  重写小波包节点（1,1）
wpt = wpjoin(wpt,[1 1]);

% Plot wavelet packet tree wpt. 
plot(wpt)

% 创建小波包树
load noisdopp; x = noisdopp;
t = wpdec(x,3,'db3');
t = wpjoin(t,[4;5]); % 合并第4、5个小波树节点
plot(t);

2.2.3 T = write(T,‘cfs’,NODE,COEFS)

T = write(T,‘cfs’,NODE,COEFS) 将
返回T：修改后的小波包树
T：输入的小波包树
NODE：小波包节点序号
COEFS：小波包系数

% 将值写入第4、5、6个小波包树节点中
sNod = read(t,'sizes',[4,5,7]);  
cfs4 = zeros(sNod(1,:));
cfs5 = zeros(sNod(2,:));
cfs7 = zeros(sNod(3,:));
t = write(t,'cfs',4,cfs4,'cfs',5,cfs5,'cfs',7,cfs7);

% Plot tree t and click the node (0,0) (see the plot function).
% 打印原始节点（0,0）,即小波包重构信号的节点
plot(t)

2.2.4 X = wprec(T)

对于一位小波包分解重构
X = wprec(T) 重组小波包
X：为小波包分解树
T：重构后的信号

wpt = write(t,'cfs',4,cfs4,'cfs',5,cfs5,'cfs',7,cfs7);% 将值写入第4、5、6个小波包树节点中
xrec = wprec(wpt); % 重构后的信号

2.2.5 原始信号小波包分解后的重构方法

1. 可由小波包分解后，利用wprcoef(T,N)函数，对最后一层的小波包进行系数重构，得到与原始数据长度一样的重构小波包系数，然后将小波包系数相加，可以得到原始信号
2. 可以使用write(T,‘cfs’,NODE,COEFS)函数，将小波包系数加入小波包分解树，然后利用wprec(T) 函数进行原始信号重构，即可得到原始信号

2.3 小波包去噪

2.3.1 [THR,SORH,KEEPAPP,CRIT] = ddencmp(IN1,IN2,X)

[THR,SORH,KEEPAPP,CRIT] = ddencmp(IN1,IN2,X)
THR是阈值，SORH是阈值函数，KEEPAPP近似系数，CRIT熵的选择（仅用于小波包）
KEEPAPP：选择0或1可对近似系数进行阈值处理
IN1：‘den’ 表示去噪， ‘cmp’ 表示压缩.
IN2：'wv’用于小波，'wp’用于小波包
X：表示待去噪的信号

下面是小波去噪代码：

[thr,sorh,keepapp] = ddencmp('den','wv',x);
[A,D] = dwt(x,'db1');
noiselev = median(abs(D))/0.6745;
thresh = sqrt(2*log(length(x)))*noiselev;

2.3.2 [XC,CXC,LXC,PERF0,PERFL2] = wdencmp(‘gbl’,X,‘wname’,N,THR,SORH,KEEPAPP)

[XC,CXC,LXC,PERF0,PERFL2] = wdencmp(‘gbl’ or ‘lvd’,X,‘wname’,N,THR,SORH,KEEPAPP)
返回值XC：
返回值CXC：XC的小波包分解结构
返回值LXC：XC的小波包分解结构
返回值PERF0：恢复L^2的范数百分比, 是用百分制表明降噪所保留的能量成分
返回值PERFL2：压缩L^2的范数百分比, 是用百分制表明压缩所保留的能量成分
‘gbl’：表示每层都采用同一个阈值进行处理
‘lvd’：对每层采用不同的阈值进行处理
X：输入的信号
‘wname’：小波基
N：小波分解的层数
THR：阈值向量的长度
SORH：阈值函数，'s’软阈值，'h’硬阈值
KEEPAPP：0对低频近似系数不进行阈值量化处理；1对低频近似系数进行阈值量化处理

[thr,sorh,keepapp] = ddencmp('den','wv',x);
xd = wdencmp('gbl',x,'db3',2,thr,sorh,keepapp);
subplot(211)
plot(x); title('Original Signal');
subplot(212)
plot(xd); title('Denoised Signal');

参考文献

1.matlab自带的help
2. https://blog.csdn.net/cqfdcw/article/details/84995904.