迭代加窗插值FFT谐波分析方法

李心一;谢志江;罗久飞

【摘 要】为了提高谐波分析精度,提出了一种基于迭代加窗插值快速傅里叶变换FFT(fast Fourier transform)的谐波分析方法,并给出了统一的谐波频率、幅值及相位的计算公式.通过主瓣拟合,将传统的基于最大旁瓣衰减窗MSDW(maximum sidelobe decay window)的插值FFT方法扩展至其他对称窗,并根据窗函数的主瓣特性选择合适的窗函数进行拟合.最后通过迭代算法计算出谐波的精确频率值.仿真结果表明:在非同步采样的条件下,该算法可精确地实现谐波和间谐波分析.与传统加窗插值FFT方法相比,所提方法不依赖窗函数的类型,针对不同的窗函数具有统一的谐波参数计算公式,通用性强,实现方式灵活. 【期刊名称】《电力系统及其自动化学报》 【年(卷),期】2019(031)002 【总页数】6页(P32-37)

【关键词】谐波分析;插值;窗函数;迭代算法;快速傅里叶变换(FFT) 【作 者】李心一;谢志江;罗久飞

【作者单位】重庆大学机械工程学院,重庆 400044;重庆大学机械工程学院,重庆 400044;重庆邮电大学先进制造工程学院,重庆 400065 【正文语种】中 文 【中图分类】TM711

大量非线性负载在电力系统中的广泛使用，由此产生的谐波导致电网中的电压和电流波形畸变，从而影响电能质量，威胁电力系统的安全稳定运行[1]。快速傅里叶变换FFT（fast Fourier transform）因其于嵌入式系统中易于实现，所以是谐波分析的主要方法之一。但是由于非同步采样等原因，采用FFT进行谐波分析会产生栅栏效应和频谱泄漏，从而导致谐波参数（频率、幅值和相位）计算不准确，影响谐波分析精度。

为了提高谐波分析精度，国内外学者提出加窗插值的方法来克服FFT的缺点，如采用Hanning窗[2]、Blackman 窗[3]、Dolph-Chebyshe窗[4]、Rife-Vincent窗[5]和Nuttall窗[6]等实现加窗插值FFT谐波分析。利用加窗插值的方法可在一定程度上抑制频谱泄漏和消除栅栏效应，从而减小谐波参数的计算误差。传统的加窗插值FFT虽能有效提高谐波分析精度，但是当采用多项余弦组合窗时，因为大多数窗函数没有简单的解析表达式，所以其插值公式较复杂，需要求解高阶方程。而迭代加窗插值不依赖于窗函数的类型，具有良好的适用性[7-9]。与普通加窗插值FFT算法相比，迭代加窗插值通用性好，但由于需要进行迭代计算逼近精确值，所以该算法在某种程度上增加了计算量。得益于集成电路技术的高速发展，如今的嵌入式处理器计算速度非常快，因此增加的计算量在高速嵌入式系统中可忽略不计。 本文在分析传统加窗插值算法的基础上，结合迭代加窗插值算法求得精确的谐波频率值，并在此基础上给出谐波幅值和相位的计算公式。该方法易于实现，校正精度高，通用性好，能很好地满足电力谐波分析的要求。 1 谐波分析

对于谐波信号yh(n)，经采样频率 fs均匀采样后，得离散序列yh(n)为

式中：H为谐波信号包含的谐波数；Ah、fh0和θh分别为信号第h次谐波的幅

值、频率和初相位。信号经加对称窗wa(n)加权截断后表示为：yh-w(n)=yh(n) wa(n)。根据卷积定理，若只考虑正频部分，则加窗信号yh-w(n)经离散傅里叶变换DFT（discrete Fourier transform）后的频谱为

式中：W（k）为对称窗 wa(n)的离散时间傅里叶变换；k为DFT后离散频谱中的谱线号；kh0为谐波信号第h次谐波的频率经DFT后的谱线号。设频率分辨率Δf=fs/N，N为截断信号的数据长度，那么信号的准确频率 fh0=kh0Δf。显然，当同步采样时谐波信号的准确频率 fh0=kh0。然而，因非同步采样等原因，在实际情况中kh0往往不是整数。因此，fh0可表示为

