您的当前位置：首页基于加窗四谱线插值FFT的谐波快速分析方法及系统[发明专利]

基于加窗四谱线插值FFT的谐波快速分析方法及系统[发明专利]

来源：微智科技网

(19)中华人民共和国国家知识产权局

(12)发明专利申请

(10)申请公布号(10)申请公布号 CN 1047960 A (43)申请公布日(43)申请公布日 2015.09.09

(21)申请号 201510326063.0(22)申请日 2015.06.15

(71)申请人中南民族大学

地址430074 湖北省武汉市洪山区民族大道

182号(72)发明人张俊敏刘开培汪立王黎

田微何顺帆姚为(74)专利代理机构北京捷诚信通专利事务所

(普通合伙) 11221

代理人王卫东(51)Int.Cl.

G01R 23/16(2006.01)

权利要求书7页说明书19页附图4页

(54)发明名称

基于加窗四谱线插值FFT的谐波快速分析方法及系统(57)摘要

本发明公开了一种基于加窗四谱线插值FFT的谐波快速分析方法及系统，涉及谐波分析领域。该方法包括以下步骤：信号预处理；确定四根谱线；计算四谱线插值算法的修正公式；计算基波参数；四谱线插值快速算法；确定谐波参数；进行误差分析。本发明从电力系统的电号加窗后的频域表达式入手，根据窗函数主瓣内任意相邻谱线相位相差π的规律，推导出加窗FFT后的真实谱线点附近最大的四根谱线之间的相位规律，提出一种加窗四谱线插值FFT快速算法。相对于双谱线和三谱线插值算法，本发明能有效提高谐波分析的精度；利用该快速算法，计算某次谐波仅需要1次模的运算，能够有效降低计算量和计算时间，显著提升计算速度。

C N 1 0 4 8 9 7 9 6 0 A CN 1047960 A

权利要求书

1/7页

1.一种基于加窗四谱线插值FFT的谐波快速分析方法，其特征在于，包括以下步骤：

S1、信号预处理：

互感器采集电号，将互感器采集到的离散电号x(n)，传输到上位机，n为采样点的序数，n为自然数；采用离散余弦窗函数w(n)，对电号x(n)进行加窗截断，得到加窗信号xw(n)：

xw(n)＝x(n)w(n) (1)

离散余弦窗函数w(n)的表达式为：

其中，N为采样点数，N为正整数，n＝0,1,2...N-1；Σ表示求和；m为窗函数的累加次数，m＝0,1,2...M-1；M为窗函数项数，M为正整数；bm为窗函数系数；

对公式(1)的加窗信号进行FFT变换后，得到加窗FFT频谱：

其中，W(·)为窗函数的频谱，k为正整数，X(k)表示第k次谐波的频谱，Ak为第k次第一次谐谐波的幅值，j表示虚数单位，e是自然对数的底数，为第k次谐波的初始相位，波为基波，fs为采样频率，f0为基波频率，Δf为离散频率间隔，且Δf＝fs/N；令：k0＝f0/Δf，k0为真实频谱的谱线位置；忽略负频率点处旁瓣的影响，公式(3)变为：

S2、确定四根谱线：

在步骤S1得到的加窗FFT频谱峰值附近区域，k0处频率点左右各两条的谱线分别为：第k1、k2、k3、k4根，k1、k2、k3、k4均为正整数，k1、k2、k3、k4的关系为：k2＝k1+1，k3＝k2+1，k4＝k3+1，这四根谱线对应的幅值分别为y1、y2、y3、y4；

记变量α＝k-k2-0.5，由于0≤k-k2≤1，则-0.5≤α≤0.5；

记变量

S3、计算四谱线插值算法的修正公式：根据公式(4)和(5)得到：

采用多项式逼近方法，计算奇函数β＝g-1(α)，表达式为：

CN 1047960 A

权利要求书

2/7页

p11，p13；…p1p为多项式逼近的奇次项系数，p是奇数；

根据公式(4)，求得电号第i次谐波的幅值Ai：

其中，i为正整数，yi为加窗FFT后第i次谐波的幅值；考虑到y2、y3是离真实谱线点最近的两根谱线，基波幅值A1为：

根据公式(7)的逼近方法，N＞1000时，窗函数系数为实系数，公式(9)表示为：A1＝N-1(y4+3y3+3y2+y1)u(α)，其中，u(α)为修正公式，且为偶函数，逼近多项式不含奇次项；四谱线修正公式为：u(α)＝(p20+p22α2+…+p2dαd) (10)公式(10)中，p20,p22…p2d为多项式逼近的偶次项系数，d为拟合的最高阶次，且d为偶数；

S4、计算基波参数：

基波幅值A1：计算基波频率f0、

f0＝k·Δf＝(k2+α+0.5)Δf (11)

A1＝N-1(y4+3y3+3y2+y1)(p20+p22α2+…+p2dαd) (12)根据公式(4)，得出基波的相位：

仿照基波参数的求取，根据公式(6)、(7)、(11)、(12)、(13)，进行各次谐波参数的分析；

考虑到其中大量窗函数的离散傅里叶分析，其表达式为：

由于N>>1，得到：

S5、加窗四谱线插值快速算法：

根据公式(5)和(12)，计算变量β和幅值A1的时候，需求出

CN 1047960 A

权利要求书

(y4+3y3+3y2+y1)：

3/7页

令变量：

根据公式(4)，得到：

根据公式(13)，得到：arg(W(k))＝-kπ (18)将公式(17)代入公式(16)，得到：

分析得出：X(k1)与X(k3)同相位，X(k2)与X(k4)同相位；且X(k1)与X(k2)，X(k2)与X(k3)，X(k3)与X(k4)之间的相位之差均为π，那么：

