长期以来,离散小波变换(Discrete Wavelet Transform)在数字信号处理、石油勘探、地震预报、医学断层诊断、编码理论、量子物理及概率论等领域中都得到了广泛的应用。各种快速傅氏变换(FFT)和离散小波变换(DWT)算法不断出现,成为数值代数方面最活跃的一个研究领域,而其意义远远超过了算法研究的范围,进而为诸多科技领域的研究打开了一个崭新的局面。本章分别对FFT 和DWT 的基本算法作了简单介绍,若需在此方面做进一步研究,可参考文献[2]。
1.1 离散小波变换DWT
1.1.1 离散小波变换DWT 及其串行算法
先对一维小波变换作一简单介绍。设f (x )为一维输入信号,记)2(2)(2/k x x j j jk -=--φφ,
)2(2)(2/k x x j j jk -=--ψψ,这里)(x φ与)(x ψ分别称为定标函数与子波函数,)}({x jk φ与
)}({x jk ψ为二个正交基函数的集合。记P 0f =f ,在第j 级上的一维离散小波变换
DWT(Discrete Wavelet Transform)通过正交投影P j f 与Q j f 将P j -1f 分解为:
∑∑+=+=-k
k
jk j k jk j k j j j d c f Q f P f P ψφ1
其中:∑
=-=-+1
1
2)(p n j n k j
k c n h c ,∑=-=-+1
1
2)(p n j n k j k c n g d )12,...,1,0,,...,2,1(-==j N k L j ,这里,{h (n )}与{g (n )}分别为低通与高通权系数,它们由基函数)}({x jk φ与)}({x jk
ψ
来确定,p 为权系数
的长度。}{0
n C 为信号的输入数据,N 为输入信号的长度,L 为所需的级数。由上式可见,每级一维DWT 与一维卷积计算很相似。所不同的是:在DWT 中,输出数据下标增加1时,权系数在输入数据的对应点下标增加2,这称为“间隔取样”。
算法 一维离散小波变换串行算法 输入:c 0
=d 0
(c 00
, c 10
,…, c N-10
)
h=(h 0, h 1,…, h L-1)
g=(g 0, g 1,…, g L-1)
输出:c i j
, d i j
(i=0, 1,…, N/2j-1
, j ≥0)
Begin (1)j=0, n=N (2)While (n ≥1) do for i=0 to n-1 do (2.1.1)c i j+1
=0, d i j+1
=0 (2.1.2)for k=0 to L-1 do
j n i k k j i j i j
n i) (k k j i j i d g d d c h c c mod )2(11mod 211,
+++++++=+=
end for
end for
j=j+1, n=n/2
end while
End
显然,算法的时间复杂度为O(N*L)。
在实际应用中,很多情况下采用紧支集小波(Compactly Supported Wavelets ),这时相应的尺度系数和小波系数都是有限长度的,不失一般性设尺度系数只有有限个非零值:
h 1,…,h N ,N 为偶数,同样取小波使其只有有限个非零值:g 1,…,g N 。为简单起见,设尺度系
数与小波函数都是实数。对有限长度的输入数据序列:M n x c n n
,,2,1,0
Λ==(其余点的值都看成0),它的离散小波变换为:
k n Z
n j n j k h c c 21-∈+∑=
k n Z
n j n j k g c d 21-∈+∑=
1,,1,0-=J j Λ
其中J 为实际中要求分解的步数,最多不超过log 2M ,其逆变换为
k n Z
k j
k k n Z
k j
k j n h c h c c 221
-∑∈-∑∈-+
=
1,,ΛJ j =
注意到尺度系数和输入系列都是有限长度的序列,上述和实际上都只有有限项。若完全
按照上述公式计算,在经过J 步分解后,所得到的J +1个序列1,,1,0,-=J j d j
k Λ和的非
零项的个数之和一般要大于M ,究竟这个项目增加到了多少?下面来分析一下上述计算过程。
j =0时计算过程为
k n M
n n k g x d 211-=∑=
不难看出,的非零值范围为:,12,,0,1,,12-??
?
???-+-
=M N k ΛΛ即有??
?
???-+=??????+--
=21212N M M N k 个非零值。的非零值范围相同。继续往下分解时,非零项出现的规律相似。分解多步后非零项的个数可能比输入序列的长度增加较多。例如,若输入序列长度为100,N =4,则有51项非零,有27项非零,有15项非零,有9项非零,有6项非零,有4项非零,有4项非零。这样分解到6步后得到的序列的非零项个数的总和为116,超过了输入序列的长度。在数据压缩等应用中,希望总的长度基本不增加,这样可以提高压缩比、减少存储量并减少实现的难度。
可以采用稍微改变计算公式的方法,使输出序列的非零项总和基本上和输入序列的非零项数相等,并且可以完全重构。这种方法也相当于把输入序列进行延长(增加非零项),因而称为延拓法。
只需考虑一步分解的情形,下面考虑第一步分解(j =1)。将输入序列作延拓,若M 为偶数,直接将其按M 为周期延拓;若M 为奇数,首先令。然后按M +1为周期延拓。作了这种延拓后再按前述公式计算,相应的变换矩阵已不再是H 和G ,事实上这时的变换矩阵类似于循环矩阵。例如,当M =8,N =4时矩阵H 变为:
2
1
4
3
432143214321
21430
0000000000000000h h h h h h h h h h h h h h h h h h h h 当M =7,N =4时矩阵H 变为:
1
4
3
321432143211430
00000000000000h h h h h h h h h h h h h h h h h
从上述的矩阵表示可以看出,两种情况下的矩阵内都有完全相同的行,这说明作了重复计算,因而从矩阵中去掉重复的那一行不会减少任何信息量,也就是说,这时我们可以对矩阵进行截短(即去掉一行),使得所得计算结果仍然可以完全恢复原输入信号。当M =8,N =4时截短后的矩阵为:
??????
????????=43
2
1
4321432
1
2143
000000000000h h h h h h h h h h h h h h h h H 当M =7,N =4时截短后的矩阵为:
??????
????????=3
2
1
4321432
1
143
0000000000h h h h h h h h h h h h h h H 这时的矩阵都只有??
?
???2M 行。分解过程成为:
向量C 1 和D 1
都只有??
????2M 个元素。重构过程为:
可以完全重构。矩阵H ,G 有等式
H *H +G *G =I
一般情况下,按上述方式保留矩阵的??
?
???2M 行,可以完全恢复原信号。
这种方法的优点是最后的序列的非0元素的个数基本上和输入序列的非0元素个数相同,特别是若输入序列长度为2的幂,则完全相同,而且可以完全重构输入信号。其代价是得到的变换系数D j
中的一些元素已不再是输入序列的离散小波变换系数,对某些应用可能是不适合的,但在数据压缩等应用领域,这种方法是可行的。
Begin
对所有处理器my_rank(my_rank=0,…, p -1)同时执行如下的算法:
(1)j =0;
(2)while (j 将数据)12/,,1,0(-=j j n N n c Λ按块分配给p 台处理器 将处理器i +1中前L -1个数据发送给处理器i 处理器i 负责1 +j n c 和)12 )1(,,2(1 11 -+=+++j j j n p N i p N i n d Λ的计算 j =j +1 end while End 这里每一步分解后数据1+j n c 和1+j n d 已经是按块存储在P 台处理器上,因此算法第一步中的数据分配除了j =0时需要数据传送外,其余各步不需要数据传送(数据已经到位)。因此,按LogP 模型,算法的总的通信时间为:2(L max(o ,g )+l),远小于计算时间O(N)。 MPI 源程序请参见所附光盘。 [1]. Verlag,1982 [2]. 陈崚.二维正交子波变换的VLSI 并行计算.电子学报,1995,23(02):95-97 基于小波变换的人脸识别 近年来,小波变换在科技界备受重视,不仅形成了一个新的数学分支,而且被广泛地应用于模式识别、信号处理、语音识别与合成、图像处理、计算机视觉等工程技术领域。小波变换具有良好的时频域局部化特性,且其可通过对高频成分采取逐步精细的时域取样步长,从而达到聚焦对象任意细节的目的,这一特性被称为小波变换的“变聚焦”特性,小波变换也因此被人们冠以“数学显微镜”的美誉。 具体到人脸识别方面,小波变换能够将人脸图像分解成具有不同分辨率、频率特征以及不同方向特性的一系列子带信号,从而更好地实现不同分辨率的人脸图像特征提取。 4.1 小波变换的研究背景 法国数学家傅立叶于1807年提出了著名的傅立叶变换,第一次引入“频率”的概念。傅立叶变换用信号的频谱特性来研究和表示信号的时频特性,通过将复杂的时间信号转换到频率域中,使很多在时域中模糊不清的问题,在频域中一目了然。在早期的信号处理领域,傅立叶变换具有重要的影响和地位。定义信号(t)f 为在(-∞,+∞)内绝对可积的一个连续函数,则(t)f 的傅立叶变换定义如下: ()()dt e t f F t j ωω-? ∞ -∞ += (4-1) 傅立叶变换的逆变换为: ()()ωωπ ωd e F t f t j ? +∞ ∞ -= 21 (4-2) 从上面两个式子可以看出,式(4-1)通过无限的时间量来实现对单个频率 的频谱计算,该式表明()F ω这一频域过程的任一频率的值都是由整个时间域上的量所决定的。可见,式(4-1)和(4-2)只是同一能量信号的两种不同表现形式。 尽管傅立叶变换可以关联信号的时频特征,从而分别从时域和频域对信号进行分析,但却无法将两者有效地结合起来,因此傅立叶变换在信号的局部化分析方面存在严重不足。但在许多实际应用中,如地震信号分析、核医学图像信号分析等,研究者们往往需要了解某个局部时段上出现了哪个频率,或是某个频率出现在哪个时段上,即信号的时频局部化特征,傅立叶变换对于此类分析无能为力。 因此需要一种如下的数学工具:可以将信号的时域和频域结合起来构成信号的时频谱,描述和分析其时频联合特征,这就是所谓的时频局部化分析方法,即时频分析法。1964年,Gabor 等人在傅立叶变换的基础上引入了一个时间局部化“窗函数”g(t),改进了傅立叶变换的不足,形成窗口化傅立叶变换,又称“Gabor 变换”。 定义“窗函数”(t)g 在有限的区间外恒等于零或很快地趋于零,用函数(t )g -τ乘以(t)f ,其效果等同于在t =τ附近打开一个窗口,即: ()()()dt e t g t f G t j f ωττω-+∞ ∞--=?, (4-3) 式(4-3)即为函数f(t)关于g(t)的Gabor 变换。由定义可知,信号(t)f 的Gabor 变换可以反映该信号在t =τ附近的频谱特性。其逆变换公式为: ()()()ττωτωπ ωd G t g e d t f f t j ,21 ? ?+∞ ∞ --- = (4-4) 可见()τω,f G 的确包含了信号(t)f 的全部信息,且Gabor 窗口位置可以随着 τ的变化而平移,符合信号时频局部化分析的要求。 虽然Gabor 变换一定程度上克服了傅立叶变换缺乏时频局部分析能力的不 10.2小波变换的基本原理 地质雷达的电磁波信号和地震波信号都是非平稳随机时变信号,长期以来,因非平稳信号处理的理论不健全,只好将其作为平稳信号来处理,其处理结果当然不满意。近年来,随着科学技术的发展和进步,国内外学术界已将注意力转向非平稳随机信号分析与处理的研究上,其中非平稳随机信号的时频表示法是研究热点之一。在这一研究中,戈勃展开、小波变换、维格纳分布与广义双线性时频分布等理论发展起来,这些方法既可以处理平稳信号过程,也可以处理非平稳随机时变信号。 小波变换是上世纪80年代中后期逐渐发展起来的一种数学分析方法。1984年法国科学家J.M OLET在分析地震波的局部特性时首先使用了小波这一术语,并用小波变换对地震信号进行处理。小波术语的含义是指一组衰减震动的波形,其振幅正负相间变化,平均值为零,是具有一定的带宽和中心频率波组。小波变换是用伸缩和平移小波形成的小波基来分解(变换)或重构(反变换)时变信号的过程。不同的小波具有不同带宽和中心频率,同一小波集中的带宽与中心频率的比是不变的,小波变换是一系列的带通滤波响应。它的数学过程与傅立叶分析是相似的,只是在傅立叶分析中的基函数是单频的调和函数,而小波分析中的基函数是小波,是一可变带宽内调和函数的组合。 小波变换在时域和频域都具有很好的局部化性质,较好地解决了时域和频域分辨率的矛盾,对于信号的低频成分采用宽时窗,对高频成分采用窄时窗。因而,小波分析特别适合处理非平稳时变信号,在语音分析和图象处理中有广泛的应用,在地震、雷达资料处理中将有良好的应用前景。 下边就小波分析的基本原理、主要作用及在雷达资料处理中的应用三方面作以介绍。 10.2.1小波分析的基本原理 小波函数的数学表达 第三章 离散傅立叶变换 一、离散傅立叶级数 计算题: 1.如果)(~n x 是一个周期为N 的周期序列,那么它也是周期为2N 的周期序列。把)(~n x 看 作周期为N 的周期序列有)(~ )(~1k X n x ?(周期为N );把)(~n x 看作周期为2N 的周期序列有)(~)(~2k X n x ?(周期为2N );试用)(k X 1~表示) (k X 2~ 。 二、离散傅立叶变换定义 填空题 2.某DFT 的表达式是∑-==10 )()(N k kl M W k x l X ,则变换后数字频域上相邻两个频率样点之间的间隔是( )。 3.某序列DFT 的表达式是∑-==1 0)()(N k kl M W k x l X ,由此可看出,该序列的时域长度是 ( ),变换后数字频域上相邻两个频率样点之间隔是( )。 4.如果希望某信号序列的离散谱是实偶的,那么该时域序列应满足条件 ( )。 5.采样频率为Hz F s 的数字系统中,系统函数表达式中1 -z 代表的物理意义是 ),其中时域数字序列)(n x 的序号n 代表的样值实际位置是( ); )(n x 的N 点DFT )k X (中,序号k 代表的样值实际位置又是( ) 。 6.用8kHz 的抽样率对模拟语音信号抽样,为进行频谱分析,计算了512点的DFT 。则频域 抽样点之间的频率间隔f ?为_______,数字角频率间隔w ?为 _______和模拟角频率间隔 ?Ω ______。 判断说明题: 7.一个信号序列,如果能做序列傅氏变换对它进行分析,也就能做DFT 对它进行分析。 ( ) 计算题 8.令)(k X 表示N 点的序列)(n x 的N 点离散傅里叶变换,)(k X 本身也是一个N 点的序列。 第三章 离散小波变换 3.1 尺度与位移的离散化方法 减小小波变换系数冗余度的做法是将小波基函数?? ? ??-= a t a t a τψψτ1)(,的 τ,a 限定在一些离散点上取值。 1. 尺度离散化:一种最通常的离散方法就是将尺度按幂级数进行离散化, 即取m m a a 0=(m 为整数,10≠a ,一般取20=a )。如果采用对数坐标,则尺度a 的离散取值如图3.1 所示。 图3.1 尺度与位移离散方法 2. 位移的离散化:当120==a 时,()τψψτ-=t t a )(,。 (1)通常对τ进行均匀离散取值,以覆盖整个时间轴。 (2)要求采样间隔τ满足Nyquist 采样定理,即采样频率大于该尺度下频率通带的2倍。 3. )(,t a τψ=? 当m 增加1时,尺度增加一倍,对应的频带减小一半(见图2.2),可见采样频率可以降低一半,即采样间隔可以增大一倍。因此,如果尺度0=m 时τ的间隔为s T ,则在尺度为m 2时,间隔可取s m T 2。此时)(,t a τψ可表示为 );(2212221 ,t T n t T n t n m s m m m s m m ψψψ记作??? ???-=??? ? ???- Z n m ∈, 为简化起见,往往把t 轴用s T 归一化,这样上式就变为 ()n t t m m n m -=-- 22)(2 ,ψψ (3.1) 4. 任意函数)(t f 的离散小波变换为 ??=R n m f dt t t f n m WT )()(),(,ψ (3.2) DWT 与CWT 不同,在尺度—位移相平面上,它对应一些如图3.1所示的离散的点,因此称之为离散小波变换。将小波变换的连续相平面离散化,显然引出两个问题: (1)离散小波变换>=<)(),(),(,t t f n m W T n m f ψ是否完全表征函数)(t f 的全部信息,或者说,能否从函数的离散小波变换系数重建原函数)(t f 。 (2)是否任意函数)(t f 都可以表示为以)(,t n m ψ为基本单元的加权和 ∑∈= Z n m n m n m t C t f ,,,)()(ψ?如果可以,系数n m C ,如何求? 上述两个问题可以归结为一个。假设条件(1)满足,可合理的选择ψ,并对τ,a 进行适当的离散(即适当的选择s T a ,0),那么一定存在与小波序列n m ,ψ对 应的n m ,~ψ序列,使得问题(1)的重建简单地表示为 ∑∈><= Z n m n m n m f t f ,,,~,)(ψψ (3.3) n m ,~ψ称为n m ,ψ的对偶,它可以由一个基本小波)(~t ψ通过位移和伸缩取得: () n t t m m n m -=--2~2)(~2,ψψ 由上式,若存在)()(2R L t g ∈,则有 ∑>><<=><>= 小波变换与傅里叶变换 的对比异同 IMB standardization office【IMB 5AB- IMBK 08- IMB 2C】小波变换详解
小波变换的基本原理
第三章 离散傅立叶变换
第三章 离散小波变换
小波变换与傅里叶变换的对比异同