折腾来折腾去

pikipity的blog

寻找一种有效估算瞬时频率的方法(EMD 求法讨论)

我试图寻找一种可以用比较少的数据点就可以较为精确的估算某个时刻瞬时频率的方法,当然了,根据“薛定尔的胖次”(对于虐猫狂人薛定谔不熟悉的可以去看这里的这本少儿读本,写的很好)在信号分析中的拓展,如果要精确计算瞬时频率就需要一个在时域上无限长的信号,当然这是不可能的,所以,如果能够找到一种计算方法,在时域上的精度和频域上的精度问题上去的平衡就可以了。

在经过一番大搜索之后,发现现在最为流行的求瞬时频率的方法就是 HHT(希尔伯特-黄转换,也就是 EMD,经验模态分析,论文在这里,洋洋洒洒90多页,大家慢慢看)。

什么是瞬时频率

第一个问题我就说不清楚,一般来说,把瞬时频率定义为

$$f(t)=\frac{1}{2\pi}\frac{d\theta(t)}{dt}$$

这个定义是 Cohen 提出的,但是这个定义仅限于 “monocomponent function",但是什么是 "monocomponent function” 呢?他没说。。。(坑爹啊 (╯‵□′)╯︵┻━┻)但是 Gabor, Bedrosian 和 Boashash 提出:

for any function to have a meaningful instantaneous frequency, the real part of its Fourier transform has to have only positive frequency.

黄锷院士(就是这个 HHT 的创始人啦)进一步验证得到,对于平均值为零的局部对称信号而言,上述定义的瞬时频率才具有物理意义。

据说,“澳门大学”(请注意,这里不是做宣传)的钱涛教授在提出 AFD 的时候定义了更为准确的瞬时频率,但是博主愚钝,到现在也没看懂他的论文,甚至都不知道从何看起,所以更不用提他提出的方法和定义了。

寻找新方法的原因

在寻找到 HHT 之前,手头已经有的方法大体可以分为两种 – Fourier Transform 和 Wavelet Transform,当然还包括一些 Fourier Transform 的变种,例如 Wigner–Ville Distribution。但是这两种方法都存在着各种各样的缺点。

Fourier Transform

我们都知道,Fourier Transform 是试图将一个信号分解为多个幅值和频率固定的三角函数线,对于其对于信号的要求是稳定、收敛,但是实际中,尤其是在脑电中,我们无法保证自己截取的信号就是符合要求的。而且 Fourier Transform 所求的是一种平均频率分布,这意味其会带来两大缺点:

  1. 能量分散到了整个频谱上,难以辨清具体频率所在。
  2. 产生了无意义的负频。

当然,以上两点都可以接受的话,其最致命的缺点就是 丧失了时间信息。尽管 STFT (Short Time Fourier Transform) 利用窗截取某个时间段的信号,可以在一定程度上获得瞬间频率,但是由于窗长度一旦选取就固定不变,导致时间精度固定,从而导致频率精度也固定,这样对频率的突变就及其不敏感,而检测突变则是信号分析中非常重要的,所以 Fourier Transform 难以分析瞬时频率。

Wavelet Transform

根据下面 Wavelet Transform 的公式,我们可以清楚的看出,为了改善 Fourier Transform 会丧失时间信息的缺点,Wavelet Transform 添加两个参数 $b$(时间刻度)和 $a$(频率刻度,可以转化为对应的频率),这样就可以同时获得时间和频率两方面的信息。同时,不再使用三角函数作为基函数,而是使用小波 $\phi$ 作为基函数,可以解决 Fourier Transform 用稳定信号分析不稳定信号的问题。

$$W(a,b;W,\phi)=\frac{1}{\sqrt{a}}\int_{-\infty}^{\infty}W(t)\phi^*(\frac{t-b}{a})dt$$

但其对 Fourier Transform 的改善还是存在缺点:

  1. 尽管随着 $a$ 的变化,频率可以在多分辨度下进行解析,但是在某一个分辨度下,其分辨精度还是固定不变的。
  2. 小波一旦选定就无法改变,没有办法保证此小波在整个波段上都合适。
  3. 小波选择问题难以解决,一般只不过是通过经验来进行选择。

