雷达原理---时频分析--3.小波变换-3.1基础知识

一、短时傅里叶变换的缺陷

窄窗口时间分辨率高、频率分辨率低；宽窗口时间分辨率低、频率分辨率高。

对于时变的非平稳信号，高频适合小窗口，低频适合大窗口。然而STFT的窗口是固定的，在一次STFT中宽度不会变化，所以STFT还是无法满足非稳态信号变化的频率的需求。

二、小波变换的优点

   小波的优点，其实就是上面傅立叶缺点的解决：

   对于加窗傅立叶变换让人头疼的就是窗口的大小问题，如果我们让窗口的大小可以改变，不就完美了吗？答案是肯定的，小波就是基于这个思路。但是不同的是，STFT是给信号加窗，分段做FFT；而小波变换并没有采用窗的思想，更没有做傅里叶变换。小波直接把傅里叶变换的基给换了—将无限长的三角函数基换成了有限长的会衰减的小波基。这样不仅能够获取频率，还可以定位到时间了。

  这个基函数会伸缩、会平移（其实是两个正交基的分解）。缩得窄，对应高频；伸得宽，对应低频。然后这个基函数不断和信号做相乘。某一个尺度（宽窄）下乘出来的结果，就可以理解成信号所包含的当前尺度对应频率成分有多少。于是，基函数会在某些尺度下，与信号相乘得到一个很大的值，因为此时二者有一种重合关系。那么我们就知道信号包含该频率的成分的多少。如前边所说，小波做的改变就在于，将无限长的三角函数基换成了有限长的会衰减的小波基。

三、小波变换和傅里叶变换的比较

1. 傅里叶变换的实质是把能量有限的信号$f(t)$分解到以$e^{j\omega t}$为正交基的空间上去；小波变换的实质是把能量有限的信号$f(t)$分解到以$W_{-j}(j=1,2,\cdots ,J)$和$V_{-j}$所构成的空间上去。
2. 傅里叶变换用到的基本函数只有$sin(\omega t),cos(\omega t),e^{j\omega t}$，具有唯一性；小波分析所用到的小波函数则不是唯一的，同一个工程问题用不同的小波函数进行分析有时结果相差甚远。小波函数的选用是小波分析应用到实际中的一个难点问题，也是分析研究的一个热点问题，目前往往是通过经验和不断的试验，将不同的分析结果进行对照分析来选择小波函数。
3. 在频域中，傅里叶变换具有较好的局部化能力，特别是对于那些频率成分比较简单的确定性信号，傅里叶变换很容易把信号表示成各频率成分的叠加和的形式，但在时域中，傅里叶变换没有局部化能力，无法从信号$f(t)$的傅里叶变换$F(\omega )$中看出$f(t)$在任一时间点附近的特性。
4. 在小波分析中，尺度$\alpha$越大相当于傅里叶变换中的$\omega$值越小。
5. 在短时傅里叶变换中，变换系数$G_f(\omega ,\tau )$主要依赖于信号在时间窗内的情况，一旦时间窗函数确定，则分辨率也就固定了。而在小波变换中，变换系数$W{T}_x(\alpha ,\tau )$虽然也是依赖于信号在时间窗内的情况，但时间尺度是随尺度$\alpha$的变化而变化的，所以小波变换具有时间局部分析的能力。
6. 若用信号通过滤波器来解释，小波变换与短时傅里叶变换不同之处在于：对短时傅里叶变换来说，带通滤波器的带宽$\Delta \omega$与中心频率$\omega$无关；相反，小波变换带通滤波器的带宽$\Delta \omega$则正比于中心频率$\omega$，即
   $$Q=\frac{\Delta \omega }{\omega }=CC为常数$$
   也就是说，这组滤波器具有品质因数恒定，即相对带宽恒定的特点，称之为等$Q$结构。

四、小波变换的基础知识（Wavelet Transform,WT）

