通俗易懂理解小波变换(Wavelet Transform)

重要说明：本文从网上资料整理而来，仅记录博主学习相关知识点的过程，侵删。

一、参考资料

小波教程原文：The Wavelet Tutorial

小波教程第二版：THE WAVELET TUTORIAL SECOND EDITION

bilibili视频：【一起学小波】01 小波理论介绍及其应用

Valens C. A really friendly guide to wavelets[J]. 1999.

Wavelets: Seeing the Forest and the Trees

《The Wavelet Tutorial》小波教程 中文翻译（上）

形象易懂讲解算法I——小波变换

【2DWT：2维离散小波变换（附Pytorch代码）】

二、相关介绍

1. 常见变换

尽管傅里叶变换(FT)可能是目前使用的最普遍的变换(尤其是在电气工程领域)，但它不是唯一的一种。还有许多其他的变换经常被工程师和数学家使用，例如：希尔伯特变换(Hilbert transform)、短时傅里叶变换(short-time Fourier transform，STFT)、维格纳分布(Wigner distributions)、拉东变换(Radon Transform)、小波变换，这些仅仅是工程师和数学家可以使用的众多变换中的一小部分。每一种变换技术都有其各自的应用领域，各有优缺点。

关于 傅里叶变换 和 短时傅里叶变换 的详细介绍，请查阅另一篇博客：通俗易懂理解傅里叶变换(Fourier Transform)

2. 正交化

什么是正交化？

为什么说小波变换的优势是能实现正交化?

简单说，如果采用正交基，变换域系数会没有冗余信息，变换前后的信号能量相等，相当于用最少的数据表达最大的信息量，有利于数值压缩。JPEG2000压缩就是使用正交小波变换。

比如典型的正交基：二维笛卡尔坐标系的(1, 0)、(0, 1)，用它们表达一个信号显然非常高效，计算简单。而如果用三个互成120°的向量表达，则会有信息冗余，有重复表达。但是并不意味着正交一定优于不正交。比如如果是做图像增强，有时候反而希望能有一些冗余信息，更利于对噪声的抑制和对某些特征的增强。

3. 正交小波(Orthogonal wavelets)

3.1 正交小波的定义

正交小波(Orthogonal wavelets)、半正交小波、双正交小波的概念，文献[1]中有解释如下：

3.2 正交小波示例：Daubechies(dbN)

Daubechies(dbN)小波是正交的，信号空间$L^2(x)$ 的一系列正交基可以由它的尺寸函数$\Phi (x)$和小波函数$\Psi (x)$推导出。Daubechies小波有一个近似阶级参数 $p$，并且它的滤波器长度为 $2p$。表格 $\mathrm{I}$ 是 $p$ 阶Daubechies小波的低通滤波器 l = { l k } \boldsymbol{l}=\{l_{k}\} l={lk​} ，其中 $1\le p\le 6$。

已知低通滤波器，可以计算高通滤波器 h = { h k } \boldsymbol{h}=\{h_k\} h={hk​} ，其推导公式如下：
$$h_k=(-1{)}^k{l}_{\mathcal{N}-k}$$
其中， $\mathcal{N}$ 是一个奇数，特别的，当$p=1$时，Daubechies(1)是一个Haar小波。

3.3 双正交小波示例：Cohen小波

Cohen小波是对称双正交小波，每一个都与尺度函数$\Phi (x)$和小波函数$\Psi (x)$以及对偶函数$\tilde{\Phi },\tilde{\Psi }$有关，它对应的有四个滤波器，分别为 $l$ ，$h$， $\tilde{l}$ 和 $\tilde{h}$。当使用滤波器 $l$ 和 $h$ 进行DWT分解信号时，可以使用对偶滤波器 $\tilde{l}$ 和 $\tilde{h}$ 对该信号进行IDWT重构。Cohen小波有两个阶级参数 $p$ 和 $\tilde{p}$。表格 $\text{II}$ 是$p$ 和 $\tilde{p}$ 阶Cohen小波的低通滤波器，其中 $2\le p=\tilde{p}\le 5$。

已知低通滤波器，可以计算高通滤波器，其推导公式如下：
$$\begin{array}{l}{h}_k=(-1{)}^k{\tilde{l}}_{\mathcal{N}-k},\\ {\tilde{h}}_k=(-1{)}^k{l}_{\mathcal{N}-k},\end{array}$$
其中， $\mathcal{N}$ 是一个奇数，特别的，当$p=1,\tilde{p}=1$时，Cohen(1,1)是一个Haar小波。

