CN107844768A

CN107844768A - 基于滑窗迭代原理的一维信号形态学滤波方法

Info

Publication number: CN107844768A
Application number: CN201711052210.5A
Authority: CN
Inventors: 丑永新; 徐本连; 顾亚; 鲁明丽
Original assignee: Changshu Institute of Technology
Current assignee: Changshu Institute of Technology
Priority date: 2017-10-30
Filing date: 2017-10-30
Publication date: 2018-03-27

Abstract

本发明公开了一种基于滑窗迭代原理的形态学滤波方法，根据缓存区存储数据长度分方法分为2个步骤。步骤1、初始化步骤：缓存区数据长度小于等于结构元素长度，缓存区只存储新采样数据点，并初始化一些中间变量；步骤2、迭代步骤：缓存区数据长度大于结构元素长度，在步骤1初始化变量的基础上，以迭代方式计算新采样数据点的滤波结果。相比现有形态学滤波方法，本发明在保持滤波精度不变的前提下，可明显提高几百倍的形态学滤波速度；同时，结构元素和缓存区长度的增加对本发明的运行时间影响不大。本发明适用于一维采样信号的实时处理环境。

Description

基于滑窗迭代原理的一维信号形态学滤波方法

技术领域

本发明涉及一种信号处理方法，具体来说是一种能够更为快速地进行一维信号形态学滤波的基于滑窗迭代原理的形态学滤波方法。

背景技术

形态学滤波方法是一种有效的一维信号处理方法，通过对信号进行膨胀、腐蚀、开、闭等操作，可进行一维信号滤波、波形分割和特征提取等。实际的一维信号处理系统中，一般设置固定长度的缓存区缓存一定长度数据用于处理。缓存区的数据一般采用滑窗法进行更新，即每采集一个新的数据点，就会剔除最早采集的数据点。然后，对缓存区所有数据进行一次形态学滤波。然而，形态学滤波过程中会进行反复的膨胀和腐蚀运算，运算量大，很难在数据更新的间隔期完成一次滤波运算。

目前，一维信号形态学滤波的大部分研究集中在结构元素形状和长度选择，以及其对滤波结果准确性的影响。这些研究大多对已采集一维信号进行离线仿真滤波，只注重处理结果准确性而忽略了其运算复杂度，研究结果很难应用于实际信号处理系统。有研究采用滑窗的方法更新缓存区数据，并对其进行形态学滤波，但由于算法运算量大，很难用于一维信号的实时处理。

因此，有必要设计一种用于一维信号处理系统信号实时处理的形态学滤波法。

发明内容

1、发明目的。

本发明是为了解决现有技术存在的更新速度慢、提高结果的准确性和速度之间的矛盾，难以在实时应用的环境下进行精确快速的使用的问题，提出了一种基于滑窗迭代原理的一维信号形态学滤波方法。

2、本发明所采用的技术方案。

本发明公开了一种基于滑窗迭代原理的一维信号形态学滤波方法，按照如下步骤进行：

步骤1：开辟缓存区，数据长度小于等于结构元素长度时，当前采集的最新数据点放在最高位，只存储新采样数据点，并初始化中间变量；

步骤2：缓存区数据长度大于结构元素长度时，在步骤1初始化变量的基础上，以缓冲区为滑窗，系统每采集一个新的数据点，就会将缓存区最低位的数据剔除，高位数据向低位移动一位，即滑窗法更新数据，并以迭代方式只计算与结构元素长度相同的数据点的滤波结果。

更进一步具体实施方式中，所述的步骤1按照如下步骤进行：

步骤11、数据更新，新数据放在最高位，对缓存区的数据进行更新，即缓存区的数据长度每次增加一位；

步骤12、对结构数据和缓存区数据进行膨胀运算后对其缓存，用于开运算或者闭运算；

步骤13、对结构数据和缓存区数据进行腐蚀运算后对其缓存，用于开运算或者闭运算；

步骤14、对所述的步骤12或13进行开运算进行进一步滤波；

步骤15、对所述的步骤12或13进行闭运算进行进一步滤波；