小波变换具有多分辨特性，也叫多尺度特性，可以由粗到精地逐步观察信号，也可以看成是用一组带通滤波器对信号作滤波。通过适当地选择尺度因子和平移因子，可得到一个伸缩窗，只要适当选择小波基，就可以使小波变换在时域和频域都具有表征信号局部特征的能力，有利于检测信号的瞬态和奇异点。

1. 连续小波变换（Continuous Wavelet Transform,CWT）

1.1 定义

将任意$L^2(R)$空间（表示为可测量的、模平方可积分的一维函数空间）中的函数$z(t)$在小波基$\phi (t)$下进行展开，称作函数$z(t)$的连续小波变换（CWT），其表达式为
$$W{T}_z(a,b)=<z(t),{\phi }_{a,b}(t)>=\frac{1}{\sqrt{a}}{\int }_{R}z(t)\phi (\frac{t-b}{a})dt$$
式中，$a>0$，称为尺度因子（在某种意义下就是频率的概念），其作用是对小波基函数作伸缩，$b$反映位移，其值可正可负，$a$和$b$都是连续的变量。
 
在不同尺度下小波的持续时间随值的增大而增宽，幅度则与$\sqrt{a}$反比减少，但波的形状保持不变。

由连续小波变换的定义可知，小波变换同傅里叶变换一样，都是一种积分变换，称$W{T}_z(a,b)$为小波变换系数。由于小波基不同于傅里叶基，因此小波变换与傅里叶变换有很大的不同，其中最重要的是，小波基具有尺度$a$和平移$b$两个参数。

根据时频分析的要求，构造的小波基函数$\phi (t)$应该满足以下条件：

1. 本身是紧支撑的，即只有小的局部非零定义域，在窗口之外函数为零；(具有有限的持续时间和突变的频率和振幅)
2. 本身是振荡的，具有波的性质，并且完全不含直流趋势成分，（在有限时间范围内平均值为0）即
   $$\Psi (0)={\int }_{-\infty }^{\infty }\phi (t)dt=0$$
   式中$\Psi (f)$是函数$\phi (t)$的傅里叶变换。
3. 包含尺度参数$a(a>0)$和平移参数$b$
   $a$增大，则时窗伸展，频宽收缩，带宽变窄，中心频率降低，而频率分辨率增高；$a$减小则带宽增加，中心频率升高，时间分辨率增高而频率分辨率降低。这恰恰符合实际问题中高频信号持续时间短、低频信号持续时间长的自然规律。因此，同固定时窗的短时傅里叶变换相比，小波变换在时频分析领域具有不可比拟的优点。

如果采用的小波满足下式，则其反变换存在
$$C_{\phi }={\int }_R\frac{|\Psi (\omega ){|}^2}{|\omega |}d\omega <\infty$$
且连续小波变换的反变换表达式为
$$x(t)=\frac{1}{{\int }_0^{\infty }\frac{|\Psi (a\omega )|}{a}d\omega }{\int }_0^{\infty }\frac{da}{a^2}{\int }_{-\infty }^{\infty }W{T}_z(a,b)\frac{1}{\sqrt{a}}\phi (\frac{t-b}{a})db$$

1.2 基本性质

连续小波变换是一种线性变换，具有以下性质：

性质1：叠加性

假设$x(t),y(t)\in L^2(R)$空间，${\lambda }_1,{\lambda }_2$为任意常数，记$x(t)$的连续小波变换为$W{T}_x(a,b)$，$y(t)$的连续小波变换为$W{T}_y(a,b)$，则$z(t)={\lambda }_{1}x(t)+{\lambda }_{2}y(t)$的连续小波变换为
$$W{T}_z(a,b)={\lambda }_{1}W{T}_x(a,b)+{\lambda }_{2}W{T}_y(a,b)$$

性质2：时移不变性

记$x(t)$的连续小波变换为$W{T}_x(a,b)$，则$x^{\prime }(t)=x(t-t_0)$的连续小波变换为
$$W{T}_{x^{\prime }}(a,b)=W{T}_x(a,b-t_0)$$
表明：延时后的信号的小波系数可将原信号的小波系数在b轴上进行相应时移即可。