3.4 其他小波

除了正交(othogonal wavelets)和双正交小波(biorthogonal wavelets)，还有各种小波和超小波(beyond-wavelets)，包括多小波(multi-wavelets)[2]、双树复小波(dual-tree complex wavelet)[3]、脊波(ridgelet)[4]、曲波(curvelet)[5]、带波(bandelet)[6]、轮廓波(contourlet)[7]等。

三、小波变换(WT)相关介绍

1. 引言

傅里叶变换给出了信号的频域信息，意味着告诉我们信号中每个频率有多少，但它没有告诉我们这些频谱成分在什么时候存在（when in time）。对于平稳(stationary)信号是不需要这个信息的。

当需要对频谱成分进行时间定位(time localization)时，需要对信号进行“时间-频率”表示的变换。短时傅里叶变换通过添加时间窗口(时间窗)来处理非平稳信号，但时间窗的宽度有时不能精准确定，所以基于小波变换的方法被提出。

2. 小波变换简介

2.1 小波变换的发展历史

“小波”的概念是由法国地质学家J.Morlet在上个世纪八十年代(1984年)在研究地下岩石油层分布时提出的，并成功应用于地质数据处理中；其后数学家Meyer创造性地构造了第一个具有一定衰减性质的光滑小波；1987年，Matllat提出了多尺度分析思想和Mallat算法（Mallat算法在离散小波变换中的地位相当于FFT在Fourier变换中的地位），成功地统一了在此之前提出的各种具体小波函数的构造。

2.2 通俗理解小波变换

小波变换是作为STFT的一种替代方法而发展起来的，但是小波变换的出发点和STFT不同。STFT是给信号加窗，分段做FFT；而小波变换直接替换了傅里叶变换的基，将无限长的三角函数基换成了有限长的会衰减的小波基。

小波变换继承和发展了STFT的局部化思想，它克服了傅里叶变换的局限性以及STFT窗口不变的缺点。小波变换主要通过伸缩和平移实现多尺度细化，突出所要处理的问题细节，有效提取局部信息。小波变换给出了一个可以调节的时频窗口，窗口的宽度随频率变化。频率增高时，时间窗口的宽度自动变窄，以提高分辨率，“采用小波分析，就像使用一架可变焦距镜头的照相机一样，可以转向任一细节部分”。

2.3 小波(Wavelet)的概念

“小波”名称，是由英文翻译而来的。wave表示为波，let是英文常见的后缀，表示小的意思，合起来翻译为“小波”。

droplet水滴

leaflet叶子，小树叶

booklet小册子

什么是小波呢？所谓小波就是小的波形，“小”即具有衰减性，“波”是指具有波动性。小波是一种将能量集中在时域的波，可以对其进行伸缩和平移。

傅里叶变换就是基于这些“波”（正弦和余弦），它从负无穷到正无穷都存在。而小波变换则是基于以下的这种“小波”：

这是一个小波基函数，比较它和前面的“波”也就可以知道为什么称它为小波（只存在于一段时间内为非零值）。

上图中的三个小波是尺度变换的，下图中的三个小波是时移变换的：

下面把“波”与“小波”作一个对比：

3. 小波基函数

小波基函数，简称小波基，又称为小波母函数。小波基函数经过尺度变换、时移 等，可得到一系列小波函数，简称小波。

设 $\Psi (\mathrm{t})$ 为可积函数，则小波基函数定义为：
$${\Psi }_{a,\tau }\left(t\right)=\frac{1}{\sqrt{\left|a\right|}}\Psi \left(\frac{t-\tau }{a}\right)$$
其中，$(a,\tau )$ 为任意实数对，参数$a$必须为非零实数。一般称参数 $a$ 为尺度因子（即STFT中的 $\omega$），负责调节小波函数的伸缩程度，与频率信息相对应。而参数 $\tau$ 为平移因子或时间中心参数（即STFT中的 $t$），负责调节小波函数的平移程度，与时间信息相对应。

小波基函数 $\Psi (\tau )$ 只有在原点附近才会有明显偏离水平轴的波动，在远离原点的地方其函数值将迅速衰减为零。所以对于任意参数 $(a,\tau )$，小波函数 ${\Psi }_{a,\tau }\left(t\right)$ 在 $t=\tau$ 附近存在明显的波动，远离 $t=\tau$的地方迅速衰减到零。小波基函数趋向于零的速度是衡量小波基函数性质好坏的一个重要标志。

