CN109194305B - 基于密度估计的数字化仪均值滤波方法 - Google Patents

Abstract

本发明公开了一种基于密度估计的数字化仪均值滤波方法，首先数字化仪进行过采样，对得到的数据样本序列进行分组得到若干子样本，然后依次将子样本添加至密度估计样本中，利用密度估计和贝叶斯信息融合方法得到数据样本序列的概率密度函数，计算得到数据样本序列的标准不确定度；按照过采样倍数对数据样本序列重新进行分组，根据标准不确定度得到重新分组后各个子样本的置信区间，进行粗大误差剔除，然后对得到的样本进行平均抽取，得到最终样本。本发明通过对采样得到的数据样本序列分组进行密度估计，基于密度估计结果设置置信区间进行粗大误差剔除，使均值滤波结果更加准确。

Description

基于密度估计的数字化仪均值滤波方法

技术领域

本发明属于数字化仪技术领域，更为具体地讲，涉及一种基于密度估计的数字化仪均值滤波方法。

背景技术

带宽、采样率和存储深度是普遍采用的评价数字化仪性能的三大核心指标。除此之外，还有一个非常重要的指标：垂直分辨率。较高的垂直分辨率意味着更加精细的波形显示和更加精确的信号测量。数字化仪的垂直分辨率取决于模数转化器(ADC)的量化位数，以及数字化仪本身的噪声和失真水平。虽然普遍用ADC的位数来描述数字化仪的垂直分辨率，但准确的参数是数字化仪整个***的有效位数(ENOB)。ENOB与信号和噪声失真比(SNR)密切相关，两者的数学关系为：

SNR(in dB)＝6.02×ENOB+1.76

根据以上关系式，SNR大约每增加6dB，ENOB增加1bit，所以提高信噪比就能提高分辨率。目前中高端数字化仪普遍具有这样的功能：在硬件性能定型的情况下，通过算法提高垂直分辨率。这些功能本质上都是对采样后的数据进行数字信号处理，目的是提高信噪比，从而在理论上提高垂直分辨率。

常见的数字化仪一般采用8位或12位ADC，在ADC的量化位数一定的情况下，数字滤波是常用的降噪手段。其中，平均是一种改善信噪比、提高有效位数(ENOB)的典型方法。一般具体是通过形成平均和高分辨率的捕获方式实现的。

如果信号是周期可重复的，数字化仪每采集n段重复波形，把它们按触发位置对齐，相加后除以n，得到1段平均后的波形。图1是信号平均的示意图。信号平均后，信噪比(dB)可以提高10lg(n)。例如4次平均，信噪比可以提高大约6dB，分辨率提高1bit。64次平均，信噪比提高18dB，分辨率提高3bit。平均算法是最简单的分辨率改善方法，但有以下局限和代价：

(1)只能用于周期信号；

(2)平均只能降低随机的、不相关的高斯白噪声，对于相关的噪声和干扰作用不大；

(3)平均算法会显著降低示波器波形更新率。例如欲将分辨率从8bit改善到12bit，需要每采集256次做一次平均，波形更新率降低不止256倍；

(4)平均可能引起波形失真。每段波形是按触发位置对齐，而触发点在每段波形上的相对位置并不固定。图2是触发抖动示意图。如图2所示，因为噪声会导致信号每次穿越触发电平的相对水平位置不一样，触发位置抖动的结果是每段波形被错位相加，平均后波形失真。当信号本身的噪声比较大时，失真现象更明显。

高分辨率捕获方式是广泛运用于数字化仪采集***的另一种平均算法。这种算法的好处是可以用于非周期信号。数字化仪以最大采样率采集波形，将一段波形的每n个采样点分成1组，对1组中n个采样点取平均得到1个数据点存入采集存储器中。最终显示的波形是样点分组平均后的数据，样点数量被抽取了n倍。图3是3次高分辨率捕获平均示意图。如图3所示，高分辨率捕获本质上也是一种数字滤波。样点抽取倍数越高，分辨率提升得越多，在数字化仪中高分辨率捕获也经常用来实现可变带宽功能。同样，高分辨率捕获方式也有自身的缺陷与不足：

(1)高采样率和高垂直分辨率不可兼得；

(2)示波器带宽和垂直分辨率也不可兼得；

(3)高分辨率捕获平均算法对抑制高斯白噪声有效果，而对于示波器ADC的INL(积分非线性)引起的噪声没有抑制作用。