HHT

就算上面两种方法可以完美解析实际信号,但是也很难从他们的算法中找到计算瞬时频率的办法,那么根据提到的定义,要找瞬时频率,就要先找瞬时相位,于是慢慢征程刚开始啊。。。

HT

Gabor 提出一种定义信号相位的方法,如果能够把信号表示成为其对应的复数形式(如下),那么相位自然也就找到了。

$$z(t)=s(t)+j\hat{s}(t)$$

上式中的 $\hat{s}(t)$ 可以通过 Hilbert Transform (HT) 来得到。下式就是 HT 的转化方法。

$$\hat{s}(t)=\frac{1}{\pi}\int_{-\infty}^{\infty}\frac{s(t)}{t-\tau}d\tau$$

为了使信号可以进行 HT,并且具有有意义的瞬时频率,必须先对信号进行处理,于是黄锷院士提出 HHT。

IMF

所谓对信号进行处理,就是将信号拆分为多个可以进行 HT 来寻找瞬时频率的信号 – Intrinsic Mode Function(IMF)。IMF 满足一下两点:

  1. 在整个信号上,局部极值(包括局部极大值和局部极小值)的个数和零点的个数必须等于或最多只能相差一个。
  2. 在任意一点,局部极大值所形成的上包络线和局部极小值所形成的下包络线的平均值等于0

这个由原始信号到 IMF 的分解流程称为 Empirical Mode Decomposition (EMD)。

EMD

EMD 的流程非常明确,在上面提供的论文中都有写,待分析的信号为 $s(t)$:

  1. 找出 $s(t)$ 的所有局部极大值和局部极小值,接着利用三次样条插值得到由局部极大值连接形成的上包络线和局部极小值连接形成的下包络线。
  2. 求出上下包络线的平均值 $$m_1(t)$$
  3. 利用下式得到第一分量:

    $$h_1(t)=s(t)-m_1(t)$$

  4. 检查 $$h_1(t)$$ 是否是 IMF,如果不是则将 $$h_1(t)$$ 作为原信号返回第一步,直到 $$h_k(t)$$ 是 IMF,就得到了一个 IMF 分量 $$c_1(t)$$,也就是

    $$c_1(t)=h_k(t)$$

  5. 求剩余分量:

    $$r_1(t)=s(t)-c_1(t)$$

  6. 将 $r_1(t)$ 当做新资料,返回步骤1,直到 $r_n(t)$ 为单调函数完成分解。

完成分解后,就将原始信号分解为了 $n$ 个 IMF 和一个趋势函数,也就是

$$s(t)=\sum_{k=1}^nc_k(t)+r_n(t)$$

EMD 的 MATLAB 实现

尽管上面的算法流程第一次看起来清晰明了,但是要注意,我们分析的是离散的数据点,而不是连续的信号,于是这就带来了很多问题,以至于上面流程中很多名词都没有严格的定义,需要我们自己来在编写程序的时候来定义,根据定义的不同,程序最后运行的结果也就会不同。一边考虑我们要分析的信号是离散信号,一边在看一遍上面的流程就会很容易发现下面几个问题:

  1. 如何定义离散信号的局部极大值与局部极小值?因为我们的离散信号是通过实际信号采样得来的,那么就会出现离散信号的局部极值与原始实际信号的局部极值不一样的情况,甚至出现多个连续一样的局部极值,处理起来非常麻烦;而且在考虑局部极值的时候还会出现边缘问题,离散信号最开始与最末尾的两个点如何处理?绝对不能舍弃,因为在上面的流程中计算局部极值是会在循环中多次出现,如果每次舍弃两个点,会丢失大量数据,而处理不当,在插值的时候边缘会产生出奇怪的震荡导致 IMF 计算不准确影响瞬时频率估计。
  2. 三次样条插值如何用 MATLAB 实现?这应该是最容易解决的问题了,因为 MATLAB 中有现成的函数,随后介绍。
  3. 如何判断 $h_k(t)$ 是 IMF?这个问题应该是第二简单的问题了,随后介绍。

