内容均来源于维基百科对db小波函数的介绍

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

一般而言的离散小波转换通常是以正交小波(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。

  • 分类方式

    多贝西小波的分类是以消失动量(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)中的每一层皆可写为下列方程式:

 

算法

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

假设输入的图片大小为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

分类:​

Logo

旨在为数千万中国开发者提供一个无缝且高效的云端环境,以支持学习、使用和贡献开源项目。

更多推荐