性质3：尺度变换

记$x(t)$的连续小波变换为$W{T}_x(a,b)$，则$x^{\prime }(t)=x(\frac{t}{\rho })$（$\rho >0$为常数）的连续小波变换为
$$W{T}_{x^{\prime }}(a,b)=\sqrt{\rho }W{T}_x(\frac{a}{\rho },\frac{b}{\rho })$$
表明：当信号在时域作某一倍数伸缩时，其小波变换在a,b两轴上也作同一倍数伸缩，且形状不变。

性质4：内积定理

假设$x(t),y(t)\in L^2(R)$空间，记$x(t)$的连续小波变换为$W{T}_x(a,b)$，$y(t)$的连续小波变换为$W{T}_y(a,b)$，则有
$$<W{T}_x(a,b),W{T}_y(a,b)>=C_{\phi }<x(t),y(t)>$$
式中$C_{\phi }$是$\phi (t)$的Cohen类时频分布，<$\cdot$>表示内积运算。

1.3 常用的连续小波基函数

1.3.1 Morlet小波

时域表达式：$\phi (t)=e^{-t^2/2}{e}^{j{w}_{0}t},w_0\ge 5$
频域表达式：$\Psi (w)=\sqrt{2\pi }{e}^{-\frac{(w-w_0{)}^2}{2}}$

Morlet小波是一种复数小波，其在时、频域都具有较好的局部性，常用于复数信号的分解及时频分析中。

%%********** Morlet 小波函数 *********%%
% 这里仅画出实部波形，且w0=5
lb = -8; ub = 8; n = 200;
[psi,xval] = morlet(lb,ub,n);
PSI=fftshift(abs(fft(psi)));
figure;
subplot(2,1,1); plot(xval,psi); title('Morlet Wavelet时域波形');
subplot(2,1,2); plot(PSI);      title('Morlet Wavelet频域波形');

1.3.2 Mexican hat小波

时域表达式：$\phi (t)=(1-t^2)e^{-t^2/2}$
频域表达式：$\Psi (w)=\sqrt{2\pi }{w}^2{e}^{-\frac{w^2}{2}}$

Mexican hat小波实际上是高斯函数的二阶导数，在时、频域都具有较好的局部性。

%%********** Mexican hat 小波函数 *********%%
lb = -8; ub = 8; n = 200;
[psi,xval] = mexihat(lb,ub,n);
PSI=fftshift(abs(fft(psi)));
figure;
subplot(2,1,1); plot(xval,psi);  title('Mexican Hat Wavelet时域波形');
subplot(2,1,2); plot(PSI);       title('Mexican Hat Wavelet频域波形');

2. 离散小波变换（Discrete Wavelet Transform,DWT）

根据连续小波变换的定义可知，在连续变化的尺度$a$和时间$b$值下，小波基函数具有很大的相关性，因此信号的连续小波变换系数的信息量是冗余的。

减小小波变换系数冗余度的方法是将小波基函数${\phi }_{a,b}(t)$的$a,b$限定在一些离散点上取值，常用的离散化方法是将尺度按幂级数进行离散化，即取$a_m=a_0^m$（$m$为整数，$a_0\ne 1$,一般取$a_0=2$），对$b$进行均匀离散取值，以覆盖整个时间轴（为了不丢失信息，要求满足采样定理），这样小波基函数变为
$${\phi }_{m,n}(t)=2^{-m/2}\phi (2^{-m}t-n)$$
因此任意函数$z(t)$的离散小波变换为
$$W{T}_z(m,n)=<z(t),{\phi }_{m,n}(t)>=2^{-m/2}{\int }_{R}z(t)\phi (2^{-m}t-n)dt$$

