小波变换[4] -- 计算问题

2022-01-14 15:20 作者:nyasyamorina 0人读过 | 我要投稿

函数近似

在第二篇专栏里提到过多分辨率框架下信号的分解和重构, 可以看到这两个步骤都只与尺度系数 {pₖ} 有关, 而于尺度函数和小波函数的具体形式无关. 但是在特定场合, 比如说可视化或者误差评估的时候, 则需要用到尺度函数的形状了, 下面介绍两种逐渐收敛得到尺度函数的方法.

第一种 是由迭代式 $%5CPhi_%7Bn%2B1%7D(x)%3D%5Csum_%7Bk%5Cin%5Cmathbb%20Z%7Dp_k%5CPhi_n(2x-k)$ 给出. 不断应用这个迭代式, 可以有 $%5CPhi_n(x)%3D%5Csum_%7Bk%5Cin%5Cmathbb%20Z%7Dc%5E%7B(n)%7D_k%5CPhi_0(2%5Enx-k)$ , 由尺度函数的标准正交性可以知道 $c%5E%7B(0)%7D_k%3D%5Cdelta_%7B0%2Ck%7D$ , 由定义得 $c%5E%7B(1)%7D_k%3Dp_k$ . 结合两式: $%5CPhi_%7Bn%2B1%7D(x)%3D%5Csum_%7Bk%5Cin%5Cmathbb%20Z%7Dp_k%5Csum_%7Bl%5Cin%5Cmathbb%20Z%7Dc%5E%7B(n)%7D_l%5CPhi_0(2%5E%7Bn%2B1%7Dx-2%5Enk-l)$ , 使用 l-2ⁿk 代替 l, 交换 k 和 l 的累加顺序后, 得到 $%5Csum_%7Bl%5Cin%5Cmathbb%20Z%7D%5CPhi_0(2%5E%7Bn%2B1%7Dx-l)%5Csum_%7Bk%5Cin%5Cmathbb%20Z%7Dp_kc%5E%7B(n)%7D_%7Bl-2%5Enk%7D$ , 亦即 $c%5E%7B(n%2B1)%7D_l%3D%5Csum_%7Bk%5Cin%5Cmathbb%20Z%7Dp_kc%5E%7B(n)%7D_%7Bl-2%5Enk%7D$ . 如果选取 Φ₀ 为 Haar, 可以有 $%5CPhi_n(x)%3Dc%5E%7B(n)%7D_%7B%5Clfloor2%5Enx%5Crfloor%7D$ .

第二种 是由尺度关系式 $%5CPhi(x)%3D%5Csum_%7Bk%5Cin%5Cmathbb%20Z%7Dp_k%5CPhi(2x-k)$ 推导二进制点的值得出. 因为 N 阶 D-小波的支撑集在 (0,2N-1) 里, 所以仅在整数 x = 1, ..., 2N-2 上有可能不为零. 记整数点上 Φ 的值为向量形式 $%5Cvec%20v%3D%5Cleft%5B%5Cbegin%7Bmatrix%7D%5CPhi(1)%26%5CPhi(2)%26%5Ccdots%26%5CPhi(2N-2)%5Cend%7Bmatrix%7D%5Cright%5D%5ET$ , 那么尺度关系式可以看作矩阵作用在向量上, 稍加推导可以得出相应的矩阵为 $%5Chat%20P_%7Ba%2Cb%7D%3Dp_%7B2a-b%7D$ , 亦即尺度关系式变为 P̂v = v. 由此可以解得 P̂ 特征值为 1 的特征向量 [实际上可以证得 P̂ 特征值为 1 的特征向量存在且只有一个, 但是不太会线代就摸了], 并且因为 D-小波满足 $%5Cint_%5Cmathbb%20R%5CPhi(x)dx%3D1$ , 所以也满足 $%5Csum_%7Bk%3D1%7D%5E%7B2N-2%7D%5CPhi(k)%3D1$ [证明在下方 * ], 即求得 v 的值. 假设已经知道了二进制点 x=2⁻ⁿk; k∈Z 的值, 那么可以通过尺度关系式 $%5CPhi(2%5E%7B-(n%2B1)%7Dl)%3D%5Csum_%7Bk%5Cin%5Cmathbb%20Z%7Dp_k%5CPhi(2%5E%7B-n%7D(l%2B2%5Enk))%3B%5C%3Bl%5Cin%5Cmathbb%20Z$ 求得下一级二进制点的值.

