第4章 连续小波变换

从小波变换的数学理论来说，它是继傅立叶(Fourier)变换之后纯粹数学和应用数学完美结合的又一光辉典范，享有“数学显微镜”的美称。
小波分析是一种变分辨率的时频分析方法。在分析低频信号时，其时间窗很大；而分析高频信号时，其时间窗较小。这恰恰符合实际问题中高频信号持续时间短，低频信号持续时间长的自然规律，因而被誉为“数学显微镜”。小波分析被广泛应用于信号处理、图像处理、语音识别、模式识别、数据压缩、故障诊断、量子物理等领域。由于小波理论的算法复杂，编制相应的软件难度较大，因而大大限制了小波分析的应用。
MATLAB小波工具箱的出现避免了程序设计中的重复性劳动，缩短了开发周期，降低了成本，因而受到工科院校师生和研究人员的青睐。本章首先介绍小波分析的基本知识，并在其基础上逐渐引出相应MATLAB中关于连续小波变换的函数以及范例。
学习目标：
(1)了解小波分析的基本概念
(2)掌握连续小波变换及其性质
(3)熟练掌握函数cwt、cwtft和icwtft
4.1 小波分析简介
小波分析是近几十年来新兴的一门学科，是时频分析的主要手段，同时也应用于信号处理、图像处理和其他工程技术领域。
4.1.1 小波分析发展概述
从小波变换的数学理论来说，它是继Fourier变换之后纯粹数学和应用数学完美结合的又一光辉典范，享有“数学显微镜”的美称。从纯粹数学的角度来说，我们研究的小波变换是调和分析(其中包括函数空间、广义函数、Fourier分析和抽象调和分析等)这一重要学科大半个世纪以来的工作结晶；从应用科学和技术科学的角度来说，小波变换又是计算机应用、信号处理、图像分析、非线性科学和工程技术近些年来在方法上的重大突破。实际上，由于小波变换在它的产生、发展、完善和应用的整个过程中都广泛受惠于计算机科学、信号和图像处理科学、应用数学和纯粹数学、物理科学和地球科学等众多科学研究领域和工程技术应用领域的专家、学者和工程师们的共同努力，所以，现在它已成为科学研究和工程技术应用中涉及面极广泛的一个热门话题。
从小波变换的发展过程来说，大致可分成以下3个阶段。
第一阶段：孤立应用时期。其主要特征是一些特殊构造的小波在某些科学研究领域的特定问题上的应用。这个时期最典型的代表性工作是法国地球物理学家Morlet和Grossmann第一次把“小波”用于分析处理地质数据，引进了以他们的名字命名的时间—尺度小波，即Grossmann-Morlet小波。
同时，著名的计算机视觉专家Marr在他的“零交叉”理论中使用的可按“尺度大小”变化的滤波算子，现在称为“墨西哥帽”的小波也是这个时期有名的工作之一，这部分工作和后来成为Mallat的正交小波构造理论支柱之“多尺度分析”或“多分辨分析”有密切联系。这个时期一个有趣的现象是各个领域的专家、学者和工程师在完全不了解别人的研究工作的状态下巧妙地、独立地构造自己需要的“小波”。
虽然如此，但通观全局可以发现，这些专家、学者和工程师所从事研究的领域广泛分布于科学和技术研究的许多方面，因此，这个现象从另一个侧面预示了小波分析理论研究和应用热潮的到来，说明了小波理论产生的历史必然性。
第二阶段：国际性研究热潮和统一构造时期。真正的小波热潮开始于1986年，当时法国数学家Meyer成功地构造出具有一定衰减性质的光滑函数，这个函数(算子)的二进尺度伸缩和二进整倍数平移产生的函数系构成著名的4-范数函数空间的标准正交基。进入这个时期之后，Lemarie和Battle 又分别独立地构造了这样的“好的”小波。再后来Meyer和计算机科学家Mallat 提出多分辨分析概念，成功地统一了此前 Strömberg、Meyer、Lemarie和Battle的各别的小波构造方法。
同时，Mallat 还在多分辨分析的基础上简洁地得到了离散小波的数值算法，即现在的Mallat 分解和合成算法，并且将此算法用于数字图像的分解与重构。几乎同时，比利时数学家 Daubechies基于多项式方式构造出具有有限支撑的正交小波基和对称双正交小波，Chui和中国学者王建忠基于样条函数构造出单正交小波函数，并讨论了具有最好局部化性质的尺度函数和小波函数的构造方法。
第三阶段：全面应用时期。从1992年开始，小波分析方法进入全面应用阶段。在前一阶段研究工作的基础上，特别是数字信号和数字图像的Mallat分解和重构算法的确定，使小波分析的应用迅速波及科学研究和工程技术应用研究的几乎所有的领域。
鉴于小波分析的“自适应性质”和“数学显微镜性质”，使其被广泛用于基础科学、应用科学，尤其是信息科学、信号分析和方方面面，例如：图像处理与传输、信号处理、模式识别、地震探测、音乐、雷达、CT成像和彩色复印等。
由于小波分析在数字信号分析方面独特的诊断效果，来自不同学科、不同背景、不同兴趣爱好的科技工作者自觉投入到小波分析理论与应用研究中，涌现出一批高水平的论文和著作，在国内外形成一次又一次研究高潮，至今方兴未艾。
4.1.2 小波分析优缺点
由于小波变换的诸多优点，将其用于各种应用领域(如参数识别)的前景是广阔的，但仍有很多缺点制约了这一方法的广泛应用，主要表现在以下几个方面。
(1)连续小波变换的时间积分区间为正负无穷大，而实际采集到的包含模态信息的信号是有限长的，由于小波函数有一定的衰减区间，在对信号的中间部分进行变换时，结果比较满意，但对数据的首尾部分进行连续小波变换时，会产生较大的误差。当数据足够长时，可通过设定辨识区间，只选取置信区间内的数据进行识别来解决这一问题。
(2)但通常情况下，数据长度是一定的，一般无法达到令人满意的长度，这就造成了模态参数无法辨识，或错误辨识。同时，这种边缘效应会对低频产生较大的影响。
(3)通过连续小波变换识别模态参数，含有模态信息的有用变换结果将在时—频平面上形成一条特征曲线，称为小波脊线。提取这一特征曲线通常采用的方法是牛顿迭代法。这一方法虽然速度快，但对时频参数敏感，求解精度完全依赖于计算时所设定的时频参数采样率，而若想求解精度提高，则运算量会大幅增加。
4.2 连续小波变换及其性质
小波变换提出了变化的时间窗，当需要精确的低频信息时，采用长的时间窗；当需要精确的高频信息时，采用短的时间窗。
4.2.1 短时傅立叶变换
由于标准傅立叶变换只在频域里有局部分析的能力，而在时域里不存在这种能力，因此DennisGabor于1946年引入了短时傅立叶变换。短时傅立叶变换的基本思想是：把信号划分成许多小的时间间隔，用傅立叶变换分析每一个时间间隔，以便确定该时间间隔存在的频率。其表达式为：
其中：*表示复共轭，g (t)是有紧支集的函数，f(t)是进入分析的信号。在这个变换中，ejωt
 起着频限的作用，g(t)起着时限的作用。
