第3章小波分析合成信号和图像

在传统的傅里叶分析中，信号完全是在频域展开的，不包含任何时域的信息，这对于某些应用来说是很恰当的，因为信号的频率信息对其非常重要。但其丢弃的时域信息可能对另外一些应用同样非常重要，所以人们对傅里叶分析进行了推广，提出了很多表征时域和频域信息的信号分析方法，如短时傅里叶变换、Gabor变换、时频分析、RandonWigner变换、小波变换等。
本章主要对小波分析进行介绍。
3.1小波变换
小波变换具有多分辨特性，也叫多尺度特性，可以由粗到精地逐步观察信号，也可以看成是用一组带通滤波器对信号滤波。通过适当地选择尺度因子和平移因子，可得到一个伸缩窗，只要适当地选择基本小波，就可以使小波变换在时域和频域都具有表征信号局部特征的能力。多分辨率分析与滤波器组相结合，丰富了小波变换的理论基础，拓宽了它的应用范围，对小波滤波器组的设计提出了更系统的方法，降低了小波变换的计算量。
小波变换的定义为：

Wf(a,b)=∫∞－∞f(t)ψa,b(t)dt=∫∞－∞f(t)a－12ψt－badt

其逆变换为：

f(t)=1Cψ∫∞－∞∫∞－∞a－2Wf(a,b)ψa,b(t)dadb

其中：

Cψ=∫∞－∞|ψ(ω)|2ωdω＜∞

式中，ψ(ω)为傅里叶变换，Cψ取有限值。

3.1.1连续小波变换
小波是通过对基本小波进行尺度伸缩和位移得到的，基本小波是一个具有特殊性质的实值函数，其振荡快速衰减，且在数学上满足积分为零的条件： 


∫∞－∞ψ(t)dt=0

Cψ=∫∞－∞ψ(s)2sds＜0


即基本小波在频域也具有很好的衰减性质。





一组小波基函数是通过尺度因子和位移因子由基本小波产生的：

ψa,b(x)=1aψx－ba

连续小波变换也称为积分小波变换，定义为：

Wf(a,b)=〈f,ψa,b(x)〉=∫∞－∞f(x)ψa,b(x)dx=1a∫∞－∞f(x)ψa,bx－badx

其逆变换为：

f(x)=1Cψ∫∞0∫∞－∞Wf(a,b)ψa,b(x)dbdaa2

二维连续小波基函数定义为：

ψabxby(x,y)=1|a|ψx－bxa,y－bya

二维连续小波变换为：

Wf(a,bx,by)=∫∞－∞∫∞－∞f(x,y)ψabxby
(x,y)dxdy

二维连续小波逆变换为： 

f(x,y)=1Cψ∫∞0∫∞－∞∫∞－∞Wf(a,bx,by)ψabxby(x,y)dbxdbydaa3

连续小波变换具有以下重要性质。
(1) 线性： 一个多分量信号的小波变换等于各个分量的小波变换之和。
(2) 平移不变性： 如果f(t)的小波变换为Wf(a,b)，则f(t－τ)的小波变换为Wf(a,b－τ)。
(3) 伸缩共变性： 若f(t)的小波变换为Wf(a,b)，则f(ct)的小波变换为1cWf(ca,cb),c＞0。
(4) 自相似性： 对应不同尺度参数a和不同平移参数b的连续小波变换之间是自相似的。
(5) 冗余性： 连续小波变换中存在信息表述的冗余度。
小波变换的冗余性事实上也是自相似性的直接反映，它主要表现在以下两个方面。
(1) 由连续小波变换恢复原信号的重构方式不是唯一的。也就是说，信号f(t)的小波变换与小波重构不存在一一对应关系，而变换与逆变换是一一对应的。
(2) 小波变换的核函数即小波函数ψa,b(t)存在许多可能的选择(例如，它们可以是非正交小波、正交小波、双正交小波，甚至允许是彼此线性相关的)。
小波变换在不同(a,b)之间的相关性增加了分析和解释小波变换结果的难度，因此小波变换的冗余度应尽可能减小，它是小波分析中的主要问题之一。
3.1.2一维离散小波变换
设ψ(t)∈L2(R)，其变换为ψ^(ω-)，当ψ^(ω)满足允许条件(完全重构条件或恒等分辨条件)：

Cψ=∫Rψ^(ω)2ωdω＜∞

时，称ψ(t)为一个基本小波或母小波。将母函数ψ(t)经伸缩和平移后得： 

ψa,b(t)=1|a|ψ(t－ba)a,b∈R;a≠0

称其为一个小波序列。其中a为伸缩因子，b为平移因子。对于任意的函数f(t)∈L2(R)的连续小波变换为： 

Wf(a,b)=〈f,ψa,b〉=|a|－1/2∫Rf(t)ψ(t－ba)dt

其重构公式(逆变换)为：

f(t)=1Cψ∫∞－∞∫∞－∞1a2Wf(a,b)ψ(t－ba)dadb

由于基小波ψ(t)生成的小波ψa,b(t)在小波变换中对被分析的信号起着观测窗的作用，所以ψ(t)还应该满足一般函数的约束条件：

∫∞－∞ψ(t)dt＜∞

故ψ^(ω)是一个连续函数。这意味着，为了满足完全重构条件式，ψ^(ω)在原点必须等于0，即：

ψ^(0)=∫∞－∞ψ(t)dt=0

为了使信号重构的实现在数值上是稳定的，除了完全重构条件外，还要求小波ψ(t)的变化满足下面的稳定性条件：

A≤∑∞-∞ψ^(2－jω)2≤B

式中0＜A≤B＜∞。
3.1.3二维离散小波变换
为了将一维离散小波变换推广到二维，只考虑尺度函数是可分离的情况，即： 

Φ(x,y)=Φ(x)Φ(y)

式中，Φ(x)是一维尺度函数，其相应的小波是ψ(x)，下列3个二维基本小波是建立二维小波变换的基础：

ψ1(x,y)=Φ(x)ψ(y),ψ2(x,y)=Φ(y)ψ(x),ψ3(x,y)=ψ(x)ψ(y)

它构成二维平方可积函数空间L2(R2)的正交归一化：

ψ1j,m,n(x,y)=2jψ1(x－2jm,y－2jn) j≥0,l=1,2,3(j,1,m,n都为整数)

1. 正变换
从一幅N×N的图像f(x,y)开始，其中上标指示尺度并且N是2的幂。对于j=0，尺度2j=20=1，也就是原图像的尺度。j值的每一次增大都使尺度加倍，而使分辨率减半。在变换的每一层次，图像都被分解为4个四分之一大小的图像，它们都是由原图与一个小波基图像内积后，再经过在行和列方向进行2倍的间隔抽样而生成的。对于第一个层次(j=1)，可写成：


f02(m,n)=〈f1(x,y),Φ(x－2m,y－2n)〉

f12(m,n)=〈f1(x,y),ψ1(x－2m,y－2n)〉