有了尺度函数之后就可以使用 $%5Cphi(x)%3D%5Csum_%7Bk%5Cin%5Cmathbb%20Z%7D(-1)%5Ek%5Coverline%7Bp_%7B1-k%7D%7D%5CPhi(2x-k)$ 求出小波函数. 用二进制点表示, 并使用 1-k 替换 k 得到 $%5Cphi(2%5E%7B-(n%2B1)%7Dl)%3D%5Csum_%7Bk%5Cin%5Cmathbb%20Z%7D(-1)%5E%7Bk-1%7D%5Coverline%7Bp_k%7D%5CPhi(2%5E%7B-n%7D(l%2B2%5En(k-1)))%3B%5C%3Bl%5Cin%5Cmathbb%20Z$ , 这个方法对于两种函数近似都是通用的.

有限采样点

当采样点为有限时, 则有可能不能进行足够的分解步骤, 这时候就需要对采样点进行延拓. 常用的延拓方法有4种, 下面介绍一下, 假设采样点的索引是从 0 到 n.

零延拓 就是简单地把在采样点外的数据当作 0, 即 $s_k%3D%5Cleft%5C%7B%5Cbegin%7Bmatrix%7Ds_k%3B%5C%3Bk%5Cin%5B0%2Cn%5D%5C%5C0%3B%5C%3B%5Cmathrm%7Belse%7D%5Cend%7Bmatrix%7D%5Cright.$ . 适用于信号端点处无关紧要或者信号本身就是有始有终的.

周期延拓 就是把信号当作周期刚好是采样区间那么长的信号, 即 $s_k%3Ds_%7Bk%5C%25n%7D$ .

线性延拓 把信号断点处的两个采样点当作线性插值的两个样本, 向外面扩展, 即 $s_k%3D%5Cleft%5C%7B%5Cbegin%7Bmatrix%7Ds_k%3B%5C%3Bk%5Cin%5B0%2Ck%5D%5C%5Cs_0%2Bk(s_1-s_0)%3B%5C%3Bk%3C0%5C%5Cs_n%2Bk(s_n-s_%7Bn-1%7D)%3B%5C%3Bk%3En%5Cend%7Bmatrix%7D%5Cright.$ . 适用于在端点处信噪比比较高的信号.

对称延拓 就是在采样点端点处放一个镜子, 让采样点对称地扩展到外面. 对称延拓有两个不同的实现, 一个是把镜子放在整数点处, 另一个是把镜子放在半整数点处, 用式子写起来怪麻烦的, 可以看看我的源码.

数列索引

可以注意到, 在前3篇专栏里出现的累加绝大部分都是 $%5Csum_%7Bk%5Cin%5Cmathbb%20Z%7D$ , 因为计算机不能实现储存无限数据, 所以所有数列都应该是有限的 [除非有推导式].

这里出现的大部分累加都可以归纳为形式 $c_k%3D%5Csum_%7Bl%5Cin%5Cmathbb%20Z%7Da_%7B%5Calpha(l)%7Db_%7B%5Cbeta(k%2Cl)%7D$ , 其中 a, b 是已知的数列, c 是需要求出的数列, α, β 是相关的索引函数, 并且索引函数都是线性的. 使用 a⁻¹(l) 替换 l, 并重新记 β(k,a⁻¹(l)) 为 β(k,l), 得到 $c_k%3D%5Csum_%7Bl%5Cin%5Cmathbb%20Z%7Da_lb_%7B%5Cbeta(k%2Cl)%7D$ , 需要注意的是, 如果前者里 α(l) 里 l 的系数绝对值大于 1, 那么后者 l 的累加范围不再是整数, 而是二倍点, 三倍点之类的.

为了方便, 记 α₀ 为数列 a 里值不为零的最小索引, α₁ 为值不为零的最大索引, 类似地, 数列 b 有 β₀, β₁. 又因为假设了索引函数是线性的, 所以重新表达为 $c_k%3D%5Csum_%7Bl%3D%5Calpha_0%7D%5E%7B%5Calpha_1%7Da_lb_%7Buk%2Bvl%2Bw%7D$ . 由 b 的索引 uk+vl+w 可以推导出数列 c 里值不为零的最小最大索引 [本来这里想贴代码, 可是通用的实在太长了]. 求得了数列 c 的索引范围后就可以先申请内存再计算数值了. 实际上, b 的索引可以进一步限制 l 的累加范围. 下面给出 $c%5E%7B(n%2B1)%7D_k%3D%5Csum_%7Bl%5Cin%5Cmathbb%20Z%7Dp_lc%5E%7B(n)%7D_%7Bk-2%5Enl%7D$ 在 julia 里的实现作为例子