步骤16、通过混合滤波器进一步滤波；

步骤17、通过交替滤波器进一步滤波；

步骤18、通过混合交替滤波器进一步滤波；

所述步骤2如下：

步骤21：数据更新，系统每采集一个新的数据点，就会将缓存区最低位的数据剔除，高位数据向低位移动一位，新采集数据点放在缓存区最高位，因此，缓存区数据为新低位剔除，高位增加一位的缓存区数据，即滑窗法；

步骤22：所用结构元素，膨胀与腐蚀运算与步骤12和步骤器相同，并采用滑窗法更新运算结果，即最低位数据剔除，高位向低位移动一位，计算结果放在最高位；

步骤23：对滑窗法更新的数据进行开运算，用于进一步滤波运算；

步骤24：对滑窗法更新的数据进行闭运算，用于进一步滤波运算；

步骤25：采用滑窗法更新结果；

步骤26：通过混合滤波器进一步滤波；

步骤27：通过交替滤波器进一步滤波；

步骤28：通过交替混合滤波器进一步滤波。

3、本发明所产生的技术效果。

(1)现有技术的形态学滤波法在实际一维信号处理系统应用中，在滑窗法更新数据点后，每次滤波都有大量重复运算，进而影响系统运行速度。本发明根据实际一维信号采集数据更新过程中滑窗造成数据位置变化的规律，对形态学滤波进行改进，使其采用迭代的方式只计算新采样数据点的滤波结果，即计算与结构元素长度相同的数据点的滤波结果，并不是将结构元素与缓存区所有数据进行运算，进而提高计算速度，并且经过论证比较，与现有形态学滤波法滤波结果的均方根误差为0，即本发明的方法与现有方法滤波效果相同，即虽然简化了90％以上的更新滤波的数据，却实现了与现有技术相同的精度。

(2)本发明改进的形态学滤波器，在一维采样信号系统中尤为适用，就平均运行时间而言，现有方法耗时是本发明方法的几百倍，现有的实时应用环境的数据大部分为一维信号数据，并且实时性要求非常高，这几百倍的速率在实时环境中应用有着巨大的差别，现有的形态学滤波法无法实现数据的实时性。本发明方法有效地提高了信号处理速度，可用于一维采样信号的实时处理。

(4)本发明改进的形态学滤波器采用迭代方法进行，与原滤波器相比，结构元素长度增加对算法运行时间影响不大。

(5)本发明改进的形态学滤波器采用迭代方法进行，与原滤波器相比，缓存区长度增加对算法运行时间影响不大。

综上，本发明提供一种适用于一维信号领域的，实时环境中的快速一维信号形态学滤波方法。将实际一维信号处理系统中采用滑窗法更新数据的方式与形态学滤波原理相结合，对现有形态学滤波方法进行改进。在每次滑窗更新缓存区数据后，采用迭代方式只计算新数据点的滤波结果。克服现有方法每更新一次数据，需要对缓存区所有数据进行处理的缺陷，减少了90％以上的数据计算量，在保持滤波结果准确性不变前提下，提高几百倍的滤波速度。

附图说明

图1是本发明基于滑窗迭代原理的形态学滤波法的流程图。

图2是本发明的基于滑窗迭代原理的形态学滤波法的初始化过程。

图3是本发明的基于滑窗迭代原理的形态学滤波法的迭代过程。

图4是本发明的仿真信号及滤波结果比较示意图。其中：左侧图为含有噪声的仿真信号；右侧图为仿真信号的滤波结果。

图5是本发明的结构元素长度对滤波时间的影响比较示意图。其中，左侧图为本发明方法对3秒仿真信号的运行时间；右侧图为现有方法对3秒仿真信号的运行时间。

图6是本发明的缓存区长度对滤波时间的影响比较示意图。其中，左侧图为本发明方法对3秒仿真信号的运行时间；右侧图为现有方法对3秒仿真信号的运行时间。

图7是是本发明的实测信号及滤波结果比较示意图。其中：左侧图为含有噪声的脉搏信号；右侧图为滤波结果。

