小波变换中的多贝西小波(DB小波函数)概述
内容均来源于维基百科对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
参考资料
- A first course in Wavelets with Fourier Analysis, A.Boggess, F.J. Narcowich, 2001
分类:
更多推荐
所有评论(0)