综上所述，以上两种平均方法对噪声的抑制和分辨率的提升均具有一定的局限性。连续捕获平均基于对信号的多次重复采集，因而受限于被测信号的重复特性；连续样点平均虽然基于对被测信号的单次采集，却在噪声的滤除过程中牺牲了模拟带宽。除此之外，受环境干扰、时钟抖动以及传输延迟等因素的影响，示波器的采样数据中还经常包含有粗大量化误差，例如一些不相干的噪声、数据失配导致的量化错误等等。对于包含粗大误差的采集样本，若不经过处理而直接平均，结果将大大偏离真实信号。

发明内容

本发明的目的在于克服现有技术的不足，提供一种基于密度估计的数字化仪均值滤波方法，通过对采样得到的数据样本序列分组进行密度估计，基于密度估计结果设置置信区间进行粗大误差剔除，使均值滤波结果更加准确。

为实现上述发明目的，本发明基于密度估计的数字化仪均值滤波方法包括以下步骤：

S1：数字化仪以V倍于实际采样率的速率进行过采样，记采样得到的数据样本序列X＝{x₁,x₂,…,x_N}，N表示采样点数量，将这N个样本数据平均分为M组，得到M个子样本，每个子样本中有Q个数据，记第m个子样本为X_m，X_m＝{x_(m-1)Q+1,x_(m-1)Q+2,…,x_mQ}，m＝1,2,…,M；

S2：初始化密度估计样本

初始化子样本序号m＝1，设置先验概率密度函数

S3：将子样本X_m添加至密度估计样本Y中，即令Y＝Y∪X_m；

S4：对当前密度估计样本Y进行密度估计得到样本Y的似然函数

S5：根据先验概率密度函数

和似然函数

利用贝叶斯信息融合方法进行概率密度函数的融合得到后验概率密度函数

S6：判断是否m＜M，如果是，进入步骤S7，否则进入步骤S8；

S7：令m＝m+1，返回步骤S3；

S8：根据

计算数据样本序列X的标准不确定度σ：

其中，

表示数据样本序列X的平均值；

S9：将数据样本序列X＝{x₁,x₂,…,x_N}根据过采样倍数V进行分组，每V个数据为一组，记所得到的子样本数量为D，记第m个子样本为X^* _d，X^* _d＝{x_(d-1)V+1,x_(d-1)V+2,…,x_dV}，d＝1,2,…,D；

S10：计算得到每个子样本X^* _d中数据样本的均值

根据步骤S8得到的标准不确定度σ设置子样本X_d′的置信区间

将子样本X^* _d中在置信区间以外的样本剔除，并加入同等数量的均值

得到样本更新后的子样本X_d′＝{x′_(d-1)V+1,x′_(d-1)V+2,…,x′_dV}；

S11：对步骤S10得到的每个子样本X_d′进行平均抽取，得到平均值

将D个平均值A_d构成的序列作为最终样本输出。

本发明基于密度估计的数字化仪均值滤波方法，首先数字化仪进行过采样，对得到的数据样本序列进行分组得到若干子样本，然后依次将子样本添加至密度估计样本中，利用密度估计和贝叶斯信息融合方法得到数据样本序列的概率密度函数，计算得到数据样本序列的标准不确定度；按照过采样倍数对数据样本序列重新进行分组，根据标准不确定度得到重新分组后各个子样本的置信区间，进行粗大误差剔除，然后对得到的样本进行平均抽取，得到最终样本。

本发明具有以下有益效果：

(1)本发明结合密度估计、贝叶斯信息融合以及过采样下的平均算法，提高了数字化仪的信噪比，提升了等效垂直分辨率的同时克服了对带宽的影响；

(2)本发明中采用贝叶斯信息融合的方法对概率密度函数进行处理，使密度函数的估计更为精确，使最终输出的测量结果更加准确；

(3)本发明在进行密度估计时提出一种频域SVD多小波密度估计，避免了阈值选择带来的不便，通过奇异值分解可以得到小波基系数的有效奇异值，滤除小波基系数中的信息量较少的分量，从而构造出概率密度函数的无偏估计量，在保证有用的信息不受损害，能够准确描述信息的同时，对密度函数进行很好的估计。

附图说明

图1是信号平均的示意图；

图2是触发抖动示意图；

图3是3次高分辨率捕获平均示意图；

图4是本发明基于密度估计的数字化仪均值滤波方法的具体实施方式流程图；