图8是本发明的结构元素长度对滤波时间的影响比较示意图。其中，左侧图为本发明方法对30秒脉搏信号的运行时间；右侧图为现有方法对30秒脉搏信号的运行时间。

图9是本发明的缓存区长度对滤波时间的影响比较示意图。其中，左侧图为本发明方法对30秒脉搏信号的运行时间；右侧图为现有方法对30秒脉搏信号的运行时间。

具体实施方式

以下结合附图，详细说明本发明的实施方式。如图1和图2所示，根据信号采集系统缓存中数据更新特点，该方法可分为初始化和迭代两个步骤。

步骤1：缓存区数据长度小于等于结构元素长度，缓存区以只存储新采样数据点，并初始化一些中间变量。具体步骤如下：

步骤11：一维信号采集系统，开辟缓存区存储采集的数据，缓存区信号长度小于结构元素长度。参照图1，当前采集的最新数据点存储在缓存区的第i位，记为x(i)，1≤i≤M。如图1所示，则与这段信号进行运算的结构元素长度应该为i。定义当前缓存区数据为X＝[x(1),x(2),…,x(i)]，所用结构元素G＝[g(1),g(2),…,g(i)]。

步骤12：膨胀运算定义为：

膨胀运算的结果记为d(i)，对其缓存可得D＝[d(1),d(2),…,d(i)]；用于开运算或者闭运算。

步骤13：腐蚀运算Θ为：

XΘG＝max{X-G}

腐蚀运算的结果记为e(i)，对其缓存可得E＝[e(1),e(2),…,e(i)]；用于开运算或者闭运算。

步骤14：由腐蚀和膨胀运算定义开闭运算，其中，开运算为：

开运算的结果记为o(i)，对其缓存可得O＝[o(1),o(2),…,o(i)]；用于进一步滤波运算。

步骤15：闭运算·为：

闭运算的结果记为c(i)，对其缓存可得C＝[c(1),c(2),…,c(i)]；用于进一步滤波运算。

步骤16：混合滤波器F_hf为：

步骤17：交替滤波器为：

步骤18：交替混合滤波器为：

F_ah＝(F_oc+F_co)/2＝(oc(i)+co(i))/2

步骤2：缓存区数据长度大于结构元素长度，在步骤1初始化变量的基础上，以缓冲区为滑窗，采用滑窗法更新数据，并以迭代方式计算新采样数据点的滤波结果。

步骤21：参照图2，当缓存区的数据长度大于结构元素长度M时，进入迭代过程。系统每采集一个新的数据点，就会将缓存区最低位的数据剔除，高位数据向低位移动一位x(1)＝x(2),x(2)＝x(3),…,x(M-1)＝x(M)，新采集数据点放在缓存区最高位x(M)。因此，缓存区数据为X＝[x(1),x(2),…,x(M)]，所用结构元素G＝[g(1),g(2),…,g(M)]；

步骤22：膨胀与腐蚀运算与步骤12和步骤13相同，所不同的是采用滑窗法更新运算结果，即D和E的最低位数据剔除，高位向低位移动一位，计算结果放在最高位d(M)和e(M)，此时D＝[d(1),d(2),…,d(M)]，E＝[e(1),e(2),…,e(M)]。

步骤23：更进一步具体实施方式中，开运算为：

步骤24：闭运算·为：

闭运算的结果记为c(i)，对其缓存可得C＝[c(1),c(2),…,c(i)]；用于进一步滤波运算；

步骤25：采用滑窗法更新结果，得到O＝[o(1),o(2),…,o(M)]，C＝[c(1),c(2),…,c(M)]。

步骤26：进一步，混合滤波器为：

步骤27：交替滤波器为：

步骤28：交替混合滤波器为：

F_ah＝(F_oc+F_co)/2＝(oc(M)+co(M))/2

实施例1：

该实施例中，对正弦波和噪声叠加的仿真信号进行滤波。其中，正弦波表达如式所示：

S(i)＝Asin(2πf₀/f_s i)

