Wavelet Transform.

本文目录：

傅里叶变换的局限性
短时傅里叶变换
连续小波变换
离散小波变换
深度学习中小波变换的应用

1. 傅里叶变换的局限性

变换 (transform) 是指将信号从一个坐标系转换到另一个坐标系的方法，通过坐标系转换能够从信号中获取原始信号中难以获取的额外信息。变换的关键是找到一组合适的基函数（描述坐标系的坐标轴）。

自然界中的大多数信号都是原始格式的时域信号（时间的函数）。信号图的自变量是时间，因变量通常是幅度。在这种坐标系下绘制的时域信号图是信号的时间-幅度表示，有时这种表示并不总是最佳的。另一种常见的坐标系是频谱，即描述信号的频率-幅度表示。

傅里叶变换 (Fourier Transform)是一种将信号从时域表示转换为频域表示的方法。FT采用复指数函数（正弦和余弦）作为新的基函数，得到了信号的频域表示（信号中每个频率的分量占比）：

\[\begin{aligned} X(f) &= \int_{-\infty}^{\infty} x(t) e^{-2\pi ift} \mathrm{d}t \\ &= \int_{-\infty}^{\infty} x(t) \left(\cos(2 \pi ft) -i \sin(2 \pi ft)\right) \mathrm{d}t \end{aligned}\]

FT的过程可以看作是一种“相关性”的计算：计算信号与频率为$f$的正弦(余弦)项的相关性（相乘并在所有时间上积分），如果信号具有频率$f$的高振幅分量，则该分量与正弦项将重合，它们的乘积将为一个（相对）较大的值，这表明信号具有频率$f$的主要分量；反之如果频率$f$不是信号的主要分量，则乘积将为一个（相对）较小的值。

然而FT不适合处理非平稳信号（unstationary signal）。如下面这个信号是一个平稳信号，因为它在任意给定时刻的频率分别为 $10$、$25$、$50$ 和 $100$ Hz：

\[x(t) = \cos(2 \pi \cdot 10 t) + \cos(2 \pi \cdot 25 t) + \cos(2 \pi \cdot 50 t) + \cos(2 \pi \cdot 100 t)\]

该函数及其FT频谱为：

非平稳信号是指频率随时间变化的信号。比如构造频率随时间减少的chirp信号：$0$ 至 $300$ 毫秒的时间间隔具有 $100$ Hz 的正弦波，$300$ 至 $600$ 毫秒的时间间隔具有 $50$ Hz 的正弦波，$600$ 至 $800$ 毫秒的时间间隔具有 $25$ Hz 的正弦波，最后 $800$ 至 $1000$ 毫秒的时间间隔具有 $10$ Hz 的正弦波。该函数及其FT频谱为：

上述两个函数的频谱差异主要体现在：

非平稳信号的频谱有很多噪声：是由于频率成分之间的突变引起的。
非平稳信号高频成分的振幅高于低频成分的振幅。这是因为高频成分持续时间更长。

除此之外，两个函数的频谱都具有四个相同位置的峰值；这是问题所在：在平稳信号中所有频率成分在整个信号持续时间内都存在，而在非平稳信号中不同的频率分量出现在不同的信号区域中！

根据海森堡测不准原理，我们无法同时准确观测一个信号的频率和时域。FT以完全抛弃了时域信息的观测为代价（即时域分辨率为$0$），获得了完全精确的频域信息。如果我们想要同时对时域和频域进行观测，比如想知道在什么时间间隔出现什么频谱成分，则不能采用FT。

2. 短时傅里叶变换 Short Term Fourier Transform

为了使FT适用于非平稳信号，我们可以把整个时域过程分解成无数个等长的小过程，每个小过程近似平稳。STFT通过在信号上滑动一个窗口函数$\omega$，并在每个窗口上计算FT，从而将信号分解为一系列频率成分随时间变化的信号。

\[STFT_X^{(\omega)}(t',f) = \int_t \left[ x(t) \cdot \omega^*(t - t') \right] \cdot e^{-2 \pi ift} dt\]

例如，对于一个包含四个不同时刻的频率分量的非平稳信号（$0$ 到 $250$ ms 的区间是一条 $300$ Hz 的简单正弦波，其余 $250$ ms 区间分别是 $200$ Hz、$100$ Hz 和 $50$ Hz 的正弦波），其对应的STFT时频图为：

上图中关于频率的对称性是因为实信号的FT始终是对称的。除此之外，STFT不仅显示出四个峰值（分别对应四个不同的频率分量），并且这四个峰值在时间轴上位于不同的时间间隔。

