本发明属于数字信号处理技术领域，尤其涉及一种基于emd(empiricalmodedecomposition,经验模式分解)的零相位滤波器设计方法。
背景技术
工程中信号绝大多数是非平稳信号，数字滤波器是工程数据采集中信号处理不可缺少的一个环节。数字滤波器处理从系统响应角度分析，就是模特性上乘以系统频率响应的模，相位特性上就是在输入相位的基础上附加了系统频率响应的相位。因此信号附加了系统频率响应的相位，改变了每个分量之间的相对相位关系，即使系统频率响应的模对信号每个分量是常数，信号在时域特性上也发生了很大的变化。因此，对于非平稳信号使用一般滤波器滤波不仅会产生相位失真，还会改变信号的瞬时频率，造成信号信息丢失。为了保证非平稳信号的准确相位信息，在相关的信号处理中需要选用零相位滤波技术。传统的零相位滤波器设计方法有二种：frr和rrf。frr滤波方法是，先将输入信号序列按顺序滤波(forwardfilter)，然后将所得结果逆转后反向通过滤波器(reversefilter)，再将所得结果逆转后输出(reverseoutput)，即可获得零相位失真输出信号序列；rrf滤波法是，先将输入信号序列反转后通过滤波器(reversefilter)，然后将所得结果逆转后再次通过滤波器(reversefilter)，这样所得的结果(forwardoutput)即为零相位失真的输出序列。
由以上看到，目前的零相位滤波技术，首先数字滤波，然后再将相位差调整为零，最终实现零相位滤波；其次该方法在第二步时需要将信号反转，也就是时间反转，再通过相同的滤波器。而时间反转在物理结构上是无法实现的。以上数字滤波器技术都是基于傅里叶变换，而该理论要求处理对象是平稳信号，而用来处理非平稳信号在理论上也存在矛盾，并不能保证非平稳信号准确的瞬时相位信息。
n.e.huang等人在1998年提出了一种新的信号处理方法，称为hilbert-huang变换(简称hht)。在hht理论中信号被定义为由有限个imf(intrisicalmodefunction,简称特征模式函数)组成。imf在每个时间点上只存在单频率。hht理论包含了平稳信号的傅里叶分解。为了得到信号的imf，n.e.huang等人提出了emd(empiricalmodedecomposition,经验模式分解)，可以将信号分解为有限个imf。emd主要依据信号自身的极值点，将信号分解为有n个imfi(i＝1,2,3……n)。该分解方法是依据信号自身的极值点，能够在信号的局部自适应的分解信号，这与基于固定基分解信号的傅里叶变换、小波变换有本质区别。因此，该方法也是公认的最适合处理非平稳信号的方法。
emd的具体实现过程见文献：huangne，shenz，longs.r，etal.theempiricalmodedecompositionandthehilbertspectrumfornonlinearandnon-stationarytimeseriesanalysis.proc.roy.soc，london.a.1998，454，903-995。由于该过程不是一个严密的数学推导过程，是由曲线拟合、边端处理、imf判断条件、筛选停止条件这些数据处理技术实现的。这些技术算法的优劣决定了emd对信号的分解效果，这其中尤其以imf判断条件对分解效果影响最大。有专家学者提出了相应的imf判断条件，但是这些判断条件不完善导致了筛选过程中会出现一些异常效果。
技术实现要素：
发明目的：为了克服现有技术中存在的不足，本发明提供一种能保证非平稳信号准确的瞬时相位信息，同时信号筛选稳定的基于emd的零相位滤波器设计方法。
技术方案：为实现上述目的，本发明采用的技术方案为：
一种基于emd的零相位滤波器设计方法，包括以下步骤：
步骤1，确定emd对信号的频率分解率。
取二个频率不相等的谐波信号x1(t)、x2(t)，emd分解x1(t)+x2(t)得到的第一阶imf1，后停止分解，此时余量为r(t)，第一阶imf1记为c1(t)，则分解效果衡量指标index如下式：
以(a2，f2，index)为坐标，绘制二维灰度图，通过二维灰度图确定分解效果分割线，根据分解效果分割线确定满足分解效果衡量指标index的emd的频率分解率β。
步骤2，构建滤波信号
设初始信号x(t)的频率范围为f0～fn，需要滤波的频率为：fn～fn，其中fn之f0。通过emd的频率分解率β根据不同的情况构建滤波信号，取x(t)＝x(t)，构建按以下规则：
如果fn×β≤fn，则取f＝fn/β，；如果fn×β＞fn，则取f＝fn，构建滤波信号：
以上amax为x(t)幅值最大值，amin为x(t)幅值最小值，t表示时间。
步骤3，利用emd分解信号实现滤波，过程如下步骤301～步骤306：
步骤301，将滤波信号y(t)与信号x(t)构建信号h(t)。
步骤302，利用emd分解信号h(t)，得到imf1与剩余r(t)。利用emd分解信号h(t)的方法如下：设当前的被筛选对象为imf'(t)，基于当前的被筛选对象imf'(t)生成的极大值包络线为mmax(t)、极小值包络线mmin(t)，均值为m(t)，则：
其中，ξ表示mmax(t)、mmin(t)均值平方和与被筛选imf'(t)平方和比值，当被筛选对象的ξ小于给定阈值时，当前的imf'(t)就为一阶imf1，剩余量为r(t)。
步骤303，取k(t)＝imf1-y(t)。
步骤304，定义ji(t)＝k(t)，如果fn×βi≤fn，则fn～fn的频段信号滤波结束，执行步骤4。如果fn×βi＞fn，则需要继续下面的过程。i为自然数，初始值为1。
步骤305，取x(t)为r(t)。取x(t)幅值最大值amax，取x(t)幅值最小值amin。如果fn×βi+1≤fn，则f＝fn/β。；如果fn×βi+1＞fn，则f＝fn×βi，利用步骤2的公式(2)构建新的滤波信号。步骤305执行结束，i＝i+1。
步骤306，重复步骤301～步骤305。
步骤4，如果步骤3结束，则得到fn～fn完整的信号j(t)。
优选的：所述步骤1中当频率比f2≥0.85时，index≥0.3，emd的频率分解率β＝f2＝0.85。
优选的：所述步骤301中构建信号h(t)：h(t)＝y(t)+x(t)。
优选的：所述步骤302中给定阈值为10-7。
优选的：所述步骤4中完整信号j(t)为：j1(t)+j2(t)+j3(t)+…ji(t)。
优选的：步骤1中信号x1的频率为1，a1幅值为10。信号x2的频率为f2，0＜f2＜1，幅值为a2，，0＜a2＜10，x1(t)、x2(t)信号间的相位差为0。
本发明相比现有技术，具有以下有益效果：
1.通过确定分解效果衡量指标index，在根据分解效果分割线确定满足分解效果衡量指标index的emd的频率分解率β，通过频率分解率β构建滤波信号，使得本发明能够保证非平稳信号准确的瞬时相位信息，实现非平稳信号的零相位滤波。
2.本发明采用均值线与被筛选对象之间的比值作为imf判断条件的方法，当这个比值小于一个恒定值时，就认为不再需要剔除均值余量，获得了合格的imf，使得信号在筛选过程中的稳定，不会产生异常情况。
附图说明
图1为emd的频率分解率示意图。
图2为非平稳信号时域图。
图3为非平稳信号频谱图。
图4为非平稳信号85k～185khz频带信号时域图。
图5为非平稳信号85k～185khz频带信号频谱图。
具体实施方式
下面结合附图和具体实施例，进一步阐明本发明，应理解这些实例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。
一种基于emd的零相位滤波器设计方法，通过确定分解效果衡量指标index，再根据分解效果分割线确定满足分解效果衡量指标index的emd的频率分解率β，通过频率分解率β构建滤波信号，使得本发明能够保证非平稳信号准确的瞬时相位信息，实现非平稳信号的零相位滤波。提出一种新的imf判断条件，以优化emd过程，更好的实现基于emd的零相位滤波技术。为了得到合格的imf，被筛选对象需要不停的剔除均值余量，就是基于极大值包络线、极小值包络线生成的均值线。本发明提出了一种均值线与被筛选对象之间的比值作为imf判断条件的方法，当这个比值小于一个恒定值时，就认为不再需要剔除均值余量，获得了合格的imf。设当前的被筛选对象为imf'(t)。设基于imf'(t)生成的极大值包络线为mmax(t)、极小值包络线mmin(t)，均值为m(t)。
这里指定ξ≤10-7。在一般的筛选过程中，也就是当被筛选对象的ξ小于给定阈值10-7时，就认为不需要再剔除均值余量了，当前的imf'(t)就为一阶imf，可以进行新的筛选，寻找下一阶imf了。本发明只需要第一阶imf，在筛选得到imf1后，就停止筛选，剩余量为r(t)。本发明所用的emd是基于新imf判断条件实现的。
具体包括以下步骤：
步骤1，首先确定emd对信号的分解能力，本发明称之为频率分解率
取二个频率不相等的谐波信号x1(t)、x2(t)，信号x1的频率f1为1，幅值a1为10。信号x2的频率为f2，0＜f2＜1，幅值为a2，，0＜a2＜10，x1(t)、x2(t)信号间的相位差为0，emd分解x1(t)+x2(t)得到的第一阶imf1，后停止分解，此时余量为r(t)，第一阶imf1记为c1(t)，则分解效果衡量指标index如下式：
以(a2，f2，index)为坐标，绘制二维灰度图，如图1所示，通过二维灰度图确定分解效果分割线，根据分解效果分割线确定满足分解效果衡量指标index的emd的频率分解率β。由图1可以得到，当频率f2/f1≥0.85时，index≥0.3，此时定义β＝f2/f1＝0.85为emd的频率分解率。也就是本发明认为index大于0.3时，就认为emd无法分解x1(t)、x2(t)这二个信号。这样定义emd的频率分解率为0.85。下面再直接执行2、3步骤。
步骤2，构建滤波信号
设初始信号x(t)的频率范围为f0～fn，需要提取滤波的频率为：fn～fn，其中fn≥f0。因为emd存在频率分解率，需要根据不同的情况构建滤波信号，因此通过emd的频率分解率β根据不同的情况构建滤波信号，取x(t)＝x(t)，初次构建按以下规则：
如果fn×β≤fn，则取f＝fn/β，；如果fn×β＞fn，则取f＝fn，构建滤波信号：
以上amax为x(t)最大幅值值，amin为x(t)最小幅值值，t表示时间。
步骤3，利用emd分解信号实现滤波，过程如下步骤301～步骤306：
步骤301，将滤波信号y(t)与信号x(t)构建信号h(t)：h(t)＝y(t)+x(t)。
步骤302，利用emd分解信号h(t)，得到imf1与剩余r(t)。利用emd分解信号h(t)的方法如下：设当前的被筛选对象为imf'(t)，基于当前的被筛选对象imf'(t)生成的极大值包络线为mmax(t)、极小值包络线mmin(t)，均值为m(t)，则：
其中，ξ表示mmax(t)、mmin(t)均值平方和与被筛选imf'(t)平方和比值，当被筛选对象的ξ小于给定阈值时，当前的imf'(t)就为一阶imf1，进行新的筛选，寻找下一阶imf。
步骤303，取k(t)＝imf1-y(t)。
步骤304，定义ji(t)＝k(t)，如果fn×βi≤fn，则fn～fn的频段信号滤波结束，执行步骤4；如果fn×βi＞fn，则需要继续下面的过程。i为自然数，初始值为1。
步骤305，取x(t)为r(t)。取x(t)幅值最大值amax，取x(t)幅值最小值amin。如果fn×βi+1≤fn，则f＝fn×βi/β。；如果fn×βi+1＞fn，则f＝fn×βi，利用步骤2的公式(2)构建新的滤波信号。步骤305执行结束，i＝i+1。
步骤306，重复步骤301～步骤305。
步骤4，如果步骤3结束，则得到fn～fn完整信号j(t)为：j1(t)+j2(t)+j3(t)+…ji(t)。
实例
图2为一非平稳信号x(t)，图3是图2的频谱图，现在要求滤波的频率范围：85k～185khz，fn＝85hz，fn＝185hz，如图4、5所示，通过以上步骤，实现零相位滤波。
第一步：确定所用emd的频率分解率，根据本发明的方法，确定所用的emd是0.85。
取x(t)＝x(t)，fn＝85hz，fn＝185hz，由于185×0.85＞85，所以取f＝185hz。
取最大幅值值amax(图4中为2)，最小幅值值amin(图4中为-2)。
第二步：利用步骤2，构建滤波谐波信号：y(t)＝20×cos(2π×185×t)。
第三步：
(1)构建h(t)＝x(t)+y(t)。
(2)利用emd分解h(t)，得到imf1，剩余量r(t)。
(3)k(t)＝imf1-y(t)
(4)取ji(t)＝k(t)，由于
fn×0.85i＝185×0.85i＞85，i初始值为1，执行(5)。
(5)fn×0.85i+1＝185×0.85i+1＞85，fn×0.85i＝157.25hz，取f＝157hz。
取x(t)＝r(t)，取x(t)的最大幅值amax，最小幅值amin，i＝i+1，利用步骤2中的公式(2)构建新的滤波信号。
(6)重复(1)～(5)，直到fn×0.855≤fn，取f＝fn＝85hz/0.85，再次执行(1)～(4)后，在(4)跳出第三步，执行第四步。
第四步：得到85～185khz频带的信号为：j1(t)+j2(t)+j3(t)+j4(t)+j5(t)(一共需要5次将此频段信号提取结束)。
以上发明设计的零相位滤波器，既滤出了指定的频段信号，也保留了信号原来的相位信息。
以上所述仅是本发明的优选实施方式，应当指出：对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。