式中，A为信号的幅值，f₀为信号的频率，f_s为信号采样频率，i＝1,2,…,N,。这些量可随机取值，本例选为A＝4，f₀＝1Hz，f_s＝1000Hz。给正弦波叠加均值为1，方差为0的白噪声，合成信号参照图3左侧示意图，对该信号进行滤波的仿真。

仿真条件是在MATLAB2016a软件下运行。

仿真结果参照图3右侧示意图，为长度为3秒的仿真信号滤波结果，与未加噪声的正弦波{S(i)}的均方根误差为0.5969，可有效地去除信号中的高频噪声。而与现有形态学滤波法滤波结果的均方根误差为0，即本发明的方法与现有方法滤波效果相同。

结构元素长度是影响形态学滤波方法运算量的重要因素之一，增加G的长度，从5个采样点增加到55个采样点，本发明方法及现有形态学滤波器对仿真信号运算时间分别参照图4左侧示意图和图4右侧示意图。对比可知，结构元素改变对本发明方法运算速度影响不大，始终保持在0.11s～0.12s之间，而对于现有形态学方法，其运算时间由20.5s增加到24.24s；就平均运行时间而言，现有方法耗时是本发明方法的192倍。这是因为相比于现有方法，本方法在每次滑窗更新数据后，没有让结构元素与缓存区所有数据进行运算，而是通过迭代的方式只更新新采集数据进行运算，从而节约了运算时间。对于3s长的信号，本发明方法完全满足信号处理的实时性，而现有形态学方法则不行。

一维信号采集系统缓冲区长度决定每次数据更新后形态学滤波的信号处理量，增加缓存区长度，从512个采样点增加到4096个采样点，本发明方法及现有形态学滤波器对仿真信号运算时间分别参照图5左侧示意图和图5右侧示意图。对比可知，缓存区长度增加对本发明方法运算速度影响不大，始终保持在0.12s～0.13s之间，而对于现有形态学方法，其运算时间由12.58s增加到38.24s；就平均运行时间而言，现有方法耗时是本发明方法的252倍。因为本方法采用迭代的方式，只有与结构元素长度相同的数据参与运算，参照图(2)，而无需对缓存区所有数据进行处理，因此节约了运算时间。对于3s长的信号，本发明方法完全满足信号处理的实时性，而现有形态学方法则不行。

实施例2：

该实施例中，将本发明的方法和现有形态学滤波法用于实测信号。

实测信号为人体脉搏信号，信号采样频率为1000Hz，采样时间为30s，参照图6左侧示意图。

仿真条件是在MATLAB2016a软件下运行。

仿真结果参照图6右侧示意图，本发明方法与现有形态学滤波法滤波结果的均方根误差为0，即本发明的方法与现有滤波效果相同。

对于结构元素长度对方法运行时间的影响，参照图7(a)图7右侧示意图，图7左侧示意图为本发明方法的运算时间，图7右侧示意图为现有方法的运算时间。对比可知，本发明方法总耗时在1s左右，而现有形态学滤波法由233.1s增长到267.9s，就平均运行时间而言，现有方法耗时是本发明方法的234倍。结构元素长度增加对本发明方法没有影响。

对于缓存区长度对方法运行时间的影响，本发明方法和现有方法运行时间分别参照图8左侧示意图和8右侧示意图。对比可知，缓存区长度增加对本发明方法运算速度影响不大，始终保持在0.92s～0.93s之间，而对于现有形态学方法，其运算时间由122.8s增加到868.2s；就平均运行时间而言，现有方法耗时是本发明方法的542倍。对于30s长的信号，本发明方法完全满足信号处理的实时性，而现有形态学方法则不行。

上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受上述实施例的限制，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。

Claims

1.基于滑窗迭代原理的一维信号形态学滤波方法，其特征在于按照如下步骤进行：

2.根据权利要求1所述的基于滑窗迭代原理的一维信号形态学滤波方法，其特征在于，所述的步骤1按照如下步骤进行：

步骤14、对所述的步骤12或13进行开运算进行进一步滤波；

步骤15、对所述的步骤12或13进行闭运算进行进一步滤波；