由于测不准原理，我们无法知道信号的精确时频表示，即我们无法知道在什么时刻存在哪些频谱分量，只能知道特定频带存在的时间间隔。STFT的分辨率与所用窗函数的宽度有关。窗函数的宽度被称为窗函数的支持度 (support)。如果窗函数较窄，则称为紧支持函数 (compactly supported)。

如果使用无限长的窗口，则STFT等价于FT，可以提供完美的频率分辨率，但没有时间信息。如果使用使用一个足够短的窗口，使信号在其中保持平稳，则窗口越窄，时间分辨率越好，平稳性假设也更佳，但频率分辨率越差。

下面给出了四个具有不同宽度的高斯窗口函数$e^{-a \left( \frac{t^2}{2} \right)}$，其对应的STFT时频图也列在下面。

如图所示，随着窗函数支持度的增大，频域分辨率逐渐提高（不同频率分量的频率间隔更明确），但时域分辨率逐渐降低（不同频率分量的时间重合度更大）。

因此使用STFT的主要困难在于窗函数的选择：窄窗口提供良好的时间分辨率，但频率分辨率较差。宽窗口提供良好的频率分辨率，但时间分辨率较差；此外，宽窗口可能违反平稳性条件。

3. 连续小波变换 Continuous Wavelet Transform

（1）连续小波变换的定义

一个对STFT自然的改进是多分辨率分析 (MRA)，即使用宽度可变的窗函数：在信号的高频部分使用窄窗提供良好的时间分辨率和较差的频率分辨率，在信号的低频部分使用宽窗提供良好的频率分辨率和较差的时间分辨率。

连续小波变换（CWT）是STFT的替代方法，旨在克服分辨率问题。STFT 和 CWT 之间有两个主要区别：

由于CWT没有对加窗信号进行傅里叶变换，因此将看到与正弦波相对应的单峰，即没有计算负频率。
在对每个频谱分量进行变换计算时，CWT窗口的宽度都会发生变化。

连续小波变换定义为：

\[CWT_x^\psi(\tau,s) = \Psi_x^\psi(\tau,s) = \frac{1}{\sqrt{|s|}} \int x(t) \psi^* \left( \frac{t - \tau}{s} \right) dt\]

变换后的信号是两个变量的函数：平移(translation)参数$\tau$和尺度(scale)参数$s$，其中$s=\frac{1}{f}$是频率的倒数。

$\psi(t)$ 是小波函数（小波基），称为母小波(mother wavelet)。其中母是指该函数是生成其他窗口函数的原型，变换过程中使用的具有不同支撑区域的函数都源自该函数；小是指该窗口函数的宽度较小（紧支撑）；波是指该函数是振荡的。

连续小波变换的计算过程如下：

将小波$\psi(t)$与原始函数$x(t)$在初始时间段位置进行比较，计算两个函数的相关程度；
每次将小波右移$\tau$，得到小波函数$\psi(t-\tau)$，计算每次平移所得的相关系数；
根据尺度参数扩展小波$\psi(\frac{t}{s})$；
重复上述步骤，计算完所有尺度下的小波变换系数。

下图为一个由 $30$ Hz、$20$ Hz、$10$ Hz 和 $5$ Hz 四个频率分量组成的非平稳信号及其连续小波变换的结果。注意到较小的尺度对应较高的频率。结果表明WT 在高频时具有良好的时间分辨率和较差的频率分辨率，在低频时具有良好的频率分辨率和较差的时间分辨率。

定义一个可接受性常数(admissibility constant) $C_\psi$，满足以下可接受性条件(admissibility condition)时，连续小波变换是可逆的：

\[C_\psi = \left\{ 2 \pi \int_{-\infty}^{\infty} \frac{|\hat{\psi}(\xi)|^2}{|\zeta|} d\xi \right\} ^{\frac{1}{2}} < \infty\]

其中$\hat{\psi}(\xi)$ 是 $\psi(t)$ 的傅里叶变换。此时逆变换为：

\[x(t) = \frac{1}{C_\psi^2} \int_s \int_\tau \left[ \Psi^\psi_x(\tau, s) \frac{1}{s^2} \psi \left( \frac{t - \tau}{s} \right) \right] d\tau \cdot ds\]

可接受性条件的一个充分条件是：

\[\hat{\psi}(0) = 0 \to \int \psi(t) \cdot dt = 0\]

要满足上述公式，小波必须是振荡的。

（2）小波级数 Wavelet Series

