数字图像处理学习笔记七小波（数学显微器）和多分辨率处理

最新推荐文章于 2024-07-18 21:05:49 发布

Belouga-

最新推荐文章于 2024-07-18 21:05:49 发布

阅读量2.7k

点赞数 1

分类专栏：数字图像处理

本文链接：https://blog.csdn.net/weixin_44330777/article/details/106717214

版权

数字图像处理专栏收录该内容

20 篇文章 10 订阅

订阅专栏

（一）傅立叶分析和小波分析之间的关系

1.1 傅里叶变换

1.2 短时傅里叶变换（Short-time Fourier Transform, STFT）

（一）傅立叶分析和小波分析之间的关系

从傅里叶变换到小波变换，并不是一个完全抽象的东西，可以理解形象一些。小波变换有着明确的物理意义，现在从它的提出时所面对的问题看起，试着整理出清晰的思路。

思路：傅里叶-->短时傅里叶变换-->小波变换，来理解为什么会出现小波这个东西、小波究竟是怎样的思路。

1.1 傅里叶变换

傅里叶变换的不足：对非平稳过程，傅里叶变换有局限性是为什么还要提出小波变换的契机。看如下一个简单的信号：

做完FFT（快速傅里叶变换）后，可以在频谱上看到清晰的四条线，信号包含四个频率成分。此时一切没有问题。再换频率随着时间变化的非平稳信号：

如上图，最上边的是频率始终不变的平稳信号。而下边两个则是频率随着时间改变的非平稳信号，它们同样包含和最上信号相同频率的四个成分。

做FFT后，对比发现这三个时域上有巨大差异的信号，频谱（幅值谱）却非常一致。尤其是下边两个非平稳信号，从频谱上无法区分它们，因为它们包含的四个频率的信号的成分确实是一样的，只是出现的先后顺序不同。可见，傅里叶变换处理非平稳信号有天生缺陷。它只能获取一段信号总体上包含哪些频率的成分，但是对各成分出现的时刻并无所知。因此时域相差很大的两个信号，可能频谱图一样。

然而平稳信号大多是人为制造出来的，自然界的大量信号几乎都是非平稳的，所以在比如生物医学信号分析等领域的论文中，基本看不到单纯傅里叶变换处理。

对于这样的非平稳信号，只知道包含哪些频率成分是不够的，还需要知道各个成分出现的时间。知道信号频率随时间变化的情况，各个时刻的瞬时频率及其幅值——即时频分析。

1.2 短时傅里叶变换（Short-time Fourier Transform, STFT）

一个简单可行的方法就是——加窗。简单说就是把整个时域过程分解成无数个等长的小过程，每个小过程近似平稳，再傅里叶变换，就知道在哪个时间点上出现了什么频率了。这就是短时傅里叶变换的思路。

时域上分成一段一段做FFT，分析出频率成分随着时间的变化，用这样的方法，可以得到一个信号的时频图：

图上既能看到10Hz, 25 Hz, 50 Hz, 100 Hz四个频域成分，还能看到出现的时间。两排峰是对称的，故只分析一排就行了。

时频分析结果到手。但是STFT依然有缺陷：窗函数的宽度如何确定？

窗太宽太窄都有问题：

窗太窄，窗内的信号太短，会导致频率分析不够精准，频率分辨率差。窗太宽，时域上又不够精细，时间分辨率低。

可以用海森堡不确定性原理来解释：类似于我们不能同时获取一个粒子的动量和位置，我们也不能同时获取信号绝对精准的时刻和频率。这也是一对不可兼得的矛盾体。我们不知道在某个瞬间哪个频率分量存在，我们知道的只能是在一个时间段内某个频带的分量存在。所以绝对意义的瞬时频率是不存在的。

实例效果：

上图对同一个信号（4个频率成分）采用不同宽度的窗做STFT，结果如右图。用窄窗，时频图在时间轴上分辨率很高，几个峰基本成矩形，而用宽窗则变成了绵延的矮山。但是频率轴上，窄窗明显不如下边两个宽窗精确。

所以窄窗口时间分辨率高、频率分辨率低，宽窗口时间分辨率低、频率分辨率高。对于时变的非稳态信号，高频适合小窗口，低频适合大窗口。然而STFT的窗口是固定的，在一次STFT中宽度不会变化，所以STFT还是无法满足非稳态信号变化的频率的需求。

