多贝西小波

多贝西小波（英语：Daubechies Wavelet），是以比利时女性物理暨数学家英格丽·多贝西（Ingrid Daubechies）的名字命名之一种小波函数，当初英格丽·多贝西发现了一种具有阶层（hierarchy）性质的小波，便将此小波以她的名字命名。多贝西小波主要应用在离散型的小波转换，是最常使用到的小波转换，通常使用在数位信号分析、信号压缩跟杂讯去除。

一般而言的离散小波转换通常是以正交小波（英语：Orthogonal_wavelet）（orthogonal wavelet）为基底，而多贝西小波也是一种正交小波。由于它很容易经由快速小波转换（fast wavelet transform（FWT））实现，所以常会放在数位信号处理的教科书中教学。

对于有限长度的小波，应用于快速小波转换（fast wavelet transform（FWT））时，会有两个实数组成的数列：一是作为高通滤波器的系数，称作小波滤波器（wavelet filter,也称为mother wavelet）；二是低通滤波器的系数，称作调整滤波器（scaling filter,也称为father wavelet）。

我们则以滤波器的长度N来形容滤波器为DN，例如：N=2的多贝西小波写作D2、N=4的多贝西小波写作D4，以此类推（N为偶数）。实际上常用的多贝西小波为D2到D20。与之对应的另一种描述方式是db N，其中N指消失矩的数量，由此，D4和db2是同一种多贝西小波。

性质

• 分类方式

  多贝西小波的分类是以消失动量（vanishing moment）的值A（亦为消失动量的个数）为依据（A称为tap），调整函式（scaling function）及小波函式（wavelet function）的平滑度（smoothness）皆会随著消失动量的值（tap）增加而增加：例如，当A=1时，多贝西小波即是哈尔小波（Haar wavelet），调整函式及小波函式都是不连续的；当A=2时，多贝西小波的调整函式及小波函式为不能平滑微分的连续函式；当A=3时，调整函式及小波函式已经是连续可微的函式了。以此类推，当A愈大时，多贝西小波的两个函式平滑度会愈来愈高。以下为多贝西小波跟不同A的调整及小波函式图：

尺度函数（蓝色）和小波函数（红色）
对应函数的频谱幅度

• 长度

  多贝西小波的长度为消失动量（vanishing moment）值A的两倍；所以当消失动量为A时，多贝西小波的小波滤波器（wavelet filter）及调整滤波器（scaling filter）长度皆为2A（N=2A）。一般而言，我们仍是以N来形容多贝西小波的长度：例如，当A=1时，有一个消失动量，多贝西小波写成D2，长度为2（也是Haar小波）；当A=2时，有两个消失动量，多贝西小波写成D4，长度为4；以此类推。但是，在MATLAB的使用上是以dbA描述多贝西小波，以下则为调整滤波器的系数及A的关系表：

正交多贝西小波系数

