CN110333389B - 基于插值dft的正弦信号频率估计方法 - Google Patents

Abstract

本发明公开了一种基于插值DFT的正弦信号频率估计方法，适用于电力系统的电能质量评估。由于实正弦信号在频谱上存在镜像的负频率分量，非相干条件下负频率分量的频谱泄露将使得DFT的直接使用产生较大的估计误差。针对该问题，提出了一种基于信号相邻谱线的两点IpDFT算法。首先，对时域信号x(n)做离散傅里叶变换，通过搜索DFT幅度谱获取峰值谱线对应的索引值，即频率的粗估计值k₀。然后比较峰值谱线左右谱线幅值定位幅度次大谱线。最后由最大的两条DFT谱线构造插值比例，计算频率粗估计值的偏差δ₀从而得到完整的信号频率。该方法可以有效抑制频谱泄露干扰，具有较高的估计精度和抗噪性能，在频谱干扰严重的场景优越性更加明显。

Description

基于插值DFT的正弦信号频率估计方法

技术领域

本发明涉及信号处理和电力系统领域，特别是涉及基于插值DFT的正弦信号频率估计方法。

背景技术

噪声背景下正弦信号的频率估计问题，不仅是现代数字信号处理中一项重要的研究内容，随着信息技术的发展，现已广泛应用于移动通信、电能质量监测等诸多工程领域。在电力系统中，系统频率是控制电网中发电负荷不平衡的关键参数，可用于确定非线性负载和其他电子元器件注入的谐波分量等。对于电网信号，只有保持标称频率值，才能维持系统的稳定性与电气设备的正常运行。因此，准确的频率估计是电网中一项必不可少的需求，而电力系统中电压波形又可建模为正弦波信号，于是对含噪正弦信号的频率估计技术研究，对分析与治理电能质量、维持电网系统稳定运行具有重大的研究意义与实用价值。

目前，国内外的文献中已经提出了许多正弦信号的频率估计方法，按照所使用的信号特征，可以将现有的估计方案分为时域估计算法和频域估计算法两大类。在基于正弦信号时域特征估计频率参数的算法中，研究和改进的方向主要包括最大似然估计算法、线性预测算法、自相关相位法等。而基于信号频域特征的正弦信号估计算法通常先对观测的采样信号进行离散傅里叶变换(DFT)，然后根据信号的频域特性提取出频率等特征。

DFT估计器由于物理意义明确、实现简单且计算效率高等诸多优点，成为正弦信号频率估计领域一种广受关注与拓展的方法，但其在非相干采样环境中伴随着的频谱泄露、栅栏效应等固有缺陷将影响其直接使用的估计精度，特别是对于频谱上包含镜像负频率成分的实正弦信号，需要额外的步骤进一步提升DFT频率估计器的性能。在基于DFT的改进频率估计算法中，加窗和插值等改进措施虽然能够在一定程度上提升DFT频率估计器的精度，但其未从根本上考虑频谱泄露的来源，而是试图忽略或者默认加窗可以将其抑制到可以忽略的程度，这将导致算法的估计结果存在一定误差，特别是在正负频率分量间隔较近时，分量间的频谱泄露干扰将使此类改进方案产生显著的估计偏差。此外，为了进一步提高估计精度，迭代DFT算法应运而生，然而迭代的引入会产生附加的计算量，造成巨大的计算负担，限制了算法的实时处理能力。

发明内容

为了解决以上问题，本发明提供基于插值DFT的正弦信号频率估计方法，将负频率分量的长程泄露完全考虑在内，可以有效解决信号时域截断带来的频谱泄露问题，具有极高的估计精度、抗噪性能和抗镜像分量干扰的能力，为达此目的，本发明提供基于插值DFT的正弦信号频率估计方法，包含以下步骤：

步骤1：采集电力系统中的待处理正弦信号，无噪环境下等间距采样得到N点离散正弦信号序列x(n)满足

其中ω₀、A和

分别表示信号的数字角频率、幅值和相位，对于样本长度为N的信号序列，数字角频率ω₀可以表示为ω₀＝2πl₀/N＝2π(k₀+δ₀)/N，其中l₀称为信号的规格化频率，物理上表征采样得到的信号周期的个数，k₀∈{0,1,2,…,N-1}和|δ₀|≤0.5分别是信号频率l₀的整数和小数部分；

步骤2：对所述的正弦信号序列x(n)做N点离散傅里叶变换，得到

X(k),k＝0,1,…,N-1，

具体表示为

步骤3：搜索信号DFT谱线序列X(k)定位幅度最大谱线获取系统频率的粗估计结果，即得到整数部分频率k₀的估计值，记为

其中|·|表示取绝对值操作；

