一种抗多频干扰的实正弦信号频率高精度估计方法

1.本发明涉及信息技术领域，具体涉及信号处理技术中的参数估计问题。


背景技术：

2.信息技术领域中普遍存在信号处理需求，由傅里叶级数可知，任何周期信号都可以由不同频率的正弦或余弦波叠加而成，因此，正弦信号参数估计问题成为了各类信号处理技术发展的基础。多频干扰实正弦信号作为一种特定的正弦信号模型，广泛涉及电力系统、仪器仪表、装置测量以及无损检测等技术领域。
3.多频干扰实正弦信号可以看成多个单频实正弦信号的叠加，通常其中能量最强(幅值最大)的单频信号为待估参数信号，根据欧拉公式，多频干扰实正弦信号可以进一步分解为多个正频率分量和负频率分量的叠加，由此将导致信号频谱分析过程中，存在着正频率和负频率成分的相互干扰，从而严重降低频谱分析的精度。为此，国内外许多学者对其展开研究，提出了许多可行的解决方法，比如采用加窗减少干扰信号的频谱泄露，降低干扰频率对待估信号的影响(张鸿博,蔡晓峰,鲁改凤.基于双窗全相位fft双谱线校正的电力谐波分析[j].仪器仪表学,2015,36,(12):2835-2841.)，但这种方法增加了等效噪声带宽，同时由于加窗后信号主瓣展宽，使得该类方法不适用于干扰频率与待估频率临近的情况。文献(djukanovic s,popovic v.efficient and accurate detection and frequency estimation of multiple sinusoids[j].ieee access,2019,7:1118-1125.)提出了一种基于频移和低通滤波思想的方法，本发明中称为ds方法，该方法采用了迭代干扰消除的算法结构，通过反复滤除负频率分量干扰，实现了较高的频率估计精度，但该方法对负频率干扰的消除并不彻底。文献(ze-long mou,ya-qing tu,peng chen,kui wang.accurate frequency estimation of multiple complex and real sinusoids based on iterative interpolation[j].digital signal processing,2021,117:103173.)提出了一种减除负频率干扰的方法，本发明中称为mzl方法，该方法同样采用迭代干扰消除的算法结构，虽然在迭代中使用了精度更高的三点插值，但在负频率干扰的抑制性能上则与ds方法类似，尤其在干扰频率与待估频率临近的情况下估计精度下降明显。
[0004]
综上，虽然已有很多关于多频实正弦信号参数估计的研究，但受负频率干扰影响，现有方法在干扰频率与待估频率临近情况下表现不佳，急需针对性研究解决。


技术实现要素：