令：

其中，C为T1的模，D为T2的模；则：

其中：Re表示取实部，Im表示取虚部；

分析得出：计算各次谐波的参数时，仅在计算幅值Ak的时候，需要进行一次求模运算；根据公式(5)计算变量β时，利用T1和T2的实部或虚部进行快速计算；同理，根据同一个主瓣相邻4根谱线的相位关系，求最大谱线时，直接利用插值前FFT运算结果的实部和虚部来寻求最大的向量；

S6、确定谐波参数：确定基波频率f0后，在范围(kf0-5,kf0+5)内重复步骤S2～S5，直

CN 1047960 A

权利要求书

4/7页

到所有谐波参数计算完毕；

S7、进行误差分析：在同样窗函数条件下，分析加窗四谱线插值快速算法的误差，并与双谱线算法的误差、三谱线插值算法的误差进行比较。

2.如权利要求1所述的基于加窗四谱线插值FFT的谐波快速分析方法，其特征在于：所述电号包括电流信号、电压信号。

3.如权利要求1所述的基于加窗四谱线插值FFT的谐波快速分析方法，其特征在于：所述真实频谱的谱线位置k0为小数。

4.如权利要求1或2或3所述的基于加窗四谱线插值FFT的谐波快速分析方法，其特征在于：

所述窗函数系数bm满足以下约束条件：

5.一种基于加窗四谱线插值FFT的谐波快速分析系统，其特征在于：包括信号预处理单元、谱线确定单元、修正公式计算单元、基波参数计算单元、四谱线插值快速计算单元、谐波参数确定单元、误差分析单元，其中：

所述信号预处理单元，用于对信号进行预处理：互感器采集电号，将互感器采集到的离散电号x(n)，传输到上位机，n为采样点的序数，n为自然数；采用离散余弦窗函数w(n)，对电号x(n)进行加窗截断，得到加窗信号xw(n)：

xw(n)＝x(n)w(n) (1)

离散余弦窗函数w(n)的表达式为：

其中，N为采样点数，N为正整数，n＝0,1,2...N-1；Σ表示求和；m为窗函数的累加次数，m＝0,1,2...M-1；M为窗函数项数，M为正整数；bm为窗函数系数；

对公式(1)的加窗信号进行FFT变换后，得到加窗FFT频谱：

其中，W(·)为窗函数的频谱，k为正整数，X(k)表示第k次谐波的频谱，Ak为第k次谐波的幅值，j表示虚数单位，e是自然对数的底数，为第k次谐波的初始相位，第一次谐波为基波，fs为采样频率，f0为基波频率，Δf为离散频率间隔，且Δf＝fs/N；令：k0＝f0/Δf，k0为真实频谱的谱线位置；忽略负频率点处旁瓣的影响，公式(3)变为：

所述谱线确定单元，用于确定四根谱线：

CN 1047960 A

权利要求书

5/7页

在信号预处理单元得到的加窗FFT频谱峰值附近区域，k0处频率点左右各两条的谱线分别为：第k1、k2、k3、k4根，k1、k2、k3、k4均为正整数，k1、k2、k3、k4的关系为：k2＝k1+1，k3＝k2+1，k4＝k3+1，这四根谱线对应的幅值分别为y1、y2、y3、y4；

记变量α＝k-k2-0.5，由于0≤k-k2≤1，则-0.5≤α≤0.5；

记变量

所述修正公式计算单元，用于计算四谱线插值算法的修正公式：根据公式(4)和(5)得到：

采用多项式逼近方法，计算奇函数β＝g-1(α)，表达式为：

p11，p13；…p1p为多项式逼近的奇次项系数，p是奇数；根据公式(4)，求得电号第i次谐波的幅值Ai：

其中，i为正整数，yi为加窗FFT后第i次谐波的幅值；考虑到y2、y3是离真实谱线点最近的两根谱线，基波幅值A1为：

所述基波参数计算单元，用于计算基波参数：计算基波频率f0、基波幅值A1：

f0＝k·Δf＝(k2+α+0.5)Δf (11)

A1＝N-1(y4+3y3+3y2+y1)(p20+p22α2+…+p2dαd) (12)根据公式(4)，得出基波的相位：

仿照基波参数的求取，根据公式(6)、(7)、(11)、(12)、(13)，进行各次谐波参数的分析；

CN 1047960 A

权利要求书

6/7页

考虑到其中大量窗函数的离散傅里叶分析，其表达式为：

由于N>>1，得到：

所述四谱线插值快速计算单元，用于进行四谱线插值快速计算：

根据公式(5)和(12)，在计算变量β和幅值A1的时候，需要求出

(y4+3y3+3y2+y1)：

令变量：

根据公式(4)，得到：

根据公式(13)，得到：arg(W(k))＝-kπ (18)将公式(17)代入公式(16)，得到：

分析得出：X(k1)与X(k3)同相位，X(k2)与X(k4)同相位；且X(k1)与X(k2)，X(k2)与X(k3)，X(k3)与X(k4)之间的相位之差均为π，那么：

CN 1047960 A

权利要求书

7/7页

令：

其中，C为T1的模，D为T2的模；则：

其中：Re表示取实部，Im表示取虚部；

分析得出：计算各次谐波的参数时，仅在计算幅值Ak的时候，需要进行一次求模运算；根据公式(5)计算变量β时，利用T1和T2的实部或虚部进行快速计算；同理，根据同一个主瓣相邻4根谱线的相位关系，求最大谱线时，直接利用插值前FFT运算结果的实部和虚部，来寻求最大的向量；