步骤4：比较峰值谱线

左右两条谱线的幅值获取插值方向i，则幅度次大的谱线为

当

时，比较

和

如果

则i＝1，反之i＝-1，当

时始终有i＝1；

步骤5：对幅度最大两根DFT谱线

和

进行相位旋转，根据谱线的旋转结果

和

计算插值比例，进而求解小数部分频率的估计值

结合整数部分频率的估计值

就可最终得到数字角频率的估计结果

进一步的，所述步骤5中计算插值比例估计小数部分频率包括如下步骤：

步骤5.1：将选取的两条DFT谱线

和

分别乘以相位旋转因子

和

并将旋转后的变换序列定义为

和

步骤5.2：计算旋转后峰值谱线的实部和虚部值，分别记为

和

其中Re[·]和Im[·]表示取实数部分和虚数部分操作；

步骤5.3：根据频率粗估计结果

插值方向i以及

和

的大小关系，定义和计算插值比例因子η：

步骤5.4：推导发现(1)式所示的插值比例可以表示为小数部分频率δ₀的一元二次函数，反解方程得到两个根，并选取绝对值较小的解作为小数部分频率的估计值

其中

步骤5.5：结合整数部分估计值，得到最终的数字角频率估计值

本申请基于插值DFT的正弦信号频率估计方法，在展开推导实正弦信号DFT峰值谱线与其相邻次大谱线的谱线值基础上，提出了一种考虑负频率分量频谱泄露的新型IpDFT算法，有效解决了现有插值DFT频率估计器普遍存在的频谱干扰问题。进一步而言，本方案同样分为粗估计与细估计两个阶段，频率的粗估计值由定位DFT幅度最大的谱线实现，而在频率细化阶段，分别建立了峰值及其相邻谱线实部或虚部之间的独立插值比例关系式，这一比例可以推导为信号小数部分频率的一元二次函数，通过求解这一比例关系即可实现估计偏差的校正。然而分析发现如果使用两种插值比例独立求解，估计的结果会出现畸变点或其他失效的情形。为了消除上述问题，增强算法的估计性能与鲁棒性，本发明通过补充的条件判断信号采用谱线实部还是虚部进行插值计算的两点IpDFT方法，与传统的正弦信号频率估计算法相比，本发明具有极高的估计精度和抗镜像分量干扰的能力，特别是频谱泄露干扰严重的正负分量间隔较近的正弦信号，该方法的优越性更加显著。

与现有的技术相比，本发明具有以下优点：1.插值过程中充分考虑了实正弦信号镜像负频率分量的长程泄露，有效抑制了频谱泄露产生的频率估计误差，增强了方法的抗镜像分量干扰能力，提升了方法的估计精度和抗噪性能。2.与传统的正弦信号插值DFT算法相比，该方法由于考虑频谱泄露问题，在提供无偏估计结果的同时，具有更加广泛的应用场景，适用于正负频率分量间隔较近的正弦信号。3.仅利用两点DFT谱线值求解信号频率，计算公式简单且计算复杂度较低。

附图说明

图1为本发明提出方法的实现流程图；

图2为本发明中不同插值比例计算的频率估计均方误差随小数部分频率δ₀的变化图，其中图2(a)为k₀＝2时的均方估计误差曲线，图2(b)为k₀＝1时的均方估计误差曲线；

图3为本发明提出方法的频率估计均方误差随整数部分频率k₀的变化图；

图4为本发明提出方法的频率估计均方误差随小数部分频率δ₀的变化图；

图5为本发明提出方法的频率估计均方误差随信噪比SNR的变化图。

具体实施方式

下面结合附图与具体实施方式对本发明作进一步详细描述：

本发明提供基于插值DFT的正弦信号频率估计方法，在通过对信号DFT幅度谱进行谱峰搜索获得频率粗估计结果的基础上，对峰值谱线左右的频率分量值进行比较获取幅度次大的谱线，然后利用这两条谱线构造插值比例，即可通过一元二次方程求解出小数部分频率值完成信号数字角频率的估计。方法将负频率分量的长程泄露完全考虑在内，可以有效解决信号时域截断带来的频谱泄露问题，具有极高的估计精度、抗噪性能和抗镜像分量干扰的能力。

本发明具体实施例如下：

无噪环境下电力系统的电压波形可以以离散时间形式表示为

其中A，ω₀∈[0,π]和

分别表示正弦信号的幅度、数字角频率和相位。对于N个样本点的观测信号序列，数字角频率ω₀可以表示为

其中l₀称为信号的规格化频率，物理上表征采样得到的信号周期的个数。k₀∈{0,1,2,…,N-1}和|δ₀|≤0.5分别是信号频率l₀的整数和小数部分。