Scaling Coefficient ${\displaystyle p_{k}}$	db1（Haar）	db2	db3	db4	db5	db6	db7	db8	db9	db10
${\displaystyle p_{0}}$	1	0.6830127	0.47046721	0.32580343	0.22641898	0.15774243	0.11009943	0.07695562	0.05385035	0.03771716
${\displaystyle p_{1}}$	1	1.1830127	1.14111692	1.01094572	0.85394354	0.69950381	0.56079128	0.44246725	0.34483430	0.26612218
${\displaystyle p_{2}}$		0.3169873	0.650365	0.8922014	1.02432694	1.06226376	1.03114849	0.95548615	0.85534906	0.74557507
${\displaystyle p_{3}}$		-0.1830127	-0.19093442	-0.03957503	0.19576696	0.44583132	0.66437248	0.82781653	0.92954571	0.97362811
${\displaystyle p_{4}}$			-0.12083221	-0.26450717	-0.34265671	-0.31998660	-0.20351382	-0.02238574	0.18836955	0.39763774
${\displaystyle p_{5}}$			0.0498175	0.0436163	-0.04560113	-0.18351806	-0.31683501	-0.40165863	-0.41475176	-0.35333620
${\displaystyle p_{6}}$				0.0465036	0.10970265	0.13788809	0.1008467	6.68194092e-4	-0.13695355	-0.27710988
${\displaystyle p_{7}}$				-0.01498699	-0.00882680	0.03892321	0.11400345	0.18207636	0.21006834	0.18012745
${\displaystyle p_{8}}$					-0.01779187	-0.04466375	-0.05378245	-0.02456390	0.043452675	0.13160299
${\displaystyle p_{9}}$					4.71742793e-3	7.83251152e-4	-0.02343994	-0.06235021	-0.09564726	-0.10096657
${\displaystyle p_{10}}$						6.75606236e-3	0.01774979	0.01977216	3.54892813e-4	-0.04165925
${\displaystyle p_{11}}$						-1.52353381e-3	6.07514995e-4	0.01236884	0.03162417	0.04696981
${\displaystyle p_{12}}$							-2.54790472e-3	-6.88771926e-3	-6.67962023e-3	5.10043697e-3
${\displaystyle p_{13}}$							5.00226853e-4	-5.54004549e-4	-6.05496058e-3	-0.01517900
${\displaystyle p_{14}}$								9.55229711e-4	2.61296728e-3	1.97332536e-3
${\displaystyle p_{15}}$								-1.66137261e-4	3.25814671e-4	2.81768659e-3
${\displaystyle p_{16}}$									-3.56329759e-4	-9.69947840e-4
${\displaystyle p_{17}}$									-5.5645514e-5	-1.64709006e-4
${\displaystyle p_{18}}$										1.32354367e-4
${\displaystyle p_{19}}$										-1.875841e-5

• 滤波器

这边列出4到10点的filter，这些在实务上已经很够用，消失点(vanish moment)也就是k/2。

g会是低频的滤波器，h会是高频的滤波器。

可以看出${\displaystyle h[n]=(-1)^{n}g[-n]}$，这边采用python的语法，h会是g的反序，且n为基数时要乘上-1。

Coiflet filter ${\displaystyle g_{k},h_{k}}$	0	1	2	3	4	5	6	7	8	9
${\displaystyle h_{4}}$	0.4829629131445	0.836516303	0.22414386	-0.129409522
${\displaystyle g_{4}}$	-0.129409522551	-0.22414386	0.83651630	-0.482962913
${\displaystyle h_{6}}$	0.3326705529500	0.80689150	0.45987750	-0.135011020	-0.08544127	0.03522629
${\displaystyle g_{6}}$	0.0352262918857	0.08544127	-0.13501102	-0.459877502	0.806891509	-0.33267055
${\displaystyle h_{8}}$	0.2303778133088	0.71484657	0.63088076	-0.027983769	-0.187034811	0.030841381	0.032883011	-0.0105974017
${\displaystyle g_{8}}$	-0.0105974017850	-0.03288301	0.03084138	0.1870348117	-0.027983769	-0.630880767	0.714846570	-0.2303778133
${\displaystyle h_{10}}$	0.1601023979741	0.60382926	0.72430852	0.1384281459	-0.242294887	-0.032244869	0.077571493	-0.0062414902	-0.012580751	0.003335725
${\displaystyle g_{10}}$	0.0033357252854	0.012580751	-0.0062414	-0.077571493	-0.032244869	0.2422948870	0.138428145	-0.7243085284	0.6038292697	-0.16010239

建立

多贝西小波具有调整函式（低通滤波）及小波函式（高通滤波）两个函式。因此，我们需先建立调整函式及小波函式的系数:

首先，调整函数在多尺度分析（multi-resolution analysis）中的每一层皆可写为下列方程式：

${\displaystyle \phi (x)=\sum _{k=0}^{N-1}a_{k}\phi (2x-k)}$,