f22(m,n)=〈f1(x,y),ψ2(x－2m,y－2n)〉

f32(m,n)=〈f1(x,y),ψ3(x－2m,y－2n)〉


后续的层次(j＞1)，依次类推。如果将内积改写卷积形式则有：


f02j+1(m,n)=［f021(x,y)×Φ(－x,－y)］(2m,2n)

f12j+1(m,n)=［f021(x,y)×ψ1(－x,－y)］(2m,2n)

f22j+1(m,n)=［f021(x,y)×ψ21(－x,－y)］(2m,2n)

f32j+1(m,n)=［f021(x,y)×ψ31(－x,－y)］(2m,2n)


因为尺度函数和小波函数都是可分离的，所以每个卷积都可分解成行和列的一维卷积。例如，在第一层，首先用h0(－x)和h1(－x)分别与图像f(x,y)的每行作卷积并丢弃奇数列(以最左列为第0列)； 接着这个(N×N)/2矩阵的每列再与h0(－x)和h1(－x)相卷积，丢弃奇数行(以最上行为第0行)，结果就是该层变换所要求的4个(N/2)×(N/2)的数组。
2. 逆变换
逆变换与上述过程相似，在每层，首先在每列的左边插入一列0来增频采样前一层的4个矩阵； 接着用h0(x)和h1(x)来卷积各行，再成对地把这几个（N/2）×N的矩阵加起来； 然后通过在每行上面插入一行0来将刚才所得的两个矩阵的增频采样为N×N； 最后用h0(x)和h1(x)与这两个矩阵的每列卷积，这两个矩阵的和就是这一层重建的结果。
3.1.4多分辨率分析
下面简要介绍一下多分辨率分析的数学理论。
定义： 空间L2(R)中的多分辨率分析是指L2(R)满足如下性质的一个空间序列{Vj}j∈Z。
(1) 单调一致性： VjVj+1，对任意j∈Z。
(2) 渐进完全性： ∩j∈ZVj=Φ，close{∪j∈ZVj}=L2(R)。 
(3) 伸缩完全性： f(t)∈Vjf(2t)∈Vj+1。
(4) 平移不变性： k∈Z，(2－j/2t)∈Vjj(2－j/2t－k)∈Vj。
(5) Riesz基存在性： 存在(t)∈V0，使得{j(2－j/2t－k)|k∈Z}构成Vj的Riesz基。关于Riesz的具体说明为，如果(t)是V0的Riesz基，则存在常数A、B，且使得：

A‖{ck}‖22≤‖∑k∈Zck(t－k)‖22≤B‖{ck}‖22

对所有双无限可平方和序列{ck}，有：

‖{ck}‖22=∑k∈Z|ck|2＜∞

成立。
满足上述几个条件的函数空间集合称为一个多分辨率分析(MultiResolution Analysis)，如果(t)生成一个多分辨率分析，那么称(t)为一个尺度函数。
多分辨率分析构造一组函数空间，这组空间是相互嵌套的，即： 

…V－2V－1V0V1V2…

那么相邻的两个函数空间的差就定义了一个由小波函数构成的空间，即：

VjWj=Vj+1

在频域中，双尺度差分方程的表现形式为：

^(2ω)=H(ω)^(ω)

如果^(ω)在ω=0连续的话，则有：

^(ω)=∑∞j=1Hω2j^(0)

说明^(ω)的性质完全由^(0)决定。
3.1.5正交小波变换
在定义了多分辨率分析以后，把函数空间L2(R)分解成为一组全嵌套的函数空间集合，从这个角度上讲，当把信号从一个集合映射到它的子集的时候，就必然会丢失信息，那么如何描述这个丢失的信息呢？也就是说，如何找到一个子空间相对于其临近的父空间的补集？这就是下面要介绍的小波函数集。
设ψ=L2(R)是多分辨率分析{Vj}j∈Z的生成元，且满足： 
(1) {(x－n)}n∈Z。
(2) (x)=2∑n∈Zhn(2x－n)。
记：

gn=(－1)nh-1－n

并定义： 

ψ(x)=2∑n∈Zgn(2x－n)

定义Wj为ψj,n的线性张成，即：

Wj=span{ψj,n(x)}n∈Z

那么，函数空间Wj满足如下性质。
(1) Wj⊥Vj，WjVj=Vj+1，Wj⊥W1,j≠1。
(2) L2(R)=j∈ZWj。
(3) {ψj,n(x)}n∈Z是Wj的标准化正交基。
(4) {ψj,n(x)}n∈Z是L2(R)的标准化正交基。
首先，第一条性质要找的多分辨率分析每个尺度的补集就由Wj给出，而且Wj还是相互正交的。
第二条性质说明Wj的闭包可以覆盖整个能量有限信号空间，也就是说任何一个物理世界的信号都可以分解到这组空间上。
第三条和第四条说明由每个函数空间内部函数的正交性质，可以构造出一组Wj或L2(R)上的标准正交基，并可以把信号分解到任意尺度。
有了这些数学基础，接下来定义正交小波变换。
由于{ψj,n(x)}j,n∈Z是L2(R)的标准化正交基，所以对于f∈L2(R)，都有：

f(x)=∑j,n∈Zgn<f,ψj,n＞ψj,n(x)

称为f(x)的小波基数。
小波分析相对于传统的时频分析的优势就在于它可以在任意的时频分辨率上将信号分解，下面就从数学上看看任意层数信号分解的原理。
设f(x)∈V0，由于V0可以做如下分解：

V0=V－1W－1=Nj=1W－jV－N

那么f{x}就可以分解为在这些函数空间上的投影之和： 

f(x)=∑Nj=1g－j(x)+f－N(x)

其中：

g－j(x)∈W－j，f－N(x)∈V－N

f(x)分别在W－j和V－N上的投影为：


g－j(x)=∑n∈Z＜f,ψ－j,n＞ψ－j,n(x)

f－N(x)=∑n∈Z＜f,－N,n＞－N,n(x)


记：


A－2－Nf={f，-N,n}n∈Z

D－2－jf={f，ψ-j,n}n∈Z


则：

{{D－2－jf}Nj=1,A－2－Nf}

完全刻画了信号f(x)∈V0，成为f(x)∈V0的正交小波变换，同理，对f(x)∈Vj也可以同样求得其正交小波变换。
在实际应用到信号处理的时候，正交小波变换取样的方法与二进制小波变换不同，它是对连续信号在小波基上进行分解。
设实际处理的有限信号为{an}0≤n＜M，则假定相应的模拟信号为：

f=∑M－1n=0an(x－n)