关于常用小波基函数的介绍，请查阅另一篇博客：小波变换之常见的小波基介绍

4. 常用的小波变换

连续小波变换、离散小波变换、二进小波变换、离散序列的小波变换、小波包

学习小波变换时会发现有一大堆形形色色的小波变换，这是因为小波变换与小波基函数 $\Psi (\mathrm{t})$ 有关，不同的小波基函数 $\Psi (\mathrm{t})$ 则对应着不同种类的小波变换，比如Haar小波、Daubechies(dbN)小波等等。

可以证明， ${\Psi }_{a,\tau }\left(t\right)$ 的时频窗口面积与参数 $a$ 和 $\tau$ 无关，仅与 $\Psi (\mathrm{t})$ 的选取有关，所以不能通过选择参数 $a$ 和 $\tau$ 使时域和频域窗口的半径同时缩小，时域和频域上的分辨率相互牵制，要想使两者的分辨率同时提高，就必须选择适当的小波基函数 $\Psi (\mathrm{t})$ 。

常用的小波变换包括：连续小波变换(CWT)，离散小波变换(DWT)，复小波变换(Complex Wavelets)，双树复小波变换(DTCWT)。

4.1 连续小波变换(CWT)

设信号$f(t)\in L^2(R)$，其连续小波变换(Continuous Wavelet Transform, CWT)定义为：
$$W_f(a,\tau )=\frac{1}{\sqrt{\mid a\mid }}{\int }_{-\infty }^{\infty }f(t)\Psi \left(\frac{t-\tau }{a}\right)dt$$
从公式中可以看出，信号 $f(t)$ 的小波变换本质上是原来的 $f(t)$ 在$t=\tau$ 附近按 ${\Psi }_{a,\tau }(t)$ 进行加权平均，体现的是以 ${\Psi }_{a,\tau }(t)$ 为标准 $f(t)$ 的快慢变化情况。这样，参数 $\tau$ 表示分析的时间中心或时间点，而参数 $a$ 体现以 $t=\tau$ 为中心的附近范围的大小。可以看出，一维信号 $f(t)$ 经过小波变换后将变成二维信号。因此，小波变换能够同时提供时间和频率信息(时频分析)，从而给出信号的 ”时间-频率“ 表示。

4.2 离散小波变换(DWT)

当尺度因子和平移因子连续变化时必然会产生冗余数据，所以便出现了离散思想。而且信号离散化后有利于计算机进一步分析和处理。连续小波变换中的尺度因子和平移因子是连续发生变化的，如果将其离散化，就会得到离散小波变换。

离散小波变换多数情况是得到一个类似于离散傅里叶变换的变换，时域和变换域都是离散有限长的，方便计算机处理。然而，离散小波变换(Discrete Wavelet Transform, DWT)指的是将连续小波变换中的尺度因子 $a$ 和平移因子 $\tau$ 离散化。特别注意：DWT并没有将信号 $f(t)$ 和小波函数 $\Psi (t)$ 中的时间变量t离散化，这与离散傅里叶变换的概念是非常不一样的。

在尺度因子 $a$ 和平移因子 $\tau$ 离散过程中，一般对 $a$ 进行幂数级离散化，对 $\tau$ 均匀离散。考虑到不同尺度下频率不同，因此不同尺度下平移因子 $\tau$ 的离散间隔不同。

对尺度因子 $a$ 进行幂级数离散化，即：

$$a=a_0^j,a_0>0,j\in Z$$

对平移因子$\tau$进行均匀离散化，即：

$$\tau =k{a}_0^j{\tau }_0$$

经过离散化后的小波函数为：

$${\Psi }_{j,k}(t)=2^{\frac{j}{2}}\Psi (2^{-j}t-k{\tau }_0),j=0,1,2,...,k\in Z$$

可以将一维离散小波变换定义为：

$$W_f(a_0^j,k{\tau }_0)=\int f(t){\Psi }_{a_0^j,k{\tau }_0}(t)dt,j=0,1,2,...,k\in Z$$

4.3 复小波变换(Complex Wavelets)

//TODO

4.4 双树复小波变换(DTCWT)

//TODO

