信号时频分析方法

  传统的信号处理主要基于信号时域和频域来进行，信号的时域信息表征了信号幅度随时间的变化关系，信号的频域信息则揭示了信号频率和能量之间的关系。电磁信号通常都是典型的非平稳信号，只进行时域角度的分析则无法获取频域的信息，而只分析了全局频谱的傅里叶变换不仅无法对非平稳信号进行准确有效的描述，还丢失了信号的时域信息。从 20 世纪 40 年代起，科研人员就开始将数学领域的理论应用到信号处理中，运用时频分析方法对电磁信号进行分析和处理，可以同时表征信号的时域和频域信息，时频分析已成为现代数字信号处理中一种常见的信号特征提取方法。

一、短时傅里叶变换

  对于平稳信号而言，傅里叶变换（Fourier Transform, FT）作为一种分析和处理的重要方法，被广泛应用于数学、信号处理、热力学分析和工程控制论等领域。但对于自然界中广泛存在的非平稳信号，傅里叶变换只能得到其全局的频谱信息而丢失在局部时间内的频率变化信息，无法准确地表征非平稳信号的特性。连续时间信号$s(t)$的傅里叶变换定义如下：
$$S(f)=\mathcal{F}[s(t)]={\int }_{-\infty }^{\infty }s(t)e^{-j2\pi ft}\mathrm{d}t$$
  1946 年 Gabor 为了解决由于傅里叶变换对信号频谱表征的全局性而无法体现出电磁信号频率随时间变化规律的问题，提出了短时傅里叶变换(Short Time Fourier Transform, STFT)的概念，其基本思想是在信号傅里叶变换前乘上一个固定长度的矩形窗函数 ，从而实现信号在局部时间里的伪平稳，然后通过矩形窗在时间轴上的移动对信号逐段进行傅里叶变换，从而得到信号的时变特性。
  连续时间信号$s(t)$的 STFT 定义如下：
$$STFT(t,f)={\int }_{-\infty }^{\infty }s(\tau )h(\tau -t)e^{-j2\pi f\tau }\mathrm{d}\tau$$
  离散时间信号$s[n]$的 STFT 定义如下：
$$STFT[m,k]=\sum _{n}s[n]\omega [n-k]e^{-j\frac{2\pi }{N}mn}m=0,1,...,N-1$$
  时间窗函数的宽度决定了信号 STFT 的时间分辨率，窗函数宽度$T_p$越小，时间分辨率越高，二者关系表达式为：
$$T_p=NT=\frac{N}{f_{sample}}$$
  频率分辨率是度量信号 STFT 效果的重要参数，是指分辨信号频谱中最小相邻谱峰间隔的能力，$\Delta f_c$越小频率分辨率越高，与时间分辨率的关系如下：
$$\Delta f_c=\frac{1}{T_p}$$
  短时傅里叶变换计算方法相对简单，是一种联合的信号线性时频分析方法。利用固定长度的矩形窗函数对信号进行时频分析，可以避免高次型非平稳信号分析方法中出现的交叉项干扰问题，适用于信号的多分量分析。

二、小波变换

  通过在时间上将非平稳信号分割成许多小的时间间隔而实现非平稳信号的局部伪平稳，然后在对每个伪平稳局部信号进行傅里叶变换的基础上乘以矩形窗函数得到信号的短时傅里叶变换。虽然短时傅里叶变换克服了傅里叶变换的频谱分布全局性的特点，但是由于其选取的矩形窗函数的宽度是固定的，无法同时兼顾较高的时间分辨率和频率分辨率，较长时间宽度通过牺牲时间分辨率得到理想的频率分辨率；较短时间宽度则通过牺牲频率分辨率来得到理想的时间分辨率。
  连续小波变换（Continue Wavelet Transform, CWT）继承和发展了短时傅里叶对非平稳信号处理采取的局部分割的思想，同时克服短时傅里叶窗函数固定宽度的缺点，为信号进行时频变换提供了一个随频率改变的时间窗口。信号连续小波变换的定义如下：
$$CWT(x,y)=\frac{1}{\sqrt{x}}{\int }_{-\infty }^{\infty }s\left(t\right)\phi \left(\frac{t-y}{b}\right)dt$$
  式中的x 和y分别表示尺度和平移量。尺度x表征小波基函数的伸缩量，与频率成反比，可以在时频图中表征频率；平移量 y 表示小波函数的平移量，对应于时频图中时间的长短。
  小波分析在信号处理、图像识别、语音合成和控制论等领域都取得了广泛的应用和巨大的成果。在通信信号处理中，小波变换主要用于分离噪声、提取弱信号、时频分析、滤波和信号识别等。