以上三个问题的解决方式不同,分解出的 IMF 也就不同,估算得到的瞬时频率的准确度也就不同,我找了多个程序,并自己编写了一个程序,但是效果都不能说是很好,其实还是蛮失望的,EMD 这个方法出现也有很多年了,结果连一个现成的比较完美的 MATLAB 程序都没有。这里使用国际比较认可的由一个法国人 G. Rilling 编写的程序(可以在这里下载到)和我在 Github 上找到的一个效果速度也还可以的 MATLAB 程序(可以在这里下载到)综合说明一下他们对上面三个问题的解决方法和我的一些见解。(我写的程序在这里,但是仅仅当做错误教程看一看就可以了,方法极度不可取,速度慢不说,分解出来的本应完全正交的 IMF 的频域还会出现重叠,边缘也出现了震荡,导致出现了很多低频杂波。)

由于上面三个问题难易程度不一,我会先讨论第二个问题,然后是第三个,最后再讨论第一个。

如何用 MATLAB 实现三次样条差值

在 MATLAB 中有三种方法可以实现三次样条差值:

  1. yy=spline(x,y,xx):我从 Github 上找到的那个程序使用的是这种插值方式。其中,x是断点的横坐标,y是断点的纵坐标,xx是插值之后重新构成的函数采样点的横坐标,yy是插值之后重新构成的函数采样点的纵坐标。当然,你也可以使用 pp=spline(x,y),这样得到的pp是插值之后的到的函数的 structure,然后可以通过 yy=ppval(pp,xx) 来得到重新构成的函数的采样点。插值效果如图

    使用```spline```进行插值

  2. pp=csape(x,y,CONDS):在我的程序中,这个函数的特点是你可以通过改变CONDS来使用不同的计算三系样条插值的方法,这些方法唯一的不同就是使用的 boundary condition 不一样,根据你选择的方法的不同,你可能需要多提供边界的一次导数或是二次导数作为额外的 boundary condition 来进行计算。其他的参数和第一个方法一样。效果如图

    使用```csape```进行插值

  3. yy=interp1(x,y,xx,'spline'):这个方法的本质还是调用第一个方法中的函数进行运算,所以其实可以合并到方法一中。G. Rilling 使用的这种方法。

将上面两种方法进行比较,得到下面这幅图

```spline```和```csape```插值效果比较

我们可以发现,只有边界处出现细微差别,其他完全一样,如果调整csape的插值方法,可以得到完全一样的结果,所以使用哪种的效果是完全一样的。

要注意我们用 MATLAB 做三次样条差值是为了寻找上下包络线,那么什么是包络线呢?维基上的解释是这样的

在几何学,某个曲线族的包络线(Envelope),是跟该曲线族的每条线都有至少一点相切的一条曲线。(曲线族即一些曲线的无穷集,它们有一些特定的关系。)

不知道你看懂没有,反正我是没有看懂(囧)。但是我至少看懂了一点,包络线不能和曲线族相交们仅仅是相切,但是用局部极值点三次样条差值得到的包络线有的时候会出现与原信号相切的情况,这里提供了一种用切点求包络线的方法,但是据博主说,计算用时会大大增加,EMD 的计算用时本来就很长,如果再让求包络线拖慢时间我认为并不值得,况且出现相切的状况只在局部斜率变化缓慢时出现,脑电上很少出现这种状况,所以我认为用极值点插值已经足够了。

然后是是否使用更高次的插值方式,很多论文中提到使用更高次的插值方式使得到的包络线更加平滑,但是我认为完全没有必要,三次样条插值已经足以满足平滑上的需求,而且我们将来使用的也只不过是插值的到的取样点,过好的平滑反而有可能会被取样点的疏松掩盖掉,三次样条差值足以。

(未完待续)



Comments