双树复小波变换(Dual-Tree Complex Wavelet Transform, DTCWT)

5. 小波变换的优势

5.1 小波变换与STFT对比

小波变换相比于STFT，有以下优势：

1. 由于小波基函数 ${\Psi }_{a,b}\left(t\right)$ 相当于窗函数，但其窗宽是可变的，较好地解决了时间分辨率和频率分辨率的矛盾，其变化规律使得小波变换具有优良的局部化特性，对分析突变信号和奇异信号非常有效，充分体现了常相对带宽频率分析和自适应分析的思想；
2. 小波变换能将各种交织在一起的由不同频率组成的混合信号分解成不同频率的信号，并对频率大小不同的信号采用相应粗细的时空域取样步长，从而能够不断聚焦到对象的任意微小细节，对时变信号的频谱分析意义重大。
3. 不要求小波变换基底是正交的，其时宽频宽乘积较小，因而展开系数的能量较为集中。

5.2 FT、STFT、小波变换对比

比较一下小波变换、短时傅里叶变换、傅里叶变换可以发现：

1. 傅里叶变换不具有局部性；
2. STFT有局部性，但有一些缺点（如前所述）；小波变换不但具有局部性，而且尺度参数a可以改变频谱结构和窗口的形状，起到“变焦”的作用，因此小波分析可能达到多分辨率分析的效果（小变波换被誉为数学显微镜）。
3. 从信号分析方法的理论发展过程可能看出：傅里叶分析特别适合分析长时间内较稳定的信号；STFT也有其一定的应用场合，但其效果取决于适当地选取窗函数；小波分析特别适合分析突变信号和奇异信号。

5.3 总结小波变换的优势

小波变换的优势总结如下：

1. 具有高的灵活性。小波能够应用于数字图像，在处理数字图像时可以看作是与离散有关的数学问题。

2. 可以表示局部特征。不管是在高频段还是低频段，小波都有相应的处理方式。

3. 减少冗余度。能够用较少的信息来表达图像的主要特征，减少不必要的信息。

4. 处理图像的精度比较高。图像是由一个个像素点组成的，通过处理像素点之间的距离和灰度量化数量级来提高精度。

6. 小波变换的应用

小波理论适用于有限(finite filters)或无限滤波器(infinite filters)，但在实际应用中很少涉及无限滤波器。
小波广泛应用于信号处理、数值分析、模式识别、计算机视觉、量子力学等领域。
一维小波变换用来处理信号，二维小波变换用来处理图像。

6.1 小波变换处理图像信号

图像信号具有非平稳特性，无法使用一种确定的数学模型来描述，而小波变换的多分辨率分析特性很好地解决了这个问题。小波变换的多分辨率特性使其既可以高效描述图像的平坦区域（低频信息、全局信息），也可以有效处理图像信号的局部突变（高频信息，即图像的边缘轮廓等部分）。小波变换在空域和频域同时具有良好的局部性，使其可以很好地聚焦到图像的任意细节。

6.2 小波变换处理突变信号

四、小波变换(Wavelet Transform)原理解析

1. 小波变换原理

1.1 小波变换的数学表达

小波变换将无限长的三角函数基换成了有限长的会衰减的小波基。

这就是为什么它叫“小波”，因为是很小的一个波。

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

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

1.2 小波变换示例

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

时域信号，如下图所示：

傅里叶变换结果，如下图所示：

小波变换结果，如下图所示：

2. 小波分解与重构

小波变换包括小波分解和小波重构过程，它可以将一个比较复杂的信号分解成若干个简单的信号来处理，经过处理之后的若干简单信号能够合成比原始信号更饱满的信号，这种变换实质上就是把信号从一个域转换到另一个域。

根据快速小波变换(FWT)在小波域进行二维变换，实现特征池化。

二维小波变换是指将输入的二维数据依次进行水平方向的行滤波和垂直方向的列滤波（其实先行后列或者先列后行不影响）。行滤波和列滤波称为一维小波变换。

$$W_{\phi }[j+1,k]=h_{\phi }[-n]\ast W_{\phi }[j,n]{|}_{n=2k,k\le 0}(6)$$

$$W_{\Psi }[j+1,k]=h_{\Psi }[-n]\ast W_{\Psi }[j,n]{|}_{n=2k,k\le 0}(7)$$