在实践中，在所有平移参数$\tau$和尺度参数$s$下计算连续小波变换计算量大、且存在数据冗余性。因此可以对连续小波变换进行离散化，即小波级数（半离散小波变换）。值得一提的是，小波级数所处理的信号仍然是连续的，只是对平移参数$\tau$和尺度参数$s$的离散化。

尺度离散化为 $\boldsymbol{s = s_0^j}$，平移离散化为 $\boldsymbol{\tau = k \cdot s_0^j \cdot \tau_0}$，其中 $\boldsymbol{s_0>1},\boldsymbol{\tau_0>0}$。比如离散尺度以 $2$ 为因子发生变化，则只计算$s=2,4,8,16…$的尺度，并且尺度越大（频率越低），平移参数$\tau$的间隔越大（相当于降低时间轴上的采样率）。

通过 $\boldsymbol{s = s_0^{\, j}}$ 和 $\boldsymbol{\tau = k \cdot s_0^{\, j} \cdot \tau_0}$，母小波可以表示为：

\[\psi_{\tau, s} = \frac{1}{\sqrt{s}} \psi \left( \frac{t - \tau}{s} \right) \\ \downarrow \\ \psi_{j, k}(t) = s_0^{\frac{-j}{2}} \psi \left( s_0^{-j} - k \tau_0 \right)\]

如果 $\boldsymbol{ \left\{ \psi_{(j, \, k)} \right\} }$ 构成正交基，则小波级数变换为：

\[\Psi^{\psi_{\, j,k}}_x = \int x(t) \, \psi^*_{j, \, k}(t) \, dt \\ x(t) = c_\psi \sum\limits_{j} \sum\limits_{k} \Psi^{\psi_{\, j,k}}_x \, \psi_{\, j,k} (t)\]

小波级数要求 $\boldsymbol{ {\psi_{(j, \, k)}} }$ 为正交、双正交或框架级数。如果 $\boldsymbol{ {\psi_{(j,k)}} }$ 不为正交，则小波级数为：

\[\Psi^{\psi_{\, j,k}}_x = \int x(t) \, \hat{\psi}^*_{(j, \, k)}(t) \, dt\]

其中 $\boldsymbol{\hat{\psi}_{j, \, k}^*(t)}$ 是对偶双正交基或对偶框架（注意 $*$ 表示共轭）。

（3）常用的连续小波

⚪ 消失矩 Vanishing Moment

消失矩用来检查母小波是否为高频的函数。消失矩越高，经过内积后被滤掉的低频成分越多。

对于小波母函数$\psi(\cdot)$，定义$k$-th动量：

\[m_k = \int t^k \psi(t) dt\]

若 $m_0=m_1=\cdots =m_{k-1} = 0$，则称$\psi(\cdot)$有$k$个消失矩。

⚪ 哈尔小波 Haar wavelet

哈尔小波的母函数定义为:

\[\psi(t) = \begin{cases} 1, & 0\leq t < \frac{1}{2} \\ -1, & \frac{1}{2} \leq t < 1 \\ 0, & \text{otherwise} \end{cases}\]

哈尔小波的消失矩为$1$，因此它并不适合用于较为平滑的函数。

⚪ 墨西哥帽小波 Mexican Hat wavelet

墨西哥帽小波的母函数定义为高斯函数的二阶导数：

\[\psi(t)= \frac{1}{\sqrt{2 \pi} \cdot \sigma^3} \left( e^{\frac{-t^2}{2 \sigma^2}} \cdot \left( \frac{t^2}{\sigma^2} - 1 \right) \right)\]

墨西哥帽小波的消失矩为$2$。

⚪ Morlet小波 Morlet wavelet

Morlet小波的母函数定义为:

\[\psi(t) = e^{iat} \cdot e^{-\frac{t^2}{2\sigma}}\]

其中$a$是调制参数，$\sigma$是影响窗口宽度的缩放参数。

4. 离散小波变换 Discrete Wavelet Transform

离散小波变换是指利用数字滤波技术获得数字信号（离散信号）的时间尺度表示，并且只选择部分缩放因子和平移参数来进行计算。DWT 系数通常是从二进制网格上的 CWT 中采样而来的，即$s = 2^j,s_0 = 2$和$t = k \cdot 2^j,t_0 = 1$。

DWT的思想是使用截止频率不同的滤波器来分析不同尺度的信号。信号通过一系列高通滤波器分析高频部分，再通过一系列低通滤波器分析低频部分。

（1）正变换（信号的分解过程）