那么，这个模拟信号在尺度20上的采样就是{an}0≤n＜M，所以二进制小波的信号本身在时(空)域就是离散的。而二进制小波变换是按Shannon定理对信号与小波函数的卷积进行取样，本质上在信号的时空域是连续的。而正交小波变换的结果每次也都是在不同尺度小波基上的展开系数，本质就是离散的。
3.1.6小波包分析
在多分辨率分析中，L2(R)=j∈ZWj表明多分辨率分析是按照不同的尺度因子j把Hilbert空间L2(R)分解为所有子空间Wj(j∈Z)的正交和的。其中，Wj为小波函数ψ(t)的闭包(小波子空间)。现在，对小波子空间Wj按照二进制分式进行频率的细分，以达到提高频率分辨率的目的。
一种自然的做法是将尺度空间Vj和小波子空间Wj用一个新的子空间Unj统一起来表征，若令：

U0j=Vj
U1j=Wjj∈Z

则Hilbert空间的正交分解Vj+1=VjWj可用Unj的分解统一为：

U0j+1=U0jU1jj∈Z(31)

定义子空间Unj是函数Un(t)的闭包空间，而Un(t)是函数U2n(t)的闭包空间，并令
Un(t)满足下面的双尺度方程： 

u2n(t)=2∑k∈Zh(k)un(2t－k)

u2n+1(t)=2∑k∈Zg(k)un(2t－k)(32)

式中，g(k)=(－1)kh(1－k)，即两系数也具有正交关系。当n=0时，以上两式直接给出：

u0(t)=∑k∈Zhku0(2t－k)

u1(t)=∑k∈Zgku0(2t－k)(33)

在多分辨率分析中，(t)和ψ(t)满足双尺度方程：

(t)=∑k∈Zhk(2t－k)
{hk}k∈Z∈l2

ψ(t)=∑k∈Zgk(2t－k){gk}k∈Z∈l2(34)

相比较，u0(t)和u1(t)分别退化为尺度函数(t)和小波基函数ψ(t)。式(32)是式(31)的等价表示。把这种等价表示推广到n∈Z+(非负整数)的情况，即得到式(32)的等价表示为：

Unj+1=UnjU2n+1jj∈Z； n∈Z+(35)

由式(32)构造的序列{un(t)}(其中n∈Z+)称为由基函数u0(t)=(t)确定的正交小波包。当n=0时，即为式(33)的情况。
由于(t)由hk唯一确定，所以又称{un(t)}n∈Z为关于序列{hk}的正交小波包。
设非负整数n的二进制表示为n=∑∞i=1εi2i－1，εi=0或1，则小波包
u^n(w)的变换由下式给出：

u^n(w-)=∏∞i=1mεi(w/2j)(36)

式中：


m0(w-)=H(w)=12∑+∞k=－∞h(k)e－jkw

m1(w-)=G(w)=12∑∞k=－∞g(k)e－jkw


设{un(t)}n∈Z是正交尺度函数(t)的正交小波包，则＜un(t－k),un(t－l)＞=δkl，即{un(t)}n∈Z构成L2(R)的规范正交基。
3.2利用小波分析合成信号和图像
小波工具箱中可利用小波分析合成信号和图像的相关算法，包括连续小波分析、小波相干性、同步压缩和数据自适应时频分析的算法等。
使用连续小波分析，可以研究光谱特征随时间演变的方式，在两个信号中识别常见的时变模式，并执行时间局域滤波。使用离散小波分析，可以分析不同分辨率的信号和图像，以检测变化点、不连续点和其他在原始数据中不容易看到的事件。可以在多个尺度上比较信号统计，并对数据进行分形分析以揭示隐藏的模式。
3.2.1可视化小波、小波包和小波滤波器
【例31】此实例展示了使用wfilters、wavefun和wpfun来获得滤波器、小波或与特定小波族对应的小波包。



>>%可以用wavefun2来分离离散二维小波

［LoD,HiD,LoR,HiR］ = wfilters('bior3.5');

subplot(2,2,1)

stem(LoD,'markerfacecolor',［0 0 1］); title('分解低通滤波器');

subplot(2,2,2)

stem(LoR,'markerfacecolor',［0 0 1］); title('重建低通滤波器');

subplot(2,2,3)

stem(HiD,'markerfacecolor',［0 0 1］); title('分解高通滤波器');

subplot(2,2,4)

stem(HiR,'markerfacecolor',［0 0 1］); title('重建高通滤波器');


运行程序，效果如图31所示。


图31分解与重建滤波器



下面代码显示可视化实值Morlet小波。



>> figure%效果如图32所示

［psi,xval］ = wavefun('morl');

plot(xval,psi,'linewidth',2)

title('$＼psi(x) = e^{-x^2/2} ＼cos{(5x)}$','Interpreter','latex',...

'fontsize',14);


得到了具有5个消失矩(sym4)的不对称Daubechies小波的前5个小波包。



>> ［wpws,x］ = wpfun('sym4',4,10);

for nn = 1:size(wpws,1)

subplot(3,2,nn)

plot(x,wpws(nn,:))

axis tight

title(［'W',num2str(nn-1)］);

end


运行程序，效果如图33所示。


图32Morlet小波




图33小波包



3.2.2时频分析和连续小波变换
连续小波变换(CWT)是一种时频变换，是分析非平稳信号的理想方法。一个信号是非平稳的意味着它的频域表示随时间而改变。许多信号都是非平稳的，如心电图、音频信号、地震数据和气候数据。
【例32】此实例说明如何利用连续小波变换的可变时频分辨率帮助获得清晰的时频表示。
1. 载入hychirp信号
加载一个有两个双曲啁啾的信号。数据以2048Hz采样，第一啁啾在0.1~0.68s内活跃，第二啁啾在0.1~0.75s内活跃。t时刻第一个啁啾的瞬时频率(Hz)为15π(0.8－t)22π； t时刻第二个啁啾的瞬时频率为5π(0.8－t)22π，并绘制信号。



>> load hychirp

plot(t,hychirp)

grid on

axis tight

xlabel('时间(s)');ylabel('幅值')

title('hychirp信号');


运行程序，效果如图34所示。
2. 时频分析： 傅里叶变换
傅里叶变换(FT)在识别信号中存在的频率分量方面非常出色，然而，FT不能识别何时出现频率分量。在放大0~200Hz的区域画出信号的幅度谱。



>> sigLen = numel(hychirp);

fchirp = fft(hychirp);

fr = Fs*(0:1/Fs:1-1/Fs);

plot(fr(1:sigLen/2),abs(fchirp(1:sigLen/2)),'x-')

xlabel('频率 (Hz)')

ylabel('幅值')

axis tight

grid on

xlim(［0 200］)


运行程序，效果如图35所示。


图34hychirp信号




图35信号的幅度谱



3. 时频分析： 短时傅里叶变换
傅里叶变换不提供时间信息，为了确定频率的变化何时发生，利用短时傅里叶变换(STFT)方法将信号分割成不同的块，并在每个块上执行FT。
选择窗口(段)大小是STFT的关键。对于使用STFT进行时频分析，选择较短的窗口有助于以牺牲频率分辨率为代价获得良好的时间分辨率； 相反，选择较大的窗口有助于以牺牲时间分辨率为代价获得良好的频率分辨率。一旦选择了窗口大小，它将在整个分析中保持固定。如果可以估计信号中期望的频率分量，那么就可以使用该信息选择一个窗口大小进行分析。

