基于加窗四谱线插值fft的谐波快速分析方法及系统的制作方法
基于加窗四谱线插值fft的谐波快速分析方法及系统的制作方法
【技术领域】
[0001] 本发明涉及谐波分析领域，具体是涉及一种基于加窗四谱线插值FFT的谐波快速 分析方法及系统。
【背景技术】
[0002] 随着大量非线性电力电子器件的使用，使得电网中谐波污染问题日益严重。谐波 问题不仅恶化电能质量，对电网的安全稳定和经济运行也造成较大影响。因此，对系统中谐 波参数进行高精度测量，对于减少谐波危害，维护电网安全稳定、高效运行是十分必要的。
[0003] FFT(Fast Fourier Transform，快速傅里叶变换）是目前用于电力系统谐波分析 最常用的算法，但采用FFT对信号进行频谱分析很难做到同步采样和整周期截断，所带来 的频谱泄露和栅栏效应会导致谐波的频率、幅值和相位等参数测量不准，无法满足测量要 求。
[0004] 为了解决这一问题，加窗插值FFT算法被广泛应用：运用各种特殊窗函数对信号 进行截断，然后结合谱线插值FFT进行谐波分析，从而提高测量精度。
[0005] 常用的窗函数有：汉宁（Hanning)窗函数、布莱克曼（Blackman)窗函数、布莱克曼 汉斯（Blackman-Harris)窗函数、纳托尔（Nuttall)窗函数、莱夫文森特（Rife-Vincent) 窗函数以及各种组合窗函数。
[0006] 在加窗基础上，D. Agrez和庞浩等人各自提出了双谱线的修正算法，Wu Jing、牛胜 锁和黄冬梅等人提出了三谱线修正算法。这些改进降低了频谱泄漏和栅栏效应的影响，提 高了谐波分析的准确性。
[0007] 然而，在工程实际使用中，双谱线和三谱线插值算法仍然无法满足高精度的谐波 分析要求，并且随着所采用谱线数目增多，求模的计算量也迅速增加，导致计算速度较慢。
【发明内容】
[0008] 本发明的目的是为了克服上述【背景技术】的不足，提供一种基于加窗四谱线插值 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)
[0013] 离散余弦窗函数w(n)的表达式为：
[0015] 其中，N为采样点数，N为正整数，n = 0, 1，2. . . N-l ; 2表示求和；m为窗函数的累 加次数，m = 0, 1，2.. .M-l ;M为窗函数项数，M为正整数；\为窗函数系数；
[0016] 对公式⑴的加窗信号进行FFT变换后，得到加窗FFT频谱：
[0018] 其中，W(_)为窗函数的频谱，k为正整数，X(k)表示第k次谐波的频谱，Ak为第k 次谐波的幅值，j表示虚数单位，e是自然对数的底数，P是为第k次谐波的初始相位，第一次 谐波为基波，f s为采样频率，f〇为基波频率，A f为离散频率间隔，且A f = f s/N ;
[0019]令：kQ= f Q/ A f，kQ为真实频谱的谱线位置；
[0020] 忽略负频率点处旁瓣的影响，公式（3)变为：
[0022] S2、确定四根谱线：
[0023] 在步骤S1得到的加窗FFT频谱峰值附近区域，匕处频率点左右各两条的谱线分 别为：第 V k2、k3、k4根，k^k2、k3、k4均为正整数，k^k2、k3、1^4的关系为：k 2= k 汴1，k3 = k2+l，k4= k 3+1，这四根谱线对应的幅值分别为yp y2、y3、y4;
[0024] 记变量 a = k-k2-0. 5,由于 0 彡 k_k2< 1，则-0? 5 彡 a 彡 〇? 5 ;
[0025]记变量，
[0026] S3、计算四谱线插值算法的修正公式：
[0027] 根据公式⑷和（5)得到：
[0029] 采用多项式逼近方法，计算奇函数0 =g4(a)，表达式为：
[0030] a ~ pnX 0+p13X 03+…+plp0p (7)
[0031] pn，p13;…p lp为多项式逼近的奇次项系数，p是奇数；
[0032]根据公式（4)，求得电网信号第i次谐波的幅值Ai:
[0034] 其中，i为正整数，yi为加窗FFT后第i次谐波的幅值；
[0035]考虑到y2、y3是离真实谱线点最近的两根谱线，基波幅值八 :为：
[0036]
[0037] 根据公式（7)的逼近方法，N> 1000时，窗函数系数为实系数，公式（9)表示为：
[0038] Aj = N (74+373+372+7!) u ( a ),
[0039] 其中，u(a)为修正公式，且为偶函数，逼近多项式不含奇次项；
[0040] 四谱线修正公式为：u(a ) = (p2Q+p22a 2+."+p2da d) (10)
[0041] 公式（10)中，p2(l，p22…p2d为多项式逼近的偶次项系数，d为拟合的最高阶次，且d 为偶数；
[0042] S4、计算基波参数：
[0043] 计算基波频率f(i、基波幅值Ay
[0044] f0= k ? A f = (k 2+ a +〇. 5) A f (11)
[0045] Af N-1 (p2CI+p22 a 2+... +p2da d) (12)
[0046] 根据公式⑷，得出基波的相位：
[0048]仿照基波参数的求取，根据公式（6)、（7)、（11)、（12)、（13)，进行各次谐波参数的 分析；
[0049] 考虑到其中大量窗函数的离散傅里叶分析，其表达式为：
[0051]由于N>>1，得到：
[0053] S5、加窗四谱线插值快速算法：
[0054] 根据公式（5)和（12)，计算变量0和幅值Ai的时候，需求出
.(y4+3y3 +3y2+yi):
[0055] 令变量：
[0056] 根据公式⑷，得到：
[0057] 根据公式（13)，得到：arg(W(k)) = -k JT (18)
[0058]将公式（17)代入公式（16)，得到：
[0059] 分析得出：X〇〇与X(k3)同相位，X(k 2)与X(k4)同相位；且XG〇与X(k2)，X(k2) 与X (k3)，X (k3)与X (k4)之间的相位之差均为31，那么：
[0065]其中，C为T1的模，D为T2的模；
[0067] 其中：Re表示取实部，Im表示取虚部；
[0068] 分析得出：计算各次谐波的参数时，仅在计算幅值Ak的时候，需要进行一次求模运 算；根据公式（5)计算变量0时，利用T1和T2的实部或虚部进行快速计算；同理，根据同 一个主瓣相邻4根谱线的相位关系，求最大谱线时，直接利用插值前FFT运算结果的实部和 虚部来寻求最大的向量；
[0069] S6、确定谐波参数：确定基波频率fQ后，在范围（kf。_5, kfQ+5)内重复步骤S2~ S5,直到所有谐波参数计算完毕；
[0070] S7、进行误差分析：在同样窗函数条件下，分析加窗四谱线插值快速算法的误差， 并与双谱线算法的误差、三谱线插值算法的误差进行比较。
[0071] 在上述技术方案的基础上，所述电网信号包括电流信号、电压信号；所述真实频谱 的谱线位置h为小数。
[0072] 在上述技术方案的基础上，所述窗函数系数bm满足以下约束条件：
[0074] 本发明还提供一种基于加窗四谱线插值FFT的谐波快速分析系统，包括信号预处 理单元、谱线确定单元、修正公式计算单元、基波参数计算单元、四谱线插值快速计算单元、 谐波参数确定单元、误差分析单元，其中：
[0075] 所述信号预处理单元，用于对信号进行预处理：
[0076] 互感器采集电网信号，将互感器采集到的离散电网信号x(n)，传输到上位机，n为 采样点的序数，n为自然数；采用离散余弦窗函数w (n)，对电网信号x (n)进行加窗截断，得 到加窗信号xw(n):
[0077] xw(n) = x(n)w(n) (1)
[0078] 离散余弦窗函数w(n)的表达式为：
[0080] 其中，N为采样点数，N为正整数，n = 0, 1，2. . . N-1 ; 2表示求和；m为窗函数的累 加次数，m = 0, 1，2. . . M-l ;M为窗函数项数，M为正整数；\为窗函数系数；
[0081] 对公式⑴的加窗信号进行FFT变换后，得到加窗FFT频谱：
[0083] 其中，W( ?)为窗函数的频谱，k为正整数，X(k)表示第k次谐波的频谱，Ak为第 k次谐波的幅值，j表示虚数单位，e是自然对数的底数，为第k次谐波的初始相位，第一 次谐波为基波，f s为采样频率，f 〇为基波频率，A f为离散频率间隔，且A f = f s/N ;
[0084] 令：kQ= f Q/ A f，kQ为真实频谱的谱线位置；
[0085] 忽略负频率点处旁瓣的影响，公式（3)变为：
[0087] 所述谱线确定单元，用于确定四根谱线：
[0088] 在信号预处理单元得到的加窗FFT频谱峰值附近区域，匕处频率点左右各两条的 谱线分别为：第1^ 1、1^2、1^、1^4根，1^1、1^2、1^、1^ 4均为正整数，1^1、1^2、1^、1^4的关系为4 2=1^1+1， k3= k2+l，k4= k3+l，这四根谱线对应的幅值分别为yi、y2、y3、y4;
[0089] 记变量 a = k-k2-0. 5,由于 0 彡 k_k2< 1，则-0? 5 彡 a 彡 〇? 5 ;
[0090] 记变量
[0091] 所述修正公式计算单元，用于计算四谱线插值算法的修正公式：
[0092] 根据公式⑷和（5)得到：
[0094] 采用多项式逼近方法，计算奇函数0 =g4(a)，表达式为：
[0095] a ~ pnX 0+p13X 03+…+plp0p (7)
[0096] pn，p13;…p lp为多项式逼近的奇次项系数，p是奇数；
[0097]根据公式（4)，求得电网信号第i次谐波的幅值Ai:
[0099] 其中，i为正整数，yi为加窗FFT后第i次谐波的幅值；
[0100] 考虑到y2、y3是离真实谱线点最近的两根谱线，基波幅值八 :为：
[0102] 根据公式（7)的逼近方法，N > 1000时，窗函数系数为实系数，公式（9)表示为：
[0103] A! = N-1 u( a )，
[0104] 其中，u(a)为修正公式，且为偶函数，逼近多项式不含奇次项；
[0105] 四谱线修正公式为：u(a) = (P2(l+P22a2+".+P2da d) (10)
[0106] 公式（10)中42(^22"12[1为多项式逼近的偶次项系数，(1为拟合的最高阶次，且(1 为偶数；
[0107] 所述基波参数计算单元，用于计算基波参数：
[0108] 计算基波频率fQ、基波幅值A1:
[0109] f0= k ? Af = (k2+a +〇. 5) Af (11)
[0110] Ai=(p20+p22a 2+…+p2da d) (12)
[0111] 根据公式⑷，得出基波的相位：
[0113]仿照基波参数的求取，根据公式（6)、（7)、（11)、（12)、（13)，进行各次谐波参数的 分析；
[0114] 考虑到其中大量窗函数的离散傅里叶分析，其表达式为：
[0116]由于N>>1，得到：
[0118] 所述四谱线插值快速计算单元，用于进行四谱线插值快速计算：
[0119] 根据公式（5)和（12)，在计算变量0和幅值心的时候，需要求出
(y4+3y 3+3y2+yi)：
[0120] 令变量：
[0121] 根据公式（4)，得到：
[0122] 根据公式（13)，得到：arg(W(k)) = -k JT (18)
[0123]将公式（17)代入公式（16)，得到：
[0124]分析得出：X〇〇与X(k3)同相位，X(k2)与X(k 4)同相位；且XG〇与X(k2)，X(k2) 与X (k3)，X (k3)与X (k4)之间的相位之差均为31，那么：
[0129]其中，C为T1的模，D为T2的模；
[0131] 其中：Re表示取实部，Im表示取虚部；
[0132] 分析得出：计算各次谐波的参数时，仅在计算幅值Ak的时候，需要进行一次求模运 算；根据公式（5)计算变量0时，利用T1和T2的实部或虚部进行快速计算；同理，根据同 一个主瓣相邻4根谱线的相位关系，求最大谱线时，直接利用插值前FFT运算结果的实部和 虚部，来寻求最大的向量；
[0133] 所述谐波参数确定单元，用于确定谐波参数：确定基波频率4后，在范围 (kfd-5, kfd+5)内，所述谱线确定单元、修正公式计算单元、基波参数计算单元、四谱线插值 快速计算单元重复进行计算，直到所有谐波参数计算完毕；
[0134] 所述误差分析单元，用于进行误差分析：在同样窗函数条件下，分析加窗四谱线插 值快速算法的误差，并与双谱线算法的误差、三谱线插值算法的误差进行比较。
[0135] 与现有技术相比，本发明的优点如下：
[0136] 本发明从电力系统的电网信号（电流信号或电压信号）加窗后的频域表达式入 手，根据窗函数主瓣内任意相邻谱线相位相差Jr的规律，推导出加窗FFT后的真实谱线点 附近最大的四根谱线之间的相位规律，提出一种加窗四谱线插值FFT快速算法。多种常用 的余弦窗函数计算实例表明，相对于双谱线和三谱线插值算法，本发明能够有效提高谐波 分析的精度；利用该快速算法，计算某次谐波仅需要1次模的运算，能够有效降低计算量和 计算时间，显著提升计算速度，具有较高的实用价值。
【附图说明】
[0137]图1是本发明实施例中基于加窗四谱线插值FFT的谐波快速分析方法的流程图。
[0138]图2是本发明实施例中相邻四条谱线的相位图。
[0139] 图3是本发明实施例中谐波信号的波形图。
[0140] 图4是Hanning窗函数的幅频特性图。
[0141] 图5是Black窗函数的幅频特性图。
[0142] 图6是Black-harris窗函数的幅频特性图。
[0143]图7是4项最大旁瓣衰减窗函数的幅频特性图。
【具体实施方式】
[0144] 下面结合附图及具体实施例对本发明作进一步的详细描述。
[0145] 参见图1所示，本发明实施例提供一种基于加窗四谱线插值FFT的谐波快速分析 方法，包括以下步骤：
[0146] S1、信号预处理：
[0147] 互感器采集电网信号，电网信号包括电流信号、电压信号，将互感器采集到的离散 电网信号x(n)，传输到上位机，n为采样点的序数，n为自然数；采用离散余弦窗函数w(n)， 对电网信号x(n)进行加窗截断，得到加窗信号\〇1):
[0148] xw (n) = x (n) w (n) (1)
[0149] 离散余弦窗函数w(n)的表达式为：
[0151] 其中，N为采样点数，N为正整数，n = 0, 1,2. ..N-1 ; 2表示求和；m为窗函数的累 加次数，m = 0, 1,2. ..M-l ;M为窗函数项数，M为正整数；\为窗函数系数，窗函数系数满 足以下约束条件：
[0153] 对公式⑴的加窗信号进行FFT变换后，得到加窗FFT频谱：
[0155] 其中，W(_)为窗函数的频谱，k为正整数，X(k)表示第k次谐波的频谱，Ak为第k 次谐波的幅值，j表示虚数单位，e是自然对数的底数，的t为第k次谐波的初始相位，第一次 谐波为基波，f s为采样频率，f 〇为基波频率，A f为离散频率间隔，且A f = f s/N。
[0156] 令：&= fV A f，&为真实频谱的谱线位置，由于非同步采样或非整周期截断，kQ 一般为小数，不为整数。
[0157] 若忽略负频率点处旁瓣的影响，公式（3)变为：
[0159] S2、确定四根谱线：
[0160] 在步骤S1得到的加窗FFT频谱峰值附近区域，匕处频率点左右各两条的谱线分 别为：第 V k2、k3、k4根，k ^ k2、k3、k4均为正整数，k ^ k2、k3、1^4的关系为：k 2= k 彳1，k3 = k2+l，k4= k 3+1，这四根谱线对应的幅值分别为yp y2、y3、y4。
[0161] 为方便计算，记变量 a = k-k2_0. 5,由于 0 < k_k2< 1，则-0? 5 < a < 〇? 5 ;
[0162] 记变量
[0163] S3、计算四谱线插值算法的修正公式：
[0164] 根据公式⑷和（5)得到：
[0166] 当N值比较大，例如：N> 1000时，公式（6)可以化简为一个函数0 = g(a)，其 反函数为a =g4(f3)。由于所采用的余弦窗系数均为实系数，其频率响应是偶对称的，因 而g( ?)和g<( ?)均为奇函数。
[0167] 采用多项式逼近方法，计算奇函数0 =g4(a)，表达式为：
[0168] a ~ pn X 0+p13X 0 3+…+plp 0 p (7)
[0169] pn，p13; "11￡)为多项式逼近的奇次项系数，p是奇数。
[0170] 根据公式（4)，求得电网信号第i次谐波的幅值Ai:
[0172] 其中，i为正整数，yi为加窗FFT后第i次谐波的幅值；
[0173] 以基波为例，考虑到y2、y3是离真实谱线点最近的两根谱线，给予较大权重，可以 得到基波幅值A 1:
[0175] 根据公式（7)的逼近方法，当N比较大，例如N > 1000时，窗函数系数为实系数， 公式（9)可表示为：
[0176] Aj = N (74+373+372+7!) u ( a ),
[0177] 其中，u(a)为修正公式，且为偶函数，逼近多项式不含奇次项。
[0178] 四谱线修正公式如下：
[0179] u ( a ) = (p2Q+p22 a 2+…+p2d a d) (10)
[0180] 公式（10)中42(^22*12[1为多项式逼近的偶次项系数，(1为拟合的最高阶次，且(1 为偶数。
[0181] S4、计算基波参数：
[0182] 考虑到y2、y3是离真实谱线点最近的两根谱线，给予较大权重，可以得到基波频率 fci、基波幅值A ::
[0183] f0= k ? A f = (k 2+ a +〇. 5) A f (11)
[0184] A! = N-1 (p2CI+p22 a 2+... +p2da d) (12)
[0185] 根据公式（4)，还可以得出基波的相位：
[0187] 仿照基波参数的求取，根据公式（6)、（7)、（11)、（12)、（13)，即可进行各次谐波参 数的分析。
[0188] 考虑到其中大量窗函数的离散傅里叶分析，其表达式为：
[0190] 由于N>>1，可以得到：
[0192] S5、加窗四谱线插值快速算法：
[0193] 根据公式（5)和（12)，在计算变量0和幅值&的时候，需要求出
、（y4+3y3+3y2+yi)。
[0194] 令变量：
[0195] 根据公式（4)可以得到：
[0196] 根据公式（13)可得到：arg(W(k)) = -k JT (18)
[0197]将公式（17)代入公式（16)可得到：
[0198] 分析得出：X〇〇与X(k3)同相位，X(k2)与X(k 4)同相位；且XG〇与X(k2)，X(k2) 与X(k3)，X(k 3)与X(k4)之间的相位之差均为31，相邻四条谱线的相位图参见图2所示，那 么：
[0201] 由图2可知：
[0202] 因此，可以再次令：
[0204]其中，C为T1的模，D为T2的模。
[0206] 其中：Re表示取实部，Im表示取虚部。
[0207] 由以上分析得出：计算各次谐波的参数时候，仅在计算幅值Ak的时候，需要进行一 次求模运算。根据公式（5)计算变量0时，可以利用T1和T2的实部或虚部进行快速计算。[0208] 同理，根据同一个主瓣相邻4根谱线的相位关系，求最大谱线时，直接利用插值前 FFT运算结果的实部和虚部，来寻求最大的向量。
[0209] S6、确定谐波参数：确定基波频率&后，在范围（kf ^5, kf^+5)内，重复步骤S2~ S5,直到所有谐波参数计算完毕。
[0210] S7、进行误差分析：在同样窗函数条件下，分析加窗四谱线插值快速算法的误差， 并与双谱线算法的误差、三谱线插值算法的误差进行比较。
[0211] Sl、S2、S3、S4、S6、S7形成一个完整的基于加窗四谱线插值FFT的高精度谐波分 析方法，能够有效提高谐波分析的精度。
[0212] 步骤S5的作用是提高计算速度，在S4、S6之间增加步骤S5,每次谐波仅需一次模 运算，计算量小，能够有效提升计算速度，在同样窗函数下，能获得更高的精度，具有较高的 实用价值。
[0213] 下面通过一个具体实施例来进行详细说明。选定一个电网信号的谐波，谐波的波 形参见图3所示。
[0214] 步骤1、信号的预处理：
[0215] 将互感器采集到的离散电网信号x(n)，信号模型参见表1所示。为了验证所提算 法的精度，进行10次谐波仿真分析。信号模型为：
[0217] 其 中：基波频率fQ为50. 5Hz，采样频率匕为5120Hz，采样点数N为1024。
[0218] 表1、电网信号的谐波参数
[0220]
[0221] 对信号x(n)进行加窗处理，选取4种常用加窗函数，其系数参见表2所示。
[0222] 表2、常用窗函数系数
[0224] 加窗FFT变换后得到信号频谱：
[0225] 步骤2、确定四根谱线：
[0226] 在45~55Hz中寻求实部（或虚部）最大四根谱线，分别为XGO、X(k 2)、X(k3)、 X (k4)〇令：
[0228] 步骤3、计算四谱线插值算法的修正公式：
[0229] 根据【具体实施方式】中步骤S3的推导，表2中的4种窗函数的修正公式a = g4(f3)，u(a)分别如下：
[0230] (1) Hanning 窗：
[0231] a = 1. 13013682 0-0. 18408680 03+〇, 07224859 05-〇. 04451832 07u( a)= 0. 53549869+0. 17622103 a 2+〇. 09310826 a 4+〇. 06644946 a 6
[0232]Hanning窗函数的幅频特性参见图4所示。
[0233] (2) Balckman 窗：
[0234] a = 1.44649012 0-0. 29326578 0 3+0. 13858447 05-〇. O9578617 07u( a)= 1. 1575129+0. 56110888 a 2+〇. 38707997 a 4+〇. 33290703 a 6
[0235]Black窗函数的幅频特性参见图5所示。
[0236] (3)Balckman-harris 窗：
[0237] a = 1. 85862073 0-0. 40299738 03+O. 19725715 05-〇. 13293118 07u( a)= 1. 24914356+0. 88776025 a 2+〇. 73879907 a 4-〇. 69429366 a 6
[0238] Black-harris窗函数的幅频特性参见图6所示。
[0239] (4) 4项最大旁瓣衰减窗：
[0240]a = 2. 37540983 0-0. 43478665 03+〇, 19124517 05-〇. 11511086 07u( a)= 1. 34456750+1. 36577423 a 2+l. 25666570 a 4+l. 17989176 a 6
[0241] 4项最大旁瓣衰减窗函数的幅频特性参见图7所示。
[0242] 步骤4、计算基波参数：
[0243] Aj=N (74+373+372+7!) u( a )
[0245] f〇= k ? A f = (k 2+ a +〇. 5) A f
[0246] 步骤5、四谱线插值快速算法：
[0247] 根据S5所述快速算法，可以对步骤2和步骤4中的两个量的求取采用快速算法：
[0248]
及Af N-1|T2|u(a)
[0249] 步骤6、确定谐波参数：确定基波频率f0后，在范围（kf^-5, kf^+5)内重复步骤2~ 5,直到所有谐波参数计算完毕。
[0250] 步骤7、进行误差分析：仿真实验的误差分析比较结果参见表3~表6所示，其中， DAi表示基波和各次谐波幅值测量值的相对误差；D f(l表示基波频率测量值的相对误差； 表示基波和各次谐波初始相位测量值的相对误差，均用百分比表示。
[0251] 表3、Hanning窗和Blackman窗的频率、幅值相对误差
[0253] 表4、Hanning窗和Blackman窗的相位相对误差
[0254]
[0255] 表5、Blackman-harris窗和四项最大旁瓣衰减窗的频率、幅值相对误差
[0257] 表6、Blackman-harris窗和四项最大旁瓣衰减窗的相位相对误差
[0259] 从表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)的表达式为：
[0265] 其中，N为采样点数，N为正整数，n = 0, 1，2. . . N-1 ; 2表示求和；m为窗函数的累 加次数，m = 0, 1,2. ..M-l ;M为窗函数项数，M为正整数；\为窗函数系数，窗函数系数满 足以下约束条件：
[0267] 对公式⑴的加窗信号进行FFT变换后，得到加窗FFT频谱：
[0269] 其中，W(_)为窗函数的频谱，k为正整数，X(k)表示第k次谐波的频谱，Ak为第k 次谐波的幅值，j表示虚数单位，e是自然对数的底数，为第k次谐波的初始相位，第一次 谐波为基波，f s为采样频率，f 〇为基波频率，A f为离散频率间隔，且A f = f s/N。
[0270] 令：ke fV A f，&为真实频谱的谱线位置，由于非同步采样或非整周期截断，h 一般不为整数。
[0271] 若忽略负频率点处旁瓣的影响，公式（3)变为：
[0273] 谱线确定单元，用于确定四根谱线：
[0274] 在信号预处理单元得到的加窗FFT频谱峰值附近区域，匕处频率点左右各两条的 谱线分别为：第1^ 1、1^2、1^、1^4根，1^1、1^2、1^、1^ 4均为正整数，1^1、1^2、1^、1^4的关系为4 2=1^1+1， k3= k2+l，k4= k3+l，这四根谱线对应的幅值分别为yi、y2、y3、y4。
[0275] 为方便计算，记变量a = k-k2-0. 5,由于0彡k_k2< 1，则-0? 5彡a彡〇? 5;
[0276] 记变量
[0277] 修正公式计算单元，用于计算四谱线插值算法的修正公式：
[0278] 根据公式⑷和（5)得到：
[0280]当N值比较大，例如：N> 1000时，公式（6)可以化简为一个函数0 = g( a)，其 反函数为a =g4(f3)。由于所采用的余弦窗系数均为实系数，其频率响应是偶对称的，因 而g( ?)和g<( ?)均为奇函数。
[0281] 采用多项式逼近方法，计算奇函数0 =g4(a)，表达式为：
[0282] a ~ pnX 0+p13X 03+…+plp0p (7)
[0283] pn，p13; "11￡)为多项式逼近的奇次项系数，p是奇数。
[0284]根据公式（4)，求得电网信号第i次谐波的幅值Ai:
[0286] 其中，i为正整数，yi为加窗FFT后第i次谐波的幅值；
[0287] 以基波为例，考虑到y2、y3是离真实谱线点最近的两根谱线，给予较大权重，可以 得到基波幅值A 1:
[0289] 根据公式（7)的逼近方法，当N比较大，例如N > 1000时，窗函数系数为实系数， 公式（9)可表示为：
[0290] Aj = N (74+373+372+7!) u ( a ),
[0291] 其中，u(a)为修正公式，且为偶函数，逼近多项式不含奇次项。
[0292] 四谱线修正公式如下：
[0293] u ( a ) = (p2〇+p22 a 2+…+p2d a d) (10)
[0294] 公式（10)中，p2(l，p22- p2$多项式逼近的偶次项系数，d为拟合的最高阶次，且d 为偶数。
[0295] 基波参数计算单元，用于计算基波参数：
[0296] 考虑到y2、y3是离真实谱线点最近的两根谱线，给予较大权重，可以得到基波频率 fci、基波幅值A ::
[0297] f0= k ? A f = (k 2+ a +〇. 5) A f (11)
[0298] Ai = N-1 (yWys+Sy^yD (p2CI+p22 a 2+... +p2da d) (12)
[0299] 根据公式（4)还可以得出基波的相位：
[0300]
[0301] 仿照基波参数的求取，根据公式（6)、（7)、（11)、（12)、（13)，即可进行各次谐波参 数的分析。
[0302] 考虑到其中大量窗函数的离散傅里叶分析，其表达式为：
[0304]由于N>>1，可以得到：
[0306] 四谱线插值快速计算单元，用于进行四谱线插值快速计算：
[0307] 由公式（5)和（12)可知，在计算变量0和幅值&的时候，需要求出
[0308] 令变量：
[0309] 根据公式（4)可以得到
[0310] 根据公式（13)可得到：arg(W(k)) = -k JT (18)
[0311] 将公式（17)代入公式（16)可得到
[0312] 分析得出：X〇〇与X(k3)同相位，X(k2)与X(k 4)同相位；且XG〇与X(k2)，X(k2) 与X(k3)，X(k 3)与X(k4)之间的相位之差均为31，相邻四条谱线的相位图参见图2所示，那 么：
[0314]
[0315] 由图2可知：
[0316] 因此，可以再次令：
[0318] 其中，C为T1的模，D为T2的模。
[0319]贝1J :
[0320] 其中：Re表示取实部，Im表示取虚部。
[0321] 由以上分析可知，计算各次谐波的参数时候，仅在计算幅值Ak的时候，需要进行一 次求模运算。根据公式（5)计算变量0时，可以利用T1和T2的实部或虚部进行快速计算。 同理，根据同一个主瓣相邻4根谱线的相位关系，求最大谱线时，直接利用插值前FFT运算 结果的实部和虚部，来寻求最大的向量。
[0322] 谐波参数确定单元，用于确定谐波参数：确定基波频率&后，在范围 (kfd-5, kfd+5)内，谱线确定单元、修正公式计算单元、基波参数计算单元、四谱线插值快速 计算单元重复进行计算，直到所有谐波参数计算完毕。
[0323] 误差分析单元，用于进行误差分析：在同样窗函数条件下，分析加窗四谱线插值快 速算法的误差，并与双谱线算法的误差、三谱线插值算法的误差进行比较。
[0324] 本领域的技术人员可以对本发明实施例进行各种修改 和变型，倘若这些修改和变 型在本发明权利要求及其等同技术的范围之内，则这些修改和变型也在本发明的保护范围 之内。
[0325] 说明书中未详细描述的内容为本领域技术人员公知的现有技术。
【主权项】
1. 一种基于加窗四谱线插值FFT的谐波快速分析方法，其特征在于，包括以下步骤： 51、 信号预处理： 互感器采集电网信号，将互感器采集到的离散电网信号x(n)，传输到上位机，η为采样 点的序数，η为自然数；采用离散余弦窗函数w (η)，对电网信号X (η)进行加窗截断，得到加 窗信号xw (η): xw (η) = X (n) w (η) (I) 离散余弦窗函数w (η)的表达式为：(2) 其中，N为采样点数，N为正整数，η = 0, 1,2. ..N-I ;Σ表示求和；m为窗函数的累加次 数，m = 0, 1，2. .. M-I ;M为窗函数项数，M为正整数；Iddi为窗函数系数； 对公式（1)的加窗信号进行FFT变换后，得到加窗FFT频谱：其中，W( ·)为窗函数的频谱，k为正整数，X(k)表示第k次谐波的频谱，Ak为第k次 谐波的幅值，j表示虚数单位，e是自然对数的底数，为第k次谐波的初始相位，第一次谐 波为基波，fs为采样频率，f 〇为基波频率，Δ f为离散频率间隔，且Λ f = f s/N ; 令：1?= fd/Δ f，1?为真实频谱的谱线位置； 忽略负频率点处旁瓣的影响，公式（3)变为：(4) 52、 确定四根谱线： 在步骤Sl得到的加窗FFT频谱峰值附近区域，1?处频率点左右各两条的谱线分别为： 第1^、1^2、1^、1^4根，1^1、1^2、1^、1^ 4均为正整数，1^1、1^2、1^、1^4的关系为 ：1^2=1^1+1，1^=1^2+1，1^4=k3+l，这四根谱线对应的幅值分别为y2、y3、y 4; 记变量 α = k-k2-0. 5,由于 0 彡 k-k2< 1，则-0· 5 彡 α 彡 〇· 5 ; 记变量（5 ) 53、 计算四谱线插值算法的修正公式： 根据公式（4)和（5)得到：采用多项式逼近方法，计算奇函数β =Ρ(α)，表达式为：(7) Pn，P13;…P 1P为多项式逼近的奇次项系数，P是奇数； 根据公式（4)，求得电网信号第i次谐波的幅值Ai:(8) 其中，i为正整数，Yi为加窗FFT后第i次谐波的幅值； 考虑到y2、y3是离真实谱线点最近的两根谱线，基波幅值A :根据公式（7)的逼近方法，N > 1000时，窗函数系数为实系数，公式（9)表示为： A1 = N (74+373+372+7!) u ( α ), 其中，u(a)为修正公式，且为偶函数，逼近多项式不含奇次项； 四谱线修正公式为：u( a ) = (ρ2(ι+ρ22 α 2+…+p2da d) (10) 公式（10)中，P2(I，P22··· P2A多项式逼近的偶次项系数，d为拟合的最高阶次，且d为偶 数； 54、 计算基波参数： 计算基波频率fci、基波幅值A1: f0= k · Af= (k 2+ a +〇. 5) Δ f (11) A1= N (74+373+372+7!) (p2〇+p22 a 2+··· +p2d a d) (12) 根据公式（4)，得出基波的相位：仿照基波参数的求取，根据公式（6)、（7)、（11)、（12)、（13)，进行各次谐波参数的分 析； 考虑到其中大量窗函数的离散傅里叶分析，其表达式为：由于N>>1，得到：55、 加窗四谱线插值快速算法： 根据公式（5)和（12)，计算变量β和幅值A1的时候，需求出根据公式（4)，得到 根据公式（13)，得到 将公式（17)代入公式（16)，得到分析得出：XG〇与X(k3)同相位，X(k2)与X(k4)同相位；且X(Ic 1)与X(k2)，X(k2)与 X(k3)，X(k3)与X(k4)之间的相位之差均为π，那么：其中，C为Tl的模，D为Τ2的模； 则：其中：Re表示取实部，Im表示取虚部； 分析得出：计算各次谐波的参数时，仅在计算幅值Ak的时候，需要进行一次求模运算； 根据公式（5)计算变量β时，利用Tl和T2的实部或虚部进行快速计算；同理，根据同一个 主瓣相邻4根谱线的相位关系，求最大谱线时，直接利用插值前FFT运算结果的实部和虚部 来寻求最大的向量； S6、确定谐波参数：确定基波频率&后，在范围（kf ^5, kf^+5)内重复步骤S2~S5,直 到所有谐波参数计算完毕； S7、进行误差分析：在同样窗函数条件下，分析加窗四谱线插值快速算法的误差，并与 双谱线算法的误差、三谱线插值算法的误差进行比较。2. 如权利要求1所述的基于加窗四谱线插值FFT的谐波快速分析方法，其特征在于： 所述电网信号包括电流信号、电压信号。3. 如权利要求1所述的基于加窗四谱线插值FFT的谐波快速分析方法，其特征在于： 所述真实频谱的谱线位置1?为小数。4. 如权利要求1或2或3所述的基于加窗四谱线插值FFT的谐波快速分析方法，其特 征在于： 所述窗函数系数\满足以下约束条件：5. -种基于加窗四谱线插值FFT的谐波快速分析系统，其特征在于：包括信号预处理 单元、谱线确定单元、修正公式计算单元、基波参数计算单元、四谱线插值快速计算单元、谐 波参数确定单元、误差分析单元，其中： 所述信号预处理单元，用于对信号进行预处理： 互感器采集电网信号，将互感器采集到的离散电网信号x(n)，传输到上位机，η为采样 点的序数，η为自然数；采用离散余弦窗函数w (η)，对电网信号X (η)进行加窗截断，得到加 窗信号xw (η): xw (η) = X (n) w (η) (I) 离散余弦窗函数w (η)的表达式为：(2) 其中，N为采样点数，N为正整数，η = 0, 1,2. ..N-I ;Σ表示求和；m为窗函数的累加次 数，m = 0, 1，2. .. M-I ;M为窗函数项数，M为正整数；Iddi为窗函数系数； 对公式（1)的加窗信号进行FFT变换后，得到加窗FFT频谱：其中，W( ·)为窗函数的频谱，k为正整数，X(k)表示第k次谐波的频谱，Ak为第k次 谐波的幅值，j表示虚数单位，e是自然对数的底数，PA-为第k次谐波的初始相位，第一次谐 波为基波，fs为采样频率，f 〇为基波频率，Δ f为离散频率间隔，且Λ f = f s/N ; 令：1?= fd/Δ f，1?为真实频谱的谱线位置； 忽略负频率点处旁瓣的影响，公式（3)变为：(4) 所述谱线确定单元，用于确定四根谱线： 在信号预处理单元得到的加窗FFT频谱峰值附近区域，1?处频率点左右各两条的谱线 分别为：第 1^1、1^2、1^、1^4根，1^1、1^ 2、1^、1^4均为正整数，1^1、1^2、1^、1^ 4的关系为：1^2=1^1+1，1^ = k2+l，k4= k 3+1，这四根谱线对应的幅值分别为yp y2、y3、y4; 记变量 α = k-k2-0. 5,由于 O 彡 k-k2< 1，则-0· 5 彡 α 彡 〇· 5 ; 记变量（5) 所述修正公式计算单元，用于计算四谱线插值算法的修正公式： 根据公式（4)和（5)得到：采用多项式逼近方法，计算奇函数β =Ρ(α)，表达式为：(7) Pn，P13;…P 1Ρ为多项式逼近的奇次项系数，P是奇数； 根据公式（4)，求得电网信号第i次谐波的幅值Ai:(8) 其中，i为正整数，Yi为加窗FFT后第i次谐波的幅值； 考虑到y2、y3是离真实谱线点最近的两根谱线，基波幅值A :根据公式（7)的逼近方法，N > 1000时，窗函数系数为实系数，公式（9)表示为： A1 = N (74+373+372+7!) u ( α ), 其中，u(a)为修正公式，且为偶函数，逼近多项式不含奇次项； 四谱线修正公式为：u( a ) = (ρ2(ι+ρ22 α 2+…+p2da d) (10) 公式（10)中，P2(I，P22··· P2A多项式逼近的偶次项系数，d为拟合的最高阶次，且d为偶 数； 所述基波参数计算单元，用于计算基波参数： 计算基波频率fci、基波幅值A1: f0= k · Af= (k 2+ a +〇. 5) Δ f (11) A1= N (74+373+372+7!) (p2〇+p22 a 2+··· +p2d a d) (12) 根据公式（4)，得出基波的相位：仿照基波参数的求取，根据公式（6)、（7)、（11)、（12)、（13)，进行各次谐波参数的分 析； 考虑到其中大量窗函数的离散傅里叶分析，其表达式为：由于N>>1，得到：所述四谱线插值快速计算单元，用于进行四谱线插值快速计算： 根据公式（5)和（12)，在计算变量β和幅值A1的时候，需要求出、（y4+3y3 +3y2+yi): 令变量：根据公式（4)，得到 根据公式（13)，得到 将公式（17)代入公式（16)，得到：分析得出：XG〇与X(k3)同相位，X(k2)与X(k4)同相位；且X(Ic 1)与X(k2)，X(k2)与 X(k3)，X(k3)与X(k4)之间的相位之差均为π，那么：令：（23) 其中，C为Tl的模，D为T2的模； 则：(24) 其中：Re表示取实部，Im表示取虚部； 分析得出：计算各次谐波的参数时，仅在计算幅值Ak的时候，需要进行一次求模运算； 根据公式（5)计算变量β时，利用Tl和T2的实部或虚部进行快速计算；同理，根据同一 个主瓣相邻4根谱线的相位关系，求最大谱线时，直接利用插值前FFT运算结果的实部和虚 部，来寻求最大的向量； 所述谐波参数确定单元，用于确定谐波参数：确定基波频率后，在范围 (kfd-5, kfd+5)内，所述谱线确定单元、修正公式计算单元、基波参数计算单元、四谱线插值 快速计算单元重复进行计算，直到所有谐波参数计算完毕； 所述误差分析单元，用于进行误差分析：在同样窗函数条件下，分析加窗四谱线插值快 速算法的误差，并与双谱线算法的误差、三谱线插值算法的误差进行比较。6. 如权利要求5所述的基于加窗四谱线插值FFT的谐波快速分析系统，其特征在于： 所述电网信号包括电流信号、电压信号。7. 如权利要求5所述的基于加窗四谱线插值FFT的谐波快速分析系统，其特征在于： 所述真实频谱的谱线位置1?为小数。8. 如权利要求5或6或7所述的基于加窗四谱线插值FFT的谐波快速分析系统，其特 征在于： 所述窗函数系数\满足以下约束条件：
【专利摘要】本发明公开了一种基于加窗四谱线插值FFT的谐波快速分析方法及系统，涉及谐波分析领域。该方法包括以下步骤：信号预处理；确定四根谱线；计算四谱线插值算法的修正公式；计算基波参数；四谱线插值快速算法；确定谐波参数；进行误差分析。本发明从电力系统的电网信号加窗后的频域表达式入手，根据窗函数主瓣内任意相邻谱线相位相差π的规律，推导出加窗FFT后的真实谱线点附近最大的四根谱线之间的相位规律，提出一种加窗四谱线插值FFT快速算法。相对于双谱线和三谱线插值算法，本发明能有效提高谐波分析的精度；利用该快速算法，计算某次谐波仅需要1次模的运算，能够有效降低计算量和计算时间，显著提升计算速度。
【IPC分类】G01R23/16
【公开号】CN104897960
【申请号】CN201510326063
【发明人】张俊敏, 刘开培, 汪立, 王黎, 田微, 何顺帆, 姚为
【申请人】中南民族大学
【公开日】2015年9月9日
【申请日】2015年6月15日
【技术领域】
[0001] 本发明涉及谐波分析领域，具体是涉及一种基于加窗四谱线插值FFT的谐波快速 分析方法及系统。
【背景技术】
[0002] 随着大量非线性电力电子器件的使用，使得电网中谐波污染问题日益严重。谐波 问题不仅恶化电能质量，对电网的安全稳定和经济运行也造成较大影响。因此，对系统中谐 波参数进行高精度测量，对于减少谐波危害，维护电网安全稳定、高效运行是十分必要的。
[0003] FFT(Fast Fourier Transform，快速傅里叶变换）是目前用于电力系统谐波分析 最常用的算法，但采用FFT对信号进行频谱分析很难做到同步采样和整周期截断，所带来 的频谱泄露和栅栏效应会导致谐波的频率、幅值和相位等参数测量不准，无法满足测量要 求。
[0004] 为了解决这一问题，加窗插值FFT算法被广泛应用：运用各种特殊窗函数对信号 进行截断，然后结合谱线插值FFT进行谐波分析，从而提高测量精度。
[0005] 常用的窗函数有：汉宁（Hanning)窗函数、布莱克曼（Blackman)窗函数、布莱克曼 汉斯（Blackman-Harris)窗函数、纳托尔（Nuttall)窗函数、莱夫文森特（Rife-Vincent) 窗函数以及各种组合窗函数。
[0006] 在加窗基础上，D. Agrez和庞浩等人各自提出了双谱线的修正算法，Wu Jing、牛胜 锁和黄冬梅等人提出了三谱线修正算法。这些改进降低了频谱泄漏和栅栏效应的影响，提 高了谐波分析的准确性。
[0007] 然而，在工程实际使用中，双谱线和三谱线插值算法仍然无法满足高精度的谐波 分析要求，并且随着所采用谱线数目增多，求模的计算量也迅速增加，导致计算速度较慢。
【发明内容】
[0008] 本发明的目的是为了克服上述【背景技术】的不足，提供一种基于加窗四谱线插值 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)
[0013] 离散余弦窗函数w(n)的表达式为：
[0015] 其中，N为采样点数，N为正整数，n = 0, 1，2. . . N-l ; 2表示求和；m为窗函数的累 加次数，m = 0, 1，2.. .M-l ;M为窗函数项数，M为正整数；\为窗函数系数；
[0016] 对公式⑴的加窗信号进行FFT变换后，得到加窗FFT频谱：
[0018] 其中，W(_)为窗函数的频谱，k为正整数，X(k)表示第k次谐波的频谱，Ak为第k 次谐波的幅值，j表示虚数单位，e是自然对数的底数，P是为第k次谐波的初始相位，第一次 谐波为基波，f s为采样频率，f〇为基波频率，A f为离散频率间隔，且A f = f s/N ;
[0019]令：kQ= f Q/ A f，kQ为真实频谱的谱线位置；
[0020] 忽略负频率点处旁瓣的影响，公式（3)变为：
[0022] S2、确定四根谱线：
[0023] 在步骤S1得到的加窗FFT频谱峰值附近区域，匕处频率点左右各两条的谱线分 别为：第 V k2、k3、k4根，k^k2、k3、k4均为正整数，k^k2、k3、1^4的关系为：k 2= k 汴1，k3 = k2+l，k4= k 3+1，这四根谱线对应的幅值分别为yp y2、y3、y4;
[0024] 记变量 a = k-k2-0. 5,由于 0 彡 k_k2< 1，则-0? 5 彡 a 彡 〇? 5 ;
[0025]记变量，
[0026] S3、计算四谱线插值算法的修正公式：
[0027] 根据公式⑷和（5)得到：
[0029] 采用多项式逼近方法，计算奇函数0 =g4(a)，表达式为：
[0030] a ~ pnX 0+p13X 03+…+plp0p (7)
[0031] pn，p13;…p lp为多项式逼近的奇次项系数，p是奇数；
[0032]根据公式（4)，求得电网信号第i次谐波的幅值Ai:
[0034] 其中，i为正整数，yi为加窗FFT后第i次谐波的幅值；
[0035]考虑到y2、y3是离真实谱线点最近的两根谱线，基波幅值八 :为：
[0036]
[0037] 根据公式（7)的逼近方法，N> 1000时，窗函数系数为实系数，公式（9)表示为：
[0038] Aj = N (74+373+372+7!) u ( a ),
[0039] 其中，u(a)为修正公式，且为偶函数，逼近多项式不含奇次项；
[0040] 四谱线修正公式为：u(a ) = (p2Q+p22a 2+."+p2da d) (10)
[0041] 公式（10)中，p2(l，p22…p2d为多项式逼近的偶次项系数，d为拟合的最高阶次，且d 为偶数；
[0042] S4、计算基波参数：
[0043] 计算基波频率f(i、基波幅值Ay
[0044] f0= k ? A f = (k 2+ a +〇. 5) A f (11)
[0045] Af N-1 (p2CI+p22 a 2+... +p2da d) (12)
[0046] 根据公式⑷，得出基波的相位：
[0048]仿照基波参数的求取，根据公式（6)、（7)、（11)、（12)、（13)，进行各次谐波参数的 分析；
[0049] 考虑到其中大量窗函数的离散傅里叶分析，其表达式为：
[0051]由于N>>1，得到：
[0053] S5、加窗四谱线插值快速算法：
[0054] 根据公式（5)和（12)，计算变量0和幅值Ai的时候，需求出
.(y4+3y3 +3y2+yi):
[0055] 令变量：
[0056] 根据公式⑷，得到：
[0057] 根据公式（13)，得到：arg(W(k)) = -k JT (18)
[0058]将公式（17)代入公式（16)，得到：
[0059] 分析得出：X〇〇与X(k3)同相位，X(k 2)与X(k4)同相位；且XG〇与X(k2)，X(k2) 与X (k3)，X (k3)与X (k4)之间的相位之差均为31，那么：
[0065]其中，C为T1的模，D为T2的模；
[0067] 其中：Re表示取实部，Im表示取虚部；
[0068] 分析得出：计算各次谐波的参数时，仅在计算幅值Ak的时候，需要进行一次求模运 算；根据公式（5)计算变量0时，利用T1和T2的实部或虚部进行快速计算；同理，根据同 一个主瓣相邻4根谱线的相位关系，求最大谱线时，直接利用插值前FFT运算结果的实部和 虚部来寻求最大的向量；
[0069] S6、确定谐波参数：确定基波频率fQ后，在范围（kf。_5, kfQ+5)内重复步骤S2~ S5,直到所有谐波参数计算完毕；
[0070] S7、进行误差分析：在同样窗函数条件下，分析加窗四谱线插值快速算法的误差， 并与双谱线算法的误差、三谱线插值算法的误差进行比较。
[0071] 在上述技术方案的基础上，所述电网信号包括电流信号、电压信号；所述真实频谱 的谱线位置h为小数。
[0072] 在上述技术方案的基础上，所述窗函数系数bm满足以下约束条件：
[0074] 本发明还提供一种基于加窗四谱线插值FFT的谐波快速分析系统，包括信号预处 理单元、谱线确定单元、修正公式计算单元、基波参数计算单元、四谱线插值快速计算单元、 谐波参数确定单元、误差分析单元，其中：
[0075] 所述信号预处理单元，用于对信号进行预处理：
[0076] 互感器采集电网信号，将互感器采集到的离散电网信号x(n)，传输到上位机，n为 采样点的序数，n为自然数；采用离散余弦窗函数w (n)，对电网信号x (n)进行加窗截断，得 到加窗信号xw(n):
[0077] xw(n) = x(n)w(n) (1)
[0078] 离散余弦窗函数w(n)的表达式为：
[0080] 其中，N为采样点数，N为正整数，n = 0, 1，2. . . N-1 ; 2表示求和；m为窗函数的累 加次数，m = 0, 1，2. . . M-l ;M为窗函数项数，M为正整数；\为窗函数系数；
[0081] 对公式⑴的加窗信号进行FFT变换后，得到加窗FFT频谱：
[0083] 其中，W( ?)为窗函数的频谱，k为正整数，X(k)表示第k次谐波的频谱，Ak为第 k次谐波的幅值，j表示虚数单位，e是自然对数的底数，为第k次谐波的初始相位，第一 次谐波为基波，f s为采样频率，f 〇为基波频率，A f为离散频率间隔，且A f = f s/N ;
[0084] 令：kQ= f Q/ A f，kQ为真实频谱的谱线位置；
[0085] 忽略负频率点处旁瓣的影响，公式（3)变为：
[0087] 所述谱线确定单元，用于确定四根谱线：
[0088] 在信号预处理单元得到的加窗FFT频谱峰值附近区域，匕处频率点左右各两条的 谱线分别为：第1^ 1、1^2、1^、1^4根，1^1、1^2、1^、1^ 4均为正整数，1^1、1^2、1^、1^4的关系为4 2=1^1+1， k3= k2+l，k4= k3+l，这四根谱线对应的幅值分别为yi、y2、y3、y4;
[0089] 记变量 a = k-k2-0. 5,由于 0 彡 k_k2< 1，则-0? 5 彡 a 彡 〇? 5 ;
[0090] 记变量
[0091] 所述修正公式计算单元，用于计算四谱线插值算法的修正公式：
[0092] 根据公式⑷和（5)得到：
[0094] 采用多项式逼近方法，计算奇函数0 =g4(a)，表达式为：
[0095] a ~ pnX 0+p13X 03+…+plp0p (7)
[0096] pn，p13;…p lp为多项式逼近的奇次项系数，p是奇数；
[0097]根据公式（4)，求得电网信号第i次谐波的幅值Ai:
[0099] 其中，i为正整数，yi为加窗FFT后第i次谐波的幅值；
[0100] 考虑到y2、y3是离真实谱线点最近的两根谱线，基波幅值八 :为：
[0102] 根据公式（7)的逼近方法，N > 1000时，窗函数系数为实系数，公式（9)表示为：
[0103] A! = N-1 u( a )，
[0104] 其中，u(a)为修正公式，且为偶函数，逼近多项式不含奇次项；
[0105] 四谱线修正公式为：u(a) = (P2(l+P22a2+".+P2da d) (10)
[0106] 公式（10)中42(^22"12[1为多项式逼近的偶次项系数，(1为拟合的最高阶次，且(1 为偶数；
[0107] 所述基波参数计算单元，用于计算基波参数：
[0108] 计算基波频率fQ、基波幅值A1:
[0109] f0= k ? Af = (k2+a +〇. 5) Af (11)
[0110] Ai=(p20+p22a 2+…+p2da d) (12)
[0111] 根据公式⑷，得出基波的相位：
[0113]仿照基波参数的求取，根据公式（6)、（7)、（11)、（12)、（13)，进行各次谐波参数的 分析；
[0114] 考虑到其中大量窗函数的离散傅里叶分析，其表达式为：
[0116]由于N>>1，得到：
[0118] 所述四谱线插值快速计算单元，用于进行四谱线插值快速计算：
[0119] 根据公式（5)和（12)，在计算变量0和幅值心的时候，需要求出
(y4+3y 3+3y2+yi)：
[0120] 令变量：
[0121] 根据公式（4)，得到：
[0122] 根据公式（13)，得到：arg(W(k)) = -k JT (18)
[0123]将公式（17)代入公式（16)，得到：
[0124]分析得出：X〇〇与X(k3)同相位，X(k2)与X(k 4)同相位；且XG〇与X(k2)，X(k2) 与X (k3)，X (k3)与X (k4)之间的相位之差均为31，那么：
[0129]其中，C为T1的模，D为T2的模；
[0131] 其中：Re表示取实部，Im表示取虚部；
[0132] 分析得出：计算各次谐波的参数时，仅在计算幅值Ak的时候，需要进行一次求模运 算；根据公式（5)计算变量0时，利用T1和T2的实部或虚部进行快速计算；同理，根据同 一个主瓣相邻4根谱线的相位关系，求最大谱线时，直接利用插值前FFT运算结果的实部和 虚部，来寻求最大的向量；
[0133] 所述谐波参数确定单元，用于确定谐波参数：确定基波频率4后，在范围 (kfd-5, kfd+5)内，所述谱线确定单元、修正公式计算单元、基波参数计算单元、四谱线插值 快速计算单元重复进行计算，直到所有谐波参数计算完毕；
[0134] 所述误差分析单元，用于进行误差分析：在同样窗函数条件下，分析加窗四谱线插 值快速算法的误差，并与双谱线算法的误差、三谱线插值算法的误差进行比较。
[0135] 与现有技术相比，本发明的优点如下：
[0136] 本发明从电力系统的电网信号（电流信号或电压信号）加窗后的频域表达式入 手，根据窗函数主瓣内任意相邻谱线相位相差Jr的规律，推导出加窗FFT后的真实谱线点 附近最大的四根谱线之间的相位规律，提出一种加窗四谱线插值FFT快速算法。多种常用 的余弦窗函数计算实例表明，相对于双谱线和三谱线插值算法，本发明能够有效提高谐波 分析的精度；利用该快速算法，计算某次谐波仅需要1次模的运算，能够有效降低计算量和 计算时间，显著提升计算速度，具有较高的实用价值。
【附图说明】
[0137]图1是本发明实施例中基于加窗四谱线插值FFT的谐波快速分析方法的流程图。
[0138]图2是本发明实施例中相邻四条谱线的相位图。
[0139] 图3是本发明实施例中谐波信号的波形图。
[0140] 图4是Hanning窗函数的幅频特性图。
[0141] 图5是Black窗函数的幅频特性图。
[0142] 图6是Black-harris窗函数的幅频特性图。
[0143]图7是4项最大旁瓣衰减窗函数的幅频特性图。
【具体实施方式】
[0144] 下面结合附图及具体实施例对本发明作进一步的详细描述。
[0145] 参见图1所示，本发明实施例提供一种基于加窗四谱线插值FFT的谐波快速分析 方法，包括以下步骤：
[0146] S1、信号预处理：
[0147] 互感器采集电网信号，电网信号包括电流信号、电压信号，将互感器采集到的离散 电网信号x(n)，传输到上位机，n为采样点的序数，n为自然数；采用离散余弦窗函数w(n)， 对电网信号x(n)进行加窗截断，得到加窗信号\〇1):
[0148] xw (n) = x (n) w (n) (1)
[0149] 离散余弦窗函数w(n)的表达式为：
[0151] 其中，N为采样点数，N为正整数，n = 0, 1,2. ..N-1 ; 2表示求和；m为窗函数的累 加次数，m = 0, 1,2. ..M-l ;M为窗函数项数，M为正整数；\为窗函数系数，窗函数系数满 足以下约束条件：
[0153] 对公式⑴的加窗信号进行FFT变换后，得到加窗FFT频谱：
[0155] 其中，W(_)为窗函数的频谱，k为正整数，X(k)表示第k次谐波的频谱，Ak为第k 次谐波的幅值，j表示虚数单位，e是自然对数的底数，的t为第k次谐波的初始相位，第一次 谐波为基波，f s为采样频率，f 〇为基波频率，A f为离散频率间隔，且A f = f s/N。
[0156] 令：&= fV A f，&为真实频谱的谱线位置，由于非同步采样或非整周期截断，kQ 一般为小数，不为整数。
[0157] 若忽略负频率点处旁瓣的影响，公式（3)变为：
[0159] S2、确定四根谱线：
[0160] 在步骤S1得到的加窗FFT频谱峰值附近区域，匕处频率点左右各两条的谱线分 别为：第 V k2、k3、k4根，k ^ k2、k3、k4均为正整数，k ^ k2、k3、1^4的关系为：k 2= k 彳1，k3 = k2+l，k4= k 3+1，这四根谱线对应的幅值分别为yp y2、y3、y4。
[0161] 为方便计算，记变量 a = k-k2_0. 5,由于 0 < k_k2< 1，则-0? 5 < a < 〇? 5 ;
[0162] 记变量
[0163] S3、计算四谱线插值算法的修正公式：
[0164] 根据公式⑷和（5)得到：
[0166] 当N值比较大，例如：N> 1000时，公式（6)可以化简为一个函数0 = g(a)，其 反函数为a =g4(f3)。由于所采用的余弦窗系数均为实系数，其频率响应是偶对称的，因 而g( ?)和g<( ?)均为奇函数。
[0167] 采用多项式逼近方法，计算奇函数0 =g4(a)，表达式为：
[0168] a ~ pn X 0+p13X 0 3+…+plp 0 p (7)
[0169] pn，p13; "11￡)为多项式逼近的奇次项系数，p是奇数。
[0170] 根据公式（4)，求得电网信号第i次谐波的幅值Ai:
[0172] 其中，i为正整数，yi为加窗FFT后第i次谐波的幅值；
[0173] 以基波为例，考虑到y2、y3是离真实谱线点最近的两根谱线，给予较大权重，可以 得到基波幅值A 1:
[0175] 根据公式（7)的逼近方法，当N比较大，例如N > 1000时，窗函数系数为实系数， 公式（9)可表示为：
[0176] Aj = N (74+373+372+7!) u ( a ),
[0177] 其中，u(a)为修正公式，且为偶函数，逼近多项式不含奇次项。
[0178] 四谱线修正公式如下：
[0179] u ( a ) = (p2Q+p22 a 2+…+p2d a d) (10)
[0180] 公式（10)中42(^22*12[1为多项式逼近的偶次项系数，(1为拟合的最高阶次，且(1 为偶数。
[0181] S4、计算基波参数：
[0182] 考虑到y2、y3是离真实谱线点最近的两根谱线，给予较大权重，可以得到基波频率 fci、基波幅值A ::
[0183] f0= k ? A f = (k 2+ a +〇. 5) A f (11)
[0184] A! = N-1 (p2CI+p22 a 2+... +p2da d) (12)
[0185] 根据公式（4)，还可以得出基波的相位：
[0187] 仿照基波参数的求取，根据公式（6)、（7)、（11)、（12)、（13)，即可进行各次谐波参 数的分析。
[0188] 考虑到其中大量窗函数的离散傅里叶分析，其表达式为：
[0190] 由于N>>1，可以得到：
[0192] S5、加窗四谱线插值快速算法：
[0193] 根据公式（5)和（12)，在计算变量0和幅值&的时候，需要求出
、（y4+3y3+3y2+yi)。
[0194] 令变量：
[0195] 根据公式（4)可以得到：
[0196] 根据公式（13)可得到：arg(W(k)) = -k JT (18)
[0197]将公式（17)代入公式（16)可得到：
[0198] 分析得出：X〇〇与X(k3)同相位，X(k2)与X(k 4)同相位；且XG〇与X(k2)，X(k2) 与X(k3)，X(k 3)与X(k4)之间的相位之差均为31，相邻四条谱线的相位图参见图2所示，那 么：
[0201] 由图2可知：
[0202] 因此，可以再次令：
[0204]其中，C为T1的模，D为T2的模。
[0206] 其中：Re表示取实部，Im表示取虚部。
[0207] 由以上分析得出：计算各次谐波的参数时候，仅在计算幅值Ak的时候，需要进行一 次求模运算。根据公式（5)计算变量0时，可以利用T1和T2的实部或虚部进行快速计算。[0208] 同理，根据同一个主瓣相邻4根谱线的相位关系，求最大谱线时，直接利用插值前 FFT运算结果的实部和虚部，来寻求最大的向量。
[0209] S6、确定谐波参数：确定基波频率&后，在范围（kf ^5, kf^+5)内，重复步骤S2~ S5,直到所有谐波参数计算完毕。
[0210] S7、进行误差分析：在同样窗函数条件下，分析加窗四谱线插值快速算法的误差， 并与双谱线算法的误差、三谱线插值算法的误差进行比较。
[0211] Sl、S2、S3、S4、S6、S7形成一个完整的基于加窗四谱线插值FFT的高精度谐波分 析方法，能够有效提高谐波分析的精度。
[0212] 步骤S5的作用是提高计算速度，在S4、S6之间增加步骤S5,每次谐波仅需一次模 运算，计算量小，能够有效提升计算速度，在同样窗函数下，能获得更高的精度，具有较高的 实用价值。
[0213] 下面通过一个具体实施例来进行详细说明。选定一个电网信号的谐波，谐波的波 形参见图3所示。
[0214] 步骤1、信号的预处理：
[0215] 将互感器采集到的离散电网信号x(n)，信号模型参见表1所示。为了验证所提算 法的精度，进行10次谐波仿真分析。信号模型为：
[0217] 其 中：基波频率fQ为50. 5Hz，采样频率匕为5120Hz，采样点数N为1024。
[0218] 表1、电网信号的谐波参数
[0220]
[0221] 对信号x(n)进行加窗处理，选取4种常用加窗函数，其系数参见表2所示。
[0222] 表2、常用窗函数系数
[0224] 加窗FFT变换后得到信号频谱：
[0225] 步骤2、确定四根谱线：
[0226] 在45~55Hz中寻求实部（或虚部）最大四根谱线，分别为XGO、X(k 2)、X(k3)、 X (k4)〇令：
[0228] 步骤3、计算四谱线插值算法的修正公式：
[0229] 根据【具体实施方式】中步骤S3的推导，表2中的4种窗函数的修正公式a = g4(f3)，u(a)分别如下：
[0230] (1) Hanning 窗：
[0231] a = 1. 13013682 0-0. 18408680 03+〇, 07224859 05-〇. 04451832 07u( a)= 0. 53549869+0. 17622103 a 2+〇. 09310826 a 4+〇. 06644946 a 6
[0232]Hanning窗函数的幅频特性参见图4所示。
[0233] (2) Balckman 窗：
[0234] a = 1.44649012 0-0. 29326578 0 3+0. 13858447 05-〇. O9578617 07u( a)= 1. 1575129+0. 56110888 a 2+〇. 38707997 a 4+〇. 33290703 a 6
[0235]Black窗函数的幅频特性参见图5所示。
[0236] (3)Balckman-harris 窗：
[0237] a = 1. 85862073 0-0. 40299738 03+O. 19725715 05-〇. 13293118 07u( a)= 1. 24914356+0. 88776025 a 2+〇. 73879907 a 4-〇. 69429366 a 6
[0238] Black-harris窗函数的幅频特性参见图6所示。
[0239] (4) 4项最大旁瓣衰减窗：
[0240]a = 2. 37540983 0-0. 43478665 03+〇, 19124517 05-〇. 11511086 07u( a)= 1. 34456750+1. 36577423 a 2+l. 25666570 a 4+l. 17989176 a 6
[0241] 4项最大旁瓣衰减窗函数的幅频特性参见图7所示。
[0242] 步骤4、计算基波参数：
[0243] Aj=N (74+373+372+7!) u( a )
[0245] f〇= k ? A f = (k 2+ a +〇. 5) A f
[0246] 步骤5、四谱线插值快速算法：
[0247] 根据S5所述快速算法，可以对步骤2和步骤4中的两个量的求取采用快速算法：
[0248]
及Af N-1|T2|u(a)
[0249] 步骤6、确定谐波参数：确定基波频率f0后，在范围（kf^-5, kf^+5)内重复步骤2~ 5,直到所有谐波参数计算完毕。
[0250] 步骤7、进行误差分析：仿真实验的误差分析比较结果参见表3~表6所示，其中， DAi表示基波和各次谐波幅值测量值的相对误差；D f(l表示基波频率测量值的相对误差； 表示基波和各次谐波初始相位测量值的相对误差，均用百分比表示。
[0251] 表3、Hanning窗和Blackman窗的频率、幅值相对误差
[0253] 表4、Hanning窗和Blackman窗的相位相对误差
[0254]
[0255] 表5、Blackman-harris窗和四项最大旁瓣衰减窗的频率、幅值相对误差
[0257] 表6、Blackman-harris窗和四项最大旁瓣衰减窗的相位相对误差
[0259] 从表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)的表达式为：
[0265] 其中，N为采样点数，N为正整数，n = 0, 1，2. . . N-1 ; 2表示求和；m为窗函数的累 加次数，m = 0, 1,2. ..M-l ;M为窗函数项数，M为正整数；\为窗函数系数，窗函数系数满 足以下约束条件：
[0267] 对公式⑴的加窗信号进行FFT变换后，得到加窗FFT频谱：
[0269] 其中，W(_)为窗函数的频谱，k为正整数，X(k)表示第k次谐波的频谱，Ak为第k 次谐波的幅值，j表示虚数单位，e是自然对数的底数，为第k次谐波的初始相位，第一次 谐波为基波，f s为采样频率，f 〇为基波频率，A f为离散频率间隔，且A f = f s/N。
[0270] 令：ke fV A f，&为真实频谱的谱线位置，由于非同步采样或非整周期截断，h 一般不为整数。
[0271] 若忽略负频率点处旁瓣的影响，公式（3)变为：
[0273] 谱线确定单元，用于确定四根谱线：
[0274] 在信号预处理单元得到的加窗FFT频谱峰值附近区域，匕处频率点左右各两条的 谱线分别为：第1^ 1、1^2、1^、1^4根，1^1、1^2、1^、1^ 4均为正整数，1^1、1^2、1^、1^4的关系为4 2=1^1+1， k3= k2+l，k4= k3+l，这四根谱线对应的幅值分别为yi、y2、y3、y4。
[0275] 为方便计算，记变量a = k-k2-0. 5,由于0彡k_k2< 1，则-0? 5彡a彡〇? 5;
[0276] 记变量
[0277] 修正公式计算单元，用于计算四谱线插值算法的修正公式：
[0278] 根据公式⑷和（5)得到：
[0280]当N值比较大，例如：N> 1000时，公式（6)可以化简为一个函数0 = g( a)，其 反函数为a =g4(f3)。由于所采用的余弦窗系数均为实系数，其频率响应是偶对称的，因 而g( ?)和g<( ?)均为奇函数。
[0281] 采用多项式逼近方法，计算奇函数0 =g4(a)，表达式为：
[0282] a ~ pnX 0+p13X 03+…+plp0p (7)
[0283] pn，p13; "11￡)为多项式逼近的奇次项系数，p是奇数。
[0284]根据公式（4)，求得电网信号第i次谐波的幅值Ai:
[0286] 其中，i为正整数，yi为加窗FFT后第i次谐波的幅值；
[0287] 以基波为例，考虑到y2、y3是离真实谱线点最近的两根谱线，给予较大权重，可以 得到基波幅值A 1:
[0289] 根据公式（7)的逼近方法，当N比较大，例如N > 1000时，窗函数系数为实系数， 公式（9)可表示为：
[0290] Aj = N (74+373+372+7!) u ( a ),
[0291] 其中，u(a)为修正公式，且为偶函数，逼近多项式不含奇次项。
[0292] 四谱线修正公式如下：
[0293] u ( a ) = (p2〇+p22 a 2+…+p2d a d) (10)
[0294] 公式（10)中，p2(l，p22- p2$多项式逼近的偶次项系数，d为拟合的最高阶次，且d 为偶数。
[0295] 基波参数计算单元，用于计算基波参数：
[0296] 考虑到y2、y3是离真实谱线点最近的两根谱线，给予较大权重，可以得到基波频率 fci、基波幅值A ::
[0297] f0= k ? A f = (k 2+ a +〇. 5) A f (11)
[0298] Ai = N-1 (yWys+Sy^yD (p2CI+p22 a 2+... +p2da d) (12)
[0299] 根据公式（4)还可以得出基波的相位：
[0300]
[0301] 仿照基波参数的求取，根据公式（6)、（7)、（11)、（12)、（13)，即可进行各次谐波参 数的分析。
[0302] 考虑到其中大量窗函数的离散傅里叶分析，其表达式为：
[0304]由于N>>1，可以得到：
[0306] 四谱线插值快速计算单元，用于进行四谱线插值快速计算：
[0307] 由公式（5)和（12)可知，在计算变量0和幅值&的时候，需要求出
[0308] 令变量：
[0309] 根据公式（4)可以得到
[0310] 根据公式（13)可得到：arg(W(k)) = -k JT (18)
[0311] 将公式（17)代入公式（16)可得到
[0312] 分析得出：X〇〇与X(k3)同相位，X(k2)与X(k 4)同相位；且XG〇与X(k2)，X(k2) 与X(k3)，X(k 3)与X(k4)之间的相位之差均为31，相邻四条谱线的相位图参见图2所示，那 么：
[0314]
[0315] 由图2可知：
[0316] 因此，可以再次令：
[0318] 其中，C为T1的模，D为T2的模。
[0319]贝1J :
[0320] 其中：Re表示取实部，Im表示取虚部。
[0321] 由以上分析可知，计算各次谐波的参数时候，仅在计算幅值Ak的时候，需要进行一 次求模运算。根据公式（5)计算变量0时，可以利用T1和T2的实部或虚部进行快速计算。 同理，根据同一个主瓣相邻4根谱线的相位关系，求最大谱线时，直接利用插值前FFT运算 结果的实部和虚部，来寻求最大的向量。
[0322] 谐波参数确定单元，用于确定谐波参数：确定基波频率&后，在范围 (kfd-5, kfd+5)内，谱线确定单元、修正公式计算单元、基波参数计算单元、四谱线插值快速 计算单元重复进行计算，直到所有谐波参数计算完毕。
[0323] 误差分析单元，用于进行误差分析：在同样窗函数条件下，分析加窗四谱线插值快 速算法的误差，并与双谱线算法的误差、三谱线插值算法的误差进行比较。
[0324] 本领域的技术人员可以对本发明实施例进行各种修改 和变型，倘若这些修改和变 型在本发明权利要求及其等同技术的范围之内，则这些修改和变型也在本发明的保护范围 之内。
[0325] 说明书中未详细描述的内容为本领域技术人员公知的现有技术。
【主权项】
1. 一种基于加窗四谱线插值FFT的谐波快速分析方法，其特征在于，包括以下步骤： 51、 信号预处理： 互感器采集电网信号，将互感器采集到的离散电网信号x(n)，传输到上位机，η为采样 点的序数，η为自然数；采用离散余弦窗函数w (η)，对电网信号X (η)进行加窗截断，得到加 窗信号xw (η): xw (η) = X (n) w (η) (I) 离散余弦窗函数w (η)的表达式为：(2) 其中，N为采样点数，N为正整数，η = 0, 1,2. ..N-I ;Σ表示求和；m为窗函数的累加次 数，m = 0, 1，2. .. M-I ;M为窗函数项数，M为正整数；Iddi为窗函数系数； 对公式（1)的加窗信号进行FFT变换后，得到加窗FFT频谱：其中，W( ·)为窗函数的频谱，k为正整数，X(k)表示第k次谐波的频谱，Ak为第k次 谐波的幅值，j表示虚数单位，e是自然对数的底数，为第k次谐波的初始相位，第一次谐 波为基波，fs为采样频率，f 〇为基波频率，Δ f为离散频率间隔，且Λ f = f s/N ; 令：1?= fd/Δ f，1?为真实频谱的谱线位置； 忽略负频率点处旁瓣的影响，公式（3)变为：(4) 52、 确定四根谱线： 在步骤Sl得到的加窗FFT频谱峰值附近区域，1?处频率点左右各两条的谱线分别为： 第1^、1^2、1^、1^4根，1^1、1^2、1^、1^ 4均为正整数，1^1、1^2、1^、1^4的关系为 ：1^2=1^1+1，1^=1^2+1，1^4=k3+l，这四根谱线对应的幅值分别为y2、y3、y 4; 记变量 α = k-k2-0. 5,由于 0 彡 k-k2< 1，则-0· 5 彡 α 彡 〇· 5 ; 记变量（5 ) 53、 计算四谱线插值算法的修正公式： 根据公式（4)和（5)得到：采用多项式逼近方法，计算奇函数β =Ρ(α)，表达式为：(7) Pn，P13;…P 1P为多项式逼近的奇次项系数，P是奇数； 根据公式（4)，求得电网信号第i次谐波的幅值Ai:(8) 其中，i为正整数，Yi为加窗FFT后第i次谐波的幅值； 考虑到y2、y3是离真实谱线点最近的两根谱线，基波幅值A :根据公式（7)的逼近方法，N > 1000时，窗函数系数为实系数，公式（9)表示为： A1 = N (74+373+372+7!) u ( α ), 其中，u(a)为修正公式，且为偶函数，逼近多项式不含奇次项； 四谱线修正公式为：u( a ) = (ρ2(ι+ρ22 α 2+…+p2da d) (10) 公式（10)中，P2(I，P22··· P2A多项式逼近的偶次项系数，d为拟合的最高阶次，且d为偶 数； 54、 计算基波参数： 计算基波频率fci、基波幅值A1: f0= k · Af= (k 2+ a +〇. 5) Δ f (11) A1= N (74+373+372+7!) (p2〇+p22 a 2+··· +p2d a d) (12) 根据公式（4)，得出基波的相位：仿照基波参数的求取，根据公式（6)、（7)、（11)、（12)、（13)，进行各次谐波参数的分 析； 考虑到其中大量窗函数的离散傅里叶分析，其表达式为：由于N>>1，得到：55、 加窗四谱线插值快速算法： 根据公式（5)和（12)，计算变量β和幅值A1的时候，需求出根据公式（4)，得到 根据公式（13)，得到 将公式（17)代入公式（16)，得到分析得出：XG〇与X(k3)同相位，X(k2)与X(k4)同相位；且X(Ic 1)与X(k2)，X(k2)与 X(k3)，X(k3)与X(k4)之间的相位之差均为π，那么：其中，C为Tl的模，D为Τ2的模； 则：其中：Re表示取实部，Im表示取虚部； 分析得出：计算各次谐波的参数时，仅在计算幅值Ak的时候，需要进行一次求模运算； 根据公式（5)计算变量β时，利用Tl和T2的实部或虚部进行快速计算；同理，根据同一个 主瓣相邻4根谱线的相位关系，求最大谱线时，直接利用插值前FFT运算结果的实部和虚部 来寻求最大的向量； S6、确定谐波参数：确定基波频率&后，在范围（kf ^5, kf^+5)内重复步骤S2~S5,直 到所有谐波参数计算完毕； S7、进行误差分析：在同样窗函数条件下，分析加窗四谱线插值快速算法的误差，并与 双谱线算法的误差、三谱线插值算法的误差进行比较。2. 如权利要求1所述的基于加窗四谱线插值FFT的谐波快速分析方法，其特征在于： 所述电网信号包括电流信号、电压信号。3. 如权利要求1所述的基于加窗四谱线插值FFT的谐波快速分析方法，其特征在于： 所述真实频谱的谱线位置1?为小数。4. 如权利要求1或2或3所述的基于加窗四谱线插值FFT的谐波快速分析方法，其特 征在于： 所述窗函数系数\满足以下约束条件：5. -种基于加窗四谱线插值FFT的谐波快速分析系统，其特征在于：包括信号预处理 单元、谱线确定单元、修正公式计算单元、基波参数计算单元、四谱线插值快速计算单元、谐 波参数确定单元、误差分析单元，其中： 所述信号预处理单元，用于对信号进行预处理： 互感器采集电网信号，将互感器采集到的离散电网信号x(n)，传输到上位机，η为采样 点的序数，η为自然数；采用离散余弦窗函数w (η)，对电网信号X (η)进行加窗截断，得到加 窗信号xw (η): xw (η) = X (n) w (η) (I) 离散余弦窗函数w (η)的表达式为：(2) 其中，N为采样点数，N为正整数，η = 0, 1,2. ..N-I ;Σ表示求和；m为窗函数的累加次 数，m = 0, 1，2. .. M-I ;M为窗函数项数，M为正整数；Iddi为窗函数系数； 对公式（1)的加窗信号进行FFT变换后，得到加窗FFT频谱：其中，W( ·)为窗函数的频谱，k为正整数，X(k)表示第k次谐波的频谱，Ak为第k次 谐波的幅值，j表示虚数单位，e是自然对数的底数，PA-为第k次谐波的初始相位，第一次谐 波为基波，fs为采样频率，f 〇为基波频率，Δ f为离散频率间隔，且Λ f = f s/N ; 令：1?= fd/Δ f，1?为真实频谱的谱线位置； 忽略负频率点处旁瓣的影响，公式（3)变为：(4) 所述谱线确定单元，用于确定四根谱线： 在信号预处理单元得到的加窗FFT频谱峰值附近区域，1?处频率点左右各两条的谱线 分别为：第 1^1、1^2、1^、1^4根，1^1、1^ 2、1^、1^4均为正整数，1^1、1^2、1^、1^ 4的关系为：1^2=1^1+1，1^ = k2+l，k4= k 3+1，这四根谱线对应的幅值分别为yp y2、y3、y4; 记变量 α = k-k2-0. 5,由于 O 彡 k-k2< 1，则-0· 5 彡 α 彡 〇· 5 ; 记变量（5) 所述修正公式计算单元，用于计算四谱线插值算法的修正公式： 根据公式（4)和（5)得到：采用多项式逼近方法，计算奇函数β =Ρ(α)，表达式为：(7) Pn，P13;…P 1Ρ为多项式逼近的奇次项系数，P是奇数； 根据公式（4)，求得电网信号第i次谐波的幅值Ai:(8) 其中，i为正整数，Yi为加窗FFT后第i次谐波的幅值； 考虑到y2、y3是离真实谱线点最近的两根谱线，基波幅值A :根据公式（7)的逼近方法，N > 1000时，窗函数系数为实系数，公式（9)表示为： A1 = N (74+373+372+7!) u ( α ), 其中，u(a)为修正公式，且为偶函数，逼近多项式不含奇次项； 四谱线修正公式为：u( a ) = (ρ2(ι+ρ22 α 2+…+p2da d) (10) 公式（10)中，P2(I，P22··· P2A多项式逼近的偶次项系数，d为拟合的最高阶次，且d为偶 数； 所述基波参数计算单元，用于计算基波参数： 计算基波频率fci、基波幅值A1: f0= k · Af= (k 2+ a +〇. 5) Δ f (11) A1= N (74+373+372+7!) (p2〇+p22 a 2+··· +p2d a d) (12) 根据公式（4)，得出基波的相位：仿照基波参数的求取，根据公式（6)、（7)、（11)、（12)、（13)，进行各次谐波参数的分 析； 考虑到其中大量窗函数的离散傅里叶分析，其表达式为：由于N>>1，得到：所述四谱线插值快速计算单元，用于进行四谱线插值快速计算： 根据公式（5)和（12)，在计算变量β和幅值A1的时候，需要求出、（y4+3y3 +3y2+yi): 令变量：根据公式（4)，得到 根据公式（13)，得到 将公式（17)代入公式（16)，得到：分析得出：XG〇与X(k3)同相位，X(k2)与X(k4)同相位；且X(Ic 1)与X(k2)，X(k2)与 X(k3)，X(k3)与X(k4)之间的相位之差均为π，那么：令：（23) 其中，C为Tl的模，D为T2的模； 则：(24) 其中：Re表示取实部，Im表示取虚部； 分析得出：计算各次谐波的参数时，仅在计算幅值Ak的时候，需要进行一次求模运算； 根据公式（5)计算变量β时，利用Tl和T2的实部或虚部进行快速计算；同理，根据同一 个主瓣相邻4根谱线的相位关系，求最大谱线时，直接利用插值前FFT运算结果的实部和虚 部，来寻求最大的向量； 所述谐波参数确定单元，用于确定谐波参数：确定基波频率后，在范围 (kfd-5, kfd+5)内，所述谱线确定单元、修正公式计算单元、基波参数计算单元、四谱线插值 快速计算单元重复进行计算，直到所有谐波参数计算完毕； 所述误差分析单元，用于进行误差分析：在同样窗函数条件下，分析加窗四谱线插值快 速算法的误差，并与双谱线算法的误差、三谱线插值算法的误差进行比较。6. 如权利要求5所述的基于加窗四谱线插值FFT的谐波快速分析系统，其特征在于： 所述电网信号包括电流信号、电压信号。7. 如权利要求5所述的基于加窗四谱线插值FFT的谐波快速分析系统，其特征在于： 所述真实频谱的谱线位置1?为小数。8. 如权利要求5或6或7所述的基于加窗四谱线插值FFT的谐波快速分析系统，其特 征在于： 所述窗函数系数\满足以下约束条件：
【专利摘要】本发明公开了一种基于加窗四谱线插值FFT的谐波快速分析方法及系统，涉及谐波分析领域。该方法包括以下步骤：信号预处理；确定四根谱线；计算四谱线插值算法的修正公式；计算基波参数；四谱线插值快速算法；确定谐波参数；进行误差分析。本发明从电力系统的电网信号加窗后的频域表达式入手，根据窗函数主瓣内任意相邻谱线相位相差π的规律，推导出加窗FFT后的真实谱线点附近最大的四根谱线之间的相位规律，提出一种加窗四谱线插值FFT快速算法。相对于双谱线和三谱线插值算法，本发明能有效提高谐波分析的精度；利用该快速算法，计算某次谐波仅需要1次模的运算，能够有效降低计算量和计算时间，显著提升计算速度。
【IPC分类】G01R23/16
【公开号】CN104897960
【申请号】CN201510326063
【发明人】张俊敏, 刘开培, 汪立, 王黎, 田微, 何顺帆, 姚为
【申请人】中南民族大学
【公开日】2015年9月9日
【申请日】2015年6月15日
最新回复(0)