在离散小波中，还有一类特殊情况，即仅在尺度上进行了二进制离散，而位移仍连续变化，称之为二进小波(Dyadic Wavelet)，其表达式为
$${\phi }_{2^k,b}(t)=2^{-k/2}\phi (\frac{t-b}{2^k})$$
离散小波变换不再具有伸缩和时移共变性质；二进小波介于连续小波和离散小波之间，仍具有连续小波变换的时移共变性质，正因为如此，它在奇异性检测、图像处理方面十分有用。

3. 正交小波变换

3.1 多尺度分析

从小波函数的条件来看，小波基不一定是正交基，但是在实际应用中希望找到正交小波基。构造正交小波基的重要方法称之为多尺度分析。

---

根据信号空间的概念，由尺度函数$\varphi (t)$可以定义其对应的小波函数$\psi (t)$，再由小波函数经过尺度伸缩与平移得到小波信号${\psi }_{m,n}(t)$，即
$$\begin{gathered} \varphi (t)⟹\psi (t)⟹{\psi }_{m,n}(t) \\ {\psi }_{m,n}(t)=2^{m/2}\psi (2^{m}t-n)m,n\in Z \end{gathered}$$

小波信号${\psi }_{m,n}(t)$设计为尺度信号${\varphi }_{m,n}(t)$的正交信号，即存在
$$<{\varphi }_{m,n}(t),{\psi }_{m,n}(t)>=\int {\varphi }_{m,n}(t)\cdot {\psi }_{m,n}(t)dt=0$$

---

多尺度分析是通过函数空间术语来严格定义的，假设$L^2(R)$空间内的子空间序列$V_m(m\in Z)$满足以下条件：

1. （嵌套性）$\cdots V_2\subset V_1\subset \cdots \subset V_{-2}\subset \cdots ;$
2. （逼近性） ⋂ m ∈ Z V m = { 0 } ， ⋃ m ∈ Z V m = L 2 ( R ) ; \bigcap_{m∈Z} V_m=\{0\}，\bigcup_{m∈Z} V_m=L^2(R); ⋂m∈Z​Vm​={0}，⋃m∈Z​Vm​=L2(R);
3. （二进制伸缩性）$f(t)\in V_m\iff f(2t)\in V_{m-1}$
4. V m = s p a n { ϕ m , n ( t ) , n ∈ Z } V_{\color{Red}m}=span\{\phi_{{\color{Red}m},n}(t),n∈Z\} Vm​=span{ϕm,n​(t),n∈Z}，即任一级子空间可由相应尺度的同一函数通过平移张成；
5. $V_{m-1}=V_m\oplus W_m$，即任一级子空间可由下一级子空间以及它的正交补空间相加而成，序列$W_m$相互之间无重叠，是正交系。

这样，我们称子空间序列$V_m,W_m$为函数空间$L^2(R)$上的一个多尺度分析，其中${\varphi }_{m,n}$称为尺度函数，m,n分别是尺度和平移参数。

---

由以上定义可得
$$V_0=V_m\oplus \sum _{i=1}^m{W}_i\text{\hspace{1em}(3.1.1)}$$
特别地，当尺度趋于无穷时，上式变为
$$L^2(R)=\sum _{i=-\infty }^{\infty }{W}_i\text{\hspace{1em}(3.1.2)}$$
存在与${\varphi }_{m,n}$相应的函数${\psi }_{m,n}$通过平移能够张成$W_m$。
W m = s p a n { ψ m , n ( t ) } (3.1.3) W_m=span\{\psi_{m,n}(t)\} \tag{3.1.3} Wm​=span{ψm,n​(t)}(3.1.3)
${\psi }_{m,n}$就是所求的小波基，它在伸缩和平移变换下都是正交的。

3.2 Mallat算法

多尺度分析在信号分析中的应用可以用式(3.1)来表达

  对于任意函数$f(t)\in V_0$，可以在下一级尺度空间$V_1$和小波空间$W_1$上进行分解，如下所示：