式中：lh为峰值频点kh0附近抽样最大幅值所对应的谱线号，为正整数；qh为频率偏差，为小数。可见，得到精确频率fh0的关键是求出频率偏差值qh。但在实际应用中，求解qh往往非常复杂。文献[4]通过Matlab求解一元超越方程；文献[5-6]通过多项式拟合的方式获得求解频谱偏差的高阶方程；Belega给出一种基于最大旁瓣衰减窗MSDW（maximum sidelobe decay window）的双谱线插值方法[10-11]，该方法无需进行复杂的多项式拟合即可求得解析解。

记lh-1为lh左边谱线，lh+1为lh右边谱线。定义δh为次大谱线幅值与最大谱线幅值之比，即

对于MSDW插值，FFT频率偏差qh则表示为

式中，u为MSDW的窗项数。当频率偏差qh确定后，第h次谐波的频率为

由式（4）～式（6）可看出，基于MSDW的插值公式比较简单，但是Daniel并

没对非MSDW的插值FFT方法进行讨论。文献[9]通过主瓣拟合技术，将基于汉宁窗（2-MSDW）插值FFT法扩展至其他经典窗。通过文献[2-6]分析可知，采用具有较快旁瓣衰减速率的窗函数进行加窗插值FFT能带来较高的计算精度，但此类窗函数往往具有较大的主瓣宽度（半主瓣宽度往往大于3个谱线间隔）。文献[9]中，当选用此类窗函数加窗插值时迭代计算收敛速度较慢。为了解决这一问题，本文将主瓣拟合和3项MSDW插值FFT相结合，在一定程度上可以减少迭代次数，提高计算速度。对于加任意对称窗的谐波信号，定义傅里叶系数为

式中，γ=±0.5。对于任意经典窗，在k∈ [-0.5，0.5]内[12]有

式中：WR(k)和W3MSDW(k)分别为任意窗和3-MSDW的归一化窗谱；参数η3MSDW=SL3MSDWSLw(n)，其中SL表示窗函数的最大扇形损失，为DFT谱线对应半频率与精确DFT谱线频率的比值[13]，计算公式为

式中，W（·）为窗函数w（n）的离散时间傅里叶变换。式（8）表明，在主瓣中间附近可以用替代W3MSDW(k)，从而可以把式（5）扩展应用于任意对称窗函数。因此为了加速迭代收敛速度，式（4）改写为

进而可以得到改进后的频率估计迭代算法，即

式中：迭代初始值l0=lh；l1=l0+qh。当|lj-lj+1|＜υ（υ的取值取决于所需的计算精度）时，停止迭代，由此可以得到精确的谐波频率 fh0=l(j)hΔf。求得精确的l(j)h后，可得到第h次谐波信号的准确幅值和相位，分别为

式中，R为窗函数的幅值恢复系数[14-15]。另外，对于具有不同主瓣宽度的窗函数，式（11）可以进一步统一为

式中，u为窗函数选择参数。当加窗插值选用窗函数主瓣较窄时（半主瓣宽度小于1.5个谱线间隔），采用矩形窗为被拟合窗函数，即u=1，式（14）为文献[8]中的迭代算法。同理，当加窗插值选用窗函数半主瓣宽度为2个谱线间隔左右时，即u=2，为文献[9]中的迭代算法。对于主瓣较宽，旁瓣衰减速率快的多项窗，其半主瓣宽度大多大于3个谱线间隔，本文采用3项MSDW为被拟合窗函数，故u=3。值得注意的是，拟合不同窗函数时，参数η会发生相应变化。 2 仿真分析与比较 2.1 谐波分析

利用文献[5]给出的信号模型进行仿真，则有

式中：f1为基波频率，取50.2 Hz；fs=2 500 Hz；N=1 024。本文算法与加5项Rife-Vincent（I）窗和4项3阶Nuttall窗双谱线插值FFT进行对比，算法测量结果的绝对误差如图1所示。 2.2 间谐波分析

间谐波是指非基波频率的整数倍谐波，间谐波会引起电压闪变和导致谐波补偿装置失效，因此精确的测量间谐波参数是十分必要的。以文献[16]给出的信号模型为例进行间谐波仿真。采样频率fs=1 250 Hz；截断信号的数据长度N=1 024，仿真结果如图2所示。