下面讨论利用插值DFT思路估计正弦信号的频率ω₀。首先对上述采集的N点正弦序列做N点离散傅里叶变换，得到变换后的谱线序列X(k),k＝0,1,…,N-1

然后根据DFT变换后的谱线序列获得信号整数部分频率k₀的估计值，即系统频率的粗估计结果。由信号处理的知识，实正弦信号反应到频域上除了自身频率分量外，还有一个镜像的负频率分量-2π(k₀+δ₀)/N。在利用DFT频谱定位最大幅度谱线位置获得频率粗估计值k₀时，如果k₀＞N/2，要将其转换到正频率区间，即

接下来则只需要利用信号的DFT谱线间的差异进行频率偏差的校正即可完成最终的系统频率估计，也就是寻求插值比例关系计算小数部分频率δ₀的估计值。

对于(2)式，为了简化表达，令

则公式(2)所示的信号DFT变换谱线序列可以表示为

对于任意索引值为k∈{0,1,…,N-1}处的谱线分量，定义DFT谱线乘以相位旋转因子e^-jπk/N后的结果为X′(k)＝X(k)·e^-jπk/N

由于X′(k)为复数，可将其表示为实部加上虚部的形式X′(k)＝X′_R(k)+jX′_I(k)，其中X′_R(k)＝Re[X′(k)]和X′_I(k)＝Im[X′(k)]分别为X′(k)的实部和虚部值，Re[·]和Im[·]分别取实部和虚部操作。具体地，X′_R(k)和X′_I(k)可以展开为

对于幅度最大的峰值谱线

旋转后的结果，则其实部和虚部具有如下近似关系

类似地，峰值谱线左右两边的相邻谱线的实部虚部同样可以近似表示为

由数字信号处理的基础知识，正弦信号的真实频率落在

相邻两条谱线

和

之间。进一步地，系统频率落在索引为

的谱线与其相邻的幅度次大谱线之间，且该条幅度次大谱线的索引值需要通过比较峰值谱线左右两边相邻谱线的幅值获得。于是，接下来需要确定DFT谱线幅度次大的谱线，即确定插值方向i

在定位到幅度最大的两根谱线后，定义变量R_Re和R_Im分别表示两条谱线旋转后的实部和虚部比值，表达式展开如下

注意到公式(13)和(14)的右边均是小数部分频率δ₀的一元二次方程，于是，以此为基础分别定义基于信号DFT谱线实部和虚部的插值比例因子η_Re和η_Im

然而需要注意的是分析公式(16)发现，如果k₀+i＝0或k₀＝0，其分母将为0，换言之，当使用的两根谱线中包含索引值为0的谱线，以上谱线虚部的比例关系可能不成立，即此时不能使用公式(16)所示的虚部插值比例关系进行小数部分频率的估计。为了使所提出的算法具有更高的估计精度与鲁棒性，除了上述问题，单独使用实部或虚部插值比例还将出现性能陡降的畸变点问题，这一现象可以通过公式(5)来解释。

在公式(5)中，存在相位因子

当相位

DFT点数N和规格化频率l₀使cos(φ)＝0，则公式(6)表示的谱线实部数值等于0，而当sin(φ)＝0时，谱线的虚部变为0，即在上述情况下，式(15)和(16)的插值比例将不再包含正弦信号的信息，而完全是高斯白噪声的描述。注意到cos(φ)＝0和sin(φ)＝0不会同时满足，结合公式(16)的成立条件，于是可以通过以下方案来避免虚部插值关系失效和奇异点的出现，即补充相关判断条件，来决定采用虚部还是实部的插值关系形成最终的插值因子。

根据公式(13)和(14)，则上述插值比例因子可以表示成为小数部分频率δ₀的一元二次方程

接下来，反解方程并选取绝对值较小的解即可得到小数部分频率的估计值

其中

最后整合整数部分估计值，得到最终的数字角频率估计值

以下结合一个实施算例的仿真结果对本发明做进一步阐述。

取采样后的待处理电网信号为

其中q(n)表示测量信号中的干扰噪声，仿真中采用均值为0、方差为

的高斯白噪声建模。信号的幅度A＝2、相位

观测的样本点数和DFT点数均为N＝128。为了更好地评估本发明方法的频率估计性能，引入克拉美罗下界(CRLB)做参照，CRLB表征信号参数估计均方误差所能达到的最低结果，即估计算法的最佳性能界限，对于本发明分析的实正弦信号模型，其频率估计CRLB的计算公式如下：

其中N为观测信号的样本点数，SNR表示正弦信号的信噪比

由纯净正弦信号的功率和噪声功率的比值计算得到。算法的实施流程如图1所示。