两种啁啾在其初始时间点的瞬时频率大约为5Hz和15Hz，下面代码使用helperPlot
Spectrogram函数绘制时间窗口大小为200ms的信号的谱图。



>> helperPlotSpectrogram(hychirp,t,Fs,200)%效果如图36所示


现在使用helperPlotSpectrogram函数绘制时间窗口大小为50ms的信号的谱图。信号中较晚出现的高频率现在被分解了，而信号开始处的低频率则没有被分解。



>> helperPlotSpectrogram(hychirp,t,Fs,50)%效果如图37所示



对于像双曲啁啾这样的非平稳信号，使用STFT是有问题的，因为没有一个单一的窗口大小可以解决这类信号的整个频率内容。


图36200ms的信号谱图




图3750ms的信号谱图



4. 时频分析： 连续小波变换
连续小波变换(CWT)的创建是为了克服STFT固有的分辨率问题。
平面的CWT平铺是非常有用的，世界的信号具有在长尺度上发生缓慢振荡，而在高频处往往是突然的或瞬态的。但是，如果高频事件持续时间长是自然的，那么使用CWT就不合适了。在没有获得任何时间分辨率的情况下，频率分辨率会更低。但通常情况下并非如此，人类的听觉系统就是这样工作的： 在低频有更好的频率局部化，在高频有更好的时间局部化。
下面绘制CWT的尺度图。尺度图是CWT的绝对值作为时间和频率的函数，图中使用对数频率轴，因为CWT中的频率是对数的。



>> cwt(hychirp,Fs)%效果如图38所示


从尺度图上可以清楚地看出信号中存在两个双曲线啁啾。使用CWT，可以准确地估计整个信号持续时间的瞬时频率，而不必担心选择一个片段长度。
图38中白色虚线标志着所谓的影响锥，影响锥显示尺度图中可能受边界效应影响的区域。要了解小波系数的大小增长有多快，可以使用函数helperPlotScalogram3d将scalogram绘制为一个三维表面。



>> helperPlotScalogram3d(hychirp,Fs)%效果如图39所示




图38CWT的尺度图




图39三维表面效果图



使用函数helperPlotScalogram来绘制信号的scalogram和瞬时频率，瞬时频率与尺度图特征可以很好地对齐。





图310尺度图和瞬时频率


>> helperPlotScalogram(hychirp,Fs);
%效果如图310所示



在以上代码中，调用到自定义的几个函数，下面为它们的源代码。



function helperPlotScalogram3d(sig,Fs)

%它可能会在未来的版本中更改或删除

figure

［cfs,f］ = cwt(sig,Fs);

sigLen = numel(sig);

t = (0:sigLen-1)/Fs;

surface(t,f,abs(cfs));

xlabel('时间(s)')

ylabel('频率(Hz)')

zlabel('幅值')

title('三维表面图')

set(gca,'yscale','log')

shading interp

view(［-40 30］)

end


function helperPlotSpectrogram(sig,t,Fs,timeres)

%它可能会在未来的版本中更改或删除

［px,fx,tx］ = pspectrum(sig,Fs,'spectrogram','TimeResolution',timeres/1000);

hp = pcolor(tx,fx,20*log10(abs(px)));

hp.EdgeAlpha = 0;

ylims = hp.Parent.YLim;

yticks = hp.Parent.YTick;

cl = colorbar;

cl.Label.String = 'Power(dB)';

axis tight

hold on

title(［'窗口大小为: ',num2str(timeres),' ms'］)

xlabel('时间(s)')

ylabel('Hz');

dt = 1/Fs;

idxbegin = round(0.1/dt);

idxend1 = round(0.68/dt);

idxend2 = round(0.75/dt);

instfreq1 = abs((15*pi)./(0.8-t).^2)./(2*pi);

instfreq2 = abs((5*pi)./(0.8-t).^2)./(2*pi);

plot(t(idxbegin:idxend1),(instfreq1(idxbegin:idxend1)),'k--');

hold on;

plot(t(idxbegin:idxend2),(instfreq2(idxbegin:idxend2)),'k--');

ylim(ylims);

hp.Parent.YTick = yticks;

hp.Parent.YTickLabels = yticks;

hold off

end



function helperPlotScalogram(sig,Fs)

%它可能会在未来的版本中更改或删除

［cfs,f］ = cwt(sig,Fs);

sigLen = numel(sig);

t = (0:sigLen-1)/Fs;

hp = pcolor(t,log2(f),abs(cfs));

hp.EdgeAlpha = 0;

ylims = hp.Parent.YLim;

yticks = hp.Parent.YTick;

cl = colorbar;

cl.Label.String = '幅值';

axis tight

hold on

title('尺度图和瞬时频率')

xlabel('秒');

ylabel('Hz');

dt = 1/2048;

idxbegin = round(0.1/dt);

idxend1 = round(0.68/dt);

idxend2 = round(0.75/dt);

instfreq1 = abs((15*pi)./(0.8-t).^2)./(2*pi);

instfreq2 = abs((5*pi)./(0.8-t).^2)./(2*pi);

plot(t(idxbegin:idxend1),log2(instfreq1(idxbegin:idxend1)),'k--');

hold on;

plot(t(idxbegin:idxend2),log2(instfreq2(idxbegin:idxend2)),'k--');

ylim(ylims);

hp.Parent.YTick = yticks;

hp.Parent.YTickLabels = 2.^yticks;

end


3.3离散小波分析
小波工具箱能够使用正交和双正交临界采样离散小波分析信号、图像和三维数据。临界采样离散小波分析又称为抽取离散小波分析，抽取离散小波分析最适用于数据压缩、去噪以及某类信号和图像的稀疏表示。在抽取离散小波分析中，尺度和平移是二进制的。
3.3.1一维小波去噪
【例33】此实例展示了如何使用离散小波分析去噪信号。



>>%创建一个参考信号

len = 2^11;

h = ［4-53-45-4.22.14.3-3.15.1-4.2］;

t = ［0.10.130.150.230.250.400.440.650.760.780.81］;

h = abs(h);

w = 0.01*［0.5 0.5 0.6 1 1 3 1 1 0.5 0.8 0.5］;

tt = linspace(0,1,len);

xref = zeros(1,len);

for j=1:11

xref = xref+(h(j)./(1+((tt-t(j))/w(j)).^4));

end

%以0.25的方差添加零均值高斯白噪声

>> rng default

x = xref + 0.5*randn(size(xref));

plot(x)%效果如图311所示

axis tight


使用具有4个消失矩的Daubechies最小非对称小波将信号降噪到3级。使用周期信号扩展模式 dwtmode('per')，利用Donoho和Johnstone的通用阈值选择规则，基于第一级小波变换系数进行软阈值处理，并将结果与参考信号一起绘制出来，以便进行比较。



