Python 离散小波变换（DWT） pywt库

最新推荐文章于 2023-09-24 14:15:16 发布

SongpingWang

最新推荐文章于 2023-09-24 14:15:16 发布

阅读量6w

点赞数 127

分类专栏： OpenCV/python 计算机视觉文章标签：计算机视觉 python

本文链接：https://blog.csdn.net/wsp_1138886114/article/details/116780542

版权

OpenCV/python 同时被 2 个专栏收录

84 篇文章 240 订阅

订阅专栏

计算机视觉

42 篇文章 30 订阅

订阅专栏

文章目录

简便安装： pip install PyWavelets 或者 conda install PyWavelets
源码安装：下载源码

tar zxvf PyWavelets-1.1.1.tar.gz PyWavelets-1.1.1
cd PyWavelets-1.1.1
sudo python3 setup.py build
sudo python3 setup.py install

使用文档：https://pywavelets.readthedocs.io/en/latest/ref/index.html

一、小波变换

尺度函数： scaling function （又称为父函数 father wavelet ）
小波函数： wavelet function（又称为母函数 mother wavelet）
连续的小波变换：CWT
离散的小波变换：DWT

小波变换的基本知识：
不同的小波基函数，是由同一个基本小波函数经缩放和平移生成的。
小波变换是将原始图像与小波基函数以及尺度函数进行内积运算, 所以一个尺度函数和一个小波基函数就可以确定一个小波变换。

离散小波变换函数

下面先列举3条关键理解：

小波分解，分解到的"不是频率域"！可以抽象理解为"小波域"，但其实没有实际内涵！傅里叶变换到频率域是有实际内涵的；
小波分解得到的"小波系数"是"没有量纲"的！它其实是"没有实际意义的数"，需要做系数重构才能从"小波域"再转回到"时域"；
“系数重构"与"重构信号"不是一个东西！系数重构就是把无量纲的小波分解系数变回到有意义的"时域”；"重构信号"就是把分解的完整恢复回去。

离散小波变换函数：不同的适用处于搭配函数.

分解与重构/恢复信号：
- 1级分解与重构原始信号函数为: dwt和dwt2 与 idwt和idwt2；
- 多级(包括1级)分解与重构原始信号函数为: wavedec和wavedec2 与 waverec和waverec2；所以wavedec可涵盖dwt。
系数重构：
- 1级分解的系数重构用函数的是: upcoef和upcoef2；
- 多级分解的系数重构用函数的是: wrcoef和wrcoef2。
(多级)系数提取：
- 多级分解低频近似系数提取：appcoef和appcoef2；
- 多解分解高频细节系数提取：detcoef和detcoef2。

说明："系数提取"只有"多级分解"才会用的到！ 1级分解是不需要"系数提取"的！因为就分成了低频和高频2个部分，直接用1维或2维分解函数的返回结果就行了！所以：多级分解的系数提取，就相当于1级分解后的返回结果的直接画图。

二、Haar 变换

2.1 一维Haar变换

设原始一维数据 $a = [a_1,a_2,a_3,a_4]$ , Haar低通滤波 $[1, 1]$ , Haar高通滤波 $[- 1, 1]$ ,

则Haar小波变换为： $[\frac{a_1+a_2}{2},\frac{a_3+a_4}{2},\frac{a_1-a_2}{2},\frac{a_3-a_4}{2}]$ ，当需要进行下2采样时计算其均值（也有保留偶数序列），直接取 $b_1,b_2]$

再次进行小波变换： $[\frac{b_1+b_2}{2},\frac{b_1-b_2}{2}]$ ，下2采样为 c1 。将 $b_3,b_4]$ 称为细节系数。

因此通过Haar 变换，一幅分辨率为4的数据就可以由：分辨率为1，以及3个细节系数表示。同样由降采样的数据和细节系数可以恢复出原始数据。从上面计算过程可以看出：矢量a 与低通滤波器卷积得到近似，与高通滤波器卷积得到细节。