其中，$\phi$ 为近似函数(approximation function)，$\Psi$为细节函数(detail function)，$W_{\phi }$ 称为近似系数(approximation coefficients)，$W_{\Psi }$ 称为细节系数(detail coefficients)。$h_{\phi }[-n]$ 和 $h_{\Psi }[-n]$ 为时间逆尺度(reversed scaling)和小波向量(and wavelet vectors)，$n$ 代表向量中的样本数量，$j$ 表示分辨率(resolution level)。对图像进行两次快速小波变换(FWT)，一次用于行，另一次用于列。最终得到每个分解层的详细子带(LH, HL, HH)和最高分解层的近似子带(LL)。其中LL表示低频特征，其余带H的表示高频特征。

在进行二维分解之后，只能用二维小波子带重构图像。使用快速小波变换的逆变换(IFWT)对图像特征进行上采样，来重构一阶的小波分解。该方法基于离散小波变换的逆变换(IDWT)。
$$\begin{array}{r}{W}_{\phi }[j,k]=h_{\phi }[-n]\ast W_{\phi }[j+1,n]+h_{\Psi }[-n]\ast W_{\Psi }[j+1,n]{|}_{n=\frac{k}{2},k\le 0}\end{array}(8)$$

2.1 小波分解本质

小波分解将图像信号分离为携带不同信息的 子带(subband, 子图像)，这些子带包含的信号各不相同，频率范围也不一样。对于不同的子带，考虑使用不同的方法来处理，这样能够凸显不同尺度下的细节信息。同时，小波分解可以对图像进行更深阶的分解，如三阶和四阶小波分解。在每一次操作时，都应在近似图像上进行分解。

2.2 小波分解过程

小波变换能够将图像的信息逐层分解，而它的分解能力来源于低通和高通滤波器。

高通滤波器(high-pass filters)：信号的高频成分可以通过，低频成分被衰减。

低通滤波器(low-pass filters)：信号的低频成分可以通过，高频成分被衰减。