1.3 小波变换

可以思考是否可以让窗口大小变起来，多做几次STFT呢——接近小波变换的思路了。但事实上小波并不是这么做的，不采用可变窗的STFT是因为这样做冗余会太严重，STFT做不到正交化，这也是它的一大缺陷。于是小波变换的出发点和STFT还是不同的。STFT是给信号加窗，分段做FFT；而小波直接把傅里叶变换的基给换了——将无限长的三角函数基换成了有限长的会衰减的小波基。这样不仅能够获取频率，还可以定位到时间。

傅里叶变换把无限长的三角函数作为基函数：

这个基函数会伸缩、会平移（其实是两个正交基的分解）。缩得窄，对应高频；伸得宽，对应低频。然后这个基函数不断和信号做相乘。某一个尺度（宽窄）下乘出来的结果，就可以理解成信号所包含的当前尺度对应频率成分有多少。于是，基函数会在某些尺度下，与信号相乘得到一个很大的值，因为此时二者有一种重合关系。借此知道信号包含该频率的成分的多少。

仔细体会可以发现这一步是在计算信号和三角函数的相关性。

分析发现这两种尺度能乘出一个大的值（相关度高），所以信号包含较多的这两个频率成分，在频谱上这两个频率会出现两个峰。

以上，就是粗浅意义上傅里叶变换的原理。

小波做的改变就在于，将无限长的三角函数基换成了有限长的会衰减的小波基。

“小波”，很小的一个波。小波变换公式：

从公式可以看出，不同于傅里叶变换，变量只有频率ω，小波变换有两个变量：尺度a（scale）和平移量 τ（translation）。尺度a控制小波函数的伸缩，平移量 τ控制小波函数的平移。尺度就对应于频率（反比），平移量 τ就对应于时间。

当伸缩、平移到这么一种重合情况时，也会相乘得到一个大的值。这时候和傅里叶变换不同的是，这不仅可以知道信号有这样频率的成分，而且知道它在时域上存在的具体位置。而当我们在每个尺度下都平移着和信号乘过一遍后，我们就知道信号在每个位置都包含哪些频率成分。非稳定信号用小波变换可以处理的很好，做时频分析。

做傅里叶变换只能得到一个频谱，做小波变换却可以得到一个时频谱。

时域信号：

傅里叶变换结果：

小波变换结果：

小波还有一些优点：

1. 我们知道对于突变信号，傅里叶变换存在吉布斯效应，我们用无限长的三角函数怎么也拟合不好突变信号

然而衰减的小波就不一样了：

2. 小波可以实现正交化，短时傅里叶变换不能。

1.4 哈尔变换（Haar）

哈尔变换的基函数是已知的最古老、也是最简单的正交小波。哈尔变换可用如下矩阵形式表示：

F是一个N x N图像矩阵，H是一个N x N哈尔变换矩阵，T是一个N x N变换结果。
H包含哈尔基函数hk(z)，z∈[0,1]，k = 0,1,2,…,N-1，其中N = 2n。要生成矩阵H，要定义整数k，k = 2p + q - 1，其中0<= p <= n-1，当p = 0时，q = 0或1；当p≠0时，1<= q <= 2n。

N x N哈尔变换矩阵的第i行包含了元素hi(z)，其中z = 0/N,1/N,2/N,…，(N-1)/N。
如2 x 2的哈尔变换矩阵为：

当N = 4时，且假设k,p,q的值如下:

时，4 x 4的变换矩阵H4为：

H2的行可用于定义一个2抽头完美重建滤波器组（见子带编码）的分析滤波器h0(n)和h1(n)，以及最简单且最古老的小波变换的缩放比例和小波向量。

1.5 Haar小波变化实验

对img进行haar小波变换：

cA,(cH,cV,cD)=dwt2(img,'haar')

小波变换之后，低频分量对应的图像：

cv2.imwrite('school_LF.png',np.uint8(cA/np.max(cA)*255))

处理效果：

小波变换之后，水平方向高频分量对应的图像：

cv2.imwrite('school_h.png',np.uint8(cH/np.max(cH)*255))

处理效果：

小波变换之后，垂直平方向高频分量对应的图像：

cv2.imwrite('school_v.png',np.uint8(cV/np.max(cV)*255))

处理效果：

小波变换之后，对角线方向高频分量对应的图像：