>> dwtmode('per'); 

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!WARNING: Change DWT Extension Mode!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!



*****************************************

**DWT Extension Mode: Periodization**

*****************************************

>> xd = wdenoise(x,3,'Wavelet','sym4',...

'DenoisingMethod','UniversalThreshold','NoiseEstimate','LevelIndependent');

plot(xd)

axis tight

hold on

plot(xref,'r-.')%效果如图312所示

legend('降噪','参考')




图311带高斯白噪声的信号




图312参考与降噪效果图



3.3.2二维抽取离散小波分析
【例34】此实例展示了如何获取输入图像的二维DWT。



%加载并显示图像。图像由垂直、水平和对角线的图案组成

>> load tartan;

imagesc(X); colormap(gray);%效果如图313所示


使用双正交B样条小波和尺度滤波器获得1级二维小波，分析滤波器中有2个消失矩，合成滤波器中有4个消失矩。提取水平、垂直和对角的小波系数和近似系数，并显示结果。



>> ［C,S］ = wavedec2(X,1,'bior2.4');

［H,V,D］ = detcoef2('all',C,S,1);

A = appcoef2(C,S,'bior2.4');

subplot(221);

imagesc(A); title('近似系数');

colormap(gray);

subplot(222);

imagesc(H); title('水平系数');

subplot(223);

imagesc(V); title('垂直系数');

subplot(224);

imagesc(D); title('对角系数');


运行程序，效果如图314所示。


图313原始图像




图314小波的四个系数



可以看到，小波细节对输入图像中的特定方向非常敏感，近似系数是对原始图像的低通近似。
3.3.3非抽取离散小波分析
【例35】此实例展示了如何获得有噪声的调频信号的非抽取(平稳)小波变换。





图315原始信号与1级、3级小波系数及4级近似值

%加载含噪多普勒信号，得到平稳小波变换到4级

>> load noisdopp

swc = swt(noisdopp,4,'sym8');

%绘制原始信号与1级、3级小波系数及4级近似值

>> subplot(4,1,1)

plot(noisdopp)

ylabel('原始信号')

subplot(4,1,2);plot(swc(1,:))

ylabel('D1')

set(gca,'ytick',［］)

subplot(4,1,3);plot(swc(3,:))

ylabel('D3')

set(gca,'ytick',［］)

subplot(4,1,4);plot(swc(5,:))

ylabel('A4')

set(gca,'ytick',［］)

load noisdopp

swc = swt(noisdopp,4,'sym8');


运行程序，效果如图315所示。

3.4匹配的追踪
3.4.1图像的稀疏表征
将原始图像分割成若干个n×n的块，这些图像块就是样本集合中的单个样本y=Rn。在固定的字典上稀疏分解y后，得到一个稀疏向量。将所有的样本进行表征，可得原始图像的稀疏矩阵，重建样本y=Rn时，通过原子集合即字典D={di}ki=1∈Rn×m(n＜m)中少量元素进行线性组合即有：

y=Dx

其中，x={x1,x2,…,xm}∈Rm是y在D上的分解系数，也称为稀疏系数。
字典矩阵中的各个列向量被称为原子(Atom)，当字典矩阵中的行数小于甚至远小于列数时，即m≤n，字典D是冗余的。图像的稀疏表征数学模型如下： 

minx‖x‖0,s.t.y=Dx

稀疏表征不仅应具有过完备性，还应具有稀疏性。对于一个完备字典D，为了可以分解出更合适且稀疏的稀疏表征，应当含有更多的原子。
稀疏编码的目标就是在满足一定的稀疏条件下，通过优化目标函数，获取信号的稀疏系数。经典的算法有匹配追踪(Matching Pursuit，MP)、正交匹配追踪(Orthogonal Matching Pursuit，OMP)、基追踪(Basis Pursuit，BP)算法等。
考虑下面一个简单例子，给定稀疏信号： 

x=(－1.2,1,0)

字典矩阵A为A=(－0.7070.800.7070.6－1)，所以有y=Ax=1.65
－0.25。
现在，给定y=1.65
－0.25和A=－0.7070.80
0.7070.6－1，如何求得x呢？
3.4.2匹配追踪
在上面的例子中A的列向量称之为Basis(基)或者Atoms(原子)，所以有如下原子： 

b1=－0.707
0.707，b2=0.8
0.6，b3=0
－1

因为A=［b1b2b3］，如果令x=［abc］，则A·x=a·b1+b·b2+c·b3。A·x是原子b1，b2，b3的线性组合： 


A·x=－0.7070.80
0.7070.6－1－1.2
1
0

=－1.2×－0.707
0.707+1×0.8
0.6+0×0
－1

=y=－1.65
0.25


从上面的方程可以看出，b1对y值的贡献最大，然后是b2，最后是b3。匹配追踪算法刚好逆方向进行计算： 首先从b1，b2，b3中选出对y贡献最大的，然后从差值(residual)中选出贡献次大的，以此类推。
而贡献值通过内积(点积)进行计算，MP算法步骤如下。
(1) 选择对y值贡献最大的原子pi=maxj〈bj,y〉。 
(2) 计算差值ri=ri－1－pi·〈ri－1,pi〉。 
(3) 选择剩余原子中与ri内积最大的。
(4) 重复步骤(2)和(3)，直到差值小于给定的阈值(稀疏度)。
3.4.3匹配追踪实现
(1) 匹配追踪字典的创建和可视化。
【例36】此实例展示了如何创建和可视化包含二级Haar小波的字典。



>> ［mpdict,~,~,longs］ = wmpdictionary(100,'lstcpt',{{'haar',2}});

%使用longs参数按级别和函数类型(缩放或小波)划分小波字典

>> for nn = 1:size(mpdict,2)

if (nn <=longs{1}(1))

plot(mpdict(:,nn),'k','linewidth',2)

grid on

xlabel('转换')

title('Haar2级缩放')

elseif (nn>longs{1}(1) && nn<=longs{1}(1)+longs{1}(2))

plot(mpdict(:,nn),'r','linewidth',2)

grid on

xlabel('转换')

title('Haar 2级小波')

else

plot(mpdict(:,nn),'b','linewidth',2)

grid on

xlabel('转换')

title('Haar 1级小波')

end

pause(0.2)

end


运行程序，效果如图316所示。
(2) 用匹配追踪法分析电力消耗。
【例37】此实例展示了如何在离散傅里叶变换基础上比较匹配追踪和非线性近似。
实现中的数据是在24小时内收集的电力消耗数据。实例表明，通过从字典中选择向量，匹配追踪通常能够比任何单一基更有效地近似向量。
① 使用DCT、正弦字典和小波字典进行匹配追踪
加载数据集并绘制数据。该数据集包含35天的电力消耗，选择第32天进行进一步分析。数据是居中和缩放的，所以实际的使用单位是不相关的。



>> load elec35_nor

x = signals(32,:);

plot(x)%效果如图317所示