小波变换计算详情：请查看

2.2 二维离散小波变换

有些背景噪声, 在一维信号中体现的是随机性，但有可能在二维信号中就显示出很强的区域性或具有较明显的特征性，二维信号的直观反映就是图像！其本质为一个二维矩阵；在往本质上说它和一维离散数据信号其实内涵都是相同的。总之：二维离散小波变换处理的是二维数值离散矩阵。

二维图像Haar变换
从水平和竖直两个方向进行低通和高通滤波（水平和竖直先后不影响），用图像表述如下图所示：
在这里插入图片描述

其中：

b: 原图信息
h1 :水平方向的细节（高频信息），
v1 表示竖直方向的细节（高频信息），
c1表示对角线方向的细节（高频信息）

A是低频信息，H是水平高频信息，V是垂直高频信息、D是对角高频信息。

假设一张图片只有4个像素（a,b,c,d），其经过2-D DWT之后得到4张子图，每个子图的详细计算过程如下：
在这里插入图片描述

2.3 pywt 小波以及详情

import pywt
for family in pywt.families():
     print("%s family: " % family + ', '.join(pywt.wavelist(family)))
========================
haar family: haar
db   family: db1, db2, db3, db4, db5, ... , db35, db36, db37, db38
sym  family: sym2, sym3, sym4, sym5, ... , sym15, sym16, sym17, sym18, sym19, sym20
coif family: coif1, coif2, coif3, coif4, coif5, ..., 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

小波函数	$\mathbf{Haar(haar)}$	Daubechies (db)	Biorthogonal (bior)	Coiflets (coif)	Symlets (sym)	Morlet (morl)	Mexican Hat (mexh)
表示形式	haar	$\mathrm{db} N$	$\mathrm{bior}Nr.Nd$	$\mathrm{coif}N$	$\mathrm{sym} N$	morl	mexh
正交性	有	有	无	有	有	无	无
双正交性	有	有	有	有	有	无	无
紧支撑性	有	有	有	有	有	无	无
连续小波变换	可以	可以	可以	可以	可以	可以	可以
离散小波变换	可以	可以	可以	可以	可以	不可以	不可以
支撑长度	1	2N-1	重构: $\textit{2Nr+1}$ 分解: $\textit{2Nd+1}$	6N-1	2N-1	有限长度	有限长度
滤波器长度	2	2N	$\mathrm{Max}(2Nr,2Nd)+2$	6N	2N	[-4, 4]	[-5, 5]
对称性	对称	近似对称	不对称	近似对称	近似对称	对称	对称
小波函数消失矩阶数	1	N	Nr-1	2N	N	-	-
尺度函数消失矩阶数	-	-		2N-1	-	-	-

小波函数	Gaus (gaus)	Dmeyer (dmey)	ReverseBior (rbioNr.Nd)	Cgau (cgau)	Cmor (cmor)	Fbsp (fbsp)	Shan (shan)
表示形式	$\mathrm{gaus} N$	dmey	$\mathrm{rbio}Nr.Nd$	$\mathrm{cgau} N$	cmor	fbsp	shan
紧支撑正交性	无	无	无	无	无	无	无
紧支撑双正交性	无	无	有	无	无	无	无
连续小波变换	可以	不可以	可以	不可以	不可以	不可以	不可以
离散小波变换	不可以	可以	可以	不可以	不可以	不可以	不可以
对称性	对称	对称	对称	对称	对称	对称	对称
小波函数消失矩阶数	-	-	-	-	-	-	-
尺度函数消失矩阶数	-	-	Nr-1	-	-	-	-

详情了解：¹