DWT 通过将信号分解为粗略近似值和细节信息，以不同的分辨率分析不同频带的信号。DWT 使用两组函数，分别称为缩放函数 (scaling function)$h[n]$（低通滤波器）和小波函数 (wavelet function) $g[n]$（高通滤波器）。

原始信号 $x[n]$ 首先经过半带高通滤波器 $g[n]$ 和低通滤波器 $h[n]$。滤波后，通过丢弃每个间隔的样本对信号进行下采样，使得信号的时间分辨率减半、频率分辨率加倍。这构成了一级分解：

\[y_{high}[k] = \sum\limits_{n} x[n] \cdot g[2k - n] \\ y_{low}[k] = \sum\limits_{n} x[n] \cdot h[2k - n]\]

上述过程（也称为子带编码subband coding）可以重复进行进一步的分解。在每一级，滤波和下采样都会使样本数量减半（时间分辨率减半）、跨越的频带减半（频率分辨率加倍）。

（2）常用的离散小波

高通滤波器和低通滤波器的脉冲响应之间应满足关系：

\[g[L - 1 - n] = (-1)^n \cdot h[n]\]

其中 $g[n]$ 为高通滤波器，$h[n]$ 为低通滤波器，$L$ 为滤波器长度（以点数表示）。满足此条件的滤波器被称为正交镜像滤波器 (Quadrature Mirror Filter, QMF)。

注意，QMF的两个滤波器有如下关系：$g$是$h$的反序，且基数索引位置要乘以$-1$。

⚪ 哈尔小波 Haar wavelet

在所有正交小波变换中，哈尔小波变换是最简单的一种变换。

哈尔小波的滤波器定义为：

\[h[n] = \begin{cases} \frac{1}{\sqrt{2}}, & n = 0, 1 \\ 0, & \text{otherwise} \end{cases} \\ g[n] = \begin{cases} \frac{1}{\sqrt{2}}, & n = 0 \\ -\frac{1}{\sqrt{2}}, & n = 1 \\ 0, & \text{otherwise} \end{cases} \\ \downarrow \\ h = [\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}}],g = [\frac{1}{\sqrt{2}},-\frac{1}{\sqrt{2}}]\]

哈尔小波是多贝西小波于$N=2$的特例，可称之为D2或db1。

⚪ 多贝西小波 Daubechies wavelet

多贝西小波是最常使用到的离散小波转换。通常以滤波器的长度$N$来形容滤波器为DN，例如$N=2$的多贝西小波写作D2、$N=4$的多贝西小波写作D4，以此类推（$N$为偶数）。

长度为$N$的多贝西小波的消失矩为$N/2$。与之对应的另一种描述方式是db N，其中$N$指消失矩的数量，由此，D4和db2是同一种多贝西小波。

关于多贝西小波等小波的设计可参考小波分析。

（3）反变换（信号的重构过程）

由于QMF构成正交基，因此DWT的重构过程非常简单：按相反顺序进行，每一层的信号都经过 $2$ 倍上采样，然后经过合成滤波器 $g[n]$ 和 $h[n]$（分别为高通和低通），最后相加。

\[x[n] = \sum\limits_{k = -\infty}^\infty \left( \, y_{high}[k] \cdot g[-n + 2k] \, \right) + \left( \, y_{low}[k] \cdot h[-n + 2k] \, \right)\]

5. 深度学习中小波变换的应用

（1）PyWavelets

paper, github, docs

PyWavelets库是一个基于 Python 的小波变换开源库，最初是一篇硕士论文的学术项目（使用小波变换对医学信号进行分析和分类）。

pip install PyWavelets

使用PyWavelets和pytorch实现一维序列的小波变换：

import pywt
import torch
import torch.nn.functional as F

def create_wavelet_filter(wave, in_size, out_size, type=torch.float):
    w = pywt.Wavelet(wave) # 创建wave类型的小波

    dec_hi = torch.tensor(w.dec_hi[::-1], dtype=type) # 高通分解滤波器，系数反转以匹配卷积操作
    dec_lo = torch.tensor(w.dec_lo[::-1], dtype=type) # 低通分解滤波器，系数反转以匹配卷积操作
    dec_filters = torch.stack([dec_lo, dec_hi], dim=0)

    dec_filters = dec_filters[:, None].repeat(in_size, 1, 1)

    rec_hi = torch.tensor(w.rec_hi[::-1], dtype=type).flip(dims=[0]) # 高通重构滤波器，系数反转以匹配卷积操作
    rec_lo = torch.tensor(w.rec_lo[::-1], dtype=type).flip(dims=[0]) # 高通重构滤波器，系数反转以匹配卷积操作
    rec_filters = torch.stack([rec_lo, rec_hi], dim=0)

    rec_filters = rec_filters[:, None].repeat(out_size, 1, 1)
    return dec_filters, rec_filters