随着时间τ的变化，g (t)所确定的“时间窗”在t轴上移动，使f(t)“逐渐”进行分析。因此，g (t)往往被称之为窗口函数，S (ω，τ)大致反映了 f(t)在时刻τ时、频率为ω的“信号成分”的相对含量。
这样信号在窗函数上的展开就可以表示为在[τ−δ，τ+δ]、[ω−ε，ω+ε]这一区域内的状态，并把这一区域称为窗口，其中，δ和ε分别称为窗口的时宽和频宽，表示时频分析中的分辨率，窗宽越小则分辨率就越高。
很显然，希望δ和ε都非常小，以便有更好的时频分析效果，但海森堡测不准原理指出δ和ε是互相制约的，两者不可能同时都任意小(事实上，
 为高斯函数时，等号成立)。
短时傅立叶变换和小波变换也是应传统的傅立叶变换不能够满足信号处理的要求而产生的。短时傅立叶变换分析的基本思想是：假定非平稳信号在分析窗函数 g(t)的一个短时间间隔内是平稳(伪平稳)的，并移动分析窗函数，使其在不同的有限时间宽度内是平稳信号，从而计算出各个不同时刻的功率谱。
但从本质上讲，短时傅立叶变换是一种单一分辨率的信号分析方法，因为它使用一个固定的短时窗函数，因而短时傅立叶变换在信号分析上还是存在着不可逾越的缺陷。
小波变换是一种信号的时间—尺度分析方法，它具有多分辨率分析的特点，而且在时、频两域都具有表征信号局部特征的能力，是一种窗口大小固定不变但其形状可改变、时间窗和频率窗都可以改变的时频局部化分析方法。即在低频部分具有较高的频率分辨率，在高频部分具有较高的时间分辨率和较低的频率分辨率，很适合于探测正常信号中夹带的瞬态反常现象并展示其成分，所以被誉为分析信号的显微镜。利用连续小波变换进行动态系统故障检测与诊断具有良好的效果。
4.2.2 一维连续小波变换
定义 设
 满足允许条件(完全重构条件或恒等分辨条件)：