xlabel('分钟');ylabel('使用')




图316小波划分效果




图317绘制的数据图



由图317可看到，电力消耗数据包含平滑的振荡，间断于使用的突然增加和减少。下面将500~1200min的时间间隔放大。



>> xlim(［500 1200］)%效果如图318所示


从图318中可以看到缓慢变化的信号在大约650min、760min和1120min时的突变。在许多现实世界的信号中，像这些数据，感兴趣和重要的信息包含在瞬态中。对这些暂态现象进行模拟是很重要的。
下面从正交匹配追踪(OMP)字典中选择35个向量构造一个信号近似。字典包括Daubechies极值相位小波和第2级的缩放向量、离散余弦变换(DCT)、正弦基、Kronecker delta基以及Daubechies最小不对称相位小波和第1级与第4级的4个消失矩的缩放向量。然后，利用OMP求出电力消耗数据的最佳项。



>> dictionary = {{'db1',2},{'db1',3},'dct','sin','RnIdent',{'sym4',4}};

［mpdict,nbvect］ = wmpdictionary(length(x),'lstcpt',dictionary);

［y,~,~,iopt］ = wmpalg('OMP',x,mpdict);

plot(x)

hold on

plot(y,'-.')

hold off

xlabel('分钟');ylabel('使用')

legend('原始信号','OMP')


运行程序，效果如图319所示。从图319可以看到，用35个系数，正交匹配追踪既近似于信号的平滑振荡部分，也近似于用电量的突变。


图318放大的效果图



图319原始信号与OMP效果图


下面代码确定OMP算法从每个子字典中选择了多少个向量。



>> basez = cumsum(nbvect);

k = 1;

for nn = 1:length(basez)

if (nn == 1)

basezind{nn} = 1:basez(nn);

else

basezind{nn} = basez(nn-1)+1:basez(nn);

end

end

dictvectors = cellfun(@(x)intersect(iopt,x),basezind, ...

'UniformOutput',false);


大多数(60%)向量来自DCT和正弦字典，考虑到电力消耗数据总体上缓慢变化的性质，这是预期的行为。
② 使用DCT、正弦字典和与完整字典匹配追踪
只使用DCT和正弦字典重复OMP。设置OMP，从DCT、正弦字典中选择35个最佳向量，构造字典并执行OMP。注意，添加小波基字典可以更准确地显示用电量的突变，小波基的优点特别明显，特别是在接近大约650min和1120min使用上行和下行峰值时。



>> dictionary2 = {'dct','sin'};

［mpdict2,nbvect2］ = wmpdictionary(length(x),'lstcpt',dictionary2);

y2 = wmpalg('OMP',x,mpdict2,'itermax',35);

plot(x)

hold on

plot(y2,'-.','linewidth',2)%效果如图320所示

hold off

title('DCT和正弦字典')

xlabel('分钟');ylabel('使用')

xlim(［500 1200］)

>> figure

plot(x)

hold on

plot(y,'-.','linewidth',2)%效果如图321所示

hold off

title('完整字典')

xlabel('分钟');ylabel('使用')

xlim(［500 1200］)




图320DCT和正弦字典匹配追踪效果图




图321完整字典匹配追踪效果图


在离散傅里叶基下得到信号的最佳35项非线性近似，并获取数据的DFT。对DFT系数进行排序，并选取35个最大的系数。实值信号的DFT是共轭对称的，所以只考虑从0(DC)到奈奎斯特(1/2周期/分钟)的频率。



>> xdft = fft(x);

［~,I］ = sort(xdft(1:length(x)/2+1),'descend');

ind = I(1:35);


检查向量ind(没有一个指标对应于0或奈奎斯特)，加上相应的复共轭，得到DFT基下的非线性近似。



>> indconj = length(xdft)-ind+2;

ind = ［ind indconj］;

xdftapp = zeros(size(xdft));

xdftapp(ind) = xdft(ind);

xrec = ifft(xdftapp);

plot(x)

hold on

plot(xrec,'-.','LineWidth',2)

hold off

xlabel('分钟');ylabel('使用');

legend('原始信号','非线性DFT近似')


运行程序，效果如图322所示。
与DCT正弦字典类似，非线性DFT近似在匹配电能消耗数据的平滑振荡方面表现良好，然而非线性DFT近似并不能准确地逼近突变。图323是放大的包含消耗突变的数据间隔。



>> plot(x)

hold on

plot(xrec,'-.','LineWidth',2)

hold off

xlabel('分钟');ylabel('使用');

legend('原始信号','非线性DFT近似')

xlim(［500 1200］)







运行程序，效果如图323所示。



图322原始信号与非线性DFT近似效果图





图323放大数据间隔效果图



3.5时频分析
对于图像，连续小波分析显示了图像的频率内容如何在图像中变化，并有助于揭示有噪声图像中的模式。为了获得更清晰的分辨率,并从信号中提取振荡模式，可以使用小波同步压缩。
使用CQT，可以对带宽进行差分采样，对于较宽的频带分量使用更多的频率采样，对于窄带分量使用更少的频率采样。利用连续小波变换（CWT）可以得到两个信号之间的小波相干性； 可以把一个非线性的或非平稳的过程分解成它固有的振荡模态； 还可以重构信号的时频局部化近似。
3.5.1用小波相干性比较信号的时频含量
【例38】此实例说明了如何利用小波相干性和小波交叉谱来识别两个时间序列的时域共振荡行为。
在实例中也比较了小波相干性和交叉谱与它们的对应傅里叶变换。在程序中使用信号处理工具箱中的mscohere函数和cpsd函数来实现分析。
许多应用程序都在两个时间序列中识别和表征公共模式，在某些情况下，两个时间序列的共同行为是一个时间序列驱动或影响另一个时间序列的结果。在其他情况下，共同模式的结果是一些未观察到的机制影响两个时间序列。
对于联合平稳时间序列，在时间或频率上表征相关行为的标准技术是互相关、(傅里叶)互谱和相干。然而，许多时间序列是非平稳的，这意味着它们的频率内容随时间而变化，对于这些时间序列，需在时频平面上有一个相关或相干的度量。
可以使用小波相干性来检测非平稳信号中常见的时域局域振荡。在一个时间序列影响另一个时间序列的情况下，可以使用小波交叉谱的相位来确定两个时间序列之间的相对滞后。
创建一个由指数加权正弦波组成的信号，该信号有两个25Hz的分量： 一个以0.2s为中心，另一个以0.5s为中心。它还有两个70Hz的分量： 一个以0.2s为中心，另一个以0.8s为中心。第一个25Hz和70Hz的分量同时出现。



>> t = 0:1/2000:1-1/2000;

dt = 1/2000;

x1 = sin(50*pi*t).*exp(-50*pi*(t-0.2).^2);

x2 = sin(50*pi*t).*exp(-100*pi*(t-0.5).^2);

x3 = 2*cos(140*pi*t).*exp(-50*pi*(t-0.2).^2);

