多贝西小波

• 分类方式

  多贝西小波的分类是以消失动量（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的关系表：

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

首先，调整函数在多尺度分析（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，让图片对高频和低频进行convolution。

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

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

对新的图片再进行高频和低频的convolution，这时变成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