此时称ψ(t)为一个基本小波或母小波。将母函数ψ(t)经伸缩和平移后得：
称其为一个小波序列。其中：a为伸缩因子，b为平移因子。对于任意的函数 f (t)∈L2
 (R)的连续小波变换为：
其重构公式(逆变换)为：
由于基小波ψ(t)生成的小波ψa，b
 (t)在小波变换中对被分析的信号起着观测窗的作用，所以ψ(t)还应该满足一般函数的约束条件：
故 是一个连续函数。这意味着，为了满足完全重构条件式，
 在原点必须等于0，即：
为了使信号重构的实现在数值上是稳定的，除了完全重构条件外，还要求小波ψ(t)的傅立叶变化满足下面的稳定性条件：
式中：
从稳定性条件可以引出一个重要的概念。
定义(对偶小波)若小波ψ(t)满足稳定性条件(4.8)式，则定义一个对偶小波
 其傅立叶变换
 由下式给出：
注意：稳定性条件(4.8)式实际上是对(4.9)式分母的约束条件，它的作用是保证对偶小波的傅立叶变换存在的稳定性。值得指出的是，一个小波的对偶小波一般不是唯一的，然而，在实际应用中，我们又总是希望它们是唯一对应的。因此，寻找具有唯一对偶小波的合适小波也就成为小波分析中最基本的问题。
连续小波变换具有以下重要性质。
(1)线性性：一个多分量信号的小波变换等于各个分量的小波变换之和。
(2)平移不变性：若f(t)的小波变换为Wf
 (a，b)，则 f (t−τ)的小波变换为：
(3)伸缩共变性：若f(t)的小波变换为Wf
 (a，b)，则f(ct)的小波变换为：
(4)自相似性：对应不同尺度参数a和不同平移参数b的连续小波变换之间是自相似的。
(5)冗余性：连续小波变换中存在信息表述的冗余度。小波变换的冗余性事实上也是自相似性的直接反映，它主要表现在以下两个方面。
①由连续小波变换恢复原信号的重构分式不是唯一的。也就是说，信号 f(t)的小波变换与小波重构不存在一一对应关系，而傅立叶变换与傅立叶反变换是一一对应的。
②小波变换的核函数即小波函数ψa，b
 (t)存在许多可能的选择(例如，它们可以是非正交小波、正交小波、双正交小波，甚至允许是彼此线性相关的)。
小波变换在不同的(a，b)之间的相关性增加了分析和解释小波变换结果的困难，因此，小波变换的冗余度应尽可能减小，它是小波分析中的主要问题之一。
4.2.3 高维连续小波变换
对f (t)∈L2
 (Rn
 )(n>1)，公式：
