(03电机)应用FFT进行电力系统谐波分析的改进算法_庞浩

synch-ronized sampling and integral period truncation in the harmonic analysis of power system with the fast Fourier transform (FFT) technique. However, not to satisfy these conditions, the results will be disturbed by the frequency leakage. Some efforts have been made including the utilization of window functions and interpolation algorithms to correct the measured frequency, phase and amplitude by FFT. In order to overcome some shortcomings of existing correction methods, an improved algorithm is present in this paper, with which the amplitudes of harmonics can be estimated from the two neighboring spectral lines. The polynomial approximation method is also employed to obtain simple formulas for frequency and amplitude correction. By these methods, the disturbance of the frequency leakage and the noise can be reduced and the accuracy of the harmonic analysis can be improved. Based on the proposed algorithm, the practical correction formulas for some typical window functions are developed. The simulation results have verified the effectiveness and practicability of the algorithm. KEY WORDS: Power system ; FFT；Harmonic analysis ； Window functions ；Interpolation；Polynomial approximation 摘要：采用快速傅立叶变换（FFT）进行电力系统谐波分析 时很难做到同步采样和整数周期截断， 由此造成的频谱泄漏 将影响到谐波分析的结果。 通过加窗以及采用插值修正算法 可以改善计算谐波频率、相位和幅值的准确度。该文针对已 有算法存在的问题， 提出了一种基于两根谱线的加权平均来 修正幅值的双峰谱线修正算法， 利用距谐波频点最近的两根 离散频谱幅值估计出待求谐波的幅值； 同时， 利用多项式逼 近方法获得了频率和幅值修正的计算公式， 这些改进能够进 一步降低泄漏和噪声干扰，提高谐波分析的准确性。基于该 改进方法， 文中推导了一些常用窗函数的实用修正公式。仿
2 算法原理
假设一个频率为 f 0 、幅值为 A 、初相位为 θ 的 单一频率信号 x(t ) ， 在经过了采样率为 f s 的模数变 换后得到如下形式的离散信号： f x (n) = A ⋅ sin( 2π 0 n + θ ) (1) fs 如果所加窗函数的时域形式为 w( n ) ， 其连续频谱为 W ( 2πf ) ，则加窗后该信号的连续傅立叶变换为 X(f ) = x( n) ⋅ w(n ) ⋅ e − j2πfn = ∑ n = −∞
52
中
国
电 机
工 程
学 报
第 23 卷
计算 α = g −1 ( β ) 可以采用多项式逼近方法。多 项式逼近是一种近似计算复杂连续函数值的数值 方法。通过控制多项式逼近的次数，可以有效地控 制逼近的精度。而且，随着硬件乘法器在微处理器 中的广泛应用，多项式逼近的计算公式易于采用程 序代码实现。当采用切比雪夫多项式逼近奇函数
y1 W (2π ⋅ ( k1 − k0 ) N ) = y2 W (2π ⋅ ( k2 − k0 ) N )
(5)
对于给定的窗函数， 由式(5)可以计算出唯一未知量 k0，从而得到修正的峰值频率。于是，相位的修正 计算公式为 θ = arg[X (k i ⋅ ∆f ) ] + π / 2 −
arg[W (2π ⋅ ( ki − k0 ) N )] (6) 式中 arg(⋅)表示复数的主辅角。 幅值的修正公式为 Ai = 2 yi W (2π ⋅ ( ki − k0 ) N ) (7) 式(6)、(7)中的 i 可以选 1 或 2[3~6]。由于该方法仅利 用幅值最大的一条谱线进行计算，即为单峰谱线修 正算法。 单峰谱线修正算法能够在一定程度上补偿短 范围泄漏造成的影响，从而改善分析结果。但是当 选择解析形式较为复杂的窗函数时， 由式(5)难以获 得 k 0 的解析解。此外，依据式 (7) 修正幅值时，直 接利用解析表达式不仅计算复杂，还会出现小数据 相除的情况。比如，当所加窗函数是矩形窗时，需 要计算 sin( π ⋅ ξ ) sin ( π ⋅ ξ ) ⋅ N ≈ (8) π ⋅ξ sin( π ⋅ ξ / N ) 式中 ξ≤1。 当ξ 接近 0 时， 式(8)的分子和分母都趋近于 0。 对于定点微处理器而言，如果直接采用除法实现该 修正算法是无法保证计算结果的精度和正确性的。 文[5] 采用了在内存中建立查找表， 然后再利用线性 插值的方法计算频率和幅值的修正系数。该方法在 选取不同的窗函数以及需要满足不同的精度要求 时，都必须重新计算查找表，设计过程比较繁琐。 而且当精度要求提高时，查找表数据的存储量也将 成倍地增加。 由于 0 ≤ k 0 − k1 ≤ 1 ，所以可以引入一个辅助参 数 α = k 0 − k1 − 0.5 。显然，α 的数值范围是 [ −0.5, 0.5] ，它是以原点为对称的。这样，将式(5)经过变 量代换和改写后，可以得到 y 2 − y1 W (2π ⋅ (−α + 0.5) N ) − W (2π ⋅ (−α − 0.5) N ) = y 2 + y1 W (2π ⋅ ( −α + 0.5) N ) + W (2π ⋅ (−α − 0.5) N ) (9) 令 β = ( y2 − y1 ) /( y2 + y1 ) ，并且当 N 较大时，式(9) 一 般 可 以 简 化 为 β = g (α ) ， 其 反 函 数 记 为 α = g −1 ( β ) 。当窗函数 w( n ) 为实系数时，其幅频响 应 W ( 2πf ) 是偶对称的，因而函数 g (⋅) 及其反函数 g −1 (⋅) 都是奇函数。
中图分类号：TM714
应用 FFT 进行电力系统谐波分析的改进算法
庞 浩， 李东霞， 俎云霄， 王赞基
（清华大学电机系，北京 100084）
AN IMPROVED ALGORITHM FOR HARMONIC ANALYSIS OF POWER SYSTEM USING FFT TECHNIQUE
PANG Hao, LI Dong-xia, ZU Yun-xiao, WANG Zan-ji (Department of Electrical Engineering, Tsinghua University, Beijing 100084, China)
ABSTRACT: There are difficulties in performing 真结果验证了该改进算法的有效性和易实现性。 关键词： 电力系统；快速傅立叶变换；谐波分析；窗函数； 插值；多项式逼近
g (⋅) 时， 所求多项式的偶次项系数将为 0， 这样就
−1
(i = 1,2) A= N
第 23 卷 第 6 期 2003 年 6 月 文章编号：0258-8013 (2003) 06-0050-05
中
国 电 机 工 程 学 报 Proceedings of the CSEE 文献标识码：A
Vol.23 No.6 Jun. 2003 ©2003 Chin.Soc.for Elec.Eng. 学科分类号：470⋅4054
∞
ቤተ መጻሕፍቲ ባይዱ
2 π( f − f 0 ) 2 π( f + f 0 ) A jθ [e W ( ) − e − jθ W ( )] 2j fs fs (2) 如果忽略负频点−f 0 处频峰的旁瓣影响， 在正频点 f 0 附近的连续频谱函数可以表达为 A 2π( f − f 0 ) X ( f ) = e jθW ( ) (3) 2j fs 对式(3)进行离散抽样， 即可得到它的离散傅立叶变 换的表达式为 X (k ⋅ ∆f ) = ( A / 2 j)e j θ W ( 2π( k ⋅ ∆f − f 0 ) / f s ) (4) 式中 离散频率间隔为 ∆f = f s N ，N 是数据截断 长度。 峰值频率 f 0 = k 0 ⋅ ∆f 很难正好位于离散谱线 频点上，也就是说， k0 一般不是整数。设峰值点左 右两侧的谱线分别为第 k1 和 k2 条谱线， 这两条谱线 也应该是峰值点附近幅值最大和次最大的谱线。显 然，k1≤ k0≤ k2= k1+1。在离散频谱中找到这两条谱 线，从而可确定 k1 和 k2。令这两条谱线幅值分别是 y1 =| X (k1 ⋅ ∆f ) | 和 y 2 =| X (k 2 ⋅ ∆f ) | ，则由式 (4) 可 知：
第6期
庞 浩等： 应用 FFT 进行电力系统谐波分析的改进算法
51
的。但可以通过选择适当的窗函数抑制长范围泄 漏，也可以根据所选择的窗函数的形式对频率、相 位和幅值进行插值修正，在一定程度上弥补短范围 泄漏造成的误差。文[3]~[8] 给出了加窗和插值修正 算法的研究成果。由于这些方法在修正幅值时只利 用了靠近被测频点的一根最高谱线幅值进行插值， 所以文中称这些方法为单峰谱线修正算法，针对单 峰谱线修正算法在实现上的困难，本文提出了双峰 谱线修正算法并结合多项式逼近获得了更为简洁 的计算公式。仿真实验验证了双峰谱线修正算法的 正确性和易实现性。