其中${\displaystyle (a_{0},\dots ,a_{N-1})}$为有限长度实数数列，称作调整系数。同时，小波函数也可以调整函数的线性组合表示：

${\displaystyle \psi (x)=\sum _{k=0}^{M-1}b_{k}\phi (2x-k)}$,

其中${\displaystyle (b_{0},\dots ,b_{M-1})}$亦为有限长度的实数数列，称做小波系数。

因为上述方程式必须是齐性的（homogeneous），在建立上，这两个函式会正规化（normalize）为和（sum）及平方和（sum of square）皆是2。

• 正交小波

  正交性质在此指调整系数就必须和位移偶数间隔后的调整系数互相垂直（内积为0），即为下式：

  ${\displaystyle \sum _{n\in \mathbb {Z} }a_{n}a_{n+2m}=2\delta _{m,0}}$.

  由于正交的特性，小波系数会满足下列条件：

  ${\displaystyle b_{n}=(-1)^{n}a_{N-1-n}}$, .

• 消失动量及多项式估计

  常用的多贝西小波为D2到D20，由于多贝西小波的消失动量为有限个，所以调整及小波系数可以表示为有限长度的多项式

  ${\displaystyle a[n]=a_{n}}$

  上式经过Z转换（Z-transform）后会变成：

  ${\displaystyle {A}(Z)=a_{0}+a_{1}Z+\dots +a_{N-1}Z^{N-1}}$,

  我们可以将上式转换为正交离散小波转换的一般表示式

  ${\displaystyle a(Z)=2^{1-A}(1+Z)^{A}p(Z)}$, ${\displaystyle order=A}$，此时，${\displaystyle p(1)=1}$、${\displaystyle p}$有实系数及${\displaystyle deg(p)=A-1}$。

  而正交的条件可写成

  ${\displaystyle a(Z)\,a(Z^{-1})+a(-Z)\,a(-Z^{-1})=4}$,或是等同于${\displaystyle (2-X)^{A}P(X)+X^{A}\,P(2-X)=2^{A}\dots }$(#),

  ${\displaystyle X}$定义为可以产生对称数列的劳伦兹多项式${\displaystyle X:=1/2\cdot (2-Z-Z^{-1})}$满足${\displaystyle X(-Z)=2-X(Z)}$

  因此${\displaystyle P(X)}$便成为对称型劳伦滋多项式，即${\displaystyle P(X(Z))=P(Z)P(Z-1)}$。因为 ${\displaystyle X(e^{iw})=1-cos(w)}$及${\displaystyle p(e^{iw})p(e^{-iw})=|p(e^{iw})|^{2}}$, ${\displaystyle {P}}$则会是区段[0,2]中的非负实数。

  方程式（#）如果除上${\displaystyle X}$的truncated power series则可求得对于每个${\displaystyle A}$的最小解

  ${\displaystyle P_{A}(X)=\sum _{k=0}^{A-1}\left({{A+k-1} \atop {A-1}}\right)2^{-k}X^{k}}$.（明显的值会是在（0,2）间的正数）

  而（#）的齐性方程式是一个对于${\displaystyle x=1}$的反对称方程式，因此可得一般解为${\displaystyle X^{A}(X-1)R((X-1)^{2})}$，此一般解有${\displaystyle R}$个多项式实系数。

  因此和为

  (sum)${\displaystyle P(X)=P_{A}(X)+X^{A}(X-1)R((X-1)^{2})}$

  ${\displaystyle {P}}$的值在区间[0,2]中并有界线（界线为${\displaystyle 4^{A-r}}$,）。为了将${\displaystyle r}$最大化的过程中会产生许多具有不等式条件的线性方程式。

  为了解出${\displaystyle P(X(Z))=P(Z)P(Z^{-1})}$的${\displaystyle P}$，这里使用Fejer-Riesz-algorithm这个方法（此为频谱分解的方法）。多项式${\displaystyle P(X)}$会因此分开成许多线性因子（linear factor）${\displaystyle P(X)=(X-\mu _{1})\dots (X-\mu _{N})}$,此时${\displaystyle N=A+1+2deg(R)}$。每一个线性因子代表可以分解成两个线性因子的一个劳伦兹多项式${\displaystyle (X(Z)-\mu )=-{\frac {1}{2}}Z+1-\mu -{\frac {1}{2}}Z^{-1}}$，任选其中一个线性因子都可设为${\displaystyle p(Z)}$。所以${\displaystyle p(Z)}$会有2^N个可能的答案。为了极端相位的目的，挑选所有根都是在单位圆上或是在单位圆内复数根的${\displaystyle P(Z)}$。

演算法

以下为示范小波转换应用于影像压缩，压缩后为原本图片的四分之一。

假设输入的图片大小为M*N，让图片对高频和低频进行卷积运算。

对M的基数进行取样，这个结果会让两个维度都变成(M/2)*N。

把低频的图片放在上面，高频的放在下面，低频的图片会长的像原本的图片，高频的图片会是只有灰色的图片。

对新的图片再进行高频和低频的卷积，这时变成M*(N/2)的大小，低频放在左边，高频放在右边，

最后可以看到[0:M/2, 0:N/2]就会是原本压缩过后的图片。

因小波转换有良好的性质，经过多次压缩还是能保有原本的资讯。也就是说可以修改以下的程式码改成回圈的方式，

进行多次小波转换，经过类似的模式再使用多次的反小波转换，还原出原本大小的图片。

import numpy as np
def subsampling(x, d):
    if d == 1:
        y = x[::2, :]
    elif d == 2:
        y = x[:, ::2]
    return y

def upsampling(x, d):
    s = x.shape
    if d == 1:
        y = np.zeros((p * s[0], s[1]))
        y[::2, :] = x
    elif d == 2:
        y = np.zeros((s[0], p * s[1]))
        y[:, ::2] = x
    return y
    
def cconv(x, h, d):
    if d == 2:
        return np.transpose(cconv(np.transpose(x), h, 1))
    y = np.zeros(x.shape)
    p = len(h)
    pc = int(round( float((p - 1) / 2 )))
    for i in range(0, p):
        y = y + h[i] * np.roll(x, i - pc, axis=0)
    return y
    
def DWT(image, h, g): # discrete wavelet transformation
    fW = image.copy()
    j = int(np.log2(image.shape[0])-1)
    A = fW[:2**(j+1):,:2**(j+1):]
    Coarse = subsampling(cconv(A,h,1),1)
    Detail = subsampling(cconv(A,g,1),1)
    A = np.concatenate( (Coarse, Detail), axis=0 )
    Coarse = subsampling(cconv(A,h,2),2)
    Detail = subsampling(cconv(A,g,2),2)
    A = np.concatenate( (Coarse, Detail), axis=1 )
    fW[:2**(j+1):,:2**(j+1):] = A
    return fW
    
def iDWT(image, fW, h, g): #image is original, fW is after DWT of that.
    f1 = fW.copy()
    j = int(np.log2(image.shape[0])-1)
    A = f1[:2**(j+1):,:2**(j+1):]
    Coarse = A[:2**j:,:]
    Detail = A[2**j:2**(j+1):,:]
    h1 = h[::-1]
    g1 = g[::-1]
    Coarse = cconv(upsampling(Coarse,1),h1,1)
    Detail = cconv(upsampling(Detail,1),g1,1)
    A = Coarse + Detail
    Coarse = A[:,:2**j:]
    Detail = A[:,2**j:2**(j+1):]
    Coarse = cconv(upsampling(Coarse,2),h1,2)
    Detail = cconv(upsampling(Detail,2),g1,2)
    A = Coarse + Detail
    f1[:2**(j+1):,:2**(j+1):] = A
    return f1

参照

• 小波分析

参考资料

1. A first course in Wavelets with Fourier Analysis, A.Boggess, F.J. Narcowich, 2001