Haar (haar)
Haar函数是小波分析中最早用到的一个具有紧支撑的正交小波函数，也是最简单的一个小波函数，它是支撑域在t∈[0,1]范围内的单个矩形波。Haar小波在时域上是不连续的，所以作为基本小波性能不是特别好。
Daubechies (db) 　多贝西小波小波函数 $Ψ (t)$ 和尺度函数 $φ (t)$ 中的支撑区为 $2 N - 1$ ， $Ψ (t)$ 的消失矩为N。 $\mathrm{db} N$ 小波具有较好的正则性，即该小波作为稀疏基所引入的光滑误差不容易被察觉，使得信号重构过程比较光滑。 $\mathrm{db} N$ 小波的特点是随着阶次（序列N）的增大消失矩阶数越大，其中消失矩越高光滑性就越好，频域的局部化能力就越强，频带的划分效果越好，但是会使时域紧支撑性减弱，同时计算量大大增加，实时性变差。另外，除N=1外， $\mathrm{db} N$ 小波不具有对称性（即非线性相位），即在对信号进行分析和重构时会产生一定的相位失真。 $\mathrm{db} N$ 没有明确的表达式（除了N=1外，N=1时即为Haar小波）
Symlets (sym) 或 $\mathbf{sym}N$
Symlet小波函数是IngridDaubechies提出的近似对称的小波函数，它是对db函数的一种改进。Symlet小波系通常表示为 $\mathrm{sym}N$ (N=2,3,…,8)。 $\mathrm{sym}N$ 小波的支撑范围为2N-1，消失矩为N，同时也具备较好的正则性。该小波与dbN小波相比，在连续性、支集长度、滤波器长度等方面与dbN小波一致，但 $\mathrm{sym}N$ 小波具有更好的对称性，即一定程度上能够减少对信号进行分析和重构时的相位失真。
Coiflets (coif)
根据R.Coifman的要求，Daubechies构造了Coiflet小波，它具有coifN $(N = 1, 2, 3, 4, 5)$ 这一系列。Coiflet的小波函数 $Ψ (t)$ 的2N阶矩为零，尺度函数 $φ (t)$ 的 $2 N - 1$ 阶矩为零。 $Ψ (t)$ 和 $φ (t)$ 的支撑长度为 $6 N - 1$ 。Coiflet的 $Ψ (t)$ 和 $φ (t)$ 具有比 $\mathrm{db}N$ 更好的对称性。
Biorthogonal (bior)
为了解决对称性和精确信号重构的不相容性，引入了双正交小波，称为对偶的两个小波分别用于信号的分解和重构。双正交小波解决了线性相位和正交性要求的矛盾。由于它有线性相位特性，所以主要应用在信号与图像的重构中。通常的用法是采用一个函数进行分解，用另外一个小波函娄进行重构。

双正交小波与正交小波的区别在于正交小波满足 $Ψ_{j,k},Ψ_{l,m}>=δ_{j,k}δ_{l,m}$ ，也就是对小波函数的伸缩和平移构成的基函数完全正交，而双正交小波满足的正交性为 $Ψ_{j,k} ,Ψ_{l,m}>=δ_{j,k}$ ，也就是对不同尺度伸缩下的小波函数之间有正交性，而同尺度之间通过平移得到的小波函数系之间没有正交性，所以用于分解与重构的小波不是同一个函数，相应的滤波器也不能由同一个小波生成。

该小波虽然不是正交小波，但却是双正交小波，具备正则性，同时也是紧支撑的，其重构支撑范围为 $2 N r + 1$ ，分解支撑范围为 $2\mathit{Nd+1}$ 。 $b i o r N r . N d$ 小波的主要特征表现在具有线性相位特性。一般来说为了获得线性相位，需要降低对于正交性的局限，为此该双正交小波降低了对于正交性的要求，保留了正交小波的一部分正交性，使小波攻得了线性相位和较短支集的特性。
Reverse biorthogonal (rbio)
“Discrete” FIR approximation of Meyer wavelet (dmey)
Dmeyer即离散的Meyer小波，它是Meyer小波基于FIR的近似，用于快速离散小波变换的计算。
Gaussian wavelets (gaus)
Gaussian小波是高斯密度函数的微分形式，它是一种非正交与非双正交的小波，没有尺度函数。
Mexican hat wavelet (mexh)
Mexican Hat函数为Gauss函数的二阶导数。因数它的形状像墨西哥帽的截面，所以我们称这个函数为墨西哥草帽函数。它在时域和频率都有很好的局部化，但不存在尺度函数，所以此小波函数不具有正交性。