cv2.imwrite('school_d.png',np.uint8(cD/np.max(cD)*255))

根据小波系数重构回去的图像：

rimg = idwt2((cA,(cH,cV,cD)), 'haar')

（二）多分辨率展开

在多分辨率分析（MRA）中，尺度函数被用于建立一个函数或一幅图像的一系列近似，相邻两近似之间的分辨率相差2倍。使用称为小波的附加函数来对相邻近似之间的差进行编码。

2.1 级数展开

一个信号或函数f(x)可展开为函数的线性组合：

k为有限和或无限和的整数下标， $\alpha _{k}$ 为展开系数， $\varphi _{k}(x)$ 为展开函数。若展开唯一，即任意给定的f(x)只有一组 $\alpha _{k}$ 与之对应。则 $\varphi _{k}(x)$ 称为基函数，集合{ $\varphi _{k}(x)$ }称为可表示这样一类函数的基。
可展开的函数形成了一个函数空间，称为展开集合的闭合跨度，表示如下：

对任意函数空间V及相应的展开集合{ $\varphi _{k}(x)$ }，都有一个表示为{ $\widetilde{\varphi _{k}(x)}$ }的对偶函数集合。展开系数αk即是对偶函数 $\widetilde{\varphi _{k}(x)}$ 和函数 f(x) 的内积。

2.2 尺度函数

考虑由实、平方可积函数φ(x)的整数平移和二值尺度组成的展开函数集合{ $\varphi _{j,k}(x)$ }

其中 $\varphi _{j,k}(x)$ 对所有的j，k∈Z（整数集）和 $\varphi (x)$ 属于 $L^{2}(R)$ 都成立（ $L^{2}(R)$ 表示度量的、平方可积的一维函数集合，R为实数集）。
整数平移k决定了 $\varphi _{j,k}(x)$ 沿x轴的位置；尺度j决定了 $\varphi _{j,k }(x)$ 的宽度；2j/2控制函数的幅度。
因为 $\varphi _{j,k}(x)$ 的形状随 j 发生变化，故 $\varphi (x)$ 称为尺度函数。
选择适当的 $\varphi (x)$ ，可使{φj，k(x)}张成 $L^{2}(R)$ 。
对任意的j，将k上张成的子空间表示为：

增加j就会增加 $V_{j}$ 的大小，进而允许子空间中包含具有更小变量或更细细节的函数。因为随着j的增大，用于表示子空间函数的 $\varphi _{j,k}(x)$ 会变窄，且x有较小变化就可以分开。

2.3 小波函数

满足MRA要求的尺度函数被定义为小波函数 $\psi (x)$ ，它与其整数平移及二值尺度一起，跨越了任意连个相邻尺度子空间 $V_{j}$ 和 $V_{j+1}$ 之间的差。

对跨越图中 $W_{j}$ 空间的所有k∈Z，定义小波集合{ $\psi _{j,k}(x)$ }
使用尺度函数，可以写出：

如果f(x)∈ $W_{j}$ ，则

尺度函数和小波函数子空间由下式联系起来

⨁表示空间并集，Vj+1中Vj的正交补集是Wj，且Vj中的所有成员与Wj中的所有成员都正交。故

所有可度量的、平方可积的函数空间可以表示为

若f(x)是V1而非 $V_{0}$ 的元素，则(1)式的展开包含使用 $V_{0}$ 尺度函数的f(x)的近似；来自 $W_{0}$ 的小波将对这种近似于实际函数之间的差进行编码。
因为小波空间存在于由相邻较高分辨率尺度函数跨越的空间中，故小波函数可表示成平移后的双倍分辨率尺度函数的加权和，表示如下：

$h\psi (n)$ 为小波函数系数， $h\varphi (n)$ 为小波向量。
利用小波跨越正交补集空间和整数小波平移是正交的条件，可得 $h\psi (n)$ )和 $h\varphi (n)$ 按下述方式相关：

（三）小波变换

3.1 一维小波变换

小波级数展开

与小波ψ(x)和尺度函数φ(x)相关的函数f(x)∈ $L^{2}(R)$ 的小波级数展开如下：

$j_{0}$ 是任意的开始尺度， $C_{j0}(k)$ 和 $d_{j}(k)$ 分别是尺度函数和小波函数下的展开系数αk的改写形式。 $C_{j0}(k)$ 称为近似和/或尺度系数， $d_{j}(k)$ 称为细节和小波系数。展开系数计算：