[0005]
为解决负频率干扰下临近频率干扰难以有效抑制的问题，本发明提出一种抗多频干扰的实正弦信号频率高精度估计方法。该方法中提出了一种新的零频谱线校正技术，可以有效抑制实信号负频率干扰影响，通过逐次迭代、逐频干扰消除的双循环运算结构，该方法可以很好地抑制临近频率干扰，实现待估信号频率的高精度估计。
[0006]
本发明的技术方案是，一种抗多频干扰的实正弦信号频率高精度估计算法，总体技术方案包括：
①
对于已知频率阶数l的多频干扰信号x(n)进行逐阶循环的粗频率估计，频
率阶数为l时的具体操作为：首先，通过补零fft运算搜索谱峰以获得信号频率估计初值其次，基于进行零频谱线校正计算信号复振幅然后，利用二分插值法得到信号频率粗估计值再次，基于更加精确的进行零频谱线校更新值；最后，基于和构造出第l阶信号并实现干扰消除。
②
采用迭代运算和逐阶干扰消除的双循环结构实现信号频率的高精度估计，迭代次数为i(最大迭代次数为q)且频率阶数为l时的具体操作为：首先，基于已知和在信号x(n)中消除l阶以外的信号以实现干扰消除；然后，采用估计精度较高的三点插值方法获得信号频率精估计值最后，基于更加准确的进行零频谱线校正更新复振幅
③
当迭代次数达到q且遍历了所有频率阶数后，返回信号频率估计结果
[0007]
所述步骤
①
中基于补零fft的频率估计初值计算步骤包括：首先，令s(n)＝x(n)并进行信号补零得到se(n)＝[s(n),0,0,
…
,0]，补零长度可根据精度和计算复杂度需要自行设定；然后，进行快速傅里叶变换得到s(m)＝fft(se(n))并进行频谱峰值搜索得到最后，频率估计初值即为
[0008]
所述步骤
①
和
②
中所有涉及零频谱线校正的计算步骤均为：根据已知条件将或h(n)＝s(n)或y(n)带入q＝round(nω
*
/π)和下式进行计算：
[0009][0010]
所述步骤
①
中利用二分插值法得到信号频率粗估计值计算步骤为：首先，进行负频率干扰消除，即然后，计算dtft变换最后，得到频率粗估计值
[0011]
所述步骤
①
中的信号干扰消除是令
[0012]
所述步骤
②
中消除l阶以外的干扰是令
[0013]
所述步骤
②
中采取的精度较高的三点插值方法的计算步骤是：首先，进行负频率干扰消除，即然后，分别计算和最后，采用下式计算信号频率的精估计值：
[0014][0015]
本发明提出的频率估计方法可以在干扰信号幅度小于待估信号幅度的50％的情况下，实现干扰信号与待估信号频率间隔在1.5倍dft频隙以上时的高精度频率估计，在对抗临近频率干扰方面明显高于其他类似方法的估计性能，同时，本发明方法也可实现各阶
干扰信号频率的准确估计。
附图说明
[0016]
为清楚地说明本发明技术方案，下面给出实施例描述所需的附图和简要说明。显然所示附图仅是本发明的一些实施例，本领域技术人员在不付出创造性劳动前提下还可根据所示附图创造其他类似附图。本发明所示附图有：
[0017]
图1本发明信号处理流程图；
[0018]
图2本发明与其他方法在不同信噪比取值下的频率估计方差比较(l＝2)；
[0019]
图3本发明与其他方法在不同信噪比取值下的频率估计方差比较(l＝4)；
[0020]
图4本发明与其他方法在不同信号长度下的频率估计方差比较；
[0021]
图5本发明与其他方法在不同信号频率间隔下的频率估计方差比较。
具体实施方式
[0022]
以下结合附图和实施例进一步说明本发明的详细技术内容，本实施例以本发明技术方案为前提给出了详细的实施方式和计算流程，但本发明的保护范围不限于下述实施例，应理解实施例仅为了说明本发明，而不是为了限制本发明的保护范围。
[0023]
根据图1所示信号处理流程，考虑输入为已知频率阶数l的多频干扰实正弦信号其中a
l
、ω
l
和θ
l
分别为第l阶信号的幅度、频率和初相位，假设信号阶数按照信号能量由强到弱排列，即a1≥a2≥，
…
，≥a
l
，则通常l＝1为应待估信号，其余则为多频干扰；z(n)为信号所叠加的零均值方差为σ2的加性高斯白噪声；信号总长度为n。
[0024]
首先，逐个阶数进行信号粗频率估计，充分利用补零fft、零频谱线校正、两点二分插值等技术手段尽可能获得准确的粗频率估计结果，具体按照l从1到l的顺序重复如下操作：
[0025]
(a)令s(n)＝x(n)，通过补零fft频谱峰值搜索获得信号频率估计初始值综合考虑计算复杂度和估计精度，优选地，可设置补零长度为n，此时补零信号s
2n
(n)＝[s(n),0,0,,0]长度为2n，s(m)＝fft(s
2n
(n))，由此即可搜索得到
[0026]
(b)进行零频谱线校正，计算从而得到信号复振幅
[0027][0028]
(c)采用二分插值估计获得信号粗频率估计值
[0029][0030]
其中
[0031]
(d)以更精确的信号频率再次进行零频谱线校正，即计算和
[0032][0033]
(e)逐步消除干扰信号，即
[0034][0035]
然后，基于所得信号粗频率估计值即可采取迭代干扰消除的方式逐次逐阶消除待估信号以外的多频干扰，按照两重迭代的嵌套运算结构，不断通过零频谱线校正和三点插值技术获得信号频率的高精度估计值，具体操作按照迭代次数i从1到q(q为最大迭代次数，通常可取q＝2)的顺序进行逐代计算，每次迭代计算中，均按照频率阶数l从1到l的顺序重复如下操作：
[0036]
(a)基于即得消除第l阶频率以外的所有干扰信号，即
[0037][0038]
(b)优选地，采用精度较高的三点插值法计算第i次迭代下的频率估计值
[0039][0040]
其中
[0041]
(c)再次以更精确的信号频率进行零频谱线校正即计算和
[0042][0043]
最后，完成双循环迭代运算后，即可返回信号频率的高精度估计结果
[0044]
为说明本发明算法的估计性能，本实施例通过数值仿真实验将本发明算法与其他几种典型算法进行类比分析，包括ds、mzl和music方法。不失一般性，假设信号幅度a1＝1，信号初相位θ
l
～u[-π,π]，非特别指明情况下均取n＝64，l＝2，q＝3，a2＝0.5，ω2＝ω1+δω，即仅有两个频率的信号，其中一个为干扰信号，频率间隔为δω，干扰信号幅度为待估信号幅度的一半。每次数值仿真结果均由2000次蒙特卡罗仿真结果平均后得到。
[0045]
图2给出了不同频率间隔下，频率估计均方误差mse随信噪比snr的变化关系，此时ω1＝0.3204rad。图2(a)中，δω＝2π
×
3/n，对应的频率间隔为3个dft频隙，此时所有估计方法都表现出较好的估计性能，其中music始终为有偏估计，而ds、mzl在高信噪比下性能略有下降。随着频率间隔减少为1.5个dft频隙，即δω＝2π
×
1.5n，如图2(b)中所示，此时ds和mzl已出现显著的估计偏差，music仍为有偏估计，而本发明方法仍然可以达到较高的估计精度。
[0046]
图3同样给出了不同频率间隔下，频率估计均方误差mse随信噪比snr的变化关系，但此时频率阶数为l＝4，干扰数量显著增加，此时ω1＝1.5413rad，对于l＝2,3,4有ω
l
＝ω1+(l-1)δω，a
l
＝0.5-0.1(l-2)，此时多频干扰仍然按照等间隔分布，但信号幅度逐渐减
弱。图3(a)中，δω＝2π
×
2/n，频率间隔为2个dft频隙，此时除本发明方法可以实现准确估计外，其余方法估计性能均出现一定的波动。而当频率间隔缩小为1.5个dft频隙，即δω＝2π
×
1.5/n，此时本发明方法在高信噪比下也出现了性能下降，但其估计精度仍然明显高于其他方法。
[0047]
图4给出了不同频率间隔下，频率估计均方误差mse随信号长度n的变化关系，此时ω1＝0.9173rad，snr＝30db。从图4(a)中可以看出，当δω＝2π
×
3/n时，本发明方法和mzl方法均能够实现准确的频率估计，而ds方法性能则表现出周期性波动，当频率间隔减小至δω＝2π
×
1.5/n，如图4(b)所示，仅有本发明方法能够实现准确的频率估计，其他估计方法性能显著下降。
[0048]
图5给出了不同信号长度下，频率估计均方误差mse随频率间隔δω的变化关系，此时ω1＝0.4612rad,snr＝25db。图4(a)中，n＝64，随着频率间隔δω从1个dft频隙增加至4个dft频隙，干扰逐步减弱，对比图4(b)中n＝128的情况可以看出，由于dft频隙间隔本身就是关于信号长度n的函数，因此n的变化并不会影响各种估计方法的相对性能表现，仅会影响估计方差的绝对值，而从图4(a)和(b)中均可以看出，本发明方法在频率间隔达到1.5倍dft频隙以上的情况下可以实现准确的频率估计，而其他几种估计方法则存在明显的估计误差。
[0049]
以上所述仅为本发明的优选实施例，并不用于限制本发明，显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。