存在几种扩展的可能性，一种可能性是选择小波 f (t)∈L2
 (Rn
 )使其为球对称，其傅立叶变换也同样球对称：
并且其相容性条件变为：
对所有的 f，g∈L2
 (gn
 )。
这里，
 ，其中
 且b∈Rn
 ，公式(4.6)也可以写为：
如果选择的小波ψ不是球对称的，但可以用旋转进行同样的扩展与平移。例如，在二维时，可定义：
这里，a>0，b∈R2
 ， ，相容条件变为：
该等式对应的重构公式为：
对于高于二维的情况，可以给出类似的结论。
4.3 连续小波变换的计算
上节针对连续小波变换理论部分作了简要概述，本节将介绍连续小波变换计算的方法，如何将其转化成可编译的算法。
4.3.1 如何计算连续小波变换
从式(4.2)可以得出，连续小波变换计算分以下5个步骤进行。
(1)选定一个小波，并与处在分析时段部分的信号相比较。
(2)计算该时刻的连续小波变换系数C。如图4-1所示，C表示该小波与处在分析时段内的信号波形相似程度，C愈大，表示两者的波形相似程度愈高。小波变换系数依赖于所选择的小波。因此，为了检测某些特定波形的信号，应该选择波形相近的小波进行分析。
图4-1 计算小波变换系数示意图
(3)如图4-2所示，调整参数b，调整信号的分析时间段，向右平移小波，重复(1)～(2)步骤，直到分析时段已经覆盖了信号的整个支撑区间。
(4)调整参数a，尺度伸缩，重复(1)～(3)步骤。
(5)重复(1)～(4)步骤，计算完所有的尺度的连续小波变换系数，如图4-3所示。
图4-2 不同分析时段下的信号小波变换系数计算
图4-3 不同尺度下的信号小波变换系数计算
由小波变换的定义式(4.2)有：
其中：
为了可编程化实现上述公式，将其离散化执行，也就是说，把一组有限长度(N)的输入信号等时间间隔采样，即：
为了计算连续小波变换系数，对每个尺度参数 a 都要计算平移系数 b。如果利用反傅立叶变换公式，则得到：
上述公式在计算机中即可实现。
4.3.2 连续小波变换的应用
上节中a对应于信号，小波变换就像用镜头相对于目标平行移动，ψa，b
 (x)代表镜头所起的作用，b相当于使镜头相对于目标平行移动，a的作用相当于镜头向目标推进或远离。当a较大时，视野宽而分析频率低，可以进行平滑部分(概貌)的观察；当a较小时，视野窄而分析频率高，可以对细节进行观察。应用MATLAB小波工具箱时主要用到以下函数。
(1)连续小波变换函数cwt，其调用格式如下：
coefs=cwt(x，scales，'wname')
coefs=cwt(x，scales，'wname'，'plot')
coefs=cwt(x，scales，'wname'，'coloration')
其中：coefs返回的是小波系数矩阵，其行是尺度大小，而列是信号长度大小；x为被分析的信号；scales是需要计算的尺度范围；wname是所用的小波基；plot表示用图像方式显示小波系数。
(2)利用FFT(快速傅立叶变换)计算的连续小波变换，其调用格式如下：
cwtstruct=cwtft(sig)
cwtstruct=cwtft(sig，Name，Value)
cwtstruct=cwtft(...，'plot')
其中：cwtstruct是返回的小波结构；sig为被分析的信号；plot表示用图像方式显示小波系数。
【例4.1】将信号y=sin(2t)
 进行连续小波变换，选取meyer小波基并画出连续小波变换的系数分析图。MATLAB代码具体设置如下：
