第45卷第1期电力系统保护与控制V ol.45 No.1 2017年1月1日Power System Protection and Control Jan. 1, 2017 DOI: 10.7667/PSPC151779
基于四谱线插值FFT的谐波分析快速算法
张俊敏1，刘开培2，汪 立3，陈文娟2
(1.中南民族大学计算机科学学院，湖北 武汉 430074; 2.武汉大学电气工程学院，湖北 武汉 430072;
3.国网天津市电力公司，天津 300000)
摘要：由于非同步采样和非整周期截断导致快速傅里叶变换(FFT)不能准确地分析出谐波参数，加窗和插值算法经常被用来改善FFT的计算精度。在信号加窗条件下，基于四谱线插值的FFT算法基础上进行快速算法研究。
该算法通过分析加窗后信号的频域表达式，利用真实谐波点附近的4根最大谱线值确定实际谱线的位置，对该次谐波进行频率、幅值和相位等参数估计。并且通过多项式拟合的方式推导出了4种典型窗函数的修正公式。根据窗函数主瓣内任意相邻谱线相位相差π的规律，提出一种快速算法，计算某次谐波开方计算量仅需要1次，大大节约了计算复杂度和计算时间。仿真实验表明，四谱线插值算法在拟合阶次较低的情况下，不仅可以获得比常用双谱线和三谱线更高的精度，还具有对偶次谐波检测精度远胜于
双、三谱线插值算法的优点。
关键词：谐波分析；窗函数；快速傅里叶变换；四谱线；插值
A rapid algorithm for harmonic analysis based on four-spectrum-line interpolation FFT
ZHANG Junmin1, LIU Kaipei2, WANG Li3, CHEN Wenjuan2
(1. College of Computer Science, South-Central University for Nationalities, Wuhan 430074, China; 2. College of Electrical
Engineering, Wuhan 430072, China; 3. State Grid Tianjin Electric Power Company, Tianjin 300000, China)
Abstract: In the case of non-synchronous sampling and in non-integral period truncation, the fast Fourier transform (FFT) can’t analyze harmonic parameters correctly. So the window functions and interpolation algorithms are made to improve the accuracy of harmonic parameter computation by FFT. This paper derives an algorithm based on four-spectrum-line interpolation. This algorithm uses the four most spectrum lines around true location to computer the accuracy position of the harmonic spectrum line through analyzing the discrete-time Fourier transform of the windowed signal. Then the true values
of frequency, amplitude and phase can be calculated. Based on four-spectrum-line interpolation algorithm, this paper gives four practical rectification formulas by using the polynomial approximation method. It is proved that the phase difference of any two adjacent spectrum is π. According to this rule, this paper induces a rapid algorithm so that square root is only calculated one time when analyzing some harmonic parameters. The algorithm can save computing complexity and computing time greatly. The simulation results show that the four-spectrum-line interpolation algorithm has higher precision than double/triple- spectrum-line interpolation algorithm especially when estimating the even harmonics and the order of polynomial is lower than two others if obtaining same precision.
Key words: harmonic analysis; window function; FFT;four-spectrum-line; interpolation
0 引言
随着大量非线性电力电子器件的使用，电网中谐波污染问题日益严重[1-4]。谐波问题不仅恶化电能质量，对电网的安全稳定和经济运行也造成较大影响。因此，对系统中谐波参数进行高精度测量，对于减少谐波危害，维护电网安全稳定、高效运行是十分必要的。
快速傅里叶变换是目前用于电力系统谐波分析最常用的算法，但采用FFT对信号进行频谱分析很难做到同步采样和整周期截断，所带来的频谱泄露和栅栏效应会导致谐波的频率、幅值和相位等参数测量不准，
无法满足测量要求。为了解决这一问题，加窗插值FFT算法被广泛应用：运用各种特殊窗函数对信号进行截断克服非同步采样带来的频谱泄漏，采用插值法克服栅栏效应而提高测量精
- 140 - 电力系统保护与控制
度[5-6]。常用窗函数例如汉宁(Hanning)窗[7]、布莱克曼(Blackman)窗[8]、布莱克曼汉斯(Blackman-Harris)
窗函数[9]
、纳托尔(Nuttall)窗函数[10-11]、莱夫文森特(Rife-Vincent)窗函数[12]以及各种组合窗[13-20]。在加窗基础上，D. Agrez 和庞浩等人各自提出了双谱线的修正算法[21,6]，Wu Jing 、牛胜锁和黄冬梅等人提出了多谱线修正算法[22-27]。插值谱线数目的增多一方面提高了谐波分析的准确性，另一方面需要大量的求模运算。模的运算在CPU 中一般是采用牛顿迭代算法实现，运算量大且所耗时间长。为了进一步提高加窗插值算法的精度，本文推导出了实用的四插值修正公式，多种常用的余弦窗函数计算实例表明，该算法相对于双谱线和三谱线插值算法，有更高的准确度。同时证明了窗函数主瓣内任意相邻谱线相位差为π，提出一种快速算法，计算某次谐波开方计算量仅需一次。
1 四谱线及其快速算法原理frequency函数计算频数
以单一频率信号进行分析，假定一个频率为0f ，幅值为A ，初始相位为ϕ的信号()x t ，以采样频率s f 采样后，得到如下形式的离散信号： s ()sin(2π/)x n A nf f ϕ=+ (1)
其中，0,1,2,,1n N =- ，N 为采样点数。
所加离散余弦窗函数的表达式为
1
2π()(1)cos()M m N m m m
w n b n N -==-∑ (2)
其中：0,1,2,,1n N =- ；M 为窗函数项数；且窗函数系数m b 满足约束条件：1
1M m m b -==∑，1
(1)0M m m m b -=-=∑。
则对离散信号()x n 进行加窗，得到()w x n =
()()x n w n ，离散傅里叶变换后得到：
j j 00()[e (e ()]2f f A
X k W k W k j f f ϕϕ-=--+∆∆ (3)
其中，s f
f N
∆=为离散频率间隔 。若忽略负频率点
处旁瓣的影响，式(3)变为
j 0()e (2f A X k W k j f
ϕ
=-∆ (4)
对信号的非同步采样或非整周期数据截断时，由于栅栏效应，信号的峰值频率00f k f =⋅∆ 很难刚好位于离散谱线的频点上，即0k 一般不为整数。设峰值频率点左右各两条的谱线分别为12k k k << 34k k << ，这四根谱线对应的幅值分别为1y ，2y ，3y ，4y 。
记20.5k k α=-- ，由于201k k ≤-≤，则0.50.5α-≤≤，另记：
4321
4321
3333y y y y y y y y β+--=+++ (5)
根据式(4)和(5)可以得到：
|( 1.5)|3|(0.5)|3|(0.5)||( 1.5)||( 1.5)|3|(0.5)|3|(0.5)||( 1.5)|
W W W W W W W W ααααβαααα-++-+------=
-++-++--+-- (6) 当N 值较大时，式(6)可以化简为()g βα= ，其反函数为1()g αβ-=。由于所采用的余弦窗系数均为实系数，其频率响应是偶对称的，因而()g 和1()g - 均为奇函数。可采用多项式逼近方法计算奇
函数1()g αβ-=，表达式为
311131p p p p p αβββ≈⨯+⨯++ (7) 式(7)中，11131,,,p p p p 为多项式逼近的奇次项系数。
求得α后，求得信号频率：
02(0.5)f k f k f α=⋅∆=++∆ (8) 信号幅值，根据式(4)可知：
02|()|i i f
A y W k f =-∆ (9)
考虑到2y ，3y 是离真实谱线点最近的两根谱线，给予较大权重。可以得到：
43212(33)
|( 1.5)|3|(0.5)|3|(0.5)||( 1.5)|
y y y y A W W W W αααα+++=
-++-++--+--
(10)
类似式(7)的逼近方法，当N 比较大，窗函数系数为实系数，式(10)可表示为：143(3A N y y -=++ 213)()y y u α+，()u 为偶函数，逼近多项式不含奇
次项。四谱线修正逼近多项式如下：
12432120222(33)()d d A N y y y y p p p αα-=++++++
(11)
式(11)中，20222,,,d p p p 为多项式逼近的偶次项系数。
根据式(4)还可以得出信号的相位：
2π
arg[()]arg[(0.5)]2
X k W ϕα=+-+ (12)
根据式(6)、式(7)、式(9)、式(11)、式(12)即可进行各次谐波参数的分析。考虑到其中大量窗函数的离散傅里叶分析，其表达式为
1
jπ0
sin()
()sin(π)e
[(1)]ππ2sin(())sin(())M k
m
m m b km W k k k m k m N N
--==--+∑
(13)
由于1N >> ，可以得到：
张俊敏，等 基于四谱线插值FFT 的谐波分析快速算法 - 141 -
1220
sin(π)|()|(1)]ππarg(())(π)M m
m m b N k W k k m W k k N -=⎧=-⎪⎪-⎨
⎪=-+⎪⎩
∑ (14) 由式(5)、式(10)可知，在计算拟合参数β和幅
值A 的时候，需要求得43214321(33)
(33)y y y y y y y y +--+++，
4321(33)y y y y +++。 令1432124321()3()3()()()3()3()()
T X k X k X k X k T X k X k X k X k =--+⎧⎨=-+-⎩，即 根据式(4)可以得到：
π
arg(())arg(())2
X k W k ϕ=-++ (15)
根据式(13)可得到：
arg(())πW k k =- (16)
将式(16)带入式(15)可得到：
π
arg(())π2
X k k ϕ=-+- (17)
因此可以得到：1()X k 与3()X k ，2()X k 与4()X k 同相位；且1()X k 与2()X k ，2()X k 与3()X k ，3()X k 与4()X k 之间相位相差均为π。如图1所示。
图1 相邻四条谱线相位图
Fig. 1 Phase diagram of adjacent four spectrum lines
那么：
143214321j(π)j j(π)j 43214321
|||()3()3()()||()|3|()|3|()||()||e 3e 3e e |33T X k X k X k X k X k X k X k X k y y y y y y y y ϕϕϕϕ++=--+=
+--=--+=+--(18)
243214321j(π)
j j(π)
j 43214321
|||()3()3()()||()|3|()|3|()||()||e
3e
3e
e
|33T X k X k X k X k X k X k X k X k y y y y y y y y ϕϕ
ϕϕ
++=-+-=
+++=-+-=
+++
(19)
由图1可知：
1122
arg()arg(())
arg()arg(())T X k T X k =⎧⎨=⎩ (20)
因此，可以再次令：
12
cos(π)sin(π)
,0cos sin T C iC C D T D iD ϕϕϕϕ=+++⎧>⎨
=+⎩(21) 则：
111222Re()Im()||Re()Im()
T T T C T D T T β==== (22)
1|2|()A N T u α-= (23)
由以上分析可知，计算各次谐波的参数时候，仅在计算幅值A 时候需求2T 模运算一次，计算拟合
参数时候可以利用1T 和2T 的实部或虚部即可。
同理由于在同一个主瓣相邻4根谱线的相位关系，在需求最大谱线时候，可以直接利用插值前FFT 运算结果的实部和虚部来寻求最大的向量。
2 典型窗函数的修正公式
选取4种常用加窗函数，其系数如表1所示。
表1 常用窗函数系数 Table 1 Coefficients of the windows
窗系数
窗类型 0b
1b
2b
3b
Hanning 0.5 0.5 Blackman 0.42 0.5 0.08 Blackman-harris 0.358 75 0.488 29 0.141 28 0.011 68 4项最大旁瓣衰减
0.312 5
0.468 75
0.187 5
0.031 25
根据算法原理节的推导，表1中的6种窗函数的修正公式1()g αβ-=，以及复制修正公式中的逼
近公式()u α分别如下：
1) Hanning 窗
357
1.130136820.184086800.072248590.04451832αββββ=-+
- 246
()0.53549869+0.176221030.093108260.06644946u αααα=+
+
2) Balckman 窗
357
1.446490120.293265780.138584470.09578617αββββ=-+
-
246
() 1.1575129+0.561108880.387079970.33290703u αααα=+
+
3) Balckman-harris 窗
357
1.858620730.402997380.197257150.13293118αββββ=-+
-
- 142 - 电力系统保护与控制
246
() 1.24914356+0.887760250.738799070.69429366u αααα=+
-
4) 4项最大旁瓣衰减窗
357
2.375409830.434786650.191245170.11511086αββββ=-+
- 246
() 1.34456750+1.365774231.25666570 1.17989176u αααα=+
+
3 算法仿真
为了验证所提算法的精度，进行10次谐波仿真分析。信号模型为
21
11
s ()sin(2π)i i i f
x n A in f ϕ==+∑ (24)
其中：基波频率1f 为50.5 Hz ；采样频率s f 为5 120
Hz ；数据的截断长度N 为1 024点。仿真所采用的信号参数如表2所示。
对截断的数据分别加如表1所示的5种窗函数，然后进行FFT 插值分析。以下研究插值算法对检测精度的影响，所涉及的插值算法有：双谱线插值，三谱线插值以及本文提出的四谱线插值。修正算法中的拟合多项式次数均取7次，算法流程图在文献中均有详细说明，此处不予赘述。
仿真结果分别由表3~表6给出。其中Ai D 表示基波和各次谐波幅值测量值的相对误差；0f D 表示基波频率测量值的相对误差；i D ϕ表示基波和各次谐波初始相位测量值的相对误差。均用百分比表示。
表2 谐波信号参数
Table 2 Parameters of the harmonic signals
谐波次数 1 2 3 4 5 6 7 8 9 10 /V i A
220 4 15 3.5 7 2 3.7 2 2.3 0.8 o /()i ϕ
25
66
120.4
50
135
80
46
43.1
-19
表3 Hanning 窗的幅值、频率和相位相对误差
Table 3 Relative errors of Hanning window in calculating amplitude, frequency and phase
比较项 双谱线 三谱线 四谱线 比较项 双谱线 三谱线 四谱线 1A D 58.1910--⨯ 62.0710-⨯ 75.0910--⨯ 1D ϕ 49.1010--⨯ 53.2510--⨯ 51.7910-⨯ 2A D 27.1710--⨯ 31.8310-⨯ 58.4110--⨯ 2D ϕ 11.8610--⨯ 26.3310-⨯ 21.5710--⨯ 3A D 41.1510-⨯ 67.2110--⨯ 76.2310-⨯ 3D ϕ 22.4810--⨯ 36.4510-⨯ 31.6410-⨯ 4A D 24.8010-⨯ 41.3210--⨯ 56.6310--⨯ 4D ϕ 16.8410-⨯ 14.8310--⨯ 33.5210-⨯ 5A D 49.0210--⨯ 53.6910--⨯ 52.0710-⨯ 5D ϕ 33.9510--⨯ 45.0710--⨯ 41.5210-⨯ 6A D 36.7810-⨯ 42.7110--⨯ 57.8010-⨯ 6D ϕ 12.3510-⨯ 37.8910--⨯ 33.7910--⨯ 7A D 31.4210--⨯ 54.7710-⨯ 52.1010--⨯ 7D ϕ 22.6510-⨯ 21.3310--⨯ 31.2810--⨯ 8A D 21.2510-⨯ 38.5710--⨯ 57.7910--⨯ 8D ϕ 13.6510--⨯ 29.5410--⨯ 35.5410--⨯ 9A D 31.0110--⨯ 67.0010--⨯ 65.5510--⨯ 9D ϕ 39.2710--⨯ 31.8510--⨯ 42.4610-⨯ 10A D 24.2110--⨯ 41.1010-⨯ 59.0610-⨯ 10D ϕ
18.5610-⨯
23.8710--⨯
31.2110--⨯
0f D
53.3210--⨯
62.2410--⨯
72.7110--⨯
表4 Blackman 窗的幅值、频率和相位相对误差
Table 4 Relative errors of Blackman window in calculating amplitude, frequency and phase
比较项 双谱线 三谱线 四谱线 比较项 双谱线 三谱线 四谱线 1A D 53.2210-⨯ 61.5010-⨯ 71.5210-⨯ 1D ϕ 44.3610--⨯ 41.2510--⨯ 68.1710-⨯ 2A D 22.9010--⨯ 51.5810-⨯ 62.9110--⨯ 2D ϕ 29.6010--⨯ 22.9110-⨯ 36.6010-⨯ 3A D 46.9810-⨯ 62.3610--⨯ 61.2810--⨯ 3D ϕ 21.4210--⨯ 33.0010-⨯ 45.6610-⨯ 4A D 21.9510-⨯ 44.8410--⨯ 52.3610-⨯ 4D ϕ 13.8610-⨯ 12.0510--⨯ 49.9710-⨯ 5A D 43.0810--⨯ 51.8610--⨯ 67.0710-⨯ 5D ϕ 31.9910--⨯ 41.8010--⨯ 43.3910-⨯ 6A D 34.5710-⨯ 42.0710--⨯ 52.9410-⨯ 6D ϕ 11.2410-⨯ 35.8410-⨯ 49.7510--⨯ 7A D 46.0110--⨯ 51.7310-⨯ 66.8210--⨯ 7D ϕ 21.5410-⨯ 35.6810-⨯ 42.1710--⨯ 8A D 21.0210-⨯ 37.3410--⨯ 52.8810--⨯ 8D ϕ 11.2810-⨯ 27.4510--⨯ 44.4210--⨯ 9A D 44.2910-⨯ 74.3010--⨯ 71.1710--⨯ 9D ϕ 34.7110-⨯ 47.2310--⨯ 43.5410-⨯ 10A D 22.1410-⨯ 3
2.5610-⨯ 54.0310--⨯ 10D ϕ
14.2010-⨯
21.5610-⨯
34.7710-⨯
0f D
68.0210--⨯
76.0110--⨯
71.3110--⨯
张俊敏，等 基于四谱线插值FFT 的谐波分析快速算法 - 143 - 表5 Blackman-harris 窗的幅值、频率和相位相对误差
Table 5 Relative errors of Blackman-harris window in calculating amplitude, frequency and phase
比较项 双谱线 三谱线 四谱线 比较项 双谱线 三谱线 四谱线 1A D 41.7610--⨯ 77.3910--⨯ 86.3710--⨯ 1D ϕ 51.7610--⨯ 78.4010--⨯ 75.3710--⨯ 2A D 32.8910--⨯ 64.1110-⨯ 62.6410-⨯ 2D ϕ 21.1510--⨯ 33.3010-⨯ 46.9310-⨯ 3A D 41.4410--⨯ 71.7310--⨯ 78.7410--⨯ 3D ϕ 32.1610-⨯ 44.1710-⨯ 61.6210-⨯ 4A D 31.0710-⨯ 69.1010--⨯ 62.2210-⨯ 4D ϕ 22.4110-⨯ 21.8910--⨯ 53.9910-⨯ 5A D 61.7210--⨯ 77.1910--⨯ 75.3010-⨯ 5D ϕ 43.2910-⨯ 56.9610-⨯ 62.9810-⨯ 6A D 31.9810--⨯ 58.9410-⨯ 74.7510-⨯ 6D ϕ 27.5210-⨯ 48.9610-⨯ 55.4310--⨯ 7A D 43.5710--⨯ 65.3910-⨯ 82.9510--⨯ 7D ϕ 21.5510-⨯ 43.1710-⨯ 53.6710--⨯ 8A D 32.5610-⨯ 45.6310-⨯ 63.2010--⨯ 8D ϕ 11.0110-⨯ 23.1010-⨯ 41.7310--⨯ 9A D 41.5410--⨯ 65.8310-⨯ 82.4410--⨯ 9D ϕ 32.3010-⨯ 45.2810-⨯ 41.5110-⨯ 10A D 21.0910-⨯ 44.5210-⨯ 66.3510-⨯ 10D ϕ
12.3410-⨯
21.1010-⨯
47.6610--⨯
0f D
74.3810--⨯
83.9210--⨯
81.2510-⨯
表6 4项最大旁瓣衰减窗的幅值、频率和相位相对误差
Table 6 Relative errors of 4-term maximum-lobe decay speed window in calculating amplitude, frequency and phase
比较项 双谱线 三谱线 四谱线 比较项 双谱线 三谱线 四谱线 1A D 72.5910-⨯ 77.8010-⨯ 81.3410--⨯ 1D ϕ 72.8710--⨯ 72.5310-⨯ 94.7710--⨯ 2A D 43.0710--⨯ 54.5910-⨯ 64.2310--⨯ 2D ϕ 32.5310-⨯ 31.9510--⨯ 54.6610-⨯ 3A D 62.4910--⨯ 74.2510--⨯ 87.6010--⨯ 3D ϕ 51.7510-⨯ 65.7910-⨯ 61.5110-⨯ 4A D 41.9610-⨯ 53.3910--⨯ 76.2810-⨯ 4D ϕ 36.6110-⨯ 34.4510--⨯ 69.1310-⨯ 5A D 76.6910--⨯ 72.0310--⨯ 72.5410-⨯ 5D ϕ 67.2310--⨯ 64.6210-⨯ 61.9810-⨯ 6A D 31.0910--⨯ 57.9810-⨯ 76.1710-⨯ 6D ϕ 34.5710-⨯ 57.5410-⨯ 68.1110--⨯ 7A D 41.0610--⨯ 72.0310--⨯ 73.3410--⨯ 7D ϕ 52.5110-⨯ 65.7310-⨯ 62.1710--⨯ 8A D 31.6610-⨯ 43.5510--⨯ 77.7410--⨯ 8D ϕ 36.3210--⨯ 51.1210-⨯ 66.8510--⨯ 9A D 72.9410--⨯ 79.2710-⨯ 72.9210--⨯ 9D ϕ 64.0310--⨯ 74.8910-⨯ 81.5010-⨯ 10A D 38.7810-⨯ 61.6910-⨯ 76.2510--⨯ 10D ϕ
35.4210-⨯
53.6710--⨯
79.2510-⨯
0f D
71.9710--⨯
86.0110--⨯
101.3310-⨯
由以上的仿真结果可以看出，本文所推导的四谱线插值FFT 计算方法，计算结果普遍好于双谱线和三谱线插值算法。同时从表3~表6数据可以看出，随着窗函数项数的增加，四谱线插值在奇数次谐波检测上的性能提升空间不大，但在检测偶数次谐波上的精度提升很大，这也是四谱线插值的主要优点。且每次谐波仅需一次求模运算，大大提升了计算速度，节约了计算量。
4 结论
在对信号进行FFT 谐波分析时，通过加窗和插值算法可以减少由于非同步采样或者对数据的非整周期截
断所引起的误差。本文从信号加窗后的频域表达式入手，推导出了一种四谱线快速插值FFT 算法。基于该算法推导出了4种典型窗函数四插值的谐波分析修正公式。在4种窗函数情况下，进行了双、三和四谱线插值算法谐波分析的精度比较，
实验结果表明本文所提出的四谱线算法的总体精度的优越性以及在检测偶次谐波参数的计算精度更高，具有较高的实用价值。并且计算每次谐波参数时仅需要一次复数求模运算，大大节约了计算量。 参考文献
[1] 黄银龙, 乐健, 毛涛, 等. 一种新型的高压直流输电系
统直流侧谐波电压测量方法[J]. 电工技术学报, 2015, 30(22): 144-152.
HUANG Yinlong, LE Jian, MAO Tao, et al. A novel measurement method of the harmonic voltage on the DC-side of HVDC transmission system[J]. Transactions of China Electrotechnical Society, 2015, 30(22): 144-152.
[2] 赵鲁, 李耀华, 葛琼璇, 等. 特定谐波消除及优化脉宽
调制单相整流器的研究[J]. 电工技术学报, 2014, 29(10): 57-64.
版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。
发表评论