图5是本实施例中频域SVD多小波密度估计方法的流程图；

图6是不同信号的奇异值特征示例图。

具体实施方式

下面结合附图对本发明的具体实施方式进行描述，以便本领域的技术人员更好地理解本发明。需要特别提醒注意的是，在以下的描述中，当已知功能和设计的详细描述也许会淡化本发明的主要内容时，这些描述在这里将被忽略。

实施例

图4是本发明基于密度估计的数字化仪均值滤波方法的具体实施方式流程图。如图4所示，本发明基于密度估计的数字化仪均值滤波方法的具体步骤包括：

S401：过采样数据分组：

数字化仪以V倍于实际采样率的速率进行过采样，V的值根据实际需要设置即可。记采样得到的数据样本序列为{x₁,x₂,…,x_N}，N表示采样点数量，将这N个样本数据平均分为M组，得到M个子样本，每个子样本中有Q个数据，记第m个子样本为X_m＝{x_(m-1)Q+1,x_(m-1)Q+2,…,x_mQ}。M的值根据实际需要设置。

对测量信号进行过采样是为了在不影响带宽的情况下对采集数据进行平均滤波，提高信噪比和等效垂直分辨率。对测量样本分组是为了实施贝叶斯信息融合得到更加精确的概率密度函数和测量不确定度，对于粗大误差的判别更加准确，进一步提升滤波效果。

S402：数据初始化：

初始化密度估计样本

初始化子样本序号m＝1，设置先验概率密度函数

S403：添加子样本：

将子样本X_m添加至密度估计样本Y中，即令Y＝Y∪X_m。

S404：密度估计：

对当前密度估计样本Y进行密度估计得到样本Y的似然函数

密度估计是概率学中的一种常用技术，一般分为参数密度估计和非参数密度估计两大类，参数密度估计主要包括最大似然估计和贝叶斯估计，非参数密度估计主要包括直方图密度估计、核密度估计、k-近邻估计以及小波密度估计。就本发明而言，可以根据被采集信号的实际情况来选择密度估计方法，本实施例中采用频域SVD多小波密度估计方法。

对于实数集合

是平方可积实函数空间，存在相互正交的尺度函数φ(x)和小波基函数ψ(x)，则任意函数

均能用小波序列函数表示。假设x₁,x₂,…,x_n是一列独立同分布的随机变量，有共同的密度函数f(x)，可以采用如下公式表示：

其中，J表示小波分解尺度，a_J,k表示尺度函数系数，d_j,k表示第j个尺度下的小波基系数，j＝1,2,…,J，k＝1,2,…,K，K表示每个分解尺度下参数数量。

把密度函数f(x)推广到多小波中，相应的多小波基和尺度函数都是由r个向量组成的向量函数，分别用Ψ(x)和Φ(x)表示，对应的小波基系数和尺度函数系数分别用a和d表示，根据单小波估计理论可以给出基于多小波的密度函数估计表达式：

用多小波基Ψ(x)和尺度函数Φ(x)对方程两边做内积可得到多小波基系数和尺度函数系数。

尺度函数系数：

其中，E表示求取均值。

多小波基系数：

则尺度函数系数和小波基系数的无偏估计为：

最终得到概率密度函数：

根据小波分解理论可以知道要处理的样本信息大部分都包含在小波基系数和尺度函数系数中，如果部分小波基系数和尺度函数系数包含的信息量较少，可以选则一定的阈值，使得这部分小于特定阈值的系数值等于零，只保留包含信息量较多的不小于特定阈値的系数。这样对数据进行处理，在保证有用的信息不受损害、能够准确描述信息的同时，还能减少计算量，对密度函数进行很好的估计。阈值有两种类型，即硬阈值和软阈值，硬阈值函数由于不连续，重构所得信号会产生振荡。软阈值函数消噪信号比较光滑，但有较大的信号失真。为了避免阈值选择带来的不便，本实施对多小波密度估计的小波基系数进行频域奇异值分解(SVD)，通过奇异值分解可以得到小波基系数的有效奇异值，滤除小波基系数中的信息量较少的分量，从而构造出概率密度函数的无偏估计量。

图5是本实施例中频域SVD多小波密度估计方法的流程图。如图5所示，本实施例中频域SVD多小波密度估计方法的具体步骤如下：

S501：多小波密度估计：

首先根据以下公式计算多小波分解尺度J：

J＝lbn-lblgn