步骤16、通过混合滤波器进一步滤波；

步骤17、通过交替滤波器进一步滤波；

步骤18、通过混合交替滤波器进一步滤波；

所述步骤2如下：

步骤25：采用滑窗法更新结果；

步骤26：通过混合滤波器进一步滤波；

步骤27：通过交替滤波器进一步滤波；

步骤28：通过交替混合滤波器进一步滤波。

3.根据权利要求1或2任一所述的基于滑窗迭代原理的一维信号形态学滤波方法，其特征在于：所述步骤1如下：

步骤11：数据更新，一维信号采集系统，开辟缓存区存储采集的数据；缓存区信号长度小于结构元素长度，当前采集的最新数据点存储在缓存区的第i位，记为x(i)，1≤i≤M，新数据放在最高位，缓存区信号进行运算的结构元素长度应该为i，即定义当前缓存区数据为X＝[x(1),x(2),…,x(i)]，所用结构元素G＝[g(1),g(2),…,g(i)]；

步骤12：X和G的膨胀运算形式为：

<mrow> <mi>X</mi> <mo>&CirclePlus;</mo> <mi>G</mi> <mo>=</mo> <mi>m</mi> <mi>a</mi> <mi>x</mi> <mo>{</mo> <mi>X</mi> <mo>+</mo> <mi>G</mi> <mo>}</mo> </mrow>

膨胀运算的结果记为d(i)，对其缓存可得D＝[d(1),d(2),…,d(i)]；用于开运算或者闭运算；

步骤13：X和G的腐蚀运算Θ为：

XΘG＝max{X-G}

腐蚀运算的结果记为e(i)，对其缓存可得E＝[e(1),e(2),…,e(i)]；用于开运算或者闭运算；

步骤14：X和G的开运算为：

开运算的结果记为o(i)，对其缓存可得O＝[o(1),o(2),…,o(i)]；用于进一步滤波运算；

步骤15：X和G闭运算·为：

<mrow> <msub> <mi>F</mi> <mi>c</mi> </msub> <mo>=</mo> <mrow> <mo>(</mo> <mi>X</mi> <mo>&CenterDot;</mo> <mi>G</mi> <mo>)</mo> </mrow> <mo>=</mo> <mrow> <mo>(</mo> <mi>X</mi> <mo>&CirclePlus;</mo> <mi>G</mi> <mi>&Theta;</mi> <mi>G</mi> <mo>)</mo> </mrow> <mo>=</mo> <mi>D</mi> <mi>&Theta;</mi> <mi>G</mi> </mrow>

步骤16：混合滤波器F_hf为：

步骤17：交替滤波器为：

步骤18：交替混合滤波器为：

F_ah＝(F_oc+F_co)/2＝(oc(i)+co(i))/2。

4.根据权利要求1或2任一所述的基于滑窗迭代原理的一维信号形态学滤波方法，其特征在于：所述步骤2如下：

步骤21：数据更新，系统每采集一个新的数据点，就会将缓存区最低位的数据剔除，高位数据向低位移动一位x(1)＝x(2),x(2)＝x(3),…,x(M-1)＝x(M)，新采集数据点放在缓存区最高位x(M)，因此，缓存区数据为X＝[x(1),x(2),…,x(M)]，即滑窗法；

步骤22：所用结构元素G＝[g(1),g(2),…,g(M)]，膨胀与腐蚀运算与步骤12和步骤器相同，并采用滑窗法更新运算结果，即D和E的最低位数据剔除，高位向低位移动一位，计算结果放在最高位d(M)和e(M)，此时D＝[d(1),d(2),…,d(M)]，E＝[e(1),e(2),…,e(M)]；

步骤23：开运算为：

步骤24：闭运算·为：

步骤25：采用滑窗法更新结果，得到O＝[o(1),o(2),…,o(M)]，C＝[c(1),c(2),…,c(M)]；

步骤26：进一步，混合滤波器为：

步骤27：交替滤波器为：

步骤28：交替混合滤波器为：

F_ah＝(F_oc+F_co)/2＝(oc(M)+co(M))/2。