我们当年分析信号时,遇到突变点常常很头疼。傅里叶变换能告诉你这个信号里有哪些频率成分,但它像个马后炮,永远说不清楚某个尖峰到底出现在哪个时刻。短时傅里叶变换试着开一扇时间窗,可窗口一固定,时间和频率的分辨率就锁死了。于是我们盯上了更聪明的框架:正交多分辨分析。
这套理论是构造小波基和实现Mallat快速算法的地基。它不像傅里叶那样把信号“一刀切”成全频段正弦波,而是像一台逐级放大的显微镜,把信号一层层拆成“骨架”和“细节”。如果你正准备入门小波分析,或者已经被各种函数空间绕得云里雾里,这篇就是第一块垫脚石。我会尽量抛开晦涩教材,用工程人的语言把“空间嵌套”“正交补”“双尺度方程”这些硬骨头讲透。
这一篇是系列的第一讲,重点放在数学原理和空间划分逻辑,先不堆代码。等这套“内功”打通了,后面看Mallat算法的C++或Python实现会轻松很多。
1. 先看直觉:从“分辨率”到“嵌套空间”的进化
1.1 为什么固定分辨率不够用
你处理过振动信号就知道,机械故障冲击往往是一瞬间的事。傅里叶变换把整段信号投影到无限延伸的正弦基上,高频分量在整个时间轴上均匀摊开,冲击发生的具体位置完全丢失。短时傅里叶变换用固定窗口切信号,窗口小了频率分辨率差,窗口大了时间定位又糊,怎么看都不对劲。
多分辨分析的核心想法很直接:既然一枚镜头做不到全焦段清晰,那就造一组“变焦镜头”。低频成分用大窗口看在宏观尺度上,高频成分用小窗口看局部细节,然后把这些结果拼起来。这套“变焦”机制在数学上的表达就是函数空间的逐级嵌套。
1.2 尺度空间:用放大倍数定义“分辨率”
我们定义一组彼此嵌套的闭子空间 (V_j),其中 (j) 是分辨率索引。(j) 越大,空间里的函数能表达的细节越丰富。数学上写成:
[ V_j = \text{span}\left{ \phi_{j,k}(t) = 2^{j/2} \phi(2^j t - k), \quad k \in \mathbb{Z} \right} ]
这里的 (\phi(t)) 叫尺度函数,是整个空间的“母版”。 (\phi_{j,k}(t)) 是母版在尺度 (j) 下平移 (k) 得到的子函数。关键在于逼近关系:
[ V_0 \subset V_1 \subset V_2 \subset \dots \subset V_j ]
想象一张 (256 \times 256) 的图像。如果只看像素的行列采样,它是 (V_8) 空间里的一个点;把它压缩成 (128 \times 128),它大约落在 (V_7) 里。压缩后图像损失的那部分高频纹理,其实就是 (V_8) 比 (V_7) 多出来的“增量信息”。随着 (j) 增大,空间不断扩张,可以分辨的函数细节也越来越细。
实操心得:初学者最难理解的就是“函数属于空间”这个抽象说法。你就把 (V_j) 看成“所有能用 (\phi(2^j t - k)) 线性组合出来的函数集合”。组合系数可以变,但基本零件是固定的。尺度越大,零件越精细。
这个逐级嵌套结构,正是“多分辨”这个名字的数学来源。接下来的问题自然变成了:当空间一大一小嵌套时,多出来的那块“增量”如何精确分离出来?这就引出了正交性。
2. 正交的威力:把“增量信号”干净地抽出来
2.1 直和分解:大空间 = 小空间 + 细节空间
前面说 (V_j) 包含 (V_{j-1}),但 (V_j) 比 (V_{j-1}) 多出来的部分不是随随便便的杂项。我们希望多出来的分量和小空间里的分量互不干扰,这样拆分后的信息才能无损合并回去。
正式定义:(W_{j-1}) 是 (V_{j-1}) 在 (V_j) 中的正交补空间,即
[ V_j = V_{j-1} \oplus W_{j-1} ]
符号 (\oplus) 表示直和,意味着 (V_j) 里任意一个函数 (f) 都能唯一分解成 (f = a + b),其中 (a \in V_{j-1}),(b \in W_{j-1}),并且 (a) 与 (b) 正交。这里“正交”指内积为零:
[ \langle a, b \rangle = \int a(t) b(t) , dt = 0 ]
这套正交特性在实际工程里异常宝贵。做信号去噪时,我们总希望把噪声阈值过滤掉而保留有效信号。如果细节分量 (W_{j-1}) 和近似分量 (V_{j-1}) 不正交,滤除噪声时就会连带损伤主体信号,怎么调阈值都别扭。正交保证了“井水不犯河水”,你动细节系数不会破坏近似分量的能量分布。
2.2 小波函数:细节空间的“骨架”
细节空间 (W_j) 同样需要一组基。这组基就由小波函数 (\psi(t)) 生成:
[ W_j = \text{span}\left{ \psi_{j,k}(t) = 2^{j/2} \psi(2^j t - k), \quad k \in \mathbb{Z} \right} ]
到这里,你可以把尺度函数 (\phi(t)) 理解成“宏观骨架”的生成器,小波函数 (\psi(t)) 理解成“微观毛刺”的生成器。两个空间正交,意味着它们代表的频率区间基本不重叠。尺度函数是低通的,小波函数是高通的,配合起来刚好把频谱“切”成两段。
避坑提示:注意这里说的是“基本不重叠”而非“严格不重叠”。现实中理想低通和高通滤波器都有无限长拖尾,工程上必须用有限长滤波器做逼近,所以常出现频谱泄露。用Daubechies小波时,正交性在数学上完美成立,但在离散序列的有限边界上会打折扣,这就引出了后面边界处理的种种讲究。
3. 双尺度方程:连接两代空间的“基因链条”
3.1 从 (\phi(t)) 到 (\phi(2t)):多分辨自相似性
前文要求所有函数空间由同一个母函数平移缩放生成,这意味着 (\phi(t)) 本身也必须能由更细尺度上的 (\phi(2t)) 拼出来。这就是双尺度方程(也叫两尺度方程或细分方程):
[ \phi(t) = \sqrt{2} \sum_{k} h[k] \phi(2t - k) ]
你一定好奇那个 (\sqrt{2}) 从哪来。原因是 (\phi_{j+1,k}(t) = 2^{(j+1)/2} \phi(2^{j+1}t - k)) 的范数是1,而 (\phi(t)) 的范数也是1。右边若直接累加 (\phi(2t)),宽度压缩一半,能量会膨胀,必须乘上 (\sqrt{2}) 保持总能量不变。这个归一化系数少写一个,后面所有系数都会等比例偏离,重建时信号幅度直接错误。
对应的,小波函数 (\psi(t)) 也满足类似的方程:
[ \psi(t) = \sqrt{2} \sum_{k} g[k] \phi(2t - k) ]
系数 (h[k]) 和 (g[k]) 本质上是离散滤波器。 (h) 是低通滤波器系数, (g) 是高通滤波器系数,二者构成正交镜像滤波器组。经典哈尔小波下:
[ h[0] = \frac{1}{\sqrt{2}},\quad h[1] = \frac{1}{\sqrt{2}},\quad g[0] = \frac{1}{\sqrt{2}},\quad g[1] = -\frac{1}{\sqrt{2}} ]
代入可得 (\phi(t) = \phi(2t) + \phi(2t-1)),这恰好描述了一个门函数的自相似性。
3.2 滤波器的“音乐节拍”:长度与正交约束
若滤波器长度是2,很容易写出上述简单系数。若换成长度4的Daubechies-2小波,系数不再是简单对称的,而要满足精确重建和正交条件。核心约束有两条:一是归一化 (\sum h[k] = \sqrt{2}),二是正交性 (\sum h[k]h[k+2m] = \delta[m])。第二条保证低通滤波器消除镜像频率分量,高通滤波器则通过交替取反构造出来。
经验之谈:读任何小波论文,先看滤波器符号约定。许多文献用 (h_0, h_1, h_2, \dots) 表示低通, (h_1) 同时表示高通中某个索引,变量撞车屡见不鲜。自己写代码时建议直接命名 (lp_coeff) 和 (hp_coeff),免得推导时被下标绕晕。
4. 动手推演:Mallat算法如何把理论变成整数运算
4.1 分解与重构:滤波器组与抽样
有了双尺度方程作为桥梁,就能把信号系数在不同尺度间搬移。设信号 (f \in V_J),它在 (V_J) 里的系数为 (c_J(k))。多分辨分析保证 (f) 还能写成:
[ f = \sum_k c_{J-1}(k) \phi_{J-1,k}(t) + \sum_k d_{J-1}(k) \psi_{J-1,k}(t) ]
通过推导能得到从 (c_J) 求 (c_{J-1}) 和 (d_{J-1}) 的递推关系,这就是Mallat快速算法的原型。离散形式长这样:
[ c_{J-1}(k) = \sum_m h[m] \cdot c_J(2k + m) ]
[ d_{J-1}(k) = \sum_m g[m] \cdot c_J(2k + m) ]
操作上就是:拿信号系数分别卷积低通和高通滤波器,然后隔一个抽一个(下采样因子2)。反过来,重构时先隔一个补一个零(上采样),再分别与滤波器卷积并求和:
[ c_{J}(k) = \sum_m h[k - 2m] c_{J-1}(m) + \sum_m g[k - 2m] d_{J-1}(m) ]
这个“分解-重构”结构因形似漏斗,常被称为“两通道滤波器组”。
4.2 4点信号手算实例:让公式落地
我们来个最简单的数字示例,用哈尔滤波器跑一遍,你看完就理解了。
输入信号:
[ c_2 = [1, 2, 3, 4] ]
使用滤波器 (h = [\frac{1}{\sqrt2}, \frac{1}{\sqrt2}]) 和 (g = [\frac{1}{\sqrt2}, -\frac{1}{\sqrt2}])。
分解第一步,先与 (h) 卷积后再下采样:
[ c_1(0) = h[0]\cdot c_2(0) + h[1]\cdot c_2(1) = \frac{1}{\sqrt2} \cdot 1 + \frac{1}{\sqrt2} \cdot 2 = \frac{3}{\sqrt2} ]
[ c_1(1) = h[0]\cdot c_2(2) + h[1]\cdot c_2(3) = \frac{1}{\sqrt2} \cdot 3 + \frac{1}{\sqrt2} \cdot 4 = \frac{7}{\sqrt2} ]
高通通道:
[ d_1(0) = g[0]\cdot c_2(0) + g[1]\cdot c_2(1) = \frac{1}{\sqrt2} \cdot 1 - \frac{1}{\sqrt2} \cdot 2 = -\frac{1}{\sqrt2} ]
[ d_1(1) = g[0]\cdot c_2(2) + g[1]\cdot c_2(3) = \frac{1}{\sqrt2} \cdot 3 - \frac{1}{\sqrt2} \cdot 4 = -\frac{1}{\sqrt2} ]
于是第一层结果是:
[ c_1 = \left[\frac{3}{\sqrt2}, \frac{7}{\sqrt2}\right], \quad d_1 = \left[-\frac{1}{\sqrt2}, -\frac{1}{\sqrt2}\right] ]
重构时,先对 (c_1) 和 (d_1) 隔点插零,再与各自的滤波器卷积叠加:
[ c_2(0) = h[0] c_1(0) + g[0] d_1(0) = \frac{1}{\sqrt2} \cdot \frac{3}{\sqrt2} + \frac{1}{\sqrt2} \cdot (-\frac{1}{\sqrt2}) = \frac{3}{2} - \frac{1}{2} = 1 ]
[ c_2(1) = h[1] c_1(0) + g[1] d_1(0) = \frac{1}{\sqrt2} \cdot \frac{3}{\sqrt2} + (-\frac{1}{\sqrt2}) \cdot (-\frac{1}{\sqrt2}) = \frac{3}{2} + \frac{1}{2} = 2 ]
可以看到原信号被“拆开再合成”,数值完美还原。这套手算过程虽然简单,但能很直观地说明“为什么滤波之后能量总和不变”。如果滤波器系数保持正交归一,分解后的近似系数和细节系数满足:
[ \sum |c_2(k)|^2 = \sum |c_1(k)|^2 + \sum |d_1(k)|^2 ]
这就是帕塞瓦尔定理在小波域的表现。实际工程中,经常用这个等式检查自己的滤波器实现是否写对了——能量不守恒,必然有系数错误。
4.3 多级分解的树状结构
实际处理信号时,我们通常会继续对 (c_1) 再分解,得到 (c_0) 和 (d_0),形成经典的小波包树状结构。每一层低频近似继续下拆,高频细节作为独立分支保存。图像压缩里真正有价值的信息大多集中在那几层低频近似里,高频细节往往大量趋近于零,这也是阈值压缩的理论依据。
5. 实战避坑:索引、归一化与边界的三个陷阱
5.1 索引约定:工程里最常见的翻车点
不同的教科书和代码库对滤波器下标、信号起始位置有截然不同的约定。有的按因果系统从0开始,有的按对称序列从负数开始。拿Matlab做小波分解,返回的近似系数和细节系数长度与Python PyWavelets不一定一致,直接套用转换公式很容易错位。
避坑方案:动手实现前,先拿单位脉冲信号(一个位置为1、其余为0的序列)过一遍完整流程。如果重构后能还原出原始脉冲位置,说明索引约定正确;如果脉冲位置偏移或出现两端毛刺,立刻检查下标公式。
5.2 归一化细节:滤波器、系数与信号幅度的三角关系
前面强调过 (\sqrt2) 的来历,但在具体实现里经常混着各种归一化版本。有的库把 (\sqrt2) 直接并进滤波器系数,比如哈尔的 h 写成 ([1,1]) 然后单独除 (\sqrt2);有的库直接在变换里做总能量归一化。你如果混着用,分解系数看起来对不上,重建出来幅值就错了。
避坑方案:每次拿到新的库,先跑一遍能量守恒检测。输入长度为8的随机序列,做一次分解重构,比较输入输出绝对误差。如果误差在 (10^{-10}) 量级,说明实现稳定;如果误差接近 (0.1),八成是归一化约定不一致。
5.3 边界处理:最容易被忽视的“伪影源头”
信号是有限长的,滤波器的卷积却天然需要边界之外的样本。最常见的做法是零填充、对称延拓、周期延拓和常数延拓。对于正交小波变换,周期延拓能保持系数数量精确为原长一半,但信号两端不连续时会产生边界跳跃伪影;对称延拓在多数图像处理任务里表现更自然,但要小心偶数长度滤波器导致的对称轴偏移。
实测心得:我当年处理地震数据时,直接用零填充做三级分解,结果重构信号开头总有一段无意义的小幅振荡。换成对称延拓后,边界伪影立刻消失。图像压缩里边界效应还会影响码率统计,所以有些库会默认对图像行做“周期化”,这对普通平缓内容没问题,对有剧烈明暗变化的边界就要慎用。
补充一个容易踩的小细节:正交小波分解后各层系数长度之和等于原信号长度加滤波器长度减一。为保证不膨胀,常规算法要求原始信号长度是 (2^j) 的倍数。如果信号长度不满足,需要预先做对称延拓或截断。许多入门者在这里栽跟头,因为底层库会自动包边界,却不告诉你已经悄悄变了数据长度。
5.4 正交变换的隐形福利:噪声能量估算
由于正交性,小波系数域的噪声能量等于时域噪声能量。这个性质在做阈值去噪时极为有用。你在时域无法直接分离噪声,但小波变换后,白噪声在不同尺度上大约均匀展开,而有效信号的系数集中在少数大值上。根据噪声标准差估计阈值,就能放心地把小系数清零,而不用担心把信号尖峰抹掉。这也是正交多分辨分析对比其他非正交基的一大优势。
如果你决定用双正交小波做图像压缩——比如JPEG2000里的CDF 9/7——就要清楚它牺牲了严格正交、换来了线性相位。这时能量守恒不再严格成立,阈值估计需要额外小心。基础概念学完之后,从正交小波转到双正交小波,你会对“为什么9/7滤波器的重构系数比分解系数光滑”有更深体会。
这套理论我自己反复读了四五遍才真正理顺,最初卡在 (V_j) 与 (W_j) 的嵌套关系上。后来用图像金字塔做类比,再把哈尔小波完整手算一遍,瞬间通了。如果你也感觉抽象,建议拿起笔把第4节那个4点示例一步步推导一次,用不了十分钟,但对理解公式的意义帮助很大。
下一篇我会从这套空间框架跳到具体的Daubechies小波构造,讲讲“消失矩”如何决定滤波器的长相,以及它是如何影响频响曲线陡峭度的。到时候我们就能亲手设计一套满足需求的滤波器组了。