所述谐波参数确定单元，用于确定谐波参数：确定基波频率f0后，在范围(kf0-5,kf0+5)内，所述谱线确定单元、修正公式计算单元、基波参数计算单元、四谱线插值快速计算单元重复进行计算，直到所有谐波参数计算完毕；

所述误差分析单元，用于进行误差分析：在同样窗函数条件下，分析加窗四谱线插值快速算法的误差，并与双谱线算法的误差、三谱线插值算法的误差进行比较。

6.如权利要求5所述的基于加窗四谱线插值FFT的谐波快速分析系统，其特征在于：所述电号包括电流信号、电压信号。

7.如权利要求5所述的基于加窗四谱线插值FFT的谐波快速分析系统，其特征在于：所述真实频谱的谱线位置k0为小数。

8.如权利要求5或6或7所述的基于加窗四谱线插值FFT的谐波快速分析系统，其特征在于：

所述窗函数系数bm满足以下约束条件：

CN 1047960 A

说明书

1/19页

基于加窗四谱线插值FFT的谐波快速分析方法及系统

技术领域

本发明涉及谐波分析领域，具体是涉及一种基于加窗四谱线插值FFT的谐波快速分析方法及系统。

[0001]

背景技术

随着大量非线性电力电子器件的使用，使得电网中谐波污染问题日益严重。谐波

问题不仅恶化电能质量，对电网的安全稳定和经济运行也造成较大影响。因此，对系统中谐波参数进行高精度测量，对于减少谐波危害，维护电网安全稳定、高效运行是十分必要的。[0003] FFT(Fast Fourier Transform，快速傅里叶变换)是目前用于电力系统谐波分析最常用的算法，但采用FFT对信号进行频谱分析很难做到同步采样和整周期截断，所带来的频谱泄露和栅栏效应会导致谐波的频率、幅值和相位等参数测量不准，无法满足测量要求。

[0004] 为了解决这一问题，加窗插值FFT算法被广泛应用：运用各种特殊窗函数对信号进行截断，然后结合谱线插值FFT进行谐波分析，从而提高测量精度。[0005] 常用的窗函数有：汉宁(Hanning)窗函数、布莱克曼(Blackman)窗函数、布莱克曼汉斯(Blackman-Harris)窗函数、纳托尔(Nuttall)窗函数、莱夫文森特(Rife-Vincent)窗函数以及各种组合窗函数。[0006] 在加窗基础上，D.Agrez和庞浩等人各自提出了双谱线的修正算法，Wu Jing、牛胜锁和黄冬梅等人提出了三谱线修正算法。这些改进降低了频谱泄漏和栅栏效应的影响，提高了谐波分析的准确性。[0007] 然而，在工程实际使用中，双谱线和三谱线插值算法仍然无法满足高精度的谐波分析要求，并且随着所采用谱线数目增多，求模的计算量也迅速增加，导致计算速度较慢。

[0002]

发明内容

本发明的目的是为了克服上述背景技术的不足，提供一种基于加窗四谱线插值FFT的谐波快速分析方法及系统，能够有效提高谐波分析的精度；利用该快速算法，计算某次谐波仅需要1次模的运算，能够有效降低计算量和计算时间，显著提升计算速度。[0009] 本发明提供一种基于加窗四谱线插值FFT的谐波快速分析方法，包括以下步骤：[0010] S1、信号预处理：

[0011] 互感器采集电号，将互感器采集到的离散电号x(n)，传输到上位机，n为采样点的序数，n为自然数；采用离散余弦窗函数w(n)，对电号x(n)进行加窗截断，得到加窗信号xw(n)：

[0012] xw(n)＝x(n)w(n) (1)

[0008] [0013] [0014]

离散余弦窗函数w(n)的表达式为：

CN 1047960 A[0015]

说明书

2/19页

其中，N为采样点数，N为正整数，n＝0,1,2...N-1；Σ表示求和；m为窗函数的累加次数，m＝0,1,2...M-1；M为窗函数项数，M为正整数；bm为窗函数系数；[0016] 对公式(1)的加窗信号进行FFT变换后，得到加窗FFT频谱：

[0017]

[0018]

其中，W(·)为窗函数的频谱，k为正整数，X(k)表示第k次谐波的频谱，Ak为第k

次谐波的幅值，j表示虚数单位，e是自然对数的底数，为第k次谐波的初始相位，第一次谐波为基波，fs为采样频率，f0为基波频率，Δf为离散频率间隔，且Δf＝fs/N；[0019] 令：k0＝f0/Δf，k0为真实频谱的谱线位置；[0020] 忽略负频率点处旁瓣的影响，公式(3)变为：

[0021]

S2、确定四根谱线：

[0023] 在步骤S1得到的加窗FFT频谱峰值附近区域，k0处频率点左右各两条的谱线分别为：第k1、k2、k3、k4根，k1、k2、k3、k4均为正整数，k1、k2、k3、k4的关系为：k2＝k1+1，k3＝k2+1，k4＝k3+1，这四根谱线对应的幅值分别为y1、y2、y3、y4；[0024] 记变量α＝k-k2-0.5，由于0≤k-k2≤1，则-0.5≤α≤0.5；

[0022] [0025] [0026]

记变量

S3、计算四谱线插值算法的修正公式：[0027] 根据公式(4)和(5)得到：

[0028] [0029] [0030]

采用多项式逼近方法，计算奇函数β＝g-1(α)，表达式为：

α≈p11×β+p13×β3+…+p1pβp (7)

[0031] p11，p13；…p1p为多项式逼近的奇次项系数，p是奇数；[0032] 根据公式(4)，求得电号第i次谐波的幅值Ai：

[0033]