其中，n表示当前密度估计样本的数据数量。

设置多小波基，得到尺度J下的多小波尺度函数Φ(x)和多小波基函数Ψ(x)，计算得到尺度函数系数a_J,k和每一分解尺度下的多小波基系数d_j,k。

S502：傅里叶变换：

对每一分解尺度下的K个多小波基系数d_j,k构成的序列进行傅里叶变换，得到长度为K的频域序列，记频率序列中第k个元素为

将频域序列重构成大小为α×β的矩阵D_j，重构方法为：采用长度为β的滑动窗口以大小为1的步长在频域序列滑动，共计得到α个子序列，将子序列作为行向量构成矩阵D_j。可见，α+β-1＝K。矩阵D_j可以表示为：

S503：奇异值分解：

对矩阵D_j的奇异值分解，得到D_j＝U_jS_jV_j，U_j为α×α阶酉矩阵，V_j为β×β阶酉矩阵，S_j为半正定α×β阶对角矩阵，其对角线元素为矩阵D_j的奇异值，对奇异值从大到小进行排序，保留前P个奇异值，将其他奇异值置为0，得到对角矩阵S_j′。P的大小根据实际需要设置即可。图6是不同信号的奇异值特征示例图。如图6所示，不同信号的奇异值曲线都存在一个拐点，将较小奇异值置为0。

S504：更新多小波基系数：

根据对角矩阵S_j′计算得到矩阵D_j′＝U_jS_j′V_j，根据步骤S502中的重构方式将矩阵D_j′恢复为长度为K的频域序列，再对频域序列进行傅里叶逆变换得到多小波基系数d′_j,k。采用以上方式，可以去除信息量较少的多小波基系数，且本实施例中在频域进行奇异值分解，可以使得到的结果更加稳定。

S505：得到似然函数：

根据尺度函数系数a_J,k和步骤S504每一分解尺度下的多小波基系数d_j,k，采用概率密度函数计算公式计算得到概率密度函数，即为当前密度估计样本的似然函数。

S405：计算后验概率密度函数：

根据先验概率密度函数

和似然函数

利用贝叶斯信息融合方法进行概率密度函数的融合得到后验概率密度函数

为了提高采用多小波密度估计得到的概率密度函数的精度，本发明采用贝叶斯理论来进行概率密度函数的融合。贝叶斯统计推断的主要特点在于以先验分布的形式融入统计推断中，先验分布是指以往的经验、资料以及直观分析与预测等先验信息。因此，即使缺少样本信息的情况下，也能根据先验分布做出统计推断。在获得现实样本后，按似然函数与先验分布乘积为核心的后验分布做出统计推断。这样的一种统计推断方法，对于进一步的提高统计推断的可靠性提供了有利保证。由于之前步骤中已经计算得到了先验概率密度函数和似然函数，本实施例中采用以下公式计算得到后验概率密度函数

S406：判断子样本是否添加完毕，即是否m＜M，如果是，进入步骤S407，否则进入步骤S408。

S407：令m＝m+1，返回步骤S403。

S408：计算标准不确定度：

根据

计算数据样本序列X的标准不确定度σ：

其中，

表示数据样本序列X的平均值。本实施例中

的计算公式如下:

S409：过采样数据重分组：

将数据样本序列X＝{x₁,x₂,…,x_N}根据过采样倍数V重新进行分组，每V个数据为一组，记所得到的子样本数量为D，记第m个子样本为X^* _d，X^* _d＝{x_(d-1)V+1,x_(d-1)V+2,…,x_dV}，d＝1,2,…,D。

S410：剔除粗大误差：

由于一次采集得到的样本数据是独立同分布的，因此其概率密度函数和不确定度都可以统一为总样本的概率密度函数和不确定度，而贝叶斯信息融合得到的概率密度函数和不确定度具有更高的精度，因此在本步骤实施过程中采用原样本序列的概率密度函数与不确定度来进行预估与剔除，其具体方法如下：

计算得到每个子样本X_d′中数据样本的均值

本实施例中均值

的计算公式如下：

根据步骤S408得到的标准不确定度σ设置子样本X^* _d的置信区间

将子样本X^* _d中在置信区间以外的样本剔除，并加入同等数量的均值

得到样本更新后的子样本X_d′＝{x′_(d-1)V+1,x′_(d-1)V+2,…,x′_dV}。

S411：基于平均的抽取：

对步骤S410得到的每个子样本X_d′进行平均抽取，得到平均值

