小波变换 | Haar 小波变换

小波变换 Haar 小波变换

连续形式

在实际实现中，尺度函数和小波函数中的系数，会被当做是滤波器参数。最终可以基于这里尺度函数和小波函数构造的基函数集合，来对输入信号进行分解（将信号投影到这些基函数上）。

• 尺度母函数 $\varphi (t)=\{\begin{array}{ll}1& 0\le t<1\\ 0& 其他\end{array}$
• 小波母函数 $\psi (t)=\{\begin{array}{ll}1& 0\le t<1/2\\ -1& 1/2\le t<1\\ 0& 其他\end{array}$

计算系数 $h[n]$ 和 $g[n]$

基于 ${\sum }_{n\in \mathbb{Z}}h[n]h[n-2m]={\delta }_m$

Haar 中，将卷积长度限定为 $L=2$，即在 n∉{0,1}n \notin \{0, 1\}n∈/{0,1} 时 $h[n]=0$。而 $h[n]=a$，$h[1]=b$，先使用待定系数。基于已知的 ${\sum }_{n\in \mathbb{Z}}h[n]h[n-2m]={\delta }_m$，我们可以对 $m$ 进行讨论，从而获得关于 Haar 中 $h$ 的更多限定：

∑n∈{0,1}h[n]h[n]=a2+b2=1m=0∑n∈{0,1}h[n]h[n−2]=0+0=0m=1 \begin{aligned} \sum_{n \in \{0, 1\}} h[n]h[n] & = a^2 + b^2 = 1 & m=0 \\ \sum_{n \in \{0, 1\}} h[n]h[n-2] & = 0 + 0 = 0 & m=1 \end{aligned} n∈{0,1}∑​h[n]h[n]n∈{0,1}∑​h[n]h[n−2]​=a2+b2=1=0+0=0​m=0m=1​

为获得最简单（对称）结构，取 $a=b$，则可以得到 $h[0]=h[1]=\frac{1}{\sqrt{2}}$。进一步，通过 Haar 中采用的构造方式，即 $g[n]=(-1{)}^{n}h[1-n]$，可以进一步推出 $g[0]=-g[1]=\frac{1}{\sqrt{2}}$。所以可以得到更具体的细化方程：

• 尺度细化方程 $\varphi (t)=2^{1/2}[h[0]\varphi (2t)+h[1]\varphi (2t-1)]=\varphi (2t)+\varphi (2t-1)$，对应的 $h[0]=2^{-1/2},h[1]=2^{-1/2}$
• 小波细化方程：$\psi (t)=2^{1/2}[g[0]\psi (2t)+g[1]\psi (2t-1)]=\psi (2t)-\psi (2t-1)$，对应的 $g[0]=2^{-1/2},g[1]=-2^{-1/2}$

基于正交展开

在计算系数的过程中，用到了正交展开，其中系数就是对应着目标函数在各个正交基上的投影。所以如下的细化公式中的系数，可以看做是函数 $\varphi (t)$ 在正交归一的基函数 $2^{1/2}\varphi (2t-n)$ 上的投影模长：

$$\varphi (t)=2^{1/2}\sum _{n\in \mathbb{Z}}h[n]\varphi (2t-n)$$

所以可以直接计算投影系数：

$$\begin{array}{rl}h[n]& =\frac{⟨\varphi (t),2^{1/2}\varphi (2t-n)⟩}{\Vert 2^{1/2}\varphi (2t-n){\Vert }^2}\\ & =⟨\varphi (t),2^{1/2}\varphi (2t-n)⟩(\Vert 2^{1/2}\varphi (2t-n){\Vert }^2=1)\\ & =2^{1/2}{\int }_0^1\varphi (2t-n)dt(\varphi (t)=1)\end{array}$$

考虑其中的放缩偏移版本的尺度函数，满足如下范围：

$$\varphi (2t-n)=\{\begin{array}{ll}1& t\in [n/2,(n+1)/2]\\ 0& 其他\end{array}$$

考虑到 $\varphi (t)$ 本身对自变量的范围限制，我们可以知道 $n$ 仅可以取值为 0 或者 1：

$$\begin{array}{rl}h[n]& =2^{1/2}{\int }_0^1\varphi (2t-n)dt\\ & =\{\begin{array}{ll}{2}^{1/2}{\int }_{1/2}^1\varphi (2t-1)dt=2^{1/2}{\int }_{1/2}^{1}1dt=\frac{1}{\sqrt{2}}& n=1\\ 2^{1/2}{\int }_0^{1/2}\varphi (2t)dt=2^{1/2}{\int }_0^{1/2}1dt=\frac{1}{\sqrt{2}}& n=0\end{array}\end{array}$$

从而可得最终的细化方程。同理 $g[n]$ 也可以计算得到。

连续时间到离散索引

应用于实际的离散采样序列（如音频或图像）时，离散小波变换的目标是将信号分解为一组多尺度、多位置的正交基函数上的展开系数。这个过程中，我们通常将原始的信号视为最细层级 $J_{\max }=\lfloor {\log }_{2}N\rfloor$（假定该尺度信号长度为 $N$）。每完成一级小波分解，就将尺度系数（低频部分）下采样一半，进入更粗一层的表示（数据分辨率降低，时间、空间尺度增大）。能够进行的最大分解层数由信号长度决定。最终，原始信号被完整表示为来自最粗层的尺度系数 $a_J$ 和所有层的小波系数 {dj}j=JJmax⁡−1\{d_{j}\}_{j=J}^{J_{\max}-1}{dj​}j=JJmax​−1​。他们分别对应信号在不同尺度的近似于细节成分。

在处理离散数据序列的时候，实际上涉及到了尺度函数和小波函数的离散形式，这里 $1/N$ 实际上是因为考虑到了正交基在离散过程中正交归一性的保持才修改的：

$${\varphi }_{j,k}(t),{\psi }_{j,k}(t)\Rightarrow {\varphi }_{j,k}^{\prime }[n]=\frac{1}{\sqrt{N}}{\varphi }_{j,k}[n],{\psi }_{j,k}^{\prime }[n]=\frac{1}{\sqrt{N}}{\psi }_{j,k}[n]$$

实现中则依照 $J_{\max }\to J$ 的方向不断拆分，其中 $a$ 表达了每个尺度的近似（低频）系数，而 $d$ 捕捉了每个尺度的细节（高频）系数。而对于 $J_{\max }$，通常直接定义 $a_{J_{\max }}=x[n]$，也就是说，原始采样信号本身就是最细尺度下的尺度系数。我们无需再从别的什么空间“计算”它，它就是我们处理的起点。

处理离散时间信号 $x[n]$，取 $J=0$，则其在 $V_J$ 和 $W_J$ 上的正交基函数展开：

$$\begin{array}{rl}x[n]& =\sum _k{a}_J[k]{\varphi }_{J,k}^{\prime }[n]+\sum _{j=J}^{J_{\max }}\sum _k{d}_j[k]{\psi }_{j,k}^{\prime }[n]\\ {\varphi }_{j,k}[n]& =\sqrt{\frac{2^j}{N}}\varphi \left(\frac{2^j}{N}n-k\right)\\ {\psi }_{j,k}[n]& =\sqrt{\frac{2^j}{N}}\psi \left(\frac{2^j}{N}n-k\right)\end{array}$$

从 [[100-收集的信息碎片/小波变换 离散小波变换|小波变换 离散小波变换]] 中可以知道如下对应分解和重构两个过程的公式：

$$\begin{array}{rl}{a}_j[k]& =\sum _{m}h[m-2k]a_{j+1}[m]\\ d_j[k]& =\sum _{m}g[m-2k]a_{j+1}[m]\\ a_{j+1}[k]& =\sum _{m}h[k-2m]a_j[m]+\sum _{m}g[k-2m]d_j[m]\end{array}$$

此时带入 $h$ 和 $g$ 的具体取值，则存在：

$$\begin{array}{rl}{a}_j[k]& =h[0]a_{j+1}[2k]+h[1]a_{j+1}[2k+1]=2^{-1/2}(a_{j+1}[2k]+a_{j+1}[2k+1])\\ d_j[k]& =g[0]a_{j+1}[2k]+g[1]a_{j+1}[2k+1]=2^{-1/2}(a_{j+1}[2k]-a_{j+1}[2k+1])\\ a_{j+1}[k]& \Leftarrow \{\begin{array}{llc}{a}_{j+1}[2m]& =h[0]a_j[m]+g[0]d_j[m]=2^{-1/2}{a}_j[m]+2^{-1/2}{d}_j[m]& k=2m\\ a_{j+1}[2m+1]& =h[1]a_j[m]+g[1]d_j[m]=2^{-1/2}{a}_j[m]-2^{-1/2}{d}_j[m]& k=2m+1\end{array}\end{array}$$

对于一个简单的例子：$x=[x_0,x_1,x_2,x_3]$，下面给出了几种不同的处理视角。

直接公式推导版

1. 预设最大层级系数（$j=2$）：$a_2=x$
2. 第一次分解（$j=2\to 1$）：$a_1[k]=2^{-1/2}(a_2[2k]+a_2[2k+1])$
   1. $a_1[0]=2^{-1/2}(a_2[0]+a_2[1])=2^{-1/2}(x_0+x_1)$
   2. $a_1[1]=2^{-1/2}(a_2[2]+a_2[3])=2^{-1/2}(x_2+x_3)$
   3. $d_1[0]=2^{-1/2}(a_2[0]-a_2[1])=2^{-1/2}(x_0-x_1)$
   4. $d_1[1]=2^{-1/2}(a_2[2]-a_2[3])=2^{-1/2}(x_2-x_3)$
3. 第二次分解（$j=1\to 0$）：$a_0[k]=2^{-1/2}(a_1[2k]+a_1[2k+1])$
   1. $a_0[0]=2^{-1/2}(a_1[0]+a_1[1])$
   2. $d_0[0]=2^{-1/2}(a_1[0]-a_1[1])$
4. 完整形式表示：$x(t)=a_0[0]{\varphi }_{0,0}(t)+d_0[0]{\psi }_{0,0}(t)+d_1[0]{\psi }_{1,0}(t)+d_1[1]{\psi }_{1,1}(t)$

矩阵系数版

对于这个计算过程，其实可以直接写成矩阵形式。Haar 分解矩阵是正交的，所以直接使用转置即可重构数据。

$$\begin{array}{rl}\left[\begin{array}{c}{a}_0[0]\\ d_0[0]\\ d_1[0]\\ d_1[1]\end{array}\right]& =\left[\begin{array}{cccc}h[0]h[0]& h[1]h[0]& h[0]h[1]& h[0]h[1]\\ h[0]g[0]& h[1]g[0]& h[0]g[1]& h[0]g[1]\\ g[0]& g[1]& 0& 0\\ 0& 0& g[0]& g[1]\end{array}\right]\left[\begin{array}{c}x[0]\\ x[1]\\ x[2]\\ x[3]\end{array}\right]\\ & =\left[\begin{array}{cccc}1/2& 1/2& 1/2& 1/2\\ 1/2& 1/2& -1/2& -1/2\\ 1/\sqrt{2}& -1/\sqrt{2}& 0& 0\\ 0& 0& 1/\sqrt{2}& -1/\sqrt{2}\end{array}\right]\left[\begin{array}{c}x[0]\\ x[1]\\ x[2]\\ x[3]\end{array}\right]\end{array}$$

下采样卷积版

这里基于公式：

$$\begin{array}{rl}{a}_j[k]& =\sum _{m}h[m-2k]a_{j+1}[m]\\ d_j[k]& =\sum _{m}g[m-2k]a_{j+1}[m]\\ a_{j+1}[k]& =\sum _{m}h[k-2m]a_j[m]+\sum _{m}g[k-2m]d_j[m](重构时，这里的a_j被不断更新)\end{array}$$

对于离散序列 a[k] 和 h[m]，其离散卷积定义为：$(a\ast h)[m]={\sum }_{k}a[k]h[m-k]$。而这里的形式其实可以看做是卷积核被翻转了。从 $2k-m$ 变为了 $m-2k$。所以在两个系数的逐层计算中，完全可以利用下采样和基于 $h$、$g$ 两种卷积核的卷积来计算。

具体实现

import mathimport numpy as npclass HaarWaveletTransform:def __init__(self):# Haar Coefficientsself.g_0 = self.h_1 = self.h_0 = 1 / math.sqrt(2)self.g_1 = -self.g_0print(f"[h_0, h_1, g_0, g_1]: {[self.h_0, self.h_1, self.g_0, self.g_1]}")# Haar Scale Functiondef phi(self, t) -> float:if 0 <= t < 1:return 1else:return 0def phi_jk(self, j, k, t) -> float:return 2 ** (j / 2) * self.phi(2**j * t - k)# Haar Wavelet Functiondef psi(self, t) -> float:if 0 <= t < 1 / 2:return 1elif 1 / 2 <= t < 1:return -1else:return 0def psi_jk(self, j, k, t) -> float:return 2 ** (j / 2) * self.psi(2**j * t - k)@propertydef coefficient_matrix(self):# Haar 分解矩阵是正交的，所以直接使用转置即可重构数据return np.array([[self.h_0 * self.h_0, self.h_1 * self.h_0, self.h_0 * self.h_1, self.h_1 * self.h_1],[self.h_0 * self.g_0, self.h_1 * self.g_0, self.h_0 * self.g_1, self.h_1 * self.g_1],[self.g_0, self.g_1, 0, 0],[0, 0, self.g_0, self.g_1],])def matrix_based_decompose(self, x):x = np.array(x).reshape(-1, 1)return self.coefficient_matrix.dot(x)def matrix_based_reconstruct(self, coef_matrix):return self.coefficient_matrix.T.dot(coef_matrix)def formula_based_decompose(self, x):"""分解（Decomposition）／分析（Analysis）阶段，又称 Forward DWT（正向离散小波变换）或 分析滤波器组（Analysis Filter Bank）。将原始信号通过低通/高通滤波+下采样，依次提取出每一层的近似系数（approximation）和细节系数（detail），构成系数序列。"""n_levels = math.ceil(math.log2(len(x)))n_elements = 2**n_levelsx.extend([0] * (n_elements - len(x)))a = x  # a_{J_{\max}} = xa_seq = []d_seq = []for j in range(n_levels, 0, -1):d_down = []a_down = []# print(f"Level {j}->{j - 1}")for k in range(0, len(a), 2):# print(f"{k} and {k + 1} (a[{j},{k}]={a[k]}, a[{j},{k + 1}]={a[k + 1]})")d_down.append(self.g_0 * a[k] + self.g_1 * a[k + 1])a_down.append(self.h_0 * a[k] + self.h_1 * a[k + 1])a = a_downd_seq.append(d_down)a_seq.append(a_down)assert len(a) == 1d_seq.reverse()  # from 0 to J_max-1a_seq.reverse()return a_seq, d_seqdef formula_based_reconstruct(self, a_seq, d_seq):"""重构（Reconstruction）／综合（Synthesis）阶段，又称 Inverse DWT（逆离散小波变换）或 综合滤波器组（Synthesis Filter Bank），将最粗尺度的近似系数和各层细节系数反向上采样、滤波并叠加，恢复出原始信号。"""n_elements = 2 ** len(d_seq)n_levels = len(d_seq)scale_factor = 1 / math.sqrt(n_elements)# reconstruct original signalreconstructed_x = []for t in range(n_elements):cur_t = t / n_elements# low frequencyx_t = a_seq[0][0] * self.phi_jk(0, 0, cur_t) * scale_factor# high frequencyfor j in range(n_levels):for k in range(2**j):# print(f"Collecting (j,k)={j, k}")x_t += d_seq[j][k] * self.psi_jk(j, k, cur_t) * scale_factorreconstructed_x.append(x_t)return reconstructed_xdef convolve(self, x, kernel, flip_kernel=True, verbose=False):"""x[m], 0 <= m < M, kernel[k], 0 <= k < Ky[n] = sum_{m=0}^{(K-1)+(M-1)} x[m] * kernel[n-m]"""M = len(x)K = len(kernel)O = M + K - 1  # 至少存在一个重叠if flip_kernel:kernel = kernel[::-1]  # => kernel[K-1-k]y = np.zeros(O)for n in range(O):for m in range(M):k = n - mif 0 <= k < K:y[n] += x[m] * kernel[k]if verbose:print(f"y[{n}]={y[n]}, x[{m}]={x[m]}, k[{k}]={kernel[k]}")return ydef downsample_1d(self, x):return x[1::2]def upsample_1d(self, x):y = np.zeros(len(x) * 2, dtype=x.dtype)y[::2] = xreturn ydef conv_based_decompose(self, x):"""a_{j}[k] =  \sum_{m} h[m-2k] a_{j+1}[m], d_{j}[k] = \sum_{m} g[m-2k] a_{j+1}[m]=> \sum_{m} h[(1+1k) - m] a_{j+1}[m] = \sum_{m} g[(1+2k) - m] a_{j+1}[m]=> 基于标准卷积，但是下采样需要取奇数项"""n_levels = math.ceil(math.log2(len(x)))n_elements = 2**n_levelsx = x + [0] * (n_elements - len(x))h_kernel = np.array([self.h_0, self.h_1])g_kernel = np.array([self.g_0, self.g_1])a_seq = []d_seq = []a = x  # a_{J_{\max}} = xfor j in range(n_levels, 0, -1):d = self.convolve(a, g_kernel, flip_kernel=True)a = self.convolve(a, h_kernel, flip_kernel=True)d = self.downsample_1d(d)a = self.downsample_1d(a)a_seq.append(a)d_seq.append(d)a_seq.reverse()d_seq.reverse()return a_seq, d_seqdef conv_based_reconstruct(self, a_0, d_seq):n_levels = len(d_seq)h_kernel = np.array([self.h_0, self.h_1])g_kernel = np.array([self.g_0, self.g_1])a = a_0for j in range(n_levels):a = self.upsample_1d(a)d = self.upsample_1d(d_seq[j])a = self.convolve(a, h_kernel, flip_kernel=False)d = self.convolve(d, g_kernel, flip_kernel=False)a = a[: 2 ** (j + 1)]d = d[: 2 ** (j + 1)]a = a + dreturn ax = [1, 4, -3, 0]
hwt = HaarWaveletTransform()print(">>>  Formula based decomposition and reconstruction")
a_seq, d_seq = hwt.formula_based_decompose(x)
print(a_seq, d_seq)
decoded_x = hwt.formula_based_reconstruct(a_seq, d_seq)
print(x, decoded_x)print(">>>  Convolution based decomposition and reconstruction")
a_seq, d_seq = hwt.conv_based_decompose(x)
print(a_seq, d_seq)
decoded_x = hwt.conv_based_reconstruct(a_seq[0], d_seq)
print(x, decoded_x)print(">>>  Matrix based decomposition and reconstruction")
coef = hwt.matrix_based_decompose(x)
print(coef.flatten().tolist())
decoded_x = hwt.matrix_based_reconstruct(coef)
print(x, decoded_x.flatten().tolist())