离散小波变换

小波系数展开将一个连续变量函数映射为一系列系数。若待展开的函数是离散的（即数字序列），则得到的系数就称为离散小波变换（DWT）。则序列f(n)的正向DWT系数（类似于上一节的 $C_{j0}(k)$ 和 $d_{j}(k)$ ）如下：

连续小波变换

离散小波变换的自然延伸是连续小波变换（CWT），连续小波变换将一个连续函数变换为两个连续变量（平移和尺度）的高冗余度函数。得到的变换易于解释并且对于时间-频率分析是有价值的。我们感兴趣的是离散图像，此处只是为了完整性。
连续平方可积函数f(x)的连续小波变换与实数值小波ψ(x)的关系定义如下：

s和τ分别为尺度参数和平移参数。
f(x)可由连续小波反变换求得：

式中Ψ(μ)是ψ(x)的傅里叶变换。

3.2 二维小波变换

在二维情况下，需要一个二维尺度函数 $\phi (x,y)$ 和三个二维小波 $\psi ^{H}(x,y)$ ， $\psi ^{V}(x,y)$ ， $\psi ^{D}(x,y)$ 。每个二维小波都是两个一维函数的乘积。排除产生一维结果的乘积，如 $\phi (x,y)\psi (x,y)$ ，4个剩下的乘积可产生尺度函数：

和可分的“方向敏感”小波：

这些小波度量函数的变化——图像的灰度变化——沿不同方向的变化： $\psi ^{H}$ 度量沿列方向的变化（如水平边缘）， $\psi ^{V}$ 响应沿行方向的变化（如垂直边缘）， $\psi ^{D}$ 度量对应对角线方向的变化。

定义尺度和平移基函数：

M*N的图像f(x,y)的离散小波变换是：

离散反小波变换：

二维DWT可以用数字滤波器和抽样来实现：

小波在图像处理中的用途，如在傅里叶域那样，基本方法是:

计算一幅图像的二维小波变换
修改变换
计算反变换

3.3 python 小波变换实验

（1）二维图像单极变换

将多通道图像变为单通道图像：

img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY).astype(np.float32)
plt.figure('二维小波一级变换')
coeffs = pywt.dwt2(img, 'haar')
cA, (cH, cV, cD) = coeffs

处理效果：

（2）二维图像多极变换

二维图像多级分解：

plt.figure('二维图像多级分解')
coeffs = pywt.wavedec2(img, 'haar', level=2)
cA2, (cH2, cV2, cD2), (cH1, cV1, cD1) = coeffs

将每个子图的像素范围都归一化到与CA2一致 CA2 [0,255* 2**level]：

AH2 = np.concatenate([cA2, cH2+510], axis=1)
VD2 = np.concatenate([cV2+510, cD2+510], axis=1)
cA1 = np.concatenate([AH2, VD2], axis=0)

AH = np.concatenate([cA1, (cH1+255)*2], axis=1)
VD = np.concatenate([(cV1+255)*2, (cD1+255)*2], axis=1)
img = np.concatenate([AH, VD], axis=0)

处理效果：

（四）小波包

4.1 小波包概述

快速小波变换将一个函数分解为一系列与对数相关的频段，其中

低频被组成窄频段
高频被组成宽频段

由于正交小波变换只对信号的低频（近似）信息做进一步分解，而对高频（细节）信息不再继续分解，使得它的频率分辨率随频率升高而降低。所以小波变换能够很好地表征一大类以低频信息为主要成分的信号，但它不能很好地分解和表示包含大量细节信息（细小边缘或纹理）的信号，如非平稳机械振动信号、遥感图象、地震信号和生物医学信号等。与之不同的是，小波包变换可以对高频部分提供更精细的分解，而且这种分解既无冗余，也无疏漏，所以对包含大量中、高频信息的信号能够进行更好的时频局部化分析。
一个三层的小波包分解如下：

分析树提供了多尺度小波变换的紧凑有效的方法，比对应的滤波器和基于子取样的方框图更容易画，并占有较少的空间，相对容易定位有效分解。

4.2 小波分解

从教程找的代码：

import pywt

print(pywt.families())  # 打印出小波族
# ['haar', 'db', 'sym', 'coif', 'bior', 'rbio', 'dmey', 'gaus', 'mexh', 'morl', 'cgau', 'shan', 'fbsp', 'cmor']