对原始图像每进行一次小波变换，会分解产生一个低频子带（LL：行低频、列低频）和三个高频子带(垂直子带LH：行低频、列高频；水平子带HL：行高频、列低频；对角子带HH：行高频、列高频）。小波变换可以进行多阶，后续小波变换基于上一阶低频子带LL进行，依次重复，可完成对图像的i阶小波变换，其中i=(1,2,3,…I）。如下图所示，A图、B图分别为i=1时的一阶小波变换分布，i=2时的二阶小波变换分布，每个子带分别包含各自对应的小波系数。

可以看到，其实每次小波变换可以看做对图像的行水平方向、列垂直方向分别进行隔点采样，空间分辨率每次变为$\frac{1}{2}$，因此第$i$阶小波变换后，其子带空间分辨率为原图的 $\frac{1}{2^i}$。

由于小波变换引入了下采样，因此四个子带map的边长是输入的一半。小波变换的逆变换可以将四个子带完美地重构原始输入。

2.2.1 一阶小波分解

一阶小波分解，实质上是采用可分离的滤波器对原始图像数据的水平方向和垂直方向分别进行一维小波变换。一阶小波分解过程如下：

首先，利用一维滤波器对原始图像进行水平方向的分解，得到原始图像在水平方向的低频分量L1和高频分量H1；然后，对上一步所得图像进行垂直方向的分解，便可以得到原始图像的全部分量。得到的四个分量分别是原始图像在水平方向和垂直方向上的低频分量LL1、水平方向上的低频和垂直方向上的高频分量LH1、水平方向上的高频和垂直方向上的低频分量HL1以及水平方向和垂直方向上的高频分量HH1。

其中LL表示低频特征，其余带H的表示高频特征。

• LL子带是原始图像的近似表示。它是由水平方向和垂直方向使用低通小波滤波器卷积后得到的小波系数。
• HL子带 (原始图像的水平方向细节子图)，用来凸显图像水平方向的奇异特性，它是使用高通小波滤波器在水平方向卷积后，再通过低通小波滤波器在垂直方向上卷积得到的小波系数；
• LH子带 (原始图像的垂直方向细节子图)，用来凸显图像垂直方向的奇异特性，它是使用低通小波滤波器在水平方向卷积后，再通过高通小波滤波器在垂直方向上卷积得到的小波系数；
• HH子带 (原始图像的对角线方向细节子图)，用来凸显图像的对角边缘特性，它是由水平和垂直两个方向使用高通小波滤波器卷积后得到的小波系数。

对于原始输入$I$，一阶小波分解记为：

$$LL1,LH1,HL1,HH1=DWT(I)$$

对应的逆变换(一阶小波重构)，可以表示为：

$$I=IDWT(LL1,LH1,HL1,HH1)$$

2.2.2 二阶小波分解

对LL1进行二阶小波分解，可以表示为：

$$LL2,LH2,HL2,HH2=DWT(LL1)$$

对应的逆变换(二阶小波重构)为：

$$LL1=IDWT(LL2,LH2,HL2,HH2)$$

2.3 小波分解样例

2.4 小波重构过程

小波变换是可逆的，小波重构与小波分解为互逆操作，小波分解得到的子图可通过组合重构原图。一阶小波重构过程：首先，对经过最终分解的图像的每一列进行一维离散小波重构，然后对所得图像的每一行再进行一维离散小波重构，便可以得到最终重建后的结果。其实现原理如下：

假设输入图像$I$的尺寸为M×N，且$M=2^m$、$N=2^n$，对其进行一阶小波分解过程如下：

1. 利用一维滤波器h和g分别对输入图像$I$进行水平方向的行滤波，丢弃奇数行，分解得到尺寸为$\frac{M}{2}\times N$的中间输出$I_L$和$I_H$；
2. 一维滤波器h和g分别对中间输出$I_L$和$I_H$进行垂直方向列的滤波，丢弃奇数列，分解得到尺寸为$\frac{M}{2}\times \frac{N}{2}$的输出$I_{LL}$、$I_{LH}$和$I_{HL}$、$I_{HH}$；

3. 基于哈尔(Haar)小波基的小波变换

3.1 Haar小波基的定义

Haar小波是最常用的小波基，公式定义如下：
$$\begin{array}{rl}\Psi (x)=& \{\begin{array}{llc}1& ,0\le x<0.5& \\ -1& ,0.5\le x<1& \\ 0& ,\text{other}& \end{array}\end{array}$$
其对应的尺度函数为：
$$\phi (x)=\{\begin{array}{ll}1& ,0\le x<1\\ 0& ,other\end{array}.$$
Haar小波具有最短的支集，支集长度为1，滤波器长度为2，具有正交性和对称性，其图示如下：

3.2 一维Haar小波变换过程

举例描述一维Haar小波变换过程。对于一维Haar小波变换来说，其一维高通滤波器FH=[1,-1]、一维低通滤波器FL=[1,1]，假设输入向量X[6]=[2,4,6,8,5,9]，对其进行一维Haar小波变换过程如下（图中蓝色填充表示滤波器的移动过程，黄色表示输入数据，绿色表示对应输出）：

1. 高通滤波，求相邻元素之间差值的平均值，存储输入数据的细节信息：

   

2. 低通滤波，求相邻元素的平均值，存储输入数据的粗略近似信息：

   

假设输入图像大小为M×N：

• 行滤波分别采用一维高通滤波器、一维低通滤波器得到对应的两个输出，输出大小均为$\frac{M}{2}\times N$；
• 列滤波对于行滤波的两个输出，同样采用一维高通滤波器、一维低通滤波器得到对应的四个输出，输出大小均为$\frac{M}{2}\times \frac{N}{2}$；

3.3 二维Haar小波变换过程

将一维Haar小波变换推广，进一步可得到二维Haar小波变换的实现过程，对应的四个滤波器分别为：

如下图所示，左上角ABCD四个元素构成一个局部区域，依次使用上述四个滤波器对该局部区域进行计算即可得到小波分解后对应子带中的一个元素，依次类推，计算公式如下：
$$\{\begin{array}{l}a=\frac{1}{2}(A+B+C+D)\\ b=\frac{1}{2}(-A-B+C+D)\\ c=\frac{1}{2}(-A+B-C+D)\\ d=\frac{1}{2}(A-B-C+D)\end{array}$$

五、相关经验

1. 小波分解与小波重构代码实现

MWCNNv2/MWCNN_code/model/common.py

1.1 ssim.py

计算结构相似度SSIM值的方法，可参考博客：【损失函数：2】Charbonnier Loss、SSIM Loss（附Pytorch实现）。

import torch
import torch.nn as nn
from torch.nn import functional as F
from torch.autograd import Variable
from PIL import Image
from torchvision import transforms
from math import exp


# 5.SSIM loss
# 生成一位高斯权重，并将其归一化
def gaussian(window_size,sigma):
    gauss = torch.Tensor([exp(-(x - window_size // 2) ** 2 / float(2 * sigma ** 2)) for x in range(window_size)])
    return gauss/torch.sum(gauss)  # 归一化


# x=gaussian(3,1.5)
# # print(x)
# x=x.unsqueeze(1)
# print(x.shape) #torch.Size([3,1])
# print(x.t().unsqueeze(0).unsqueeze(0).shape) # torch.Size([1,1,1, 3])

# 生成滑动窗口权重，创建高斯核：通过一维高斯向量进行矩阵乘法得到
def create_window(window_size,channel=1):
    _1D_window = gaussian(window_size, 1.5).unsqueeze(1)  # window_size,1
    # mm:矩阵乘法 t:转置矩阵 ->1,1,window_size,_window_size
    _2D_window = _1D_window.mm(_1D_window.t()).float().unsqueeze(0).unsqueeze(0)
    # expand:扩大张量的尺寸，比如3,1->3,4则意味将输入张量的列复制四份，
    # 1,1,window_size,_window_size->channel,1,window_size,_window_size
    window = Variable(_2D_window.expand(channel, 1, window_size, window_size).contiguous())
    return window


def _ssim(img1, img2, window, window_size, channel, size_average=True):
    mu1 = F.conv2d(img1, window, padding=window_size // 2, groups=channel)
    mu2 = F.conv2d(img2, window, padding=window_size // 2, groups=channel)

    mu1_sq = mu1.pow(2)
    mu2_sq = mu2.pow(2)
    mu1_mu2 = mu1 * mu2

    sigma1_sq = F.conv2d(img1 * img1, window, padding=window_size // 2, groups=channel) - mu1_sq
    sigma2_sq = F.conv2d(img2 * img2, window, padding=window_size // 2, groups=channel) - mu2_sq
    sigma12 = F.conv2d(img1 * img2, window, padding=window_size // 2, groups=channel) - mu1_mu2

    C1 = 0.01 ** 2
    C2 = 0.03 ** 2

    ssim_map = ((2 * mu1_mu2 + C1) * (2 * sigma12 + C2)) / ((mu1_sq + mu2_sq + C1) * (sigma1_sq + sigma2_sq + C2))

    if size_average:
        return ssim_map.mean()
    else:
        return ssim_map.mean(1).mean(1).mean(1)


# 构造损失函数用于网络训练或者普通计算SSIM值
class SSIM(torch.nn.Module):
    def __init__(self, window_size=11, size_average=True):
        super(SSIM, self).__init__()
        self.window_size = window_size
        self.size_average = size_average
        self.channel = 1
        self.window = create_window(window_size, self.channel)

    def forward(self, img1, img2):
        (_, channel, _, _) = img1.size()

        if channel == self.channel and self.window.data.type() == img1.data.type():
            window = self.window
        else:
            window = create_window(self.window_size, channel)

            if img1.is_cuda:
                window = window.cuda(img1.get_device())
            window = window.type_as(img1)

            self.window = window
            self.channel = channel

        return _ssim(img1, img2, window, self.window_size, channel, self.size_average)


# 普通计算SSIM
def ssim(img1, img2, window_size=11, size_average=True):
    (_, channel, _, _) = img1.size()
    window = create_window(window_size, channel)

    if img1.is_cuda:
        window = window.cuda(img1.get_device())
    window = window.type_as(img1)

    return _ssim(img1, img2, window, window_size, channel, size_average)

1.2 搭建DWT&IDWT

import torch
import torch.nn as nn
from PIL import Image
from torchvision import transforms
import matplotlib.pyplot as plt
from ssim import ssim


def dwt_init(x):
    x01 = x[:, :, 0::2, :] / 2
    x02 = x[:, :, 1::2, :] / 2
    x1 = x01[:, :, :, 0::2]
    x2 = x02[:, :, :, 0::2]
    x3 = x01[:, :, :, 1::2]
    x4 = x02[:, :, :, 1::2]
    x_LL = x1 + x2 + x3 + x4
    x_HL = -x1 - x2 + x3 + x4
    x_LH = -x1 + x2 - x3 + x4
    x_HH = x1 - x2 - x3 + x4

    return torch.cat((x_LL, x_HL, x_LH, x_HH), 1)


def iwt_init(x):
    r = 2
    in_batch, in_channel, in_height, in_width = x.size()
    # print([in_batch, in_channel, in_height, in_width])
    out_batch, out_channel, out_height, out_width = in_batch, int(
        in_channel / (r ** 2)), r * in_height, r * in_width
    x1 = x[:, 0:out_channel, :, :] / 2
    x2 = x[:, out_channel:out_channel * 2, :, :] / 2
    x3 = x[:, out_channel * 2:out_channel * 3, :, :] / 2
    x4 = x[:, out_channel * 3:out_channel * 4, :, :] / 2

    h = torch.zeros([out_batch, out_channel, out_height, out_width]).float().cuda()

    h[:, :, 0::2, 0::2] = x1 - x2 - x3 + x4
    h[:, :, 1::2, 0::2] = x1 - x2 + x3 - x4
    h[:, :, 0::2, 1::2] = x1 + x2 - x3 - x4
    h[:, :, 1::2, 1::2] = x1 + x2 + x3 + x4

    return h


class DWT(nn.Module):
    def __init__(self):
        super(DWT, self).__init__()
        self.requires_grad = False

    def forward(self, x):
        return dwt_init(x)


class IWT(nn.Module):
    def __init__(self):
        super(IWT, self).__init__()
        self.requires_grad = False

    def forward(self, x):
        return iwt_init(x)

1.3 小波分解

lenna标准图像：http://www.lenna.org/len_std.jpg

# 小波分解
dwt_module=DWT()
x=Image.open('./lena_std.tif')
# x=Image.open('./mountain.png')
x=transforms.ToTensor()(x)
x=torch.unsqueeze(x,0)
x=transforms.Resize(size=(256,256))(x)
subbands=dwt_module(x)

title=['LL','HL','LH','HH']

plt.figure()
for i in range(4):
    plt.subplot(2,2,i+1)
    temp=torch.permute(subbands[0,3*i:3*(i+1),:,:],dims=[1,2,0])
    plt.imshow(temp)
    plt.title(title[i])
    plt.axis('off')
plt.show()

小波分解结果

1.4 小波重构

# 小波重构
title=['Original Image','Reconstruction Image']
reconstruction_img=IWT()(subbands).cpu()
ssim_value=ssim(x,reconstruction_img)  # 计算原图与重构图之间的结构相似度
print("SSIM Value:",ssim_value) # tensor(1.)
show_list=[torch.permute(x[0],dims=[1,2,0]),torch.permute(reconstruction_img[0],dims=[1,2,0])]

plt.figure()
for i in range(2):
    plt.subplot(1,2,i+1)
    plt.imshow(show_list[i])
    plt.title(title[i])
    plt.axis('off')
plt.show()

输出结果

SSIM Value: tensor(1.0000)

小波重构结果

从视觉效果或者结构相似度值来看，小波变换整个过程是封闭的、无损的。

2. 基于PyTorch的2D小波变换

原始论文：Cotter F. Uses of complex wavelets in deep convolutional neural networks[D]. , 2020.

github代码：2D Wavelet Transforms in Pytorch

文档：Welcome to Pytorch Wavelets’s documentation!

pip install pywavelets

3. 复小波在DCNN中的使用

论文：《Uses of complex wavelets in deep convolutional neural networks》

七、参考文献

[1] 董新洲，贺家李，葛耀中.小波变换：第2讲离散小波变换[J].1999，27(2)：57-60.

[2] Keinert F. Wavelets and multiwavelets[M]. CRC Press, 2003.

[3] Selesnick I W, Baraniuk R G, Kingsbury N C. The dual-tree complex wavelet transform[J]. IEEE signal processing magazine, 2005, 22(6): 123-151.

[4] Ma J, Fenn M. Combined complex ridgelet shrinkage and total variation minimization[J]. SIAM Journal on Scientific Computing, 2006, 28(3): 984-1000.

[5] Candes E, Demanet L, Donoho D, et al. Fast discrete curvelet transforms[J]. multiscale modeling & simulation, 2006, 5(3): 861-899.

[6] Le Pennec E, Mallat S. Image compression with geometrical wavelets[C]//Proceedings 2000 international conference on image processing (cat. No. 00CH37101). IEEE, 2000, 1: 661-664.

[7] Do M N, Vetterli M. The contourlet transform: an efficient directional multiresolution image representation[J]. IEEE Transactions on image processing, 2005, 14(12): 2091-2106.