将D个平均值A_d构成的序列作为最终样本输出。

尽管上面对本发明说明性的具体实施方式进行了描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

Claims (5)

1.一种基于密度估计的数字化仪均值滤波方法，其特征在于，包括以下步骤：

S1：数字化仪以V倍于实际采样率的速率进行过采样，记采样得到的数据样本序列X＝{x₁,x₂,…,x_N}，N表示采样点数量，将这N个样本数据平均分为M组，得到M个子样本，每个子样本中有Q个数据，记第m个子样本为X_m，X_m＝{x_(m-1)Q+1,x_(m-1)Q+2,…,x_mQ}，m＝1,2,…,M；

S2：初始化密度估计样本

初始化子样本序号m＝1，设置先验概率密度函数

S3：将子样本X_m添加至密度估计样本Y中，即令Y＝Y∪X_m；

S4：对当前密度估计样本Y进行密度估计得到样本Y的似然函数

S5：根据先验概率密度函数

和似然函数

利用贝叶斯信息融合方法进行概率密度函数的融合得到后验概率密度函数

S6：判断是否m＜M，如果是，进入步骤S7，否则进入步骤S8；

S7：令m＝m+1，返回步骤S3；

S8：根据

计算数据样本序列X的标准不确定度σ：

其中，

表示数据样本序列X的均值；

S9：将数据样本序列X＝{x₁,x₂,…,x_N}根据过采样倍数V进行分组，每V个数据为一组，记所得到的子样本数量为D，记第m个子样本为X^* _d，X^* _d＝{x_(d-1)V+1,x_(d-1)V+2,…,x_dV}，d＝1,2,…,D；

S10：计算得到每个子样本X^* _d中数据样本的均值

根据步骤S8得到的标准不确定度σ设置子样本X′_d的置信区间

将子样本X^* _d中在置信区间以外的样本剔除，并加入同等数量的均值

得到样本更新后的子样本X′_d＝{x′_(d-1)V+1,x′_(d-1)V+2,…,x′_dV}；

S11：对步骤S10得到的每个子样本X′_d进行平均抽取，得到平均值

将D个平均值A_d构成的序列作为最终样本输出。

2.根据权利要求1所述的数字化仪均值滤波方法，其特征在于，所述步骤S3中密度估计采用频域SVD多小波密度估计方法，包括以下步骤：

S3.1：根据以下公式计算多小波分解尺度J：

J＝lbn-lblgn

其中，n表示当前密度估计样本Y的数据数量；

设置多小波基，得到尺度J下的多小波尺度函数Φ(x)和多小波基函数Ψ(x)，计算得到尺度函数系数a_J,k和每一分解尺度下的多小波基系数d_j,k；

S3.2：分别对每一分解尺度下的K个多小波基系数d_j,k构成的序列进行傅里叶变换，得到长度为K的频域序列，记频率序列中第k个元素为

将频域序列重构成大小为α×β的矩阵D_j，重构方法为：采用长度为β的滑动窗口以大小为1的步长在频域序列滑动，共计得到α个子序列，将子序列作为行向量构成矩阵D_j；

S3.3：对矩阵D_j的奇异值分解，得到D_j＝U_jS_jV_j，U_j为α×α阶酉矩阵，V_j为β×β阶酉矩阵，S_j为半正定α×β阶对角矩阵，其对角线元素为矩阵D_j的奇异值，对奇异值从大到小进行排序，保留前P个奇异值，将其他奇异值置为0，得到对角矩阵S′_j；

S3.4：根据对角矩阵S′_j计算得到矩阵D′_j＝U_jS′_jV_j，根据步骤S3.2中的重构方式将矩阵D′_j恢复为长度为K的频域序列，再对频域序列进行傅里叶逆变换得到多小波基系数d′_j,k；

S3.5：根据尺度函数系数a_J,k和步骤S3.4每一分解尺度下的多小波基系数d_j,k，采用概率密度函数计算公式计算得到概率密度函数，即为当前密度估计样本的似然函数。

3.根据权利要求1所述的数字化仪均值滤波方法，其特征在于，所述步骤S4中概率密度函数的融合采用以下计算公式完成：

4.根据权利要求1所述的数字化仪均值滤波方法，其特征在于，所述步骤S8中均值

的计算公式如下：

5.根据权利要求1所述的数字化仪均值滤波方法，其特征在于，所述步骤S10中均值

的计算公式如下：

Images