$\mathrm{mexh}(x) = c*\exp(-\frac{x^2}{2})*(1-x^2)\\~ \\ c = \frac{2}{\sqrt{3}*\pi^{1/4}}$
Morlet wavelet (morl)
Morlet小波是高斯包络下的单频率正弦函数，没有尺度函数，是非正交分解。
$\mathrm{morl}(x) = \exp(-\frac{x^2}{2})*\cos(5x)$
Complex Gaussian wavelets (cgau)
属于一类复小波，没有尺度函数。
Shannon wavelets (shan)
Frequency B-Spline wavelets (fbsp)
样条函数（splinefunction）指一类分段（片）光滑、并且在各段交接处也有一定光滑性的函数，简称样条。
Complex Morlet wavelets (cmor)
Morlet小波是一种单频复正弦调制高斯波，也是最常用的复值小波该小波，在时频两域均具有良好的分辨率，将此小波加以改造特别适用于地震资料的分析。

效果展示：https://blog.csdn.net/HanFeiKei/article/details/83584963

三、代码演示

文档查看 ²
代码演示 ³

import pywt


db3 = pywt.Wavelet('db3')  # 创建一个小波对象
print(db3)
"""
Family name:    Daubechies
Short name:     db
Filters length: 6             #滤波器长度
Orthogonal:     True          #正交
Biorthogonal:   True          #双正交
Symmetry:       asymmetric    #对称性，不对称
DWT:            True          #离散小波变换
CWT:            False         #连续小波变换
"""