for family in pywt.families():  # 打印出每个小波族的每个小波函数
    print('%s family: ' % (family) + ','.join(pywt.wavelist(family)))
# haar family: haar
# db family: db1,db2,db3,db4,db5,db6,db7,db8,db9,db10,db11,db12,db13,db14,db15,db16,db17,db18,db19,db20,db21,db22,db23,db24,db25,db26,db27,db28,db29,db30,db31,db32,db33,db34,db35,db36,db37,db38
# sym family: sym2,sym3,sym4,sym5,sym6,sym7,sym8,sym9,sym10,sym11,sym12,sym13,sym14,sym15,sym16,sym17,sym18,sym19,sym20
# coif family: coif1,coif2,coif3,coif4,coif5,coif6,coif7,coif8,coif9,coif10,coif11,coif12,coif13,coif14,coif15,coif16,coif17
# bior family: bior1.1,bior1.3,bior1.5,bior2.2,bior2.4,bior2.6,bior2.8,bior3.1,bior3.3,bior3.5,bior3.7,bior3.9,bior4.4,bior5.5,bior6.8
# rbio family: rbio1.1,rbio1.3,rbio1.5,rbio2.2,rbio2.4,rbio2.6,rbio2.8,rbio3.1,rbio3.3,rbio3.5,rbio3.7,rbio3.9,rbio4.4,rbio5.5,rbio6.8
# dmey family: dmey
# gaus family: gaus1,gaus2,gaus3,gaus4,gaus5,gaus6,gaus7,gaus8
# mexh family: mexh
# morl family: morl
# cgau family: cgau1,cgau2,cgau3,cgau4,cgau5,cgau6,cgau7,cgau8
# shan family: shan
# fbsp family: fbsp
# cmor family: cmor

db3 = pywt.Wavelet('db3')  # 创建一个小波对象
print(db3)


# Filters length: 6        #滤波器长度
# Orthogonal:     True    #正交
# Biorthogonal:   True    #双正交
# Symmetry:       asymmetric    #对称性，不对称
# DWT:            True    #离散小波变换
# CWT:            False    #连续小波变换

def print_array(arr):
    print('[%s]' % ','.join(['%.14f' % x for x in arr]))


# 离散小波变换的小波滤波系数
# dec_lo Decomposition filter values 分解滤波值， rec 重构滤波值
# db3.filter_bank 返回4 个属性
print(db3.filter_bank == (db3.dec_lo, db3.dec_hi, db3.rec_lo, db3.rec_hi))  # True
print(db3.dec_len)
print(db3.rec_len)  # 6

# DWT 与 IDWT
# 使用db2 小波函数做dwt
x = [3, 7, 1, 1, -2, 5, 4, 6]
cA, cD = pywt.dwt(x, 'db2')  # 得到近似值和细节系数
print(cA)  # [5.65685425 7.39923721 0.22414387 3.33677403 7.77817459]
print(cD)  # [-2.44948974 -1.60368225 -4.44140056 -0.41361256  1.22474487]

# IDWT
print(pywt.idwt(cA, cD, 'db2'))  # [ 3.  7.  1.  1. -2.  5.  4.  6.]

# 传入小波对象，设置模式
w = pywt.Wavelet('sym3')
cA, cD = pywt.dwt(x, wavelet=w, mode='constant')
print(cA)  # [ 4.38354585  3.80302657  7.31813271 -0.58565539  4.09727044  7.81994027]
print(cD)  # [-1.33068221 -2.78795192 -3.16825651 -0.67715519 -0.09722957 -0.07045258]

print(pywt.Modes.modes)
# ['zero', 'constant', 'symmetric', 'periodic', 'smooth', 'periodization', 'reflect', 'antisymmetric', 'antireflect']

print(pywt.idwt([1, 2, 0, 1], None, 'db3', 'symmetric'))
print(pywt.idwt([1, 2, 0, 1], [0, 0, 0, 0], 'db3', 'symmetric'))
# [ 0.83431373 -0.23479575  0.16178801  0.87734409]
# [ 0.83431373 -0.23479575  0.16178801  0.87734409]

# 小波包 wavelet packets
X = [1, 2, 3, 4, 5, 6, 7, 8]
wp = pywt.WaveletPacket(data=X, wavelet='db3', mode='symmetric', maxlevel=3)
print(wp)
print(wp.data)  # [1 2 3 4 5 6 7 8 9]
print(repr(wp.path))
print(wp.level)  # 0    #分解级别为0
print(wp['ad'].maxlevel)  # 3