$$f(t)=p_{1}f(t)+q_{1}f(t)$$
式中$p_{1}f(t)={\sum }_k{C}_{1k}{\varphi }_{1,k}$，$q_{1}f(t)={\sum }_k{D}_{1k}{\psi }_{1,k}$。
 
  $p_{1}f(t)$是逼近部分，$q_{1}f(t)$是细节部分。然后将逼近部分进一步分解，如此重复就可得到任意尺度上的逼近部分和细节部分，迭代公式为
$$p_{m-1}f(t)=p_{m}f(t)+q_{m}f(t)=\sum _k{C}_{mk}{\varphi }_{m,k}+D_{mk}{\psi }_{m,k}$$
式中$C_m=H{C}_{m-1}$，$D_m=G{C}_{m-1}$，H是低通滤波器。每一次分解，$p_{m}f(t)$的采样都比原来稀疏两倍，采样率越来越粗，波形越来越光滑；G是H的镜像高通滤波器，带宽每次也以两倍缩减。
 
  经过m次分解，得到
$$f(t)=p_{m}f(t)+\sum _{j=1}^m{q}_{j}f(t)$$
式中$p_{m}f(t)$给出的是函数$f(t)$的低频全局信息，第二次是逐次分解中分离的从$V_0$到$V_{m-1}$各个尺度上$f(t)$的相应局部细节信息。上述信号分解过程称之为Mallat算法。

其过程如图所示：

3.3 常用的正交小波基函数

常用的正交小波基有Haar小波、Meyer小波、Daubechies小波。

1. Haar小波

%%********** Haar 小波函数 **********%%
[phi,psi,xval]=wavefun('haar',20);
figure;
subplot(2,1,1); plot(xval,psi); title('Harr Wavelet Function');
subplot(2,1,2); plot(xval,phi); title('Harr Scaling Function');

2. Meyer小波
   

3. Daubechies 小波
   

4. 小波包（Wavelet Packet）

4.1 原理

在多尺度分析中，每一步分解都只是对尺度子空间$V_i$进行，对子空间$W_i$则不再分解，因此高频跨度较宽。实际中，在很多问题中，我们只是对某些特定时间和频域段(点)的信号感兴趣，只要提取这些位置的信息即可。为此，我们希望在感兴趣的频率点上最大可能地提高频率分辨率，在感兴趣的时间点上最大程度地提高时间分辨率，此时正交小波不再满足这种要求。

这需要对子空间$W_i$也一样进行分解，以提高频率分辨率，称为小波包分解。

4.2 分解过程

其过程如图所示：

对小波子空间$W_i$可以像尺度子空间$V_i$一样利用镜像滤波器进行正交分解，对其频带进行分割。将其子空间用符号$U$统一表示，子空间分解的公式为
$$U_i^m=U_{i-1}^{2m}\oplus U_{i-1}^{2m+1}\text{\hspace{1em}(4.2.1)}$$
式中下标表示尺度层次，上标区别同一层上不同的子空间。
 
按照式(4.2.2)，小波子空间$W_i=U_i^1$可以如下分解：
$$\begin{array}{rl}{W}_i& =U_{i-1}^2\oplus U_{i-1}^3\\ W_i& =U_{i-2}^4\oplus U_{i-2}^5\oplus U_{i-2}^6\oplus U_{i-2}^7\\ W_i& =U_{i-k}^{2^k}\oplus U_{i-k}^{2^k+1}\oplus \cdots \oplus U_{i-k}^{2^{k+1}-1}\end{array}\text{\hspace{1em}(4.2.2)}$$
每层小波包将原频带一分为二，$k$层小波包可将原频带分割为$2^k$个子频带，从而实现频带细分，提高了频率分辨率。对整个小波包来说，它是一个按二进制组织的包含从宽到窄各个频带的带通滤波器组，各种应用都能够从中找到符合需求的最优组合。可以通过构造代价函数来计算如何实现最优（利益最大化或代价最小化）。

 
 

参考资料：Matlab 时频分析技术及其应用/葛哲学，陈仲生编著.–北京:人民邮电出版社，2006.1