t=0:.2:2;
y=sin(2*t);
len=length(y);
cw1=cwt(y，1:32，'meyr'，'plot');
title('ContinuousTransform，absolutecoefficients.')
ylabel('Scale')
小波系数如图4-4所示。
图4-4 小波系数
【例4.2】作信号vonkoch的连续小波变换mesh图，选取小波基为db9。MATLAB程序设置如下：
Load vonkoch
vonkoch=vonkoch(1:510);
coefs=cwt(vonkoch，1:32，'db9');
colormap([0 0 0])
mesh(abs(coefs));
结果如图4-5所示。
图4-5 连续小波变换mesh图
注意：示例中使用的函数就是MATLAB小波工具箱中的连续小波变换函数，即coefs=cwt(x，scales，'wname')。其中：coefs为小波系数；x为被分析的信号；scales为小波变换尺度；wname为选取的小波。
【例 4.3】利用 Haar小波函数，计算一维信号 x(长度为1 000，除了在 500 点外，其他都为0)的连续小波变换。选取不同尺度参数，比较其信号的连续小波变换的系数。
MATLAB程序设置如下。
首先，画出Haar小波基函数的图像：
[～，psi，xval]=wavefun('haar'，10);
plot(xval，psi);
axis([0 1 -3.5 3.5]);
title('Haar Wavelet');
Haar函数如图4-6所示。
图4-6 Haar函数
接下来，输入一维信号：
x=zeros(1000，1);
x(500)=1;
画出连续小波变换系数图：
CWTcoeffs=cwt(x，1:128，'haar'，'plot');
colormap jet;
colorbar
连续小波变换系数图如图4-7所示。
比较尺度系数为10、50和90的小波系数：
subplot(311)
plot(CWTcoeffs(10，:)); title('Scale 10');
subplot(312)
plot(CWTcoeffs(50，:)); title('Scale 50');
subplot(313)
plot(CWTcoeffs(90，:)); title('Scale 90');
尺度比较如图4-8所示。
图4-7 连续小波变换系数图
图4-8 尺度比较
从图4-8中不难发现，当尺度越大时，视野越宽，反之亦然。
【例4.4】利用cwtft 函数计算一组长度为1 000的信号 y=y1
 +y2
 ，其中
