矩形窗三点插值傅里叶变换高精度频率估计方法.pdf

引用格式：引用格式：熊德智,肖宇,胡军华,等.矩形窗三点插值傅里叶变换高精度频率估计方法J.中国测试，2023,49(9):57-62.XIONGDezhi,XIAO Yu,HU Junhua,et al.Rectangular window three-point interpolation Fourier transform high-precision frequencyestimationmethodJ.ChinaMeasurement&Test,2023,49(9):57-62.DOI:10.11857/j.issn.1674-5124.2021120030矩形窗三点插值傅里叶变换高精度频率估计方法熊德智1,2，肖宇1,2，胡军华1,2，黄瑞1,2，温和2,3，温乐鹏2,3(1.国网湖南省电力有限公司，湖南长沙410004;2.智能电气量测与应用技术湖南省重点实验室，湖南长沙410004;3.湖南大学电气与信息工程学院，湖南长沙410082)摘要:针对电力系统信号频率估计的高准确性要求，该文提出一种基于矩形窗的三点复频域插值频率估计方法。所提出的方法首先采用矩形窗对信号加权处理，然后进行离散傅里叶变换，最后选择离散频谱中幅度最大的三根谱线进行复频域插值，得到频率估计结果。该方法同时考虑离散频谱中的正、负频谱，消除短程和长程频谱泄漏对频率估计误差的影响。采用主瓣最窄的矩形窗进行加窗，能最大程度地减少白噪声对频率估计误差的影响。仿真和实验结果表明：在不同周期和噪声强度情况下，所提出的方法均可以实现对信号频率的准确估计。在与已有的加窗插值傅里叶变换方法相比，该文所提出的插值算法抗噪性更好、频率估计误差更小，适用于对电力系统信号频率的准确估计。关键词:频率估计;傅里叶变换;频谱插值;窗函数中图分类号:TM933文献标志码:A文章编号:16745124(2023)09005706Rectangular window three-point interpolation Fourier transform high-precisionfrequency estimation methodXIONGDezhi1,2,XIAOYu1,2,HUJunhua1,2,HUANGRui1,2,WENHe2,3,WENLepeng2,3(1.StateGridHunanElectricPowerCorporationLimited,Changsha410004,China;2.HunanProvinceKeyLaboratoryofIntelligentElectricalMeasurementandApplicationTechnology,Changsha410004,China;3.CollegeofElectricalandInformationEngineering,HunanUniversity,Changsha410082,China)Abstract:Aimingatthehighaccuracyrequirementofpowersystemsignalfrequencyestimation,thispaperproposes a three-point complex frequency domain interpolation frequency estimation method based onrectangularwindow.Theproposedmethodfirstusesarectangularwindowtoweightthesignal,thenperformsdiscreteFouriertransform,andfinallyselectsthethreespectrallineswiththelargestamplitudeinthediscretespectrumforcomplexfrequencydomaininterpolationtoobtainthefrequencyestimationresult.Thismethodsimultaneouslyconsidersthepositiveandnegativespectruminthediscretespectrumtoeliminatetheinfluenceofshort-rangeandlong-rangespectralleakageonthefrequencyestimationerror.Therectangularwindowwiththenarrowestmainlobeisusedforwindowing,whichcanminimizetheinfluenceofwhitenoiseonthefrequencyestimationerror.Simulationandexperimentalresultsshowthattheproposedmethodcanachieveaccurate estimation of signal frequency under different periods and noise intensities.Compared with the收稿日期:2021-12-06；收到修改稿日期:2022-02-09作者简介:熊德智（1979-），男，湖南长沙市人，高级工程师，主要从事电能计量与电气测量研究。通信作者:温和（1982-），男，湖南长沙市人，教授，主要从事智能电气与量测应用技术的研究。第49卷第9期中国测试Vol.49No.92023年9月CHINAMEASUREMENT&TESTSeptember,2023existingFouriertransformmethodofwindowedinterpolation,theinterpolationalgorithmproposedinthispaperhasbetteranti-noiseandsmallerfrequencyestimationerror.Soitissuitableforaccurateestimationofpowersystemsignalfrequency.Keywords:frequencyestimation;Fouriertransform;spectralinterpolation;windowfunction0 引言电力系统中频率是电能质量的重要参数，发电和负载之间动态不平衡会导致信号频率发生变化，准确的频率估计是电网稳定和电气设备正常运行的先决条件1-2。因此，从混有白噪声的正弦波中准确估计信号频率对电力系统是十分必要的。为了实现对信号频率的准确估计，国内外学者提出了很多改进或全新的频率估计算法，根据基本原理，主要分为时域和频域两类。基于时域的算法利用连续的时域信号值进行计算，典型算法主要有Prony3、MUSIC4、ESPRIT5等。此类算法只有在已知模型阶数的条件下才可实现对频率的准确估计，并且存在大量复杂的矩阵运算，计算复杂度限制了此类算法的在线实时应用。基于频域的算法主要有快速傅里叶变换（fastFouriertransform，FFT），其具有较好的计算速度，得到了广泛应用。但是，由于电力系统频率往往存在波动，非同步采样下的离散傅里叶变换（discreteFouriertransform，DFT）存在固有的频谱泄漏和栅栏效应，导致频率估计存在较大误差。值得注意的是，DFT 的频谱泄漏可分为由正频率引起的短程泄漏和由负频率引起的长程泄漏。为了提高频率估计精度，通常采用时域加窗和频域插值的方法来减少频谱泄漏和栅栏效应的影响6-12。加窗插值 DFT（interpolatedDFT，IpDFT）被广泛应用于实正弦信号的频率估计。两点 IpDFT 算法采样最大和次大两根峰值谱线，根据峰值谱线的幅度比值计算频率。文献 6 提出一种基于最大旁瓣衰减窗（maximumsidelobedecaywindow，MDW）的两点 IpDFT，在不考虑负频率的情况下，提取两根最大谱线进行频率估计。三点 IpDFT 算法使用三个最大峰值谱线进行频率估计。文献 7 利用最大旁瓣衰减窗口推导出三点 IpDFT 频率估计公式，实验证明在白噪声干扰下，能够准确估计信号频率。但是，这些算法都忽略了来自负频率的长程泄漏，在获取信号周期数较少时，频率估计的精度不够高。最近，已有学者提出同时考虑短程和长程泄漏的IpDFT 方法。文献8 基于MDW 的三点IpDFT 算法，降低了短程和长程泄漏对频率估计的影响。文献9 利用短程和长程泄漏的校正因子来提高算法的估计精度。文献 10 考虑负频率对频率估计的影响，提出一种频率估计校正模型，提高了频谱分析准确性。为了进一步提高频率估计精度，本文提出一种基于矩形窗的三点复频域插值方法，实现较短观测时间窗口的频率准确估计。该方法采用矩形窗对被采样信号进行截断，选择傅里叶变换后离散频谱中幅度最大的三根谱线进行复频域插值，以消除短程和长程泄漏对频率估计误差的影响。因矩形窗能具有最窄等效噪声增益带宽，可最大程度地减少白噪声对频率估计误差的影响。仿真和实验结果表明，所提出的方法可以实现准确的频率估计。1 设计方案假设以采样频率 fs对实正弦信号 y(t)进行采样，则 N 点离散采样序列 y(n)可表示为：y(n)=x(n)+e(n)=Acos(2f0n/fs+)+e(n)(1)式中：A、0=2f0/fs和 离散信号的幅值、频率和初始相位；e(n)方差为 2的加性高斯白噪声；f0连续时间信号的频率；fs信号的采样频率。用一个长度为 N 的矩形窗 w(n)对信号进行加窗，得到长度为 N 的加窗序列 s(n)：s(n)=y(n)w(n)(2)n=0,1,2,N1其中。s(n)的 N 点 DFT 可表示为：Y(k)=X(k)+E(k)=N1n=0Acos(0n+)WknN+E(k)WkN=ej(2Nk)0=2l0/N=2(k0+0)/N(3)式中：X(k)和 E(k)离散信号 x(n)和噪声 e(n)的DFT；l0正弦信号的周期数；0,1,N/210.5,0.5k0和 0归一化频率的整数部分和小数部分，且 k0和 0的取值范围分别为和。58中国测试2023年9月通过欧拉公式，X(k)可表示为：X(k)=12Aej1ejN01ej(2Nk0)+12Aej1ejN01ej(2Nk+0)(4)2(k0+0)/Nk0=argmax|X(k)|由式(4)可以看出，实正弦信号频谱包括，根据给定序列通过 X(k)估计频率0，可以通过寻找 DFT 最大谱峰值的位置实现粗略估计，即。v=Aej(1ejN)/2v=Aej(1ejN)/2本文采用三点复频域插值 DFT 消除短程和长程频谱泄漏，进而实现精确估计信号频率。定义，则 式(4)可写为：X(k)=v1WkN+v1WkN=(v+v)(v+v)WkN1(+)WkN+(WkN)2(5)=ej0+=2cos0=1其中，根据欧拉公式可得和。假设已知离散信号任意两点的 X(k1)和X(k2)，则通过式(5)可得如下线性方程组：X(k1)12cos(0)Wk1N+(Wk1N)2=(v+v)(v+v)Wk1NX(k2)12cos(0)Wk2N+(Wk2N)2=(v+v)(v+v)Wk2N(6)同理，任取一个不同于 X(k1)和 X(k2)的 DFT值 X(k3)，与 X(k2)构成一个新的方程，如下：v+v=Wk2NX(k3)Wk2NWk3N(12cos0Wk3N+(Wk2N)2)Wk3NX(k2)Wk3NWk3N(12cos0Wk2N+(Wk2N)2)(7)根据式(6)和式(7)，能够推导关于 0的求解方程：1Wk2NWk1NWk3NX(k1)(12cos0Wk1N+(Wk1N)2)1Wk2NWk3NWk2NX(k2)(12cos0Wk2N+(Wk2N)2)=1Wk2NWk3NWk2NX(k3)(12cos0Wk3N+(Wk3N)2)1Wk2NWk3NWk3NX(k2)(12cos0Wk2N+(Wk2N)2)(8)0在实际应用中，采样信号往往会被噪声干扰，从而使 X(ki)成为 Y(ki)的采样值，定义通过 Y(k)序列估计的信号频率为，频率估计公式为：0=arccos a=arccos(1Y(k3)2Y(k2)3Y(k1)+4Y(k2)1(Y(k2)Y(k1)2(Y(k2)Y(k3)(9)aaa其中 是 的估计值，真实值 由下式得出：a=1X(k3)2X(k2)3X(k1)+4X(k2)1(X(k2)X(k1)2(X(k2)X(k3)(10)式(9)中各系数计算如下：1=(Wk2NWk1N)Wk2N(1+(Wk3N)2)2=(Wk2NWk1N)Wk3N(1+(Wk2N)2)3=(Wk2NWk3N)Wk2N(1+(Wk1N)2)4=(Wk2NWk3N)Wk1N(1+(Wk2N)2)(11)1=2(Wk2NWk3N)Wk1NWk2N2=2(Wk2NWk1N)Wk2NWk3N(12)由式(9)可知，本文所提出的频率估计方法，可以通过任意三点 DFT 值实现。在实际应用中，考虑到噪声对信号的干扰，优先选择幅值最大的谱线 Y(k0)及其相邻谱线 Y(k01)和 Y(k0+1)，从而减小频谱泄漏对 DFT 样本值造成的影响。2 噪声对频率估计的影响本节推导所提出算法的理论方差，从而用来评估噪声对所提出的频率估计算法的影响。a a首先分析噪声对 a 的影响，把式(3)代入 的表达式，则 公式可改写为：a=1X(k3)2X(k2)3X(k1)+4X(k2)+11(X(k2)X(k1)2(X(k2)X(k3)+2(13)式中：1=1Q(k3)2Q(k2)3Q(k1)+4Q(k2)2=1(Q(k2)Q(k1)2(Q(k2)Q(k3)(14)a=a+用 表示噪声引起的估计误差，则，的计算公式为：=11+22(15)式中：1=11(X(k2)X(k1)2(X(k2)X(k3)2=1X(k3)2X(k2)3X(k1)+4X(k2)(1(X(k2)X(k1)2(X(k2)X(k3)2(16)第49卷第9期熊德智，等：矩形窗三点插值傅里叶变换高精度频率估计方法59cos1 acos1a然后，利用泰勒级数将按照展开，如下式所示：=01cos20+o1cos20201cos20(17)E()=0需要注意的是，由于噪声 Q(k)的均值为 0，所以式中 是无偏的，即。经过理论推导，可得到方差计算公式：var(0)=E(00)2E2(00)=1sin20E(2)=N2sin202121+23+(42)2+22(12)2+21+22+21212+31+(42)(12)(18)3 仿真分析和实验验证3.1 所提出算法性能(0,20,1000.5,4.5为了评估噪声干扰下所提出算法性能，首先要验证解析式(9)和式(18)的有效性。本文通过设置不同 l0和 SNR 值进行仿真来评估所提出算法的性能。在仿真实验中，设置采样点数N=512，采样频率fs=1024Hz，信号幅值 A=1，相位 在范围随机变化。图 1 为在 l0=1.8 的情况下，SNR 在dB范围内逐渐变化时的仿真和理论结果。图 2 为在SNR=50dB 的情况下，周期数 l0在范围内逐渐变化时的仿真和理论结果。0102030405060708090 100SNR/dB1801601401201008060方差Practicetheory图 1 l0=1.8 时，方差随 SNR 变化如图 1 所示，所提出算法频率估计方差在 SNR变化过程中，实验方差和理论方差几乎重合，随着SNR 的增加，方差逐步降低。在图 2 中，实验方差和理论方差也十分吻合。所提出算法频率估计方差在l0变化过程中一直很稳定，仿真结果也很接近 CRLB。从图 2 可以看出，在 SNR=50dB 下，l0=1 时频率估计的方差约为115，这表明所提出的算法在一个周期内就可以实现准确稳定的频率估计。3.2 与其他算法性能对比(0,2在本节中，分析了白噪声和频率变化对频率估计精度的影响。在实验中，初始相位在范围内随机选择。对于每个 SNR 值，仿真 1000 次以评估统计特性。方差用于评估所提出的算法和其他算法的性能，计算公式如下：var(0)=10lg1NNi=1(0(i)0)2(19)0(i)其中代表第 i 次独立实验的频率估计值。本节将所提出的算法与其他基于 DFT 的频率估计算法进行比较，这些算法有基于不同加窗（Minl3,BH4）的两点 IpDFT 算法13、基于不同加窗（MDW2,MDW3,MDW4,Blackman）的三点 IpDFT 算法14。各窗函数系数如表 1 所示。表 1 窗函数系数Windowa0a1a2a3HMDW20.500000.500002MDW30.3750000.500000.125003MDW40.312500.468750.187500.0312503Blackman0.420000.500000.080003Minl30.42438010.49734060.07827933BH40.402170.497030.0989200.001884图 3 展示了信号加入 SNR=50dB 噪声后，周期数 l0从 0.54.5 变化时，所提出算法与其他算法频率估计的方差。从图中可以看出，在 l0整个变化范0.51.01.52.02.53.03.54.04.5l0/cycle1301201101009080706050方差Practicetheory图 2 SNR=50 dB 时，方差随 l0的变化60中国测试2023年9月围内，所提出的算法均保持稳定且具有很小的方差，方差基本稳定在115 水平，非常接近 CRLB。基于Minl3、BH4 加窗的 I2pDFT 算法在 l0较小时还能保持很高的估计精度，在随着 l0不断增加，算法的方差会出现较大的波动。基于 MDW2、MDW3、MDW4 加窗的 I3pDFT 算法的方差虽然不会出现明显的波动，但整体上性能劣于所提出的算法。0.51.01.52.02.53.03.54.04.5l0/cycle130120110100908070506040方差ProposedMDW2BlackmanMDW3MDW4I2pDFTMinl3I2pDFTBH4CLRB图 3 SNR=50 dB 时，方差随 l0变化情况图 4 为在所获取周期数 l0=1.25，SNR 从 0100dB 变化时，所提出的算法与其他算法频率估计的方差。从图可以看出，在 SNR 整个变化范围内，所提出的算法频率估计精度明显优于其他算法，并且频率估计方差随着 SNR 线性变化，随着 SNR 的增加，频率估计精度进一步提高。在所有对比算法中，所提出的频率估计算法性能是最好的。0102030405060807090 100SNR/dB180160140120100608040方差ProposedMDW2BlackmanMDW3MDW4I2pDFTMinl3I2pDFTBH4CLRB图 4 l0=1.25 时，方差随 SNR 变化情况3.3 实测信号分析本文分别用所提出的算法和经典两点插值与三点插值算法，分析福禄克源实测的信号数据，对比结果验证所提出算法的鲁棒性。实测信号如图 5 所示，福禄克源采样频率为8000Hz，信号主频为50.5Hz，包含微弱的奇数次谐波分量，采样点数设置为N=512。050 100 150 200 250 300400350450 500N420311234幅值/V图 5 实测信号为了验证所提出算法的鲁棒性，设计实验对实测信号分别加入 0100dB 的白噪声，分析白噪声对所提出算法性能的影响。在实验中，SNR 值以 5dB递增，在每个 SNR 下，进行 1000 次循环以增加结果的统计特性。对比算法为经典的两点 IpDFT 和三点 IpDFT，实验结果如图 6 所示。0102030405060807090 100SNR/dB180140100160120806020400方差ProposedI2pDFTI3pDFTCLRB图 6 不同 SNR 下方差变化情况图 6 为所提出算法和对比算法，在实测信号加入了 0100dB 白噪声时，频率估计精度的变化情况。由图可以看出，所提出算法频率估计方差在 SNR 变化过程中，随着 SNR 的增加，方差呈线性降低，并且所提出算法的频率估计精度始终优于对比算法。4 结束语本文提出了一种新的频率估计算法，采用三根谱线进行频域插值实现频率估计。所提出的算法考虑了正负频率的频谱叠加，克服了长短频程泄漏的影响，频率估计精度高。本文提出的方法采用矩形窗，无需额外加窗运算，且因矩形窗具有最小的等第49卷第9期熊德智，等：矩形窗三点插值傅里叶变换高精度频率估计方法61效噪声带宽增益，可确保所提出的方法的抗噪性。仿真和实验结果表明，与已有加窗插值 DFT 方法相比，所提出的算法稳定性高，频率估计误差小。参考文献WENH,ZHANGJH,MENGZ,etal.HarmonicestimationusingsymmetricalinterpolationFFTbasedontriangularself-convolution windowJ.IEEE Transactions on IndustrialInformatics,2015,11(1):16-26.1舒骁骁,祝君剑,朱亮,等.计及噪声影响的高准确度迭代滤波电网频率测量方法J.中国测试,2020,46(7):54-59,132.2HAUERJF,DEMEURECJ.InitialresultsinPronyanalysisofpowersystemresponsesignalsJ.IEEETransactionsonPowerSystems,1990,5(1):80-89.3SCHMIDTR.MultipleemitterlocationandsignalparameterestimationJ.IEEETransactionsonAntennasandPropagation,1986,34(3):276-280.4WANG X C,TANG F,WANG X R,et al.Estimation ofelectromechanicalmodesunderambientconditionviarandomdecrement technique and TLS-ESPRIT algorithmC/2014International Conference on Power System Technology(Powercon),2014.5BELEGAD,PETRID.Accuracyanalysisofthemulticyclesynchrophasor estimator provided by the interpolated DFTalgorithmJ.IEEE Transactions on Instrumentation andMeasurement,2013,62(5):942-953.6WENH,LICC,TANGL.Novelthree-pointinterpolation7DFT method for frequency measurement of sine-waveJ.IEEE Transactions on Industrial Informatics,2017,13(5):2333-2338.BORKOWSKIJ,KANIAD,MROCZKAJ.Interpolated-DFT-basedfastandaccuratefrequencyestimationforthecontrolofpowerJ.IEEETransactionsonIndustrialElectronics,2014,61(12):7026-7034.8BELEGAD,PETRID,DALLETD.Frequencyestimationofasinusoidalsignalviaathree-pointinterpolatedDFTmethodwithhighimagecomponentinterferencerejectioncapabilityJ.DigitalSignalProcessing,2014,24:162-169.9张金平,李建立,段晨.计及负频率影响的新能源发电低频间谐波检测方法J.电测与仪表,2020,57(2):95-100.10李恺,向鑫,卜文彬,等.基于频谱分辨率自适应的双插值DFT 谐波分析方法J.电测与仪表,2021,58(5):92-97.11陈向群,高云鹏,李典卿,等.基于 Blackman-Nuttall 窗改进FFT 校正的新型三相动态谐波电能表设计J.电测与仪表,2020,57(9):132-139.12BELEGAD,DALLETD,PETRID.Statisticaldescriptionofthesine-wavefrequencyestimatorprovidedbytheinterpolatedDFTmethodJ.Measurement,2012,45(1):109-117.13WEN H,LI C C,YAO W X.Power system frequencyestimation of sine-wave corrupted with noise by windowedthree-pointinterpolatedDFTJ.IEEETransactionsonSmartGrid,2018,9(5):5163-5172.14(编辑:商丹丹)(上接第 30 页)参考文献国务院印发计量发展规划(2021-2035 年)鼓舞全国计量工作者 J.中国计量,2022(3):16-29.1方向.全面提升计量科技水平推进中国式现代化进程J.中国计量,2023(5):6-7.2赵卢雷,沈伯平.新时代背景下实施创新驱动发展战略的若干思考J.改革与战略,2020,36(2):68-79.3周洪宇,李宇阳.新时代新征程实施科教兴国战略:新内涵、重大意义与推进路径J.人民教育,2022(22):24-27.4姚冠男.中国式教育现代化之路探讨基于十八大以来教育改革实践及二十大报告相关论述J.内蒙古师范大学学报(教育科学版),2023,36(3):1-10.5国务院关于印发计量发展规划(20212035 年)的通知 J.中华人民共和国国务院公报,2022(5):52-62.6铁路计量发展规划(20212035 年)J.铁道技术监督,2022,50(5):2-5+10.7杨林,李敏.改革开放以来中国企业技术创新主体地位演变历程、机理及态势J.江南大学学报(人文社会科学版),2023,22(3):27-43.8湖南省人民政府关于推进计量事业高质量发展的实施意见J.湖南省人民政府公报,2022(20):14-19.9刘国传.深入学习贯彻党的十九届六中全会精神大力实施计量优先发展战略J.中国计量,2022(4):28-29.10加快构建国家现代先进测量体系国务院印发计量发展规划(2021-2035 年)J.中国质量监管,2022(2):10.11杨奕,郭瑞萍.高质量发展促进共同富裕的逻辑结构与实现路径J.中共石家庄市委党校学报,2023,25(3):29-33.12王战军,赵敏.新时代我国研究生教育的新使命、新举措J.现代教育管理,2023(4):44-53.13方向.全面开创计量科技创新战略发展新局面J.中国计量,2022(12):16-17.142022 年全国计量工作要点 J.中国计量,2022(3):30-34.15市场监管总局全力推进计量“硬能力”建设夯实经济社会高质量发展“硬支撑”J.中国计量,2022(5):14-15.16秦宜智.实施计量发展规划(2021-2035 年)全面开启加快计量发展新征程J.中国计量,2022(6):5-9.17陈丹丹,尚晓旭,张晶.构建质量基础设施协同发展生态体系J.中国质量与标准导报,2022(6):63-66.18陈岳飞,徐学林,朱美娜等.国家现代先进测量体系的基本内涵、生成逻辑及实现路径J.中国测试,2021,47(2):1-5.19(编辑:徐柳)62中国测试2023年9月