x4 = 2*sin(140*pi*t).*exp(-80*pi*(t-0.8).^2);

x = x1+x2+x3+x4;

plot(t,x)

grid on;

title('叠加信号')


运行程序，效果如图324所示。



%获取并显示CWT

>> cwt(x,2000);%效果如图325所示

title('使用默认Morse小波的解析CWT');

xlabel('时间(s)');ylabel('频率(Hz)')




图324叠加信号图




图325CWT效果图



通过将CWT系数归零，移除发生在0.07~0.3s内的25Hz分量。使用逆CWT(icwt)来重建信号的近似。



>> ［cfs,f］ = cwt(x,2000);

T1 = .07;T2 = .33;

F1 = 19;F2 = 34;

cfs(f > F1 & f < F2, t> T1 & t < T2) = 0;

xrec = icwt(cfs);

%显示重构信号的CWT，初始25Hz部分被移除

>> cwt(xrec,2000)%效果如图326所示

>> title('量级图');

>> xlabel('时间(s)');ylabel('频率(Hz)')

%绘制原始信号和重建图

>> subplot(2,1,1);

plot(t,x);

grid on;

title('原始信号');

subplot(2,1,2);

plot(t,xrec)

grid on;

title('第一个25Hz分量被移除的信号');


运行程序，效果如图327所示。


图326重构信号的CWT图




图327原始信号与重建信号





图328原始信号与重构信号的比较效果图

将重构信号与没有以0.2s为中心的25Hz分量的原始信号进行比较。




>> y = x2+x3+x4;

figure;

plot(t,xrec)

hold on

plot(t,y,'r--')

grid on;

legend('逆CWT近似','移除25Hz分量后的原始信号');

hold off


运行程序，效果如图328所示。
3.5.2时变一致性
傅里叶域相干性是一种在0到1尺度上测量两个平稳过程之间的线性相关关系的成熟技术。由于小波在时间和尺度(频率)上提供关于数据的局部信息，基于小波的相干性允许测量作为频率函数的时变相关性。换句话说，其是一种适用于非平稳过程的相干测量。
为了说明这一点，下面实例检测两名人体受试者的近红外光谱(NIRS)数据。近红外光谱是利用含氧和脱氧血红蛋白的不同吸收特性来测量大脑活动的，两组受试者的记录部位均为额叶上皮层，数据采样率均为10Hz。
【例39】检测近红外光谱数据的时变一致性。



%实验中，在被测试的一项任务上进行交替合作和竞争，任务的周期约为7.5s

>> load NIRSData;

figure

plot(tm,NIRSData(:,1))

hold on

plot(tm,NIRSData(:,2),'r-.')

legend('任务1','任务2','Location','NorthWest')

xlabel('秒')

title('NIRS数据')

grid on;hold off;


运行程序，效果如图329所示。
检查时域数据，如果不清楚在单个时间序列中存在什么振荡，或者两个数据集共有什么振荡，可用小波分析来回答这两个问题。
为了得到小波相干性作为时间和频率的函数，可以使用wcoherence函数输出小波相干性、跨谱、尺度。在本例中，利用函数helperPlotCoherence打包了一些用于绘制wcoherence函数输出的有用参数。



>> ［wcoh,~,f,coi］ = wcoherence(NIRSData(:,1),NIRSData(:,2),10,'numscales',16);

helperPlotCoherence(wcoh,tm,f,coi,'Seconds','Hz');%效果如图330所示

title('小波一致性');xlabel('秒');ylabel('Hz');



在图330中可以看到，在整个数据收集期间，一个高度一致的区域在1Hz左右，这是由两位受试者的心脏节律引起的。此外，还可以看到0.13Hz附近有强相干区域，这表明这项任务在受试者的大脑中引起了连贯的振荡。如果用周期而不是用频率来观察小波相干性会显得更自然。可以输入采样间隔，通过采样间隔，wcoherence函数提供了尺度到周期的转换。


图329NIRS数据效果图




图330小波一致性效果图







>> ［wcoh,~,P,coi］ = wcoherence(NIRSData(:,1),NIRSData(:,2),seconds(1/10),...

'numscales',16);

helperPlotCoherence(wcoh,tm,seconds(P),seconds(coi),'Time (secs)','Periods (Seconds)');

title('小波一致性');xlabel('时间(s)');ylabel('周期(s)');




图331采用频率实现小波一致性效果图

运行程序，效果如图331所示。

同样，注意记录中出现了与受试者心脏活动相对应的连续振荡，周期约为1s。与任务相关的活动也很明显，其周期大约为8s。总之，这个例子展示了如何利用小波相干性来寻找两个时间序列的时间局域相干振荡行为。对于非平稳信号，同时提供时间和频率(周期)信息的相干性度量更有用。
3.6离散多分辨率分析实现
离散小波变换(DWT)，包括最大重叠离散小波变换(MODWT)，将信号和图像分析成逐步细化倍频带。DODWT能够检测原始数据中不可见的模式，可以使用小波来获得信号的多尺度方差估计，或测量两个信号之间的多尺度相关性，还可以重建仅保留所需特征的信号(1D)和图像(2D)近似，并比较信号中跨频带的能量分布。小波包提供了一系列的变换，将信号和图像的频率内容分割成逐步细化的等宽间隔。
在3.1.4节中已介绍了多分辨率分析的相关概念，下面直接通过一个例子来演示离散多分辨率分析的实现。
【例310】离散多分辨率分析的实现。
这个例子展示了如何执行和解释基本的信号多分辨率分析(MRA)，例子中使用模拟数据和真实数据来回答诸如： 多分辨率分析意味着什么？执行多分辨率分析能获得关于信号的哪些见解？不同的MRA技术有哪些优点和缺点？ 
现实世界的信号是不同成分的混合，通常，只对这些组件的子集感兴趣，多分辨率分析允许通过将信号分成不同分辨率的组件来缩小分析范围。
以不同分辨率提取信号成分相当于在不同的时间尺度上分解数据的变化，或者等效地在不同的频带(不同的振荡速率)。因此，可以同时在不同的尺度或频带可视化信号变异性。
下面利用小波MRA分析和绘制合成信号，信号被分析在8个分辨率或级别。信号由三个主要成分组成： 频率为60个周期/s的时域振荡，频率为200个周期/s的时域振荡，以及一个趋势项。这里的趋势项也是正弦的，但是频率是每秒1/2个周期，所以它在1s的间隔内只完成1/2个周期。60周期/s或60Hz振荡发生在0.1~0.3s，而200Hz振荡发生在0.7~1s。



>> Fs = 1e3;

t = 0:1/Fs:1-1/Fs;

comp1 = cos(2*pi*200*t).*(t>0.7);%信号1

comp2 = cos(2*pi*60*t).*(t>=0.1 & t<0.3);%信号2

trend = sin(2*pi*1/2*t);%信号3

rng default

wgnNoise = 0.4*randn(size(t));

