多分辨率分析和连续小波变换多分辨率分析虽然时间和频率分辨率问题是由一个物理现象(海森堡测不准原理)造成的,还是可以用另外一种叫做多分辨率分析的方法来分析信号。如同它的名字一样,用不同的分辨率分析信号的不同频率分量。每个谱分量都不是由像STFT里那样用相同的方式解决。 多分辨率分析被用来在处理高频信号中获得一个好的时间分辨率和较差的频率分辨率,低频信号中获得较好的频率分别率和较高的时间分辨率。这个方法在处理高频信号持续时间较短,低频信号持续时间较长时才有用。幸运的是,实际应用中遇到的大多数信号都满足这一点。举例来说,下面这幅图就是这样的一种。在整个信号周期内,低频分量一直存在,高频信号只在中间很短的一段时间内出现了。
连续小波变换为了解决多分辨率的问题,小波变换作为替换短时傅立叶变换的一种方法被提出来。小波分析与短时傅立叶分析采用相同的处理方法,也是用窗口的形式将一个函数与信号相乘(小波变换中这个用来相乘的函数就是小波函数),变换结果被分成在时域内不同的片段。尽管如此,连续小波变换与短时傅立叶变换还是有两个主要的不同点:
离散小波变换的定义如下:
如上式,小波变换其实是一个含有两个自变量的函数,tau和s,分别作为平移和缩放参数。psi(t)是变换函数,叫做母小波。之所以叫做母小波,是基于小波分析的两个特征,如下所示: 这里“小波”的意思就是小的波形。满足这个条件的窗函数应该是有限宽度的(即紧支撑的),函数是振荡的。“母”这个字眼揭示了变换中用到的含有不同支撑域的函数,都可以追溯到一个主要函数——母小波。换句话说,母小波就是产生其他窗函数的原型。 “变换”这个词和在快速傅立叶变换中的变换是一样应用的。随着窗口在信号上平移,和窗口的位置相关。这个说法,在变换域内明显反映了时间信息。不过,不像之前说的短时傅立叶变换,我们没有频率参数。作为替代,我们有定义为frequency的缩放参数。频率这个词主要在STFT中使用。 尺度在小波分析中用到的这个缩放参数和地图中的缩放是类似的。在地图中,尺度高就意味着没有细节,只是一个整体的视图,比例低则意味着可以看到更多的细节。类似的,在频率中,低频(高尺度)一般是对信号的全景信息(一般在整个信号周期持续),高频(低尺度)则能给出信号中更多更详细的信息。下面四幅图展示了不同比例下的余弦信号:
幸运的是,在实际应用中,高频信号经常持续时间很短,不像图中展示的那样,他们经常以一种短时突变或尖峰的形式出现。而低频信号则在信号的整个周期内存在。 作为一个数学运算,缩放是对信号的放大或者缩小。大的尺度意味着对信号的扩大,小的尺度意味着对信号的缩小。图中所有的信号都是用同一个正弦信号变换而得,即他们是同一个函数的不同缩放版本。在上图中,s=0.05是最小的尺度,s=1是最大的尺度。 在数学函数中,如果给出一个函数f(t),如果s>1,f(st)就是f(t)的缩小版本,如果s<1,f(st)就是f(t)的放大版本。 尽管如此,在小波变换的定义中,缩放变量是作为分母存在的,因此,上面这些说明的相反一面就表明:s>1时是放大了信号,而s<1时则是压缩了信号。这种对缩放的解释将贯穿本文的始终。 连续小波变换的计算过程对上面公式的解释将在本节中进行详细说明。以x(t)作为被分析的信号。选用的小波作为信号处理中用到的所有窗函数的原型。应用的所有窗都是母小波的放大(或缩小)和平移版本。有很多函数可以满足这个条件。Morlet小波和墨西哥帽小波是其中最有代表性的,本章中后面部分中所举的例子也会用这两个小波进行小波分析。 一旦选择了母小波,就可以从s=1开始计算了,连续小波变换就是计算对应所有值的s,或者小于1,或者大于1。不过,与要分析的信号有关,一般不需要完整的变换。对所有的应用来说,信号是有带宽限制的,因此,在有限时间内做变换就经常能够满足要求了。在这篇文章里的后续部分,只用到s在有限时间内的值。 为方便起见,计算过程将会始于s=1,然后s值逐渐增大,即分析将会从高频开始,然后逐步到低频。s的第一个值是对大多数缩小的小波的反映。随着s增大,小波逐渐被放大。 小波要被放在信号的最初点,即t=0时刻。用尺度为1的小波函数与信号相乘,然后在所有时间内做积分。积分结果再乘上这个常量——1/sqrt(s)。后面这次相乘是为了使能量归一化处理而作的,其目的是使变换后的信号在任意比例上都有相同的能量。最终结果就是变换后的值,即在t=0时刻和s=1的情况下的连续小波变换。换句话说,就是在时间-尺度平面内,tau=0,s=1时刻信号的响应。 然后平移s=1时的小波值至t=tau时刻的位置,得到在时间-尺度平面内,t=tau,s=1时刻信号的响应。 重复这个过程,知道小波到达了信号的末端。这是,在时间-尺度平面内,就会得到一系列的点。 然后,将s增大一点。注意到,这是小波变换,所以tau和s都必须连续的增加。不过,如果是由计算机来进行这个变换过程,两个参数都以一个很小的步长增加。这是由于采样造成的。 对所有s值,重复上面这个过程。每一个根据给定的s计算的结果都对应时间-尺度平面内的一行。当对所应所有的s都计算完成后,对信号的连续小波变换就完成了。 下面的图说明了计算过程的每一步:
在图3.3中,显示了在四个不同时刻tau时的信号和小波函数。信号是图3.1所示的信号的裁剪版本,对应着信号的高频部分。可以看到它多么紧凑(蓝色的窗口)。它的宽度应该与信号中最高频分量出现的时间一致。图中显示了to=2,to=40,to=90和to=140时刻小波的位置。在每一个位置,都将它与信号相乘。很明显,只有在小波的支撑域内,乘积不为0,其余部分全为0。通过在时间轴上平移小波,信号被定位在时间轴上。进一步,通过改变s的值,信号又被定为在频率范围内。 如果信号对当前的s值有一个谱分量(这个例子中s的值是1),在信号的谱分量出现的时刻,信号与小波的乘积会相当大。如果对于当前的s值,谱分量不存在,那么,乘积将会很小或为0。图3.3中,在s=1和t=100ms附近,信号中存在一个窗口宽度的频谱分量。 图3.3中,在100ms时,对信号做连续小波变换后将会产生大的结果,在其他时候则值很小。另一方面,当尺度很高时,连续小波变换将对在整个信号周期内得到一个很大的值,因为低频信号在整个周期内都存在。
图3.4和3.5是s=4和s=5时的相同处理过程。注意到窗口宽度的改变是如何使频率分辨率降低的。随着窗口宽度的增大,变换将会夹杂一些低频分量。 结果,对每一个尺度和时间(间隔),都会得到时间——尺度平面内的一个点。同样比例时计算出来的结果作为时间——尺度平面的行,不同的比例计算出来的结果时间——尺度平面的列。 现在,看下面这个例子,看一下小波变换到底是什么样的。考虑图3.6中的非平稳信号。这个例子和在讲STFT时给出的例子很像,除了频率不同之外。如图所示,图中的信号是由四个频率分量组成的,分别为30Hz,20Hz,10Hz和5Hz。
图3.7是信号的连续小波变换。图中的轴分别是平移和缩放,不是时间和频率。不过,平移就对应着时间,因为它表明了母小波的位置。对母小波的平移可以看成是自从t=0时刻开始时间的流逝。而尺度,则是另外一种完全不同的情况了。还记得式3.1吗,s对应频率的反比。换句话说,考虑到分辨率问题时,无论何时我们谈到小波变换的性质,s的反比将会在图中展示时域信号的小波变换。
注意到图3.7中小的缩放反映了高频信息,频率随着尺度增加而降低。尺度为0附近的图形实际反映了信号最高频部分,最高尺度处的图形反映了信号的最低配部分。信号中最早出现的分量是30Hz(最高频分量)的,这幅图里则对应着最低尺度的部分。然后是20Hz的分量,第二高频的。5Hz的分量在轴上最后出现,对应着对大的尺度,正如我们所期望的。
现在,回忆一下分辨率特征:不像STFT那样,在信号的整个周期内分辨率都是固定的,小波变换的结果对高频信号具有高时间分辨率和低频率分辨率,对低频信号具有高频率分辨率和低时间分辨率。为了更好的揭示分辨率特征,图3.8用了一个不同的视角,显示了和图3.7一样的小波变换结果。在图3.8中,比例低的部分(比例低意味着比例的值更确切)反映了低频率分辨率。类似的,高比例对应着高的频率分辨率。 图3.7和3.8是归一化的信号,所以应该相比较进行评估。在平移坐标轴上个点的尖峰对应着,在尺度轴上个点对的频率(平移坐标轴和缩放坐标轴的单位不是秒和,他们仅仅是用来计算的采样点数)。 时间和频率分辨率这一节我们会更进一步的分析小波变换的分辨率特征。还记得,正是由于分辨率的问题,才使得我们快速傅立叶变换转到小波变换上。 图3.9经常被用来解释怎样诠释时间和频率分辨率。图3.9中的每个方块都反映了在时频平面内的小波变换结果。所有的方块都没有0区域,这说明在时频平面内,不能知道某个特定点对应的值。时频平面内每个方块中的点都代表了小波变换的一个结果。
更深入的看一下图3.9,虽然方块的宽高改变了,但是其面积却是常数。即每一个方块都代表时频平面内相同的部分,只不过时间和频率不相同。注意到在最低频处,方块的高度比较短(意味着频率分辨率最高,更容易知道确切的频率),但是最高频处的方块高度最高(意味着频率分辨率最低,更不易知道确切的频率)。在高频处,方块的宽度减小,时间分辨率提高了,随着方块高度降低,频率分辨率越来越低。 在总结这一节之前,我们应该再提一下快速傅立叶变换中的时频图是怎么样的。回想一下快速傅立叶变换中的时频分辨率是由窗口宽度决定的,一旦确定就不变了,因此时频分辨率也就固定了。因此它的时频平面是由方块组成的。 不考虑方块的维数,STFT和WT的结果都满足海森堡测不准原理。总结一下,方块的面积对STFT和CWT来说都是固定的,不同的窗函数将产生不同的面积。但是,不管怎样,方块的面积都至多在1/4\pi左右。根据海森堡测不准原理,我们不能将方块无限缩小。但在另一方面,对于一个给定的小波,方块的维数可以改变,但是其面积保持不变。这就是小波变换的结果。 小波理论:数学逼近这一节将会描述小波分析理论的主要概念,这些概念也可以被看成是大部分信号分析方法的准则。傅立叶定义的傅立叶变换是用一些基础函数来分析和重构一个函数。向量空间中的每一个向量都是向量基的线性组合,如把一些常数和向量相乘,然后计算点积。对信号的分析就包括估计这些常数(变换系数,傅立叶系数,小波系数等等)。合成或者说重构即计算线性组合方程。 所有有关这个题目的定义和理论都可以在凯瑟的书中找到,他的书对学习小波有很好的指导作用,也可以作为理解在小波理论中这些基础函数是如何起作用的入门级读物。因此,这些信息会出现在本节中。 向量基注:所有的公式中包含了很多希腊字母。这些字母都明确的用他们的名字写出,如tau,psi,phi等等。对大写字母,第一个字母都是大写的,如Tau,Psi,Phi等等。附注中的字母都会有下划线_。次方符号用^来表示。还要知道,用黑体写出的这些字母或字母名称都是代表性的向量,一些重点也以黑体字表示。但是其意思必须在上下文中明确。 向量空间V的基是一系列线性不相关向量,因为任何v都可以看作是向量基的线性组合。在一个空间中可能有超过一个向量基。但是,所有的空间都有相同的向量基数量,这个数量就是向量空间的维度。如二维空间的向量基有两个向量。
式3.2显示了向量v是如何由向量基b_k和相应系数nu^k线性组合而成的。 用向量的形式给出的这个原理可以用一个公式来表示,只需要把b_k替换为phi_k(t),把v替换为f(t),式3.2就变成了:
复指数函数(正弦和余弦)是傅立叶变换用到的基本公式,这些函数都是正交的,这些是进行重构必要的特性。 设f(t)和g(t)是L^2[a,b]空间内的两个函数(L^2[a,b]表示在整数区间[a,b]平方可积的一系列函数)。两个函数的内积定义如下式:
根据以上对内积的定义,连续小波变换可以被看成是测试信号与基本函数psi_(tau,s)(t)的内积:
其中
这个连续小波变换的定义说明了小波分析就是用来测量基本函数(即小波函数)与信号本身的相似性,这里的相似是指相似的频率分量。计算出来的连续小波变换系数说明了在当前尺度下原始信号与小波信号的近似程度。 这一点更进一步的阐明了前面关于确定尺度下原始信号与小波的相关性讨论。如果在当前比例下,信号中含有一个主要的频率分量,那么小波(基本函数)将会与信号中该频率分量出现的地方很相似。因此,这一点处的连续小波变换的系数将会是一个很大的数字。 内积,正交和正交归一化如果两个向量v和w的内积为0,则说它们是正交的:
类似的,如果两个函数的内积也为0,则可以说两个函数是正交的:
如果一个向量序列互相对偶正交,并且长度都为1,那么就说它们是正交归一化的。如下式:
类似的,一个函数序列phi_k(t),k=1,2,3…如果满足下面的公式也可说是正交归一化的:
且
相当于
其中delta_(kl)是克罗内克δ表示,定义如下:
如上所说,基本函数(或向量)可能不只一系列。在他们中间,正交函数基(或向量基)极其重要,因为它们在查找这些分析系数时表现出良好的特征。利用正交归一化性质,正交归一化基使得人们可以用一个更简单和直接的方法计算这些系数。
函数f(t)可以通过替换mu_k系数由式3.2a来重构。即:
在双正交基应用的场合正交归一化基不适用,而双正交基却是从正交归一化基总结出来的。这里的“双正交”意思是两个不同的基,彼此正交,但是二者都不是正交序列。 在坐标系统应用的场合,双正交基也不适用了。坐标系统是构成小波理论的一部分,感兴趣的读者可以去读前文提到的凯瑟的书。 和第二章讲解快速傅立叶变换一样的顺序,我们将会举一系列连续小波变换的例子。例子里给出的图都是从一个计算连续小波变换的程序而来。 在我们结束这一节之前,我将说一下两个应用最广的母小波。墨西哥帽小波被定义为高斯函数的衍生版本:
即
Morlet小波定义为:
其中a为调制参数,sigma为影响窗宽度的尺度参数。 CWT例子下面给出的所有例子均为现实生活中的非平稳信号。这些信号都来自包括正常人与阿尔茨海默氏症(Alzheimer)患者的事件相关电位(ERP)数据库。因为这些不是如简单正弦信号一样的测试信号,要解释它们并不容易,这里只是为了说明现实信号的连续小波变换(CWT)到底是怎样的。 图3.11所示为正常人的ERP信号:
图3.12为图3.11所示信号的CWT结果。坐标轴上的数字可以不必太关心。这些数字仅是表明,计算CWT时,在平移-尺度平面上有350个平移单位和60个尺度单位。需要特别说明的是,这里的计算并不是真正的连续小波变换。因为显然,上述结果是在有限个点上计算出来的。它仅仅是CWT的一种离散形式,这点稍后会再解释。同时还要说明的是,这更不是离散小波变换(DWT)。DWT的内容会在更后面介绍。
图3.13与图3.12一样同样是图3.11所示信号的CWT,只是观看的角度与图3.12有所不同。
图3.14是Alzheimer 症患者的ERP信号:
图3.15为图3.14所示信号的CWT结果。
图3.16与图3.15一样是图3.14所示信号的CWT结果是与图3.15一样的结果,只是观察的角度有所不同。
小波合成如果满足式3.18所示的条件,则CWT为可逆变换。幸运的是,这并不是一个非常苛刻的条件。只要满足式3.18所示的条件,即便基函数并不是归一化正交基。由小波系数计算原始信号值的小波重构过程可用如下公式计算:
其中 C_psi为与所用小波有关的常数。这个与小波重构过程有关的常数称为容许性常数。式3.18所示的小波重构条件称为容许性条件。
其中psi^hat(xi)为psi(t)的傅立叶变换。式3.18还表明:psi^hat(0) = 0,即:
由前面的讨论可知,式3.19并不是一个非常苛刻的条件,因为许多小波函数的积分值为0。满足式3.19所示条件的小波必定是振荡的。 小波级数:CWT的离散化如今,人们大量使用计算机来完成大数据量的运算。显然,无论是傅立叶变换(FT),短时傅立叶变换(STFT)还是连续小波变换(CWT),都需要用解析式、积分等方式来计算。于是在用计算机实现的过程中就会遇到离散化的问题。如果FT与STFT一样,最直观的做法是直接在时-频平面上进行采样。更直观地,对时-频平面进行均匀采样是最自然的选择。但是,在小波变换中,变化的尺度可以用来降低采样率。 在高尺度部分(即低频部分),根据奈奎斯特定理,采样率可以降低。换句话说,在时间-尺度平面上,如果可以用采样率N_1对尺度s_1进行采样,那么同样可以用采样率N_2对尺度s_2进行采样。其中s_1< s_2(对应频率f1>f2),并且N_2 < N_1。N_2 与 N_1之间的关系为:
或者用频率表示,可写为:
这意味着,在低频部分可以用较低的采样率进行采样,从而节省相当可观的运算量。 需要说明的是,如果仅考虑信号的分解,那么离散化的过程可以不受任何条件的限制。 如果不需要信号的合成,离散化的过程甚至都不需要满足奈奎斯特定理。但是如果还需要对信号进行重构,那么对离散化及采样频率的限制就变得非常重要。奈奎斯特采样频率是能够保证连续信号能够从离散信号完全重构的最小频率。正是因为这个原因,前面提到的基矢量才特别重要。 前面已经提到,小波psi(tau,s)如果满足式3.18所示的容许性条件,则能够利用式3.17完全恢复原始信号。对连续变换而言这是正确的。可问题是,如果我们在时间-尺度平面上进行了离散化,还能重构吗?回答是能够,但是必须满足一定的条件。 尺度参数 s首先以对数方式进行离散化。然后再在对应的尺度参数上对时间参数进行离散化。即不同的尺度上使用了不同的采样频率。这也就是说,采用了如图3.17所示的二进采样栅格来对时间-尺度平面进行采样。
考虑整个时间-尺度平面,连续小波变换的计算要在整个平面上逐点进行,因此,连续小波变换系数的数目为无穷多。离散化的过程首先考虑尺度轴。虽然尺度轴上的点数为无穷多,但利用对数规则,仅需要用到很少的一部分。对数的基可以根据需要选择。最常见的是选2,因为这样非常方便。如果选了以2为底的对数,即仅有2,4,8,16,32,64,…,等有限的一些尺度需要计算。当然,对数的基也可选为3,那样的话仅有3,9,27,81,243,…,等有限的一些尺度需要计算。在对尺度轴进行离散化之后再对时间轴进行离散化。如果选定对数的基为2,那么离散的尺度以2为因子变化。于是不同尺度上的时间采样率也同样以2为因子变化。 这里需要说明的是,图3.17中,在最小的尺度上(s=2),时间轴上仅采样到32点数据。在下一个尺度上,即s=4,时间轴上的采样率降低了2倍,因为尺度参数增加了2倍,于是在s=4这个尺度上,仅采样到16点数据。同理,再下一个尺度s=8上,仅有有8个采样点。 虽然常称为时间-尺度平面,实际上更准确的叫法是平移-尺度平面。因为变换域中的时间实际上对应着小波在时间上的平移。对小波级数而言,时间实际上仍然是连续的。 与傅立叶变换(FT)、傅立叶级数(FS)和离散傅立叶变换(DFT)之间的关系相同,同样也有连续小波变换(CWT)、半离散小波变换(即小波级数, WS)和离散小波变换(DWT)。 如果用数学公式来描述上述离散化过程,尺度参数离散为s = s_0^j,平移参数离散为tau = k * s_0^j * tau_0,其中 s_0>1,tau_0>0。 由此可以看出平移参数的离散化是如何依赖于尺度离散化参数s_0。 连续小波函数为:
将s = s_0^j,tau = k * s_0^j * tau_0代入上式,则小波函数变为:
如果{psi_(j,k)}为一组正交基,则小波级数变换变为:
或
小波级数需要{psi(j,k)}满足正交,或者双正交,或者框架条件。如果{psi(j,k)}不是归一化正交,则式3.24变为:
其中hat{ psi_{j,k}^*(t)}为二重双正交或者二重框架(此处*表示共轭)。 如果{ psi_(j,k) }为正交或者双正交,上述变换为非冗余的。如果{psi_(j,k) }为框架,则变换是冗余的。当然,从另外一方面讲,构造框架比发现正交或者双正交基要容易得多。 下面的类比可能会帮助更清晰地理解上述概念。将整个处理过程当作观察一个具体的物体,人眼首先是粗略地对物体进行观察,其程度与眼睛距物体的远近而有所不同。这对应着调整尺度参数s_0^(-j)。如果是近距离详细观察物体,j为负,并且绝对值很大(低尺度,高频率,对信号进行详细分析)。以很小的幅度缓慢地移动眼睛,对应着很小的tau = k * s_0^j * tau_0。如果j为负数并且绝对值很大时,tau则对应着时间上的微小变化(高采样频率),而s_0^-j则变化很大(低尺度,高频率,采样频率很高)。尺度参数也可以看作是放大。 采样频率到底最低到多少仍然能保证信号的完全重构?这是对上述过程进行最优化时必须回答的主要问题。从编程的角度来说,最方便的数值是 s_0=2,tau=1。显然,当采样频率尽可能低的时候,可能的正交小波数目也自然减少了。 前面给出的连续小波变换的例子实际上就是给定信号的小波级数。其参数的选择与信号有关。由于并不需要信号的重构,采样频率有时候远偏离于标准值。在不同的例子中,s_0的取值范围从2到10,tau_0的范围从2到8。 在结束小波级数讨论的时候,还有一个问题需要讨论。虽然离散化的小波变换,即小波级数可以利用计算机来运算,但运算时间从数秒到数小时不等。这主要依赖于信号的长度及所需的分辨率。计算信号的小波变换还有一种令人着迷的快速算法,即离散小波变换(DWT)。关于DWT的内容稍后介绍。 |