图1 谐波参数误差曲线Fig.1 Error curves of harmonic parameters

从仿真结果来看，对于间谐波次数（1，4，6，8，11），在大多数情况下，基于本文方法的间谐波参数测量精度都很高，其中基于非对称Rife-Vincent（I）5窗

的间谐波频率绝对误差最小可达10-11。 2.3 其他经典窗

在加窗插值FFT中通过频率偏差qh与谱线的幅值比值之间的关系可以求得qh，从而可以得到谐波的频率精确值。谱线的幅值比值又与窗函数关系密切，因此在用加窗插值FFT进行谐波分析时选择合适的窗函数对分析精度有较大的影响。信号同样采用文献[16]的间谐波信号模型，采样频率fs=1 250 Hz，采样点N=1 024。图3所示为利用常用窗函数进行谐波参数估计的结果。从图中可以看出，对于其他经典窗，本文算法同样具有良好的精度和适用性。值得注意的是，由于5-MSDW具有最快的旁瓣衰减速率。因此，基于5-MSDW的插值计算精度更高。 图2 间谐波参数误差曲线Fig.2 Error curves of inter-harmonic parameters 2.4 基波频率波动的谐波分析

图3 其他经典窗间谐波参数误差曲线Fig.3 Error curves of inter-harmonic parameters with other classical windows

在谐波分析中，基波频率波动会影响谐波参数估计的结果。对于式（15）的复杂谐波信号，基波频率在 50～51 Hz之间波动，fs=2 500 Hz，N=1 024。本文算法与基于Nuttall窗和RV（I）-5窗插值算法的结果如图4～图7所示。 图4 基波频率波动时幅值绝对误差（Nuttall窗）Fig.4 Amplitude absolute errors when the fundamental frequencey fluctuates（Nuttall window） 图5 基波频率波动时相位绝对误差（Nuttall窗）Fig.5 Phase absolute errors when the fundamental frequency fluctuates（Nuttall window）

图6 基波频率波动时幅值绝对误差（RV（I）-5窗）Fig.6 Amplitude absolute errors of the fundamental frequencey fluctuates（RV（I）-5 window） 图7 基波频率波动时相位绝对误差（RV（I）-5窗）Fig.7 Phase absolute errors of the fundamental frequency fluctuates（RV（I）-5 window）

由图4～图7可见，对于基波频率波动的情况，基于本文方法的谐波幅值和相位受基波频率波动的影响较小，误差较为平缓。尤其是对于幅值，与另外两种方法相比，本文方法的误差曲线非常平缓。由此可见，本文方法能有效抑制基波频率波动对谐波参数估计的影响。 3 结论

（1）本文在分析基于MSDW插值FFT方法的基础上，扩展了基于迭代加窗插值FFT谐波分析方法，并给出了谐波参数计算公式。在传统的加窗插值FFT谐波分析算法中，针对不同的窗函数频率偏差qh计算公式无法统一，尤其是对于主瓣较宽，旁瓣衰减速率快，窗项数较多的窗函数，虽然具有良好的频谱泄漏抑制能力，但其大多数没有简单的解析表达式，因此需要进行多项式拟合来求得qh。Belega等在文献[10-11]中给出了基于MSDW插值的频率、幅值和相位的计算公式，其公式较为简单，但对非MSDW并没有给出统一的计算公式。

（2）与传统加窗插值FFT算法相比，本文方法具有更好的通用性，针对不同的窗函数可以选用不同的主瓣拟合方法，并且有统一的迭代和谐波参数计算公式。虽然需要迭代来求得精确值，在一定程度上增加了计算负担，但若选用合适的拟合窗函数，仅需要几步迭代即可求得精确值，在高速的嵌入式系统中影响并不大。 （3）通过理论分析和仿真结果可以看出，在大多数情况下，本文方法都能取得较高的计算精度，并且实现方式灵活，适合嵌入式系统的应用。

本文方法最大的特点是适用性强，因此非常适合工程实际应用。虽然采用3-MSDW代替汉宁窗在一定程度上能提高计算速度，但如何进一步加快迭代求解速度则还有待于研究。

【相关文献】