其中，i为正整数，yi为加窗FFT后第i次谐波的幅值；[0035] 考虑到y2、y3是离真实谱线点最近的两根谱线，基波幅值A1为：

[0034] [0036]

CN 1047960 A

说明书

3/19页

根据公式(7)的逼近方法，N＞1000时，窗函数系数为实系数，公式(9)表示为：

-1

[0038] A1＝N(y4+3y3+3y2+y1)u(α)，[0039] 其中，u(α)为修正公式，且为偶函数，逼近多项式不含奇次项；[0040] 四谱线修正公式为：u(α)＝(p20+p22α2+…+p2dαd) (10)[0041] 公式(10)中，p20,p22…p2d为多项式逼近的偶次项系数，d为拟合的最高阶次，且d为偶数；[0042] S4、计算基波参数：[0043] 计算基波频率f0、基波幅值A1：

[0044] f0＝k·Δf＝(k2+α+0.5)Δf (11)

-12d

[0045] A1＝N(y4+3y3+3y2+y1)(p20+p22α+…+p2dα) (12)[0046] 根据公式(4)，得出基波的相位：

[0037] [0047]

[0048] 仿照基波参数的求取，根据公式(6)、(7)、(11)、(12)、(13)，进行各次谐波参数的考虑到其中大量窗函数的离散傅里叶分析，其表达式为：

分析；

[0049] [0050]

[0051] 由于N>>1，得到：

[0052]

S5、加窗四谱线插值快速算法：

[0054] 根据公式(5)和(12)，计算变量β和幅值A1的时候，需求出

[0053]

(y4+3y3+3y2+y1)：

[0055] 令变量：

CN 1047960 A[0056] [0057] [0058] [0059] [0060]

说明书

4/19页

根据公式(4)，得到：

根据公式(13)，得到：arg(W(k))＝-kπ (18)将公式(17)代入公式(16)，得到：

分析得出：X(k1)与X(k3)同相位，X(k2)与X(k4)同相位；且X(k1)与X(k2)，X(k2)与X(k3)，X(k3)与X(k4)之间的相位之差均为π，那么：

[0061]

[0062]

[0063]

[00] [0065] [0066] [0067]

令：

其中，C为T1的模，D为T2的模；则：

其中：Re表示取实部，Im表示取虚部；[0068] 分析得出：计算各次谐波的参数时，仅在计算幅值Ak的时候，需要进行一次求模运算；根据公式(5)计算变量β时，利用T1和T2的实部或虚部进行快速计算；同理，根据同一个主瓣相邻4根谱线的相位关系，求最大谱线时，直接利用插值前FFT运算结果的实部和虚部来寻求最大的向量；[0069] S6、确定谐波参数：确定基波频率f0后，在范围(kf0-5,kf0+5)内重复步骤S2～S5，直到所有谐波参数计算完毕；[0070] S7、进行误差分析：在同样窗函数条件下，分析加窗四谱线插值快速算法的误差，并与双谱线算法的误差、三谱线插值算法的误差进行比较。

CN 1047960 A[0071]

说明书

5/19页

在上述技术方案的基础上，所述电号包括电流信号、电压信号；所述真实频谱的谱线位置k0为小数。

[0072] 在上述技术方案的基础上，所述窗函数系数bm满足以下约束条件：

[0073]

本发明还提供一种基于加窗四谱线插值FFT的谐波快速分析系统，包括信号预处理单元、谱线确定单元、修正公式计算单元、基波参数计算单元、四谱线插值快速计算单元、谐波参数确定单元、误差分析单元，其中：[0075] 所述信号预处理单元，用于对信号进行预处理：[0076] 互感器采集电号，将互感器采集到的离散电号x(n)，传输到上位机，n为采样点的序数，n为自然数；采用离散余弦窗函数w(n)，对电号x(n)进行加窗截断，得到加窗信号xw(n)：

[0077] xw(n)＝x(n)w(n) (1)[0078] 离散余弦窗函数w(n)的表达式为：

[0074] [0079]

其中，N为采样点数，N为正整数，n＝0,1,2...N-1；Σ表示求和；m为窗函数的累

加次数，m＝0,1,2...M-1；M为窗函数项数，M为正整数；bm为窗函数系数；[0081] 对公式(1)的加窗信号进行FFT变换后，得到加窗FFT频谱：

[0080] [0082]

[0083]

其中，W(·)为窗函数的频谱，k为正整数，X(k)表示第k次谐波的频谱，Ak为第

为第k次谐波的初始相位，第一

k次谐波的幅值，j表示虚数单位，e是自然对数的底数，

次谐波为基波，fs为采样频率，f0为基波频率，Δf为离散频率间隔，且Δf＝fs/N；[0084] 令：k0＝f0/Δf，k0为真实频谱的谱线位置；[0085] 忽略负频率点处旁瓣的影响，公式(3)变为：

[0086]

[0087] [0088]

所述谱线确定单元，用于确定四根谱线：

在信号预处理单元得到的加窗FFT频谱峰值附近区域，k0处频率点左右各两条的谱线分别为：第k1、k2、k3、k4根，k1、k2、k3、k4均为正整数，k1、k2、k3、k4的关系为：k2＝k1+1，k3＝k2+1，k4＝k3+1，这四根谱线对应的幅值分别为y1、y2、y3、y4；[00] 记变量α＝k-k2-0.5，由于0≤k-k2≤1，则-0.5≤α≤0.5；

CN 1047960 A

说明书

6/19页

[0090] [0091]

记变量

所述修正公式计算单元，用于计算四谱线插值算法的修正公式：[0092] 根据公式(4)和(5)得到：

[0093]

采用多项式逼近方法，计算奇函数β＝g-1(α)，表达式为：