x = comp1+comp2+trend+wgnNoise;

mra = modwtmra(modwt(x,8));

helperMRAPlot(x,mra,t,'wavelet','小波MRA ',［2 3 4 9］)




图332合成的信号效果图


运行程序，效果如图332所示。

在不解释图332上的符号是什么意思的情况下，利用对信号的知识，试着理解小波MRA展示的是什么。如果从最上层的图开始，然后继续向下，直到到达原始数据的图，会看到组件变得越来越平滑。如果更喜欢从频率的角度来考虑数据，那么分量中包含的频率会越来越低。回想一下，原始信号有三个主要成分，一个200Hz的高频振荡，一个60Hz的低频振荡，以及一个趋势项，所有这些都被加性噪声破坏了。
如果观察D～2图，可以看到时域高频分量是孤立的，和研究这个重要的信号特征本质上是孤立的。D～3和D～4两幅图包含低频振荡，这是多分辨率分析的一个重要方面，即重要的信号成分可能不会被孤立在一个MRA成分中，但它们很少位于两个以上。最后，看到S8图包含趋势项。为了方便起见，这些组件中的轴的颜色已被更改，以在MRA中突出显示它们。如果希望可视化此图或后续的图，而不突出显示，请将最后的数字输入省略到helppermraplot中。
经验模态分解(EMD)是一种数据自适应的多分辨率技术，EMD递归地从数据中提取不同的分辨率，而不使用固定的函数或过滤器。EMD认为一个信号是一个快振荡叠加在一个慢振荡上，在快速振荡被提取后，处理过程将剩余的较慢分量作为新信号，并再次将其视为一个快速振荡叠加在一个较慢振荡上，这个过程会继续，直到达到某个停止条件。虽然EMD不使用小波等固定函数来提取信息，但EMD方法在概念上与小波方法非常相似，小波方法将信号分离为细节和近似，然后再将近似分离为细节和近似。EMD中的MRA成分称为本征模态函数(IMF)。



>> %绘制同一信号的EMD分析图

［imf_emd,resid_emd］ = emd(x);

helperMRAPlot(x,imf_emd,t,'emd','经验模态分解',［1 2 3 6］)%效果如图333所示


虽然MRA分量的数量不同，但EMD和小波MRA产生的信号图像是相似的，这不是偶然的。在EMD分解中，高频振荡主要局限于第一本构模态函数(IMF 1)，低频振荡主要局限于IMF 2，但在IMF 3中也能看到一些影响，IMF 6的趋势分量与小波技术提取的趋势分量非常相似。
另一种自适应多分辨率分析技术是变分模态分解(VMD)，像EMD一样，VMD试图在不使用固定函数进行分析的情况下从信号中提取固有模态函数或振荡模态。但是EMD和VMD以非常不同的方式决定模式，EMD递归地在时域信号上工作，逐步提取低频的IMF； VMD在频域识别信号峰值，并同时提取所有模式。



>> ［imf_vmd,resid_vmd］ = vmd(x);

helperMRAPlot(x,imf_vmd,t,'vmd','变分模态分解',［2 4 5］)%效果如图334所示




图333经验模态分解图




图334变分模态分解图



需要注意的关键是，与小波分解和EMD分解类似，VMD将三个感兴趣的成分分离为完全独立的模式或少量相邻的模式，这三种技术都允许在与原始信号相同的时间尺度上可视化信号分量。
3.7小波分析用于降噪
在小波分析中，应用最广泛的无疑是信号处理和图像处理，而在这两个领域中，应用最多的就是信号(图像)的降噪和压缩。
小波分析用于降噪的过程如下。
(1) 分解过程： 选定一种小波，对信号进行N层小波(小波包)分解。
(2) 作用阈值过程： 对分解得到的各层系数选择一个阈值，并对细节系数做软阈值处理。
(3) 重建过程： 降噪处理后的系数通过小波(小波包)重建恢复原始信号。
【例311】此实例展示了如何使用小波去噪信号和图像。
因为小波将数据中的特征定位到不同的尺度，所以可以在去除噪声的同时保留重要的信号或图像特征。小波去噪或小波阈值化背后的基本思想是，小波变换导致许多真实世界的信号和图像的稀疏表示，这意味着小波变换将信号和图像特征集中在几个大的小波系数中，若小波系数的数值很小，则通常是噪声，这时可以在不影响信号或图像质量的情况下“缩小”这些系数或去除它们。对系数设置阈值后，使用小波逆变换重构数据。
1. 消除干扰信号
为了说明小波降噪，创建一个有噪声的“颠簸”信号，在这种情况下，包含原始信号和含噪声信号。



>> rng default;

［X,XN］ = wnoise('bumps',10,sqrt(6));

subplot(211)

plot(X); title('原始信号');

AX = gca;

AX.YLim = ［0 12］;

subplot(212)

plot(XN); title('噪声信号');

AX = gca;

AX.YLim = ［0 12］;


运行程序，效果如图335所示。
使用函数wdenoise默认设置将信号降噪到4级，降噪使用抽取小波变换。将结果与原始信号一起绘制出来。



>> xd = wdenoise(XN,4);

figure;

plot(X,'k-.')

hold on;

plot(xd)

legend('原始信号','降噪信号','Location','NorthEastOutside')

axis tight;

hold off;


运行程序，效果如图336所示。

也可以使用非抽取小波变换去噪信号，使用非抽取小波变换将信号再次降噪到4级，并将结果与原始信号一起绘制出来。



>> xdMODWT = wden(XN,'modwtsqtwolog','s','mln',4,'sym4');

figure;

plot(X,'b-.')

hold on;

plot(xdMODWT)

legend('原始信号','降噪信号')

axis tight;hold off;




图335原始信号与噪声信号




图336信号降噪到4级效果图



运行程序，效果如图337所示。

从图336及图337可以看出，在这两种情况下，小波降噪消除了相当多的噪声，同时保留了信号中的尖锐特征。这是基于傅里叶变换去噪不具备的优势，在基于傅里叶的降噪或滤波中，会应用一个低通滤波器来去除噪声，然而当数据具有高频特征时，如信号中的峰值或图像中的边缘，低通滤波器会将这些平滑。
2. 降噪图像
也可以使用小波降噪图像。在图像中，边缘是图像亮度变化迅速的地方，图像降噪时保持边缘对感知质量至关重要。传统的低通滤波在去除噪声的同时，往往会使图像边缘平滑，从而影响图像质量。小波降噪能够在去除噪声的同时保留感知上的重要特征。



%加载一个有噪声的图像。使用默认设置的wdenoise2函数降噪图像，默认情况下，wdenoise2函数使

%用双正交小波bior4.4。若想显示原始图像和降噪图像，不要提供任何输出参数

>> load('jump.mat')

wdenoise2(jump)

subplot(121);title('原始图像');

subplot(122);title('降噪图像')


运行程序，效果如图338所示。


图337信号再降噪到4级效果图




图338图像降噪效果