技术特征：
1.一种抗多频干扰的实正弦信号频率高精度估计方法，其特征在于，所述方法步骤为：
①
对最大频率阶数为l的多频干扰信号x(n)进行逐阶循环的粗频率估计，第l阶频率下的具体操作为：首先，通过补零fft搜索谱峰以获得信号频率估计初值其次，基于进行零频谱线校正计算信号复振幅然后，利用二分插值法得到信号频率粗估计值再次，基于进行零频谱线校更新值；最后，基于和构造出第l阶信号并实现干扰消除。
②
采用迭代运算和逐阶干扰消除的双循环结构实现信号频率的高精度估计，迭代次数为i(最大迭代次数为q)且频率阶数为l时的具体操作为：首先，基于和在信号x(n)中消除第l阶频率以外的信号干扰；然后，采用三点插值法获得信号频率精估计值最后，基于进行零频谱线校正更新复振幅
③
当迭代次数达到q且遍历了所有频率阶数l后，返回信号频率估计结果2.根据权利要求1所述算法，其特征在于，所述步骤
①
和
②
中所有涉及零频谱线校正的计算步骤均为本发明独创，根据不同环节所使用的零频谱线校正，可分别将计算步骤均为本发明独创，根据不同环节所使用的零频谱线校正，可分别将或h(n)＝s(n)或y(n)带入q＝round(nω
*
/π)和以下公式计算得到信号复振幅/π)和以下公式计算得到信号复振幅该公式的等价变型包括但不限于：

技术总结
本发明提供了一种抗多频干扰的实正弦信号频率高精度估计方法，属于信息技术领域。所述算法主要针对已知频率阶数的多频干扰实正弦信号，首先，通过补零FFT频谱峰值搜索、二分插值、零频谱线校正等技术逐阶消除频率干扰并获得信号频率和复振幅的粗估计值；然后，基于已有频率和复振幅估计结果消除多频干扰，并利用三点插值、零频谱线校正等技术实现频率和复振幅估计结果的不断更新；最后，通过逐阶逐次的双循环迭代运算，获得信号频率的高精度估计结果。本发明创新提出且在计算中多次使用的零频谱线校正技术可以高效地抑制实正弦信号下的正负频率干扰，使得本发明方法能够在干扰信号与待估信号频率间隔极小的情况下获得很强的抗干扰能力，并实现高精度的频率估计，尤其在高信噪比下性能优势明显。在高信噪比下性能优势明显。在高信噪比下性能优势明显。


技术研发人员：王魁 陈辉国 毛育文 李学新 景琦 杨雨涵 张皓 黄炀
受保护的技术使用者：中国人民解放军陆军勤务学院
技术研发日：2023.05.09
技术公布日：2023/9/12