[0095] α≈p11×β+p13×β+…+p1pβ (7)[0096] p11，p13；…p1p为多项式逼近的奇次项系数，p是奇数；[0097] 根据公式(4)，求得电号第i次谐波的幅值Ai：

[0094] [0098]

其中，i为正整数，yi为加窗FFT后第i次谐波的幅值；[0100] 考虑到y2、y3是离真实谱线点最近的两根谱线，基波幅值A1为：

[0099] [0101]

根据公式(7)的逼近方法，N＞1000时，窗函数系数为实系数，公式(9)表示为：

-1

[0103] A1＝N(y4+3y3+3y2+y1)u(α)，[0104] 其中，u(α)为修正公式，且为偶函数，逼近多项式不含奇次项；[0105] 四谱线修正公式为：u(α)＝(p20+p22α2+…+p2dαd) (10)[0106] 公式(10)中，p20,p22…p2d为多项式逼近的偶次项系数，d为拟合的最高阶次，且d为偶数；

[0107] 所述基波参数计算单元，用于计算基波参数：[0108] 计算基波频率f0、基波幅值A1：

[0109] f0＝k·Δf＝(k2+α+0.5)Δf (11)

-12d

[0110] A1＝N(y4+3y3+3y2+y1)(p20+p22α+…+p2dα) (12)[0111] 根据公式(4)，得出基波的相位：

[0102] [0112]

[0113] 仿照基波参数的求取，根据公式(6)、(7)、(11)、(12)、(13)，进行各次谐波参数的考虑到其中大量窗函数的离散傅里叶分析，其表达式为：

分析；

[0114]

CN 1047960 A

说明书

7/19页

[0115]

[0116] 由于N>>1，得到：

[0117]

所述四谱线插值快速计算单元，用于进行四谱线插值快速计算：

[0119] 根据公式(5)和(12)，在计算变量β和幅值A1的时候，需要求出

[0118]

(y4+3y3+3y2+y1)：

[0120] 令变量：

[0121] [0122] [0123] [0124] [0125]

根据公式(4)，得到：

根据公式(13)，得到：arg(W(k))＝-kπ (18)将公式(17)代入公式(16)，得到：

分析得出：X(k1)与X(k3)同相位，X(k2)与X(k4)同相位；且X(k1)与X(k2)，X(k2)

与X(k3)，X(k3)与X(k4)之间的相位之差均为π，那么：

[0126]

CN 1047960 A

说明书

8/19页

[0127]

[0128] [0129] [0130] [0131]

令：

其中，C为T1的模，D为T2的模；则：

其中：Re表示取实部，Im表示取虚部；[0132] 分析得出：计算各次谐波的参数时，仅在计算幅值Ak的时候，需要进行一次求模运算；根据公式(5)计算变量β时，利用T1和T2的实部或虚部进行快速计算；同理，根据同一个主瓣相邻4根谱线的相位关系，求最大谱线时，直接利用插值前FFT运算结果的实部和虚部，来寻求最大的向量；

[0133] 所述谐波参数确定单元，用于确定谐波参数：确定基波频率f0后，在范围(kf0-5,kf0+5)内，所述谱线确定单元、修正公式计算单元、基波参数计算单元、四谱线插值快速计算单元重复进行计算，直到所有谐波参数计算完毕；[0134] 所述误差分析单元，用于进行误差分析：在同样窗函数条件下，分析加窗四谱线插值快速算法的误差，并与双谱线算法的误差、三谱线插值算法的误差进行比较。[0135] 与现有技术相比，本发明的优点如下：

[0136] 本发明从电力系统的电号(电流信号或电压信号)加窗后的频域表达式入手，根据窗函数主瓣内任意相邻谱线相位相差π的规律，推导出加窗FFT后的真实谱线点附近最大的四根谱线之间的相位规律，提出一种加窗四谱线插值FFT快速算法。多种常用的余弦窗函数计算实例表明，相对于双谱线和三谱线插值算法，本发明能够有效提高谐波分析的精度；利用该快速算法，计算某次谐波仅需要1次模的运算，能够有效降低计算量和计算时间，显著提升计算速度，具有较高的实用价值。附图说明

[0137] [0138] [0139] [0140] [0141] [0142] [0143]

图1是本发明实施例中基于加窗四谱线插值FFT的谐波快速分析方法的流程图。图2是本发明实施例中相邻四条谱线的相位图。

图3是本发明实施例中谐波信号的波形图。图4是Hanning窗函数的幅频特性图。图5是Black窗函数的幅频特性图。

图6是Black-harris窗函数的幅频特性图。图7是4项最大旁瓣衰减窗函数的幅频特性图。

具体实施方式

[0144] 下面结合附图及具体实施例对本发明作进一步的详细描述。[0145] 参见图1所示，本发明实施例提供一种基于加窗四谱线插值FFT的谐波快速分析

CN 1047960 A

说明书

9/19页

方法，包括以下步骤：[0146] S1、信号预处理：

[0147] 互感器采集电号，电号包括电流信号、电压信号，将互感器采集到的离散电号x(n)，传输到上位机，n为采样点的序数，n为自然数；采用离散余弦窗函数w(n)，对电号x(n)进行加窗截断，得到加窗信号xw(n)：[0148] xw(n)＝x(n)w(n) (1)[0149] 离散余弦窗函数w(n)的表达式为：

[0150]

其中，N为采样点数，N为正整数，n＝0,1,2...N-1；Σ表示求和；m为窗函数的累

加次数，m＝0,1,2...M-1；M为窗函数项数，M为正整数；bm为窗函数系数，窗函数系数bm满足以下约束条件：

[0151] [0152] [0153] [0154]

对公式(1)的加窗信号进行FFT变换后，得到加窗FFT频谱：