三、魏格纳-维利分布

  1932年美国物理学家Eugene Wigner提出了Wigner-Ville分布（WVD）的概念，Vile于1948年将WVD应用到信号处理领域中。连续时间信号 的Wigner-Ville分布可以通过对信号自身的瞬时自相关函数进行傅里叶变换得到，信号的瞬时自相关表达式如下：
$$r_s(t,\tau )=s(t+\frac{\tau }{2})s^{\ast }(t-\frac{\tau }{2})$$
  其WVD的定义式如下：
$$WVD\left(t,f\right)={\int }_{-\infty }^{\infty }s\left(t+\frac{\tau }{2}\right)s^{\ast }\left(t-\frac{\tau }{2}\right)e^{-j2\pi f\tau }\mathrm{d}\tau$$
  信号的魏格纳-维利分布具有较高的时间分辨率和频率分辨率，并且有良好的时频聚集性，克服了 STFT 的缺点，适合分析非稳态的随机信号。

四、希尔伯特-黄变换

  1998 年，Norden E. Huang 等人提出了经验模态分解（Empirical Mode Decomposition, EMD）方法，并引入了希尔伯特谱的概念和分析方法，这一方法后来被命名为 Hilbert-Huang Transform，简称 HHT，即希尔伯特-黄变换。HHT 处理非平稳信号的基本思路是：首先利用经验模态分解方法将给定的信号分解为若干固有模态函数（Intrinsic Mode Function, IMF），然后对每一个固有模态函数进行希尔伯特变换得到相应的希尔伯特谱，最后汇总所有固有模态函数的希尔伯特谱就得到原始信号的希尔伯特谱。
  希尔伯特-黄变换主要包括两个过程，经验模态分解和希尔伯特变换，经验模态分解方法过程如下：
  （1）分别提取待分析信号$s(t)$的所有局部极大值和极小值点，再分别拟合极大值点、极小值点形成极大值和极小值包络线，得到极值包络线均值。
$$m_1=\frac{e_{max}(t)+e_{min}(t)}{2}$$
  （2）信号$s(t)$减去极值包络线的均值$m_1$得到新信号$y_1(t)$
$$y_1(t)=s(t)-m_1$$
得到的新信号$y_1(t)$就是固有模态函数的分量形式，然后判断该分量形式是否满足条件。如果$y_1(t)$不满足条件，则将$y_1(t)$作为原始数据重复步骤（1）和（2），不断筛选一直到$y_1(t)$满足固有模态函数的分量条件，记满足条件的$y_1(t)=c_1(t)$，则$c_1(t)$是信号$s(t)$的第一个固有模态函数分量。

  （3）将剩余信号$r_1(t)=s(t)-c_1(t)$作为原始信号，同上面原理分解 $j$次得到 $j$ 阶固有模态函数，余量为$r_n(t)$。如果得到的余量$r_n(t)$是一个单调信号或$m_n$的值足够小时，经验模态分解过程结束。综合上述结果，将各阶固有模态函数及余项$r_n(t)$重构获得如下的原始信号：
$$s(t)=\sum _{i=1}^n{c}_i(t)-c_n(t)$$
  对原信号进行经验模态分解后，再对分解得到的固有模态函数分量进行希尔伯特变换就可以得到每个固有模态函数的瞬时频率，然后将所有分量的瞬时频谱综合即可获得原始信号的希尔伯特谱。
  各阶固有模态函数分量进行希尔伯特变换：
$$H[c(t)]=\frac{1}{\pi }PV{\int }_{-\infty }^{\infty }\frac{s(\tau )}{t-\tau }\mathrm{d}\tau$$
式中 PV 表示柯西主值，根据瞬时频率公式的定义：
$$\omega (t)=\frac{d\left[\arctan \left(\frac{\hat{s}(t)}{s(t)}\right)\right]}{dt}$$
可以将原始信号表示为：
$$s(t)=\text{Re}\left(\sum _{i=1}^n{a}_i\left(t\right)e^{j{\varphi }_i\left(t\right)}\right)=\text{Re}\left(\sum _{i=1}^n{a}_i\left(t\right)e^{j\int {\omega }_i\left(t\right)dt}\right)$$
展开后就是希尔伯特谱，可以写成：
$$H\left(\omega ,t\right)=\mathrm{Re}\left(\sum _{i=1}^n{a}_i\left(t\right)e^{j\int {\omega }_i\left(t\right)dt}\right)$$
  希尔伯特谱是一种完整的时间、频率和能量的联合谱，其中可以同时看到不同时间段的频率变化情况和能量随着时间与频率的变化而变化的情况。