基于ST,的EMD,算法在FPGA,上的设计与实现
康世勋，孔德杰，冯进良，马晨阳
（长春理工大学 光电工程学院，长春 130022）
数字信号处理是许多领域理论和应用的重要组成部分，一些传统的信号处理方法如傅里叶变换是处理线性系统和平稳信号的有效方法，但是现实世界中的大多数系统都是非线性，且信号是非平稳的[1]。
经验模态分解（EMD）算法是希尔伯特-黄变换（Hilbert-Huang Transform，HHT）的主要组成部分。其在结构动力分析、医学检测、地震检测以及设备故障诊断[2-5]等工程领域有着重要的意义。其筛选过程需要大量的计算，利用软件实现EMD 的高速计算是困难的。但是在许多应用中，特别是与安全相关的应用，如轴承故障检测和转子条断裂检测，高速分析和反应是必要的。基于上述问题，提出一种提高EMD 计算速度的解决方案。
为了提高EMD 的计算速度，现场可编程逻辑门阵列（FPGA）凭借其低功耗、高能耗比、并行化处理和可重配置的优点得到了研究人员的重视，并且其并行化与流水线处理方式能够很好地处理EMD 算法的高计算量。例如周浩等人[6]提出了基于FPGA 的EMD 算法的实现，提高了EMD 算法的计算速度及其灵活性。
然而，在实现EMD 算法的实现步骤中，因为CSI 需要求取二阶导数，并且需要执行大量的乘除法运算，为了进一步提高EMD 算法在FPGA 上的采样率与计算速度，提出了使用计算更为简单的锯齿变换（ST）来代替CSI 生成包络。以ST为基础的EMD 算法可以更好地适应硬件环境并可以提高采样率与计算速度，使得在FPGA 上实现的EMD 算法可以更快的速度处理各个频率的信号。
EMD 算法是由Huang 等人[1]提出的，它可以将任何一种非线性和非平稳信号分解成若干个本征模函数（Intrinsic Mode Function，IMF）。
1.1 IMF 的特性
由EMD 算法产生的IMF 需要满足以下两个条件：
（1）在函数的整个时间范围内，局部极值点与函数零点的个数必须相等或最多只差一个。
（2）在任意时间点，由局部极大值定义的上包络线和局部极小值定义的下包络线的均值为0。
1.2 EMD 算法的步骤
对于任意的信号x（t），EMD 算法可以分为以下几个步骤：
（1）求取输入信号x(t)的极大值与极小值；
（2）分别使用CSI 或ST（本文使用的是ST）连接输入信号的极大值和极小值，生成上包络线u(t)和下包络线l(t)；
（3）计算局部均值，即mean(t) =u(t) +l(t)/2；
（4）计算hj(t)，hj(t)=x(t)- meant，其 中j为生成IMF 的索引；
（5）如果hj(t) 满足成为IMF 的条件，则hj(t)是一个IMF，c1(t)=hj(t)；
若hj(t) 不满足，则令x(t)=hj(t)，重新从步骤（1）开始；
（6）计算第一个残差信号，将残差信号记作r1(t)=x(t)-c1(t)，将r1(t) 作为新的输入信号，重复步骤（1）～（4）得到下一个IMF。当最终的IMF为单调分量即信号不可再分解时，无法得到进一步的IMF，EMD 算法的分解结束，经EMD 分解后，最初的信号x(t)可以写成：
其中，得到的IMF个数为n个；
第k个IMF记为ck(t)；
最终残差记为rn(t)。
整个EMD 算法的流程图如图1 所示。
图1 EMD 算法流程图
由Lu[7]提出了基于锯齿变换（ST）的EMD 算法更适合硬件环境的实现，基于ST 的EMD 算法不同于三次样条插值（CSI）在原始数据上进行处理，它是将原始数据函数转换到锯齿空间，然后在锯齿空间中，连接函数的极大值与极小值构建上下包络线。IMF 为变换后的函数与上下包络均值的差值，最后再将IMF 变换到原始的函数空间，得到IMF 分量。这种方法对函数进行一次数据处理，得到唯一的IMF 分量，且不需要进行EMD 算法的筛选过程[7]，速度更快，可预测性更强，且对硬件环境更为友好。
2.1 将原始函数变换到锯齿空间
将原始函数变换到锯齿空间中需要其对应的锯齿函数，锯齿函数可以通过用直线连接函数的极值点实现，如图2 所示为原始函数及其对应的锯齿函数。在锯齿函数的每一段上，其极点值是与原函数的极大极小值一一对应的，而原始函数数据的变化是与锯齿函数的直线段一一对应的。
图2 原始函数与对应的锯齿波函数
锯齿变换可以分为以下几个步骤：
（1）设原始信号有m个极值点，记为：
（2）设有k个极大值组成上包络线U(t)，记作：
设有l个极小值组成下包络线L(t)，记作：
其中，k+l=m。
（3）将原函数的锯齿变换定义为：
除了函数上的点，原函数坐标空间也需要进行变换为对应的锯齿空间，在原函数空间中，点的坐标为(t,y)；
在锯齿空间中，点的坐标为(u,s)，将坐标上的锯齿变换定义为：
从式（5）和式（6）中可以看出锯齿变换并不对纵坐标（函数的值）做变换，而是对横坐标（时间）进行压缩或者拓展，使函数变为分段的锯齿函数。
2.2 在锯齿空间中求函数上下包络线及均值
锯齿函数在极大值U(ti)和极小值L(ti)之间连续的线性变换。
求上包络线可以通过连接连续的极大值取出来，记作：
下包络线同样可以通过连接连续的极小值取出来，记作：
极大值和极小值的和等于极值的个数；
残差是包络的平均值，记作：
IMF 是锯齿函数与均值的差值，记作：
函数c(u)具有上下包络对称性的，满足IMF的两个特性，上包络线和下包络线的均值为0，且在锯齿空间中IMF 的极大值与极小值是关于零值对称的。
图3 原始函数的锯齿变换及在锯齿空间内的上下包络线和均值
图4 锯齿空间中的IMF 及其上下包络线
在锯齿空间中，IMF 和上下包络线都是分段线性的函数，只需要在极值点处操作，不需要重复的操作。
2.3 锯齿空间到原函数空间的转换
要将在锯齿空间中得到的IMF、上下包络线和残差转换到原始函数空间得到cdata(t)、rdata(t)、Udata(t)、Ldata(t)中，只需要将锯齿空间中的结果转换为原始空间中的结果。而锯齿空间对函数的值（纵坐标）没有进行改变，所以只需要在锯齿空间中找到对应横坐标即可，公式如下：
由于锯齿反变换并不会破坏IMF 的包络对称性，所以变换回原空间的函数依然满足IMF 的性质。
基于ST 的EMD 算法将被分解为四个模块在FPGA 上实现，分别是极值寻取模块、ST 生成包络模块、输出计算模块、门控时钟模块。
3.1 算法在FPGA 上的整体实现
本文提出的实现EMD 算法的方法采用自底向上方法实现。首先介绍了极值寻取模块、包络模块和输出计算模块等基本模块，并使用Verilog-HDL 在Vivado 软件上实现，然后将这些模块整合到一起，形成顶层模块。
在本文的设计中，采样后的输入信号被发送到极值寻取模块的输入缓冲区中，用于临时存储输入数据，之后用于计算IMF。如果出现极大值或者极小值的话，那么极值寻取模块会把极大值或极小值输出。此外，该模块根据输出的极大值和极小值情况，分别将极大值或极小值能使信号置为高电平。包络模块的输入控制器加载上包络线和下包络线的输入数据，然后通过包络线生成模块计算上下包络线并输出，之后借助输出计模块计算信号的IMF 和残差信号。在所提出的设计中，极值寻取模块、包络线输入控制器模块和输出计算模块以比例时钟频率工作，而其余模块，即包络线生成模块以输入频率工作。使用这两种不同频率的目的是使系统进行实时处理。当极值寻取模块忙于扫描输入时在求极大值和极小值的同时，包络模块利用循环缓冲区中已经存储的极值来计算上包络和下包络，从而使系统速度更快。整个算法的结构框图如图5 所示。
图5 算法的整体结构图
3.2 极值寻取模块
EMD 算法的第一步是找到输入信号的极大值与极小值，在本文提出的设计中，有一个单独的模块用来实现这一功能，其结构框图如图6 所示。这一模块主要由一个输入缓冲区（由三个寄存器组成）、两个比较器、两个计数器、与门和三态缓冲器组成。程序开始运行时，在每个时钟周期，将输入信号发送到由寄存器X、Y、Z组成的输入缓冲区中。在两个比较器的作下，寄存器X中的数据与寄存器Y中的数据作比较，寄存器Z中的数据也与寄存器Y中的数据作比较。两个比较器的输出G_1 与L_1 通过一个与门（图6 中A1），两个比较器的输出G_2 与L_2 通过另一个与门（图6 中A2）。如果在寄存器Y中的数据是极大值，寄存器Y中的数据比在寄存器X和寄存器Z中的数据更大，与门A1 将会输出高电平，输出极大值的使能端将为高电平。如果寄存器Y中的数据小于寄存器X和寄存器Z中的数据，那么寄存器Y中的数据是极小值，与门A2将输出高电平，极小值的使能端将为高电平。极值寻取模块也有两个独立的计数器。当找到极大值时，极大值计数器清除为0，开始新的计数。当找到极小值时，极小值计数器清除为0，开始新的计数。
图6 极值寻取模块
3.3 使用ST 的生成包络线模块
使用ST 的生成包络线模块的输入是极值寻取模块所得到的极值。此模块的实现借助了有限状态机（FSM），模块内部包括减法、除法、乘法和加法器结构，并且为了进一步提高计算速度，模块内的乘除法器采用了流水线结构。计算上包络线和下包络线的结构是一样的。
计算包络线的模块的算法原理是基于第二节提到的ST，如果A0、A1是连续的极大值或者极小值，C是这两个值之间的时间跨度，那么根据ST 算法，包络线的输出是：
其中，j是计数器从0 到C- 1 的输出。模块的输入与ready 信号是由包络模块的输入控制器给出的，下文将详细介绍。
当ready 信号的值为高电平时，包络模块的FSM 控制器通过减法器获取A0、A1的差值，然后把差值和C的值输入到除法器上。FSM 控制器将除法器输出的结果与j作为乘法器的输入（j的初始值设为0），将乘法器的输出结果与A0相加，并输出结果。与此同时FSM 控制器将包络线生成模块的输出使能信号置为高电平，表示模块的输出结果在输出端可用。然后模块内的计数器将j的值加1，并将再次输入到乘法器中，与除法器的输出结果相乘，乘法器的结果再次与A0相加，然后将加法器的输出结果放在模块的输出处。此后一直重复这个过程，直到j的值等于C- 1。当j的值等于C- 1时，FSM 将完成信号置为高电平并将输出使能信号置为低电平，表示所有点的计算已经完成。
使用ST 计算包络线模块的框图如图7 所示。
图7 使用ST 计算包络线模块
包络模块的输入控制器（以下简称控制器A）的功能是将极值模块的输出加载到生成包络模块。因为使用了单独的模块寻找信号的极值，因此为了检测出是否找到了极大值与极小值，并将其输入到包络模块，设计了控制器A。
控制器A监测极值输入的使能信号，如果信号的极大值被找到，则极大值输入使能信号为高电平，极大值计数寄存器的值加1，控制器A将信号的极大值作为A0输入到上包络生成模块，极大值使能信号恢复为低电平；
紧接着当连续的下一个极值被找到时，极大值输入使能信号为高电平，控制器A将该极值作为A1输入到上包络生成模块，并且极大值计数寄存器的值为2，ready 信号为高电平，表示上包络生成模块已准备好计算，并将极值模块中极大值的计数器的输出作为C输入到上包络生成模块，极大值输入使能信号为低电平，上包络输入模块开始运行，极值数量寄存器的值减1。同样的，下包络线生成模块的输入也由相同结构的控制器给出，控制器A的功能如图8 所示。
图8 控制器A流程图
3.4 输出计算模块
输出计算模块的主要任务是计算上下包络线的平均值，并将均值从输入信号中减去得到IMF 信号。算法模块由加法器、除法器、减法器和控制器组成。控制器用来监测包络生成模块中完成信号的值，如果包络生成模块的完成信号为高电平，则表明包络线计算完成，将包络线的值加载到输出计算模块，计算模块的框图如图9 所示。
图9 输出计算模块
3.5 门控时钟模块
由于在本文的设计中，每个模块不是同时工作的，所以使用门控时钟来降低整个系统的功耗。在这个电路中，使能信号来自另一个模块。由于使能信号的不同，门控信号可能有故障。因此，为了减少故障，门控时钟使用了D锁存器来输入使能信号。在门控时钟的帮助下可以降低各电路的动态功耗。
在本文的设计中，输入信号的一些点可能没有极大值或极小值可以用来生成上、下包络线。在这种情况下，包络生成模块中的FSM 控制器将等待下一次可以生成包络线的点的到来，模块中的其他部分的时钟信号将停止翻转，从而降低动态功耗。门控时钟的电路图如图10 所示。
图10 门控时钟电路
将上述各个模块及总体设计在Vivado 上编译、综合、布局布线后，得到的FPGA 资源使用报告如表1 所示
表1 FPGA 主要资源表
4.1 仿真的结果
将第三节的设计编译、综合、布线之后再在Vivado 上进行波形仿真，仿真的输入波形是由Matlab 产生的。与DSP 相关的EMD 算法大部分适用于低频信号，本文提出的基于FPGA 的EMD算法可用于高频信号的分解，仿真波形是对频率为100 kHz、500 kHz、和5 000 kHz 的混合信号进行验证，其仿真结果在Matlab 如图11 所示，Vivado 中仿真波形如图12 所示。经仿真后发现本文提出的EMD 算法可以通过其不同的模量分解出不同的频率分量。
图11 利用Matlab 绘制EMD 输入信号及IMF
图12 Vivado 中IMF 的仿真波形
4.2 两种设计之间的对比
在本文中，提出了一种不同于传统的生成包络线的方法，本文提出的设计使用ST 生成包络线。该方法在准确度方面相对于传统的EMD 算法较低（ST 生成的包络线与CSI 生成的包络线的相关度为0.969），但是在复杂度方面大大降低复杂性的对比度，如表2 所示，且对硬件环境更为友好。在采样率方面，相对于在161 kHz 下工作的传统的EMD 算法，本文提出的设计可以在25 MHz 上进行采样，并且计算速度也有提高，其对比如表3 所示。
表2 ST 与CSI 的计算步骤复杂度对比
表3 使用ST与使用CSI的EMD算法的采样率与计算速度对比
为了提高EMD 算法的实时性能与计算速度，本文提出了基于ST 的EMD 算法，并在FPGA上进行了实现。相对于传统的基于CSI 的EMD算法，本文的设计在准确度方面略微降低，但在采样率和计算速度上大幅度提高，提高了算法的实时性，使得其在故障检测等实时性要求较高的场合能够得到更好的应用。并且与其他与DSP 相关的设计大部分用于低频信号的分解相比，本文的设计也可以用于处理高频信号。在以后的工作中，考虑将锯齿变换得到的结果进行平滑处理，进一步优化，提高其准确度。
猜你喜欢极小值包络线锯齿猛犸杀手——锯齿虎天天爱科学(2021年10期)2021-10-11一道抽象函数题的解法思考与改编*中学数学研究(江西)(2020年5期)2020-07-03构造可导解析函数常见类型例析*中学数学研究(江西)(2019年11期)2019-12-31由椭圆张角为直角的弦所在直线形成的“包络”中学数学研究(广东)(2019年21期)2019-12-16抛体的包络线方程的推导新教育论坛(2019年35期)2019-09-10极小值原理及应用科技风(2018年19期)2018-05-14一种用于故障隔离的参数区间包络线计算方法现代机械(2018年1期)2018-04-17基于庞特里亚金极小值原理的多运载体有限时间编队控制自动化学报(2017年1期)2017-03-11一种免测量的锯齿磨削方法*精密制造与自动化(2015年4期)2015-11-27EMD 包络线拟合算法改进及在泡沫尺寸趋势提取中的应用服装学报(2015年6期)2015-01-15