[0155]

其中，W(·)为窗函数的频谱，k为正整数，X(k)表示第k次谐波的频谱，Ak为第k

次谐波的幅值，j表示虚数单位，e是自然对数的底数，为第k次谐波的初始相位，第一次谐波为基波，fs为采样频率，f0为基波频率，Δf为离散频率间隔，且Δf＝fs/N。令：k0＝f0/Δf，k0为真实频谱的谱线位置，由于非同步采样或非整周期截断，k0

一般为小数，不为整数。

[0157] 若忽略负频率点处旁瓣的影响，公式(3)变为：

[0156] [0158]

S2、确定四根谱线：

[0160] 在步骤S1得到的加窗FFT频谱峰值附近区域，k0处频率点左右各两条的谱线分别为：第k1、k2、k3、k4根，k1、k2、k3、k4均为正整数，k1、k2、k3、k4的关系为：k2＝k1+1，k3＝k2+1，k4＝k3+1，这四根谱线对应的幅值分别为y1、y2、y3、y4。[0161] 为方便计算，记变量α＝k-k2-0.5，由于0≤k-k2≤1，则-0.5≤α≤0.5；

[0159] [0162] [0163]

记变量

S3、计算四谱线插值算法的修正公式：

CN 1047960 A[01] [0165]

说明书

10/19页

根据公式(4)和(5)得到：

当N值比较大，例如：N＞1000时，公式(6)可以化简为一个函数β＝g(α)，其反函数为α＝g-1(β)。由于所采用的余弦窗系数均为实系数，其频率响应是偶对称的，因而g(·)和g-1(·)均为奇函数。[0167] 采用多项式逼近方法，计算奇函数β＝g-1(α)，表达式为：

[0168] α≈p11×β+p13×β+…+p1pβ (7)[0169] p11，p13；…p1p为多项式逼近的奇次项系数，p是奇数。[0170] 根据公式(4)，求得电号第i次谐波的幅值Ai：

[0166] [0171]

其中，i为正整数，yi为加窗FFT后第i次谐波的幅值；[0173] 以基波为例，考虑到y2、y3是离真实谱线点最近的两根谱线，给予较大权重，可以得到基波幅值A1：

[0172] [0174]

根据公式(7)的逼近方法，当N比较大，例如N＞1000时，窗函数系数为实系数，

公式(9)可表示为：

-1

[0176] A1＝N(y4+3y3+3y2+y1)u(α)，[0177] 其中，u(α)为修正公式，且为偶函数，逼近多项式不含奇次项。[0178] 四谱线修正公式如下：

[0179] u(α)＝(p20+p22α+…+p2dα) (10)[0180] 公式(10)中，p20,p22…p2d为多项式逼近的偶次项系数，d为拟合的最高阶次，且d为偶数。[0181] S4、计算基波参数：[0182] 考虑到y2、y3是离真实谱线点最近的两根谱线，给予较大权重，可以得到基波频率f0、基波幅值A1：

[0183] f0＝k·Δf＝(k2+α+0.5)Δf (11)

-12d

[0184] A1＝N(y4+3y3+3y2+y1)(p20+p22α+…+p2dα) (12)[0185] 根据公式(4)，还可以得出基波的相位：

[0175] [0186]

仿照基波参数的求取，根据公式(6)、(7)、(11)、(12)、(13)，即可进行各次谐波参数的分析。

[0188] 考虑到其中大量窗函数的离散傅里叶分析，其表达式为：

[0187]

CN 1047960 A

说明书

11/19页

[01]

[0190] 由于N>>1，可以得到：

[0191]

S5、加窗四谱线插值快速算法：

[0193] 根据公式(5)和(12)，在计算变量β和幅值A1的时候，需要求出

[0192]

(y4+3y3+3y2+y1)。

[0194] 令变量：

[0195] [0196] [0197] [0198]

根据公式(4)可以得到：

根据公式(13)可得到：arg(W(k))＝-kπ (18)将公式(17)代入公式(16)可得到：

分析得出：X(k1)与X(k3)同相位，X(k2)与X(k4)同相位；且X(k1)与X(k2)，X(k2)

与X(k3)，X(k3)与X(k4)之间的相位之差均为π，相邻四条谱线的相位图参见图2所示，那么：

[0199]

[0200]

CN 1047960 A

说明书

12/19页

[0201] [0202] [0203]

由图2可知：因此，可以再次令：

[0204] [0205] [0206]

其中，C为T1的模，D为T2的模。则：

其中：Re表示取实部，Im表示取虚部。[0207] 由以上分析得出：计算各次谐波的参数时候，仅在计算幅值Ak的时候，需要进行一次求模运算。根据公式(5)计算变量β时，可以利用T1和T2的实部或虚部进行快速计算。[0208] 同理，根据同一个主瓣相邻4根谱线的相位关系，求最大谱线时，直接利用插值前FFT运算结果的实部和虚部，来寻求最大的向量。[0209] S6、确定谐波参数：确定基波频率f0后，在范围(kf0-5,kf0+5)内，重复步骤S2～S5，直到所有谐波参数计算完毕。[0210] S7、进行误差分析：在同样窗函数条件下，分析加窗四谱线插值快速算法的误差，并与双谱线算法的误差、三谱线插值算法的误差进行比较。[0211] S1、S2、S3、S4、S6、S7形成一个完整的基于加窗四谱线插值FFT的高精度谐波分析方法，能够有效提高谐波分析的精度。[0212] 步骤S5的作用是提高计算速度，在S4、S6之间增加步骤S5，每次谐波仅需一次模运算，计算量小，能够有效提升计算速度，在同样窗函数下，能获得更高的精度，具有较高的实用价值。