在(0≤t<1)时，计算信号y连续小波变换的结果。
MATLAB程序设置如下：
N=1000;
t=linspace(0，1，N);
y=sin(2*pi*4*t).*(t<=0.5)+sin(2*pi*8*t).*(t>0.5);
cwtS1=cwtft(y，'plot');
结果如图4-9所示。
注：图4-9左下方显示为重构信号的关闭按键，即可显示或隐藏所要重构的信号。
图4-9 连续小波变换后生成的结构图
4.3.3 连续小波界面式应用实例
启动小波工具箱有两种方法：第一方法是在MATLAB窗口中键入wavemenu命令；第二种方法是，单击APPS，然后找到wavelet单击即可，具体如图4-10所示。
图4-10 小波工具箱启动
下面向读者介绍有关一维连续小波分析工具分析信号的具体操作步骤。
【例4.5】利用一维连续小波分析工具分析正弦曲线噪声信号。具体操作步骤如下。
(1)打开一维连续小波分析工具。首先，在窗口中键入wavemenu，弹出小波工具箱主界面，如图4-11所示。
图4-11 小波工具箱主界面
(2)在图4-11上单击Continuous Wavelet 1-D按钮，弹出一维连续小波的主界面，如图4-12所示。
图4-12 一维连续小波的主界面
(3)下载信号源。选择File→Load Signal命令，在弹出的对话框中选择noissin.mat 文件，它的路径为toolbox/wavelet/wavedemo，单击OK按钮。这样信号就被下载到此工具中，如图4-13所示。
图4-13 加载信号
(4)执行连续小波变换。在图4-13所示的界面中选择db4小波，并且尺度设置为从1到48。
(5)单击Analyze按钮。经过短暂的计算后，工具显示系数图，对应尺度a=24的系数图和最大尺度图，如图4-14所示。
图4-14 小波变换系数
(6)观看小波系数行。右击系数图，如图4-15所示，然后单击New Coefficients Line按钮观察Coefficients Line图的变化。
图4-15 小波系数行观测
(7)观看最大系数行。右击系数图，然后单击Refresh Maxima Lines按钮，如图4-16所示。
图4-16 最大系数行观测
在系数图中，鼠标移动的位置信息显示在工具屏幕下的Info框中，如图4-17所示。
图4-17 鼠标移动的位置信息
(8)尺度转换到伪频率。在工具屏幕的右下方选中Frequencies单选按钮，再次右击系数图，如图4-18所示。
图4-18 频率显示
(9)在Selected Axes 选项组中选中系数选项，如图4-19所示。
图4-19 图形数目显示
(10)放大细节。按鼠标左键画矩形框，如图4-20所示。
图4-18 频率显示
(9)在Selected Axes 选项组中选中系数选项，如图4-19所示。
图4-19 图形数目显示
(10)放大细节。按鼠标左键画矩形框，如图4-20所示。
图4-20 放大细节
(11)单击X+按钮，如图4-20所示，只水平放大，如图4-21所示。
图4-21 水平放大
(12)观察初始化系数图与当前系数图。在Coloration Mode 选项组中分别选中init +by scale+abs和current+by scale+abs 选项，结果如图4-22和图4-23所示。
图4-22 初始化模式系数显示
下面，向读者介绍有关一维连续小波分析工具(利用FFT解法)分析信号的实例。
【例4.6】利用一维连续小波分析工具分析正弦曲线噪声信号。具体操作步骤如下。
(1)打开一维连续小波分析工具。首先，在窗口中键入wavemenu，弹出小波工具箱主界面，如图4-11所示。
(2)在图4-11上单击Continuous Wavelet 1-D按钮，弹出一维连续小波的主界面。
图4-23 当前模式系数显示
(3)下载信号源。选择File→Load Signal命令，在弹出的对话框中选择noissin.mat 文件，它的路径为toolbox/wavelet/wavedemo，单击OK按钮。这样信号就被下载到此工具中，如图4-24所示。
图4-24 加载信号
(4)执行连续小波变换。在图4-24所示的界面中选择morl小波，并且参数设置为4。单击Analyze按钮。经过短暂的计算后，如图4-25所示。
图4-25 信号分析
(5)单击Close按钮，会出现更新合成图像(Update the syethesized image)对话框，单击yes按钮结束返回到小波分析工具箱。
注意：经过连续小波变换处理的信号其实是经过“复小波变换”成的，因此，在图4-25中出现了实部变换和虚部变换，还有相位等。
4.3.4 连续小波反变换的应用
公式(4.6)介绍了连续小波变换的反(逆)变换，即已知尺度和小波系数，即可重构信号。MATLAB小波工具箱主要用到以下函数：
xrec=icwtft(cwtstruct)
xrec=icwtft(cwtstruct，'plot')
xrec=icwtft(cwtstruct，'signal'，SIG，'plot')
其中：cwtstruct是返小波结构；plot表示用图像方式显示小波系数。
【例4.7】将例4.4中的小波变换后的系数按反变换重构信号。MATLAB程序设置如下：
sigrec=icwtft(cwtS1，'signal'，y，'plot');
结果如图4-26所示。
图4-26 重构信号
4.4 本章小结
本章主要结合连续小波分析的基本概念和基本原理，详细讨论了MATLAB小波工具箱中的相应函数，并结合一些实例来说明函数的用法。
读者通过本章的学习应该能够了解小波分析的基本方法理论和实用效果，并且通过操作掌握MATLAB提供的函数以从事一些基本操作，如对一维信号进行连续小波变换或反变换重构信号等。另外，读者也可以通过简单地利用连续小波的相关函数命令处理一些常见信号的基本问题。
本章的编写重点不在于让晦涩难懂的数学理论成为初学小波者的负担，而是结合理论与实际应用把简单的计算结果呈现给广大科技工作者和小波爱好者，使得读者学习完本章以后，能够借着MATLAB平台，加深小波的实用性，既不会脱离理论去搞实践，又不会因只实践不回馈理论。