在第一组模拟中，通过改变信号的小数部分频率δ₀来验证本发明提出的基于插值DFT的正弦信号频率估计方法不存在单独使用DFT谱线实部或虚部插值比例所具有的畸变点和虚部插值关系失效的问题。图2(a)和图2(b)分别画出了SNR＝40dB下整数部分频率k₀＝2和k₀＝1时，使用实部插值比例因子η_Re、虚部插值比例因子η_Im以及本发明提出的实部虚部插值关系相结合方法的频率估计均方误差(MSE)随δ₀的变化情况。由图中本发明提出的估计方法不存在关系式失效和性能陡降的畸变点现象，算法具有更高的鲁棒性和估计精度。

为了进一步验证与分析所提出频率估计算法的性能特点与优势，接下来从不同参数维度观察本发明方法的性能表现。图3分析了信号整数部分k₀发生变化时，算法的频率估计均方误差，这里设置信噪比SNR＝40dB、小数部分频率δ₀＝0.25，如图所示，算法的均方误差在k₀发生改变时，始终逼近克拉美罗下界CRLB，具有良好的估计性能。图4则给出了固定k₀＝1时，δ₀从-0.5到0.5之间按照0.01步长递增的情况下，算法的频率估计结果，观察发现算法的估计精度始终接近最优的估计上限，即在信号正负频率分量相距较近的频谱泄露干扰严重情形，算法仍然可以提供准确的估计结果。

最后为了更为深入地分析本发明所提出方法的抗噪性能，固定信号的频率k₀＝1、δ₀＝0.25，观察所提出算法的频率估计表现随着信噪比SNR的变化情况，如图5所示，本方法的频率估计均方误差随着信噪比的增加呈近似线性下降趋势，算法的估计精度随着信噪比的增加而提高，因此，算法有效地解决了频谱泄露对DFT估计结果的影响。

综上所述，本发明所提出的基于插值DFT的正弦信号频率估计方法性能优异，且有效抑制了频谱泄露产生的频率估计误差，增强了方法的抗镜像分量干扰能力，具有运算简单、估计精度高和抗噪性能强等特点。

本发明的技术方案可应用于电能质量监测等领域的正弦信号频率估计场景。

以上所述，仅是本发明的较佳实施例而已，并非是对本发明作任何其他形式的限制，而依据本发明的技术实质所作的任何修改或等同变化，仍属于本发明所要求保护的范围。

Claims (1)

1.基于插值DFT的正弦信号频率估计方法，其特征在于：包含以下步骤：

步骤1：采集电力系统中的待处理正弦信号，无噪环境下等间距采样得到N点离散正弦信号序列x(n)满足

其中ω₀、A和

分别表示信号的数字角频率、幅值和相位，对于样本长度为N的信号序列，数字角频率ω₀可以表示为ω₀＝2πl₀/N＝2π(k₀+δ₀)/N，其中l₀称为信号的规格化频率，物理上表征采样得到的信号周期的个数，k₀∈{0,1,2,…,N-1}和|δ₀|≤0.5分别是信号的规格化频率l₀的整数和小数部分；

步骤2：对所述的正弦信号序列x(n)做N点离散傅里叶变换，得到X(k),k＝0,1,…,N-1，具体表示为

步骤3：搜索信号DFT谱线序列X(k)定位幅度最大谱线获取数字角频率的粗估计结果，即得到整数部分频率k₀的估计值，记为

其中|·|表示取绝对值操作；

步骤4：比较峰值谱线

左右两条谱线的幅值获取插值方向i，则幅度次大的谱线为

当

时，比较

和

如果

则i＝1，反之i＝-1，当

时始终有i＝1；

步骤5：对幅度最大两根DFT谱线

和

进行相位旋转，根据谱线的旋转结果

和

计算插值比例，进而求解小数部分频率的估计值

结合整数部分频率的估计值

就可最终得到数字角频率的估计结果

所述步骤5中计算插值比例估计小数部分频率包括如下步骤：

步骤5.1：将选取的两条DFT谱线

和

分别乘以相位旋转因子

和

并将旋转后的变换序列定义为

和

步骤5.2：计算旋转后峰值谱线的实部和虚部值，分别记为

和

其中Re[·]和Im[·]表示取实数部分和虚数部分操作；

步骤5.3：根据频率粗估计结果

插值方向i以及

和

的大小关系，定义和计算插值比例因子η：

步骤5.4：推导发现(1)式所示的插值比例可以表示为小数部分频率δ₀的一元二次函数，反解方程得到两个根，并选取绝对值较小的解作为小数部分频率的估计值

其中

步骤5.5：结合整数部分估计值，得到最终的数字角频率估计值

Images