def wavelet_transform(x, filters):
    b, c, l = x.shape
    pad = filters.shape[2] // 2 - 1
    x = F.conv1d(x, filters, stride=2, groups=c, padding=pad)
    x = x.reshape(b, c, 2, l // 2)
    return x

def inverse_wavelet_transform(x, filters):
    b, c, _, l_half = x.shape
    pad = filters.shape[2] // 2 - 1
    x = x.reshape(b, c * 2, l_half)
    x = F.conv_transpose1d(x, filters, stride=2, groups=c, padding=pad)
    return x

wt_type  = 'db1'
"""
>>> import pywt
>>> pywt.families()
['haar', 'db', 'sym', 'coif', 'bior', 'rbio', 'dmey', 'gaus', 'mexh', 'morl', 'cgau', 'shan', 'fbsp', 'cmor']
>>> pywt.families(short=False)
['Haar', 'Daubechies', 'Symlets', 'Coiflets', 'Biorthogonal', 'Reverse biorthogonal', 'Discrete Meyer (FIR Approximation)', 'Gaussian', 'Mexican hat wavelet', 'Morlet wavelet', 'Complex Gaussian wavelets', 'Shannon wavelets', 'Frequency B-Spline wavelets', 'Complex Morlet wavelets']
"""
wt_level = 1

wt_filter, iwt_filter = create_wavelet_filter(wt_type, in_channels, in_channels, torch.float)

x_ll_in_levels = []
x_h_in_levels = []
shapes_in_levels = []
curr_x_ll = x

for i in range(wt_level):
    curr_shape = curr_x_ll.shape
    shapes_in_levels.append(curr_shape)
    if curr_shape[2] % 2 > 0:
        curr_pads = (0, curr_shape[2] % 2)
        curr_x_ll = F.pad(curr_x_ll, curr_pads)

    curr_x = wavelet_transform(curr_x_ll, wt_filter)
    curr_x_ll = curr_x[:,:,0,:]
    curr_x_h  = curr_x[:,:,1,:]

    x_ll_in_levels.append(curr_x_ll)
    x_h_in_levels.append(curr_x_h)

next_x_ll = 0

for i in range(wt_levels-1, -1, -1):
    curr_x_ll = x_ll_in_levels.pop()
    curr_x_h = x_h_in_levels.pop()
    curr_shape = shapes_in_levels.pop()

    curr_x_ll = curr_x_ll + next_x_ll

    curr_x = torch.cat([curr_x_ll.unsqueeze(2), curr_x_h.unsqueeze(2)], dim=2)
    next_x_ll = inverse_wavelet_transform(curr_x, iwt_filter)

    next_x_ll = next_x_ll[:, :, :curr_shape[2]]

x_tag = next_x_ll
x = x + x_tag

（2）Pytorch Wavelets

paper, github, docs

Pytorch Wavelets库支持使用 pytorch 计算离散小波和双树复小波变换。相比于PyWavelets库，Pytorch Wavelets库的封装更好。

pip install pytorch_wavelets

使用Pytorch Wavelets实现一维序列的小波变换：

import torch
from pytorch_wavelets import DWT1DForward, DWT1DInverse  # or simply DWT1D, IDWT1D
dwt = DWT1DForward(wave='db6', J=3)
X = torch.randn(10, 5, 100) # (batch, channels, sizes)
yl, yh = dwt(X)
print(yl.shape)
>>> torch.Size([10, 5, 22])
print(yh[0].shape)
>>> torch.Size([10, 5, 55])
print(yh[1].shape)
>>> torch.Size([10, 5, 33])
print(yh[2].shape)
>>> torch.Size([10, 5, 22])
idwt = DWT1DInverse(wave='db6')
x = idwt((yl, yh))

（3）小波卷积 FWavelet Convolution (WTConv)

paper：Wavelet Convolutions for Large Receptive Fields

小波卷积先对输入进行小波变换，将输入分解为不同频率的分量，然后在这些分量上分别进行小卷积核的深度可分离卷积，最后通过逆小波变换将结果重构。该操作不仅能够将卷积分离到不同的频率分量上，还能够让小卷积核在更大的输入区域内操作，从而扩大其感受野。

⚪ Reference

The Wavelet Tutorial