[1] 曾巧燕，杨洪耕，王泽，等（Zeng Qiaoyan，Yang Honggeng，Wang Ze，et al）.全相位谱法在闪变检测中的应用（Application of all-phase spectrum method to flicker detection）[J].电力系统及其自动化学报（Proceedings of the CSU-EPSA），2017，29（4）：66-71. [2] Zhang Fusheng，Geng Zhongxing，Yuan Wei.The algorithm of interpolating windowed FFT for harmonic analysis of electric power system[J].IEEE Trans on Power Delivery，2001，16（2）：160-1.

[3] 许珉，杨阳，章梦哲，等（Xu Min，Yang Yang，Zhang Mengzhe，et al）.一种加三项余弦窗的加窗插值FFT算法（An interpolated FFT algorithm based on three-item cosine window）[J].电力系统保护与控制（Power System Protection and Control），2010，38（17）：11-15. [4] 孙仲民，黄俊，杨健维，等（Sun Zhongmin，Huang Jun，Yang Jianwei，et al）.基于切比雪夫窗的电力系统谐波/间谐波高精度分析方法（A high accuracy method for harmonic and inter-harmonic in power systems based on Dolph-Chebyshev windows）[J].电力系统自动化（Automation of Electric Power Systems），2015，39（7）：117-123.

[5] 曾博，滕召胜，高云鹏，等（Zeng Bo，Teng Zhaosheng，Gao Yunpeng，et al）.基于Rife-Vincent窗的高准确度电力谐波相量计算方法（An accurate approach for power harmonic phasor calculation based on Rife-Vincent window）[J].电工技术学报（Transactions of China Electrotechnical Society），2009，24（8）：154-159.

[6] 卿柏元，滕召胜，高云鹏，等（Qing Boyuan，Teng Zhaosheng，Gao Yunpeng，et al）.基于Nuttall窗双谱线插值FFT的电力谐波分析方法（An approach for electrical harmonic analysis bases on Nuttall window double-spectrum-line interpolation FFT）[J].中国电机工程学报（Proceedings of the CSEE），2008，28（25）：153-158.

[7] Liu Yanhui，Nie Zaiping，Zhao Zhiqin，et al.Generalization of iterative Fourier interpolation algorithms for single frequency estimation[J].Digital Signal Processing，2011，21（1）：141-149.

[8] Aboutanios E，Mulgrew B.Iterative frequency estimation by interpolation on Fourier coefficients[J].IEEE Trans on Signal Processing，2005，53（4）：1237-1242.

[9] Luo Jiufei，Xie Zhijiang，Xie Ming.Frequency estimation of the weighted real tones or resolved multiple tones by iterative interpolation DFT algorithm[J].Digital Signal Processing，2015，41：118-129.

[10]Belega D，Dallet D.Frequency estimation via weighted multipoint interpolated DFT[J].IET Science，Measurement&Technology，2008，2（1）：1-8.

[11]Belega D，Dallet D.Multifrequency signal analysis by Interpolated DFT method with maximum sidelobe decay windows[J].Measurement，2009，42（3）：420-426.

[12]Luo Jiufei，Xie Zhijiang，Xie Ming.Interpolated DFT algorithms with zero padding for classic windows[J].Mechanical Systems and Signal Processing，2016，s70-71：1011-1025. [13] 谢长贵，陈平（Xie Changgui，Chen Ping）.一种DFT的频率校正算法（DFT

interpolation algorithm for frequency estimation）[J].云南大学学报：自然科学版（Journal

of Yunnan University：Natural Science），2017，39（2）：200-206.

[14] 焦新涛，丁康（Jiao Xintao，Ding Kang）.加窗频谱分析的恢复系数及其求法（Resetting module and solutions in windowing spectrum analysis）[J].汕头大学学报：自然科学版（Journal of Shantou University：Natural Science），2003，18（3）：26-30，38. [15]Harris F J.On the use of windows for harmonic analysis with the discrete Fourier transform[J].Proceedings of the IEEE，1978，66（1）：51-83.

[16]Qian Hao，Zhao Rongxiang，Chen Tong.Interharmonics analysis based on

interpolating windowed FFT algorithm[J].IEEE Trans on Power Delivery，2007，22（2）：10-1069.