# 访问小波包的子节点
# 第一层：
print(wp['a'].data)
# [ 4.52111203  1.54666942  2.57019338  5.3986205   8.19182134 11.27067814
#  12.65348525]        # 当设置分解的 maxlevel 时，分解得到的data

# [ 4.52111203  1.54666942  2.57019338  5.3986205   8.20681003 11.18125264] 设置为2 时
print(wp['a'].path)  # a
# 第2 层
print(wp['aa'].data)
# [ 3.63890166  6.00349136  2.89780988  6.80941869 15.41549196]
print(wp['ad'].data)
# [ 1.25531439 -0.60300027  0.36403471  0.59368086 -0.53821027]
print(wp['aa'].path)  # aa
print(wp['ad'].path)  # ad

# 第3 层时：
print(wp['aaa'].data)
# [ 6.7736584   5.78857317  5.69392399 10.98672847 19.92241106]

# print(wp['aaaa'].data)  #超过最大层时，会报错
# 获取特定层数的所有节点
print([node.path for node in wp.get_level(3, 'natural')])  # 第3层有8个
# ['aaa', 'aad', 'ada', 'add', 'daa', 'dad', 'dda', 'ddd']

# 依据频带频率进行划分
print([node.path for node in wp.get_level(3, 'freq')])
# ['aaa', 'aad', 'add', 'ada', 'dda', 'ddd', 'dad', 'daa']

# 从小波包中 重建数据
X = [1, 2, 3, 4, 5, 6, 7, 8]
wp = pywt.WaveletPacket(data=X, wavelet='db1', mode='symmetric', maxlevel=3)
print(wp['ad'].data)  # [-2,-2]
new_wp = pywt.WaveletPacket(data=None, wavelet='db1', mode='symmetric')
new_wp['a'] = wp['a']
new_wp['aa'] = wp['aa'].data
new_wp['ad'] = [-2, -2]  # wp['ad'].data
new_wp['d'] = wp['d']
print(new_wp.reconstruct(update=False))
# new_wp['a'] = wp['a']  直接使用高低频也可进行重构
# new_wp['d'] = wp['d']
print(new_wp)  #: None
print(new_wp.reconstruct(update=True))  # 更新设置为True时。
print(new_wp)
# : [1. 2. 3. 4. 5. 6. 7. 8.]

# 获取叶子结点
print([node.path for node in new_wp.get_leaf_nodes(decompose=False)])
# ['aa', 'ad', 'd']
print([node.path for node in new_wp.get_leaf_nodes(decompose=True)])
# ['aaa', 'aad', 'ada', 'add', 'daa', 'dad', 'dda', 'ddd']

# 从小波包树中移除结点
dummy = wp.get_level(2)
for i in wp.get_leaf_nodes(False):
    print(i.path, i.data)
# aa [ 5. 13.]
# ad [-2. -2.]
# da [-1. -1.]
# dd [-1.11022302e-16  0.00000000e+00]
node = wp['ad']
print(node)  # ad: [-2. -2.]
del wp['ad']  # 删除结点
for i in wp.get_leaf_nodes(False):
    print(i.path, i.data)
# aa [ 5. 13.]
# da [-1. -1.]
# dd [-1.11022302e-16  0.00000000e+00]

print(wp.reconstruct())  # 进行重建
# [2. 3. 2. 3. 6. 7. 6. 7.]

wp['ad'].data = node.data  # 还原已删除的结点
print(wp.reconstruct())
# [1. 2. 3. 4. 5. 6. 7. 8.]

print(wp['a'])
print(wp.a)

filename = r'D:\ml_datasets\PHM\c6\c_6_001.csv'
data = pd.read_csv(filename)
data = data.iloc[100000:110000, 3]
cA1, cD1 = pywt.dwt(data, 'db3')  # 得到近似值和细节系数

wap = pywt.WaveletPacket(data=data, wavelet='db3')
dataa = wap['a'].data
print(wap['a'].data)
print(len(wap['a'].data))
plt.figure(num='ca')
plt.plot(dataa)
plt.figure(num='data')
plt.plot(dataa)
plt.show()

运行出了点错误，改好了后再和大家进一步分享吧。