def dwt_and_idwt():
    '''
    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]
    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, cD)
    print(pywt.Modes.modes)
    # [ 4.38354585  3.80302657  7.31813271 -0.58565539  4.09727044  7.81994027]
    # [-1.33068221 -2.78795192 -3.16825651 -0.67715519 -0.09722957 -0.07045258]
    # ['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]


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

    # 访问小波包的子节点
    # 第一层：
    print(wp['a'].data)       # [ 4.52111203  1.54666942  2.57019338  5.3986205   8.20681003 11.18125264]
    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)
    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'] = 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)])
    print([node.path for node in new_wp.get_leaf_nodes(decompose=True)])
    # ['aa', 'ad', 'd']
    # ['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.]

    assert wp.a == wp["a"]
    print(wp["a"])
    # a: [ 2.12132034  4.94974747  7.77817459 10.60660172]


if __name__ == '__main__':
    dwt_and_idwt()
    wavelet_packets()

haar 演示

import numpy as np
import pywt
import cv2


raw_img = cv2.imread("./len_std.jpg")
coeffs = pywt.dwt2(raw_img, 'haar')
cA, (cH, cV, cD) = coeffs           # 低频分量，(水平高频，垂直高频，对角线高频)
re_raw_img = pywt.idwt2(coeffs, 'haar')
cv2.imshow("re_raw_img", re_raw_img.astype(np.uint8))
cv2.waitKey(0)
cv2.destroyAllWindows()

#==============================================
import numpy as np
import pywt
from PIL import Image


img = Image.open("len_std.jpg")
imgarr = np.array(img)
coeffs = pywt.dwt2(imgarr, 'haar')
re_img = pywt.idwt2(coeffs, 'haar')

import numpy as np
import pywt
import cv2
import matplotlib.pyplot as plt


def haar_img():
    img_u8 = cv2.imread("len_std.jpg")
    img_f32 = cv2.cvtColor(img_u8, cv2.COLOR_BGR2GRAY).astype(np.float32)

    plt.figure('二维小波一级变换')
    coeffs = pywt.dwt2(img_f32, 'haar')
    cA, (cH, cV, cD) = coeffs

    # 将各个子图进行拼接，最后得到一张图
    AH = np.concatenate([cA, cH], axis=1)
    VD = np.concatenate([cV, cD], axis=1)
    img = np.concatenate([AH, VD], axis=0)
    return img

if __name__ == '__main__':
    img = haar_img()

    plt.imshow(img, 'gray')
    plt.title('img')
    plt.show()

左图为img，右图为三个高频+255后显示的图片。查看更多示例
在这里插入图片描述

Daubechies （db2）演示

import numpy as np
import pywt
from matplotlib import pyplot as plt
from pywt._doc_utils import wavedec2_keys, draw_2d_wp_basis

x = pywt.data.camera().astype(np.float32)
shape = x.shape
max_lev = 3       # 要绘制多少级分解
label_levels = 3  # 图上要显式标注多少层

fig, axes = plt.subplots(2, 4, figsize=[14, 8])
for level in range(0, max_lev + 1):
    if level == 0:
        # 显示原始图像
        axes[0, 0].set_axis_off()
        axes[1, 0].imshow(x, cmap=plt.cm.gray)
        axes[1, 0].set_title('Image')
        axes[1, 0].set_axis_off()
        continue

    # 绘制标准DWT基的子带边界
    # draw_2d_wp_basis(shape, wavedec2_keys(level), ax=axes[0, level], label_levels=label_levels)
    # axes[0, level].set_title('{} level\ndecomposition'.format(level))

    # 计算二维 DWT
    c = pywt.wavedec2(x, 'db2', mode='periodization', level=level)

    # 独立规范化每个系数数组以获得更好的可见性
    c[0] /= np.abs(c[0]).max()
    for detail_level in range(level):
        c[detail_level + 1] = [d / np.abs(d).max() for d in c[detail_level + 1]]

    # 显示归一化系数 (normalized coefficients)
    arr, slices = pywt.coeffs_to_array(c)
    axes[1, level].imshow(arr, cmap=plt.cm.gray)
    axes[1, level].set_title('Coefficients\n({} level)'.format(level))
    axes[1, level].set_axis_off()
plt.tight_layout()
plt.show()

在这里插入图片描述

Biorthogonal （bior）演示

import numpy as np
import matplotlib.pyplot as plt
import pywt
import pywt.data


def test_pywt():
    original = pywt.data.camera()

    # 图像的小波变换，图形的逼近和细节
    titles = ['Approximation', ' Horizontal detail', 'Vertical detail', 'Diagonal detail']
    coeffs2 = pywt.dwt2(original, 'bior1.3')
    LL, (LH, HL, HH) = coeffs2
    plt.imshow(original)
    plt.colorbar(shrink=0.8)
    fig = plt.figure(figsize=(12, 3))
    for i, a in enumerate([LL, LH, HL, HH]):
        ax = fig.add_subplot(1, 4, i + 1)
        ax.imshow(a, interpolation="nearest", cmap=plt.cm.gray)
        ax.set_title(titles[i], fontsize=10)
        ax.set_xticks([])
        ax.set_yticks([])

    fig.tight_layout()
    plt.show()


if __name__ == '__main__':
    test_pywt()

在这里插入图片描述

鸣谢：
一维离散小波变换：
https://www.jianshu.com/p/fbbc8573f7ed#
https://www.jianshu.com/p/56733f6c0a10

SongpingWang

关注

127
点赞
踩
650

收藏

觉得还不错? 一键收藏
打赏
4
评论
Python 离散小波变换（DWT） pywt库

文章目录一、小波变换离散小波变换函数二、Haar 变换2.1 一维Haar变换2.2 二维离散小波变换三、代码演示简便安装：pip install PyWavelets 或者 conda install PyWavelets源码安装：下载源码tar zxvf PyWavelets-1.1.1.tar.gz PyWavelets-1.1.1cd PyWavelets-1.1.1sudo python3 setup.py buildsudo python3 setup.py install使用文
复制链接

扫一扫