[0213] 下面通过一个具体实施例来进行详细说明。选定一个电号的谐波，谐波的波形参见图3所示。[0214] 步骤1、信号的预处理：

[0215] 将互感器采集到的离散电号x(n)，信号模型参见表1所示。为了验证所提算法的精度，进行10次谐波仿真分析。信号模型为：

[0216]

其中：基波频率f0为50.5Hz，采样频率fs为5120Hz，采样点数N为1024。[0218] 表1、电号的谐波参数

[0217] [0219]

[0220]

CN 1047960 A

说明书

13/19页

对信号x(n)进行加窗处理，选取4种常用加窗函数，其系数参见表2所示。[0222] 表2、常用窗函数系数

[0221] [0223]

[0224] [0225]

加窗FFT变换后得到信号频谱：

步骤2、确定四根谱线：

[0226] 在45～55Hz中寻求实部(或虚部)最大四根谱线，分别为X(k1)、X(k2)、X(k3)、X(k4)。令：

[0227] [0228]

那么

步骤3、计算四谱线插值算法的修正公式：

[0229] 根据具体实施方式中步骤S3的推导，表2中的4种窗函数的修正公式α＝g-1(β)，u(α)分别如下：[0230] (1)Hanning窗：

357

[0231] α＝1.13013682β-0.18408680β+0.07224859β-0.04451832βu(α)＝0.53549869+0.17622103α2+0.09310826α4+0.0946α6[0232] Hanning窗函数的幅频特性参见图4所示。[0233] (2)Balckman窗：

357

[0234] α＝1.449012β-0.29326578β+0.13858447β-0.09578617βu(α)＝1.1575129+0.56110888α2+0.38707997α4+0.33290703α6[0235] Black窗函数的幅频特性参见图5所示。[0236] (3)Balckman-harris窗：

357

[0237] α＝1.85862073β-0.40299738β+0.19725715β-0.13293118βu(α)＝1.24914356+0.88776025α2+0.73879907α4-0.69429366α6[0238] Black-harris窗函数的幅频特性参见图6所示。[0239] (4)4项最大旁瓣衰减窗：

357

[0240] α＝2.37540983β-0.43478665β+0.19124517β-0.11511086βu(α)＝

CN 1047960 A

说明书

14/19页

1.34456750+1.36577423α2+1.25666570α4+1.179176α6

[0241] 4项最大旁瓣衰减窗函数的幅频特性参见图7所示。[0242] 步骤4、计算基波参数：

-1

[0243] A1＝N(y4+3y3+3y2+y1)u(α)

[0244]

f0＝k·Δf＝(k2+α+0.5)Δf

[0246] 步骤5、四谱线插值快速算法：[0247] 根据S5所述快速算法，可以对步骤2和步骤4中的两个量的求取采用快速算法：

[0245] [0248] [0249] [0250]

及A1＝N-1|T2|u(α)

步骤6、确定谐波参数：确定基波频率f0后，在范围(kf0-5,kf0+5)内重复步骤2～5，直到所有谐波参数计算完毕。

步骤7、进行误差分析：仿真实验的误差分析比较结果参见表3～表6所示，其中，

DAi表示基波和各次谐波幅值测量值的相对误差；Df0表示基波频率测量值的相对误差；表示基波和各次谐波初始相位测量值的相对误差，均用百分比表示。[0251] 表3、Hanning窗和Blackman窗的频率、幅值相对误差

[0252]

[0253] [0254]

表4、Hanning窗和Blackman窗的相位相对误差

CN 1047960 A

说明书

15/19页

[0255] [0256]

表5、Blackman-harris窗和四项最大旁瓣衰减窗的频率、幅值相对误差

[0257] [0258]

表6、Blackman-harris窗和四项最大旁瓣衰减窗的相位相对误差

CN 1047960 A[0259]

说明书

16/19页

从表3～表6数据可以看出，本发明实施例所推导的加窗四谱线插值FFT快速算法，计算结果普遍好于双谱线和三谱线插值算法，并且计算每次谐波仅需一次模运算，能够有效节约计算量和计算时间，显著提升计算速度。

[0260] 本发明实施例还提供一种基于加窗四谱线插值FFT的谐波快速分析系统，包括信号预处理单元、谱线确定单元、修正公式计算单元、基波参数计算单元、四谱线插值快速计算单元、谐波参数确定单元、误差分析单元，其中：[0261] 信号预处理单元，用于对信号进行预处理：互感器采集电号，电号包括电流信号、电压信号，将互感器采集到的离散电号x(n)，传输到上位机，n为采样点的序数，n为自然数；采用离散余弦窗函数w(n)，对电号x(n)进行加窗截断，得到加窗信号xw(n)：

[0262] xw(n)＝x(n)w(n) (1)[0263] 离散余弦窗函数w(n)的表达式为：

[02]

其中，N为采样点数，N为正整数，n＝0,1,2...N-1；Σ表示求和；m为窗函数的累

加次数，m＝0,1,2...M-1；M为窗函数项数，M为正整数；bm为窗函数系数，窗函数系数bm满足以下约束条件：

[0265] [0266] [0267] [0268]

对公式(1)的加窗信号进行FFT变换后，得到加窗FFT频谱：

[0269]

其中，W(·)为窗函数的频谱，k为正整数，X(k)表示第k次谐波的频谱，Ak为第k

次谐波的幅值，j表示虚数单位，e是自然对数的底数，为第k次谐波的初始相位，第一次谐波为基波，fs为采样频率，f0为基波频率，Δf为离散频率间隔，且Δf＝fs/N。

[0270] 令：k0＝f0/Δf，k0为真实频谱的谱线位置，由于非同步采样或非整周期截断，k0一般不为整数。

[0271] 若忽略负频率点处旁瓣的影响，公式(3)变为：

[0272]

谱线确定单元，用于确定四根谱线：

[0274] 在信号预处理单元得到的加窗FFT频谱峰值附近区域，k0处频率点左右各两条的谱线分别为：第k1、k2、k3、k4根，k1、k2、k3、k4均为正整数，k1、k2、k3、k4的关系为：k2＝k1+1，k3＝k2+1，k4＝k3+1，这四根谱线对应的幅值分别为y1、y2、y3、y4。

[0273]

CN 1047960 A[0275] [0276] [0277]

说明书

17/19页

为方便计算，记变量α＝k-k2-0.5，由于0≤k-k2≤1，则-0.5≤α≤0.5；记变量

修正公式计算单元，用于计算四谱线插值算法的修正公式：[0278] 根据公式(4)和(5)得到：

[0279]

当N值比较大，例如：N＞1000时，公式(6)可以化简为一个函数β＝g(α)，其

反函数为α＝g-1(β)。由于所采用的余弦窗系数均为实系数，其频率响应是偶对称的，因而g(·)和g-1(·)均为奇函数。[0281] 采用多项式逼近方法，计算奇函数β＝g-1(α)，表达式为：

[0282] α≈p11×β+p13×β+…+p1pβ (7)[0283] p11，p13；…p1p为多项式逼近的奇次项系数，p是奇数。[0284] 根据公式(4)，求得电号第i次谐波的幅值Ai：

[0280] [0285]

其中，i为正整数，yi为加窗FFT后第i次谐波的幅值；[0287] 以基波为例，考虑到y2、y3是离真实谱线点最近的两根谱线，给予较大权重，可以得到基波幅值A1：

[0286] [0288]

根据公式(7)的逼近方法，当N比较大，例如N＞1000时，窗函数系数为实系数，

公式(9)可表示为：

-1

[0290] A1＝N(y4+3y3+3y2+y1)u(α)，[0291] 其中，u(α)为修正公式，且为偶函数，逼近多项式不含奇次项。[0292] 四谱线修正公式如下：

[02]

u(α)＝(p20+p22α2+…+p2dαd) (10)

[0294] 公式(10)中，p20,p22…p2d为多项式逼近的偶次项系数，d为拟合的最高阶次，且d为偶数。

[0295] 基波参数计算单元，用于计算基波参数：[0296] 考虑到y2、y3是离真实谱线点最近的两根谱线，给予较大权重，可以得到基波频率f0、基波幅值A1：

[0297] f0＝k·Δf＝(k2+α+0.5)Δf (11)

-12d

[0298] A1＝N(y4+3y3+3y2+y1)(p20+p22α+…+p2dα) (12)[0299] 根据公式(4)还可以得出基波的相位：

[0293] [0300]

CN 1047960 A

说明书

18/19页

仿照基波参数的求取，根据公式(6)、(7)、(11)、(12)、(13)，即可进行各次谐波参数的分析。

[0302] 考虑到其中大量窗函数的离散傅里叶分析，其表达式为：

[0301] [0303]

[0304] 由于N>>1，可以得到：

[0305]

四谱线插值快速计算单元，用于进行四谱线插值快速计算：[0307] 由公式(5)和(12)可知，在计算变量β和幅值A1的时候，需要求出

[0306]

(y4+3y3+3y2+y1)。

[0308] 令变量：

[0309] [0310] [0311] [0312]

根据公式(4)可以得到：

根据公式(13)可得到：arg(W(k))＝-kπ (18)将公式(17)代入公式(16)可得到：

分析得出：X(k1)与X(k3)同相位，X(k2)与X(k4)同相位；且X(k1)与X(k2)，X(k2)与X(k3)，X(k3)与X(k4)之间的相位之差均为π，相邻四条谱线的相位图参见图2所示，那么：

[0313]

[0314]

CN 1047960 A

说明书

19/19页

[0315] [0316] [0317]

由图2可知：因此，可以再次令：

[0318] [0319] [0320]

其中，C为T1的模，D为T2的模。则：

其中：Re表示取实部，Im表示取虚部。[0321] 由以上分析可知，计算各次谐波的参数时候，仅在计算幅值Ak的时候，需要进行一次求模运算。根据公式(5)计算变量β时，可以利用T1和T2的实部或虚部进行快速计算。同理，根据同一个主瓣相邻4根谱线的相位关系，求最大谱线时，直接利用插值前FFT运算结果的实部和虚部，来寻求最大的向量。[0322] 谐波参数确定单元，用于确定谐波参数：确定基波频率f0后，在范围(kf0-5,kf0+5)内，谱线确定单元、修正公式计算单元、基波参数计算单元、四谱线插值快速计算单元重复进行计算，直到所有谐波参数计算完毕。[0323] 误差分析单元，用于进行误差分析：在同样窗函数条件下，分析加窗四谱线插值快速算法的误差，并与双谱线算法的误差、三谱线插值算法的误差进行比较。[0324] 本领域的技术人员可以对本发明实施例进行各种修改和变型，倘若这些修改和变型在本发明权利要求及其等同技术的范围之内，则这些修改和变型也在本发明的保护范围之内。

[0325] 说明书中未详细描述的内容为本领域技术人员公知的现有技术。

CN 1047960 A

说明书附图

1/4页

图1

CN 1047960 A

说明书附图

2/4页

图2

图3

CN 1047960 A

说明书附图

3/4页

图4

图5

CN 1047960 A

说明书附图

4/4页

图6

图7

因篇幅问题不能全部显示，请点此查看更多更全内容

查看全文