基于互乘法窗函数的三谱线插值fft谐波分析方法及系统的制作方法
基于互乘法窗函数的三谱线插值fft谐波分析方法及系统的制作方法
【技术领域】
[0001] 本发明涉及谐波分析领域，具体是涉及一种基于互乘法窗函数的三谱线插值FFT 谐波分析方法及系统。
【背景技术】
[0002] 随着大量非线性电力电子器件的使用，使得电网中谐波污染问题日益严重。谐波 问题不仅恶化电能质量，对电网的安全稳定和经济运行也造成较大影响。因此，对系统中谐 波参数进行高精度测量，对于减少谐波危害，维护电网安全稳定、高效运行是十分必要的。
[0003] 谐波检测的关键问题是：非同步采样下，如何解决频谱泄露和栅栏效应。目前，常 见的检测谐波和间谐波的方法主要有：傅里叶变换，机器学习法，小波变换，功率谱估计，希 尔伯特-黄变换。
[0004] 在傅里叶变换分析谐波中，所用方法是：运用各种特殊窗函数，对信号进行截断， 然后结合谱线插值FFT (Fast Fourier Transform，快速傅里叶变换）进行谐波分析。
[0005] 常用的窗函数有：汉宁（Hanning)窗函数、布莱克曼（Blackman)窗函数、布莱克曼 汉斯（Blackman-Harris)窗函数、纳托尔（Nuttall)窗函数、莱夫文森特（Rife-Vincent) 窗函数以及各种组合窗函数。
[0006] 在加窗基础上，D. Agrez和庞浩等人各自提出了双谱线的修正算法，Wu Jing、牛胜 锁和黄冬梅等人提出了三谱线修正算法。这些改进降低了频谱泄漏和栅栏效应的影响，提 高了谐波分析的准确性。
[0007] 然而，在工程实际使用中，常用的窗函数插值算法仍然无法满足高精度的谐波分 析要求。
【发明内容】
[0008] 本发明的目的是为了克服上述【背景技术】的不足，提供一种基于互乘法窗函数的三 谱线插值FFT谐波分析方法及系统，与常规窗函数比，具有更高的精度。
[0009] 本发明提供一种基于互乘法窗函数的三谱线插值FFT谐波分析方法，包括以下步 骤：
[0010] S1、构造互乘法窗函数：
[0011] 乘法窗函数的通用公式是由多个窗函数乘积产生的，乘法窗函数W(n)的通用公 式为：
[0013] 其中，n为采样点的序数，n为自然数；m为参加构造乘法窗的基本窗函数的个数， m为正整数为第i个基本窗函数表达式，i = 1，2…m七为参加第i个基本窗函数 的个数，称为Wi (n)的子阶次
，C为互乘法窗函数的总阶次；根据Cl、c 2、(33的 不同取值，形成不同的组合（c^cv-cj，即构造出不同形式的乘法窗函数；
[0014] 当^〇1)表达式相同时，即所采用的基本窗相同时，w(n)为自乘法窗函数，当wjn) 表达式不同时，即所采用的基本窗不同时，w(n)为互乘法窗函数；
[0015] 三个常用窗函数的表达式分别为％ (n)，w2 (n)，w3(n)，总阶次c = Wq，给c 一 个最大值，令：c = 3 ;
[0016] 根据Ci、c2、c3的不同取值，形成不同的组合（cffcj，构造出以下9种乘法窗函 数：
[0017] 当Ci= 3、c 2= 0、c 3= 0时，构造出互乘300窗函数；
[0018] 当(^=2、(32= l、c3=0时，构造出互乘210窗函数；
[0019] 当(^=2、(32=0、(33= 1时，构造出互乘201窗函数；
[0020] 当Ci= 1、c 2= 2、c 3= 0时，构造出互乘120窗函数；
[0021] 当~二l、c2=0、c3=2时，构造出互乘102窗函数；
[0022] 当Ci= 0、c 2= 3、c 3= 0时，构造出互乘030窗函数；
[0023] 当(^= l、c2= l、c3= 1时，构造出互乘111窗函数；
[0024] 当Ci= 0、c 2= 2、c 3= 1时，构造出互乘021窗函数；
[0025] 当Ci= 0、c 2= 0、c 3= 3时，构造出互乘003窗函数；
[0026] 其中，互乘300窗函数、互乘030窗函数、互乘003窗函数，这三种窗函数属于自乘 法窗函数，其余6种为构造出的互乘法窗函数；
[0027] S2、信号预处理：
[0028] 互感器采集电网信号，将互感器采集到的电网信号x(n)，传输到上位机；对电网 信号x(n)进行加互乘法窗函数w(n)进行截断，得到加窗信号\〇1):
[0029] xw(n) = x(n)w(n) (2)
[0030] 对公式⑵的加窗信号进行FFT变换后，得到加窗FFT频谱：
[0032] 其中，W( ?)为窗函数的频谱，k为正整数，X(k)表示第k次谐波的频谱，Ak为第 k次谐波的幅值，j表示虚数单位，e是自然对数的底数，外;为第k次谐波的初始相位，第一 次谐波为基波，fs为采样频率，f〇为基波频率，A f为离散频率间隔，且A f = f s/N ;
[0033] 令：k0= f。/Af，k。为真实频谱的谱线位置，
[0034]忽略负频率点处旁瓣的影响，公式（3)变为：
[0036] S3、确定三根谱线：
[0037] 在S2得到的加窗FFT频谱峰值附近区域，&处频率点较大的三根谱线分别为：第 V k2、k3根，k i、k2、k3均为正整数，k k 3的关系为：k2 = k i+1，k3= k 2+1，这三根谱线对 应的幅值分别为yp y2、
[0038] 记变量 a = k_k2,则-0? 5 彡 a 彡 〇? 5 ;
[0039] 另记变量
[0040] S4、计算三谱线插值算法的修正公式：
[0041] 根据公式⑷和（5)得到：
[0043] 采用多项式逼近方法计算奇函数0 =g4(a)，表达式为：
[0044] a ~ pnX 0+p13X 03+…+plp0p (7)
[0045] pn，p13; -p 1￡)为多项式逼近的奇次项系数，p是奇数；
[0046] 根据公式⑷，求得电网信号第i次谐波的幅值Ai:
[0048] 其中，i为正整数，yi为加窗FFT后第i次谐波的幅值；
[0049] 考虑到y2是离真实谱线点最近谱线，得到：
[0051] N > 1000时，窗函数系数为实系数，公式（9)表示为：
[0052] Aj= N_1 (y3+2y2+y1)u( a )
[0053] 其中，u(a)为修正公式，且为偶函数，逼近多项式不含奇次项；
[0054] 三谱线修正逼近多项式如下：
[0055] u ( a ) = (p2〇+p22 a 2+…+p2d a d) (10)
[0056] 公式（10)中，p2(l，py p2$多项式逼近的偶次项系数，d为拟合的最高阶次，且d 为偶数；
[0057] S5、计算基波参数：
[0058] 计算基波频率fQ、基波幅值Ay
[0059] f0= k? Af = (k2+a)Af (11)
[0060] A! = N-1 (p2CI+p22 a 2+... +p2d a d) (12)
[0061]根据公式（4)，计算基波的相位：
[0063] 仿照基波参数的求取，根据公式（6)、（7)、（9)、（11)、（12)、（13)进行各次谐波参 数的分析；
[0064] S6、确定谐波参数：确定基波频率fQ后，在范围（kf。_5, kfQ+5)内，重复步骤S3~ S5,直到所有谐波参数计算完毕；
[0065] S7、进行误差分析：分析基于互乘法窗函数的三谱线插值FFT算法的误差，并与常 规窗函数插值算法精度进行比较。
[0066] 在上述技术方案的基础上，所述电网信号包括电流信号、电压信号。
[0067] 在上述技术方案的基础上，所述真实频谱的谱线位置&为小数。
[0068] 本发明还提供一种基于互乘法窗函数的三谱线插值FFT谐波分析系统，包括互乘 法窗函数构造单元、信号预处理单元、谱线确定单元、修正公式计算单元、基波参数计算单 元、谐波参数确定单元、误差分析单元，其中：
[0069] 所述互乘法函数构造单元，用于构造互乘法窗函数：
[0070] 乘法窗函数的通用公式是由多个窗函数乘积产生的，乘法窗函数w(n)的通用公 式为：
[0072] 其中，n为采样点的序数，n为自然数；m为参加构造乘法窗的基本窗函数的个数， m为正整数为第i个基本窗函数表达式，i = 1，2…m七为参加第i个基本窗函数 的个数，称为K (n)的子阶次；
，c为互乘法窗函数的总阶次；根据Cl、c 2、(：3的 不同取值，形成不同的组合（c^cv-cj，即构造出不同形式的乘法窗函数；
[0073] 当^〇1)表达式相同时，即所采用的基本窗相同时，w(n)为自乘法窗函数，当wjn) 表达式不同时，即所采用的基本窗不同时，w(n)为互乘法窗函数；
[0074] 三个常用窗函数的表达式分别为％ (n)，w2 (n)，w3(n)，总阶次c = Ci+c2+c3，给c 一 个最大值，令：c = 3 ;
[0075] 根据Cl、c2、c3的不同取值，形成不同的组合（ Clcvcm)，构造出以下9种乘法窗函 数 ：
[0076] 当Ci= 3、c 2= 0、c 3= 0时，构造出互乘300窗函数；
[0077] 当Ci= 2、c 2= 1、c 3= 0时，构造出互乘210窗函数；
[0078] 当Ci= 2、c 2= 0、c 3= 1时，构造出互乘201窗函数；
[0079] 当Ci= 1、c 2= 2、c 3= 0时，构造出互乘120窗函数；
[0080] 当Ci= 1、c 2= 0、c 3= 2时，构造出互乘102窗函数；
[0081] 当Ci= 0、c 2= 3、c 3= 0时，构造出互乘030窗函数；
[0082] 当(^= l、c2= l、c3= 1时，构造出互乘111窗函数；
[0083] 当Ci= 0、c 2= 2、c 3= 1时，构造出互乘021窗函数；
[0084] 当Ci= 0、c 2= 0、c 3= 3时，构造出互乘003窗函数；
[0085] 其中，互乘300窗函数、互乘030窗函数、互乘003窗函数，这三种窗函数属于自乘 法窗函数，其余6种为构造出的互乘法窗函数；
[0086] 所述信号预处理单元，用于对信号进行预处理：
[0087] 互感器采集电网信号，将互感器采集到的电网信号x(n)，传输到上位机；对电网 信号x(n)进行加互乘法窗函数w(n)进行截断，得到加窗信号\〇1):
[0088] xw(n) = x(n)w(n) (2)
[0089] 对公式⑵的加窗信号进行FFT变换后，得到加窗FFT频谱：
[0091] 其中，W( ?)为窗函数的频谱，k为正整数，X(k)表示第k次谐波的频谱，Ak为第 k次谐波的幅值，j表示虚数单位，e是自然对数的底数，？^为第k次谐波的初始相位，第一 次谐波为基波，f s为采样频率，f〇为基波频率，A f为离散频率间隔，且A f = f s/N ;
[0092] 令：k0= f。/ A f，k。为真实频谱的谱线位置，
[0093] 忽略负频率点处旁瓣的影响，公式（3)变为：
[0095] 所述谱线确定单元，用于确定三根谱线：
[0096] 在S2得到的加窗FFT频谱峰值附近区域，h处频率点较大的三根谱线分别为：第 根，k 为正整数，k广k 3的关系为：k 2= k汴1，k3= k 2+1，这三根谱线对应 的幅值分别为y2、y3;
[0097] 记变量 a = k_k2，则 _0? 5 < a < 〇? 5 ;
[0098] 另记变量
[0099] 所述修正公式计算单元，用于计算三谱线插值算法的修正公式：
[0100] 根据公式⑷和（5)得到：
[0102] 采用多项式逼近方法计算奇函数0 =g4(a)，表达式为：
[0103] a ~ pnX 0+p13X 03+…+plp0p (7)
[0104] pn，p13;…p lp为多项式逼近的奇次项系数，p是奇数；
[0105] 根据公式（4)，求得电网信号第i次谐波的幅值Ai:
[0107] 其中，i为正整数，yi为加窗FFT后第i次谐波的幅值；
[0108] 考虑到y2是离真实谱线点最近谱线，得到：
[0110] N > 1000时，窗函数系数为实系数，公式（9)表示为：
[0111] Aj= N_1 (y3+2y2+y1)u( a )
[0112] 其中，u(a)为修正公式，且为偶函数，逼近多项式不含奇次项；
[0113] 三谱线修正逼近多项式如下：
[0114] u ( a ) = (p2Q+p22 a 2+…+p2d a d) (10)
[0115] 公式（10)中，p2Q，p22…p 2d为多项式逼近的偶次项系数，d为拟合的最高阶次，且d 为偶数；
[0116] 所述基波参数计算单元，用于计算基波参数：
[0117] 计算基波频率fQ、基波幅值Ay
[0118] f0= k? Af = (k2+a)Af (11)
[0119] Ai= N 1 (y3+2y2+y) (p20+p22 a 2+…+p2d a d) (12)
[0120]根据公式（4)，计算基波的相位：
[0122] 仿照基波参数的求取，根据公式（6)、（7)、（9)、（11)、（12)、（13)进行各次谐波参 数的分析；
[0123] 所述谐波参数确定单元，用于确定谐波参数：确定基波频率4后，在范围 (kfd-5, kf^+5)内，所述谱线确定单元、修正公式计算单元、基波参数计算单元重复进行计 算，直到所有谐波参数计算完毕；
[0124] 所述误差分析单元，用于进行误差分析：分析基于互乘法窗函数的三谱线插值 FFT算法的误差，并与常规窗函数插值算法精度进行比较。
[0125] 在上述技术方案的基础上，所述电网信号包括电流信号、电压信号。
[0126] 在上述技术方案的基础上，所述真实频谱的谱线位置k0为小数。
[0127] 与现有技术相比，本发明的优点如下：
[0128] 本发明提出一种新的窗函数的构造方式，利用常规窗函数进行不同的乘法组合， 首次实现构造互乘法窗函数，将互乘法窗函数应用到三谱线插值FFT算法中，进行谐波分 析，计算谐波参数。多种常用的余弦窗函数计算实例表明，本发明构造出的互乘法窗函数相 对于常规窗函数插值算法，有更高的准确度，实现了谐波的高精度测量，具有较高的实用价 值。
【附图说明】
[0129] 图1是本发明实施例中基于互乘法窗函数的三谱线插值FFT谐波分析方法的流程 图。
[0130] 图2是谐波信号的波形图。
[0131] 图3是互乘210窗函数的幅频特性图。
[0132] 图4是互乘201窗函数的幅频特性图。
[0133] 图5是互乘120窗函数的幅频特性图。
[0134] 图6是互乘102窗函数的幅频特性图。
[0135] 图7是互乘111窗函数的幅频特性图。
[0136] 图8是互乘021窗函数的幅频特性图。
【具体实施方式】
[0137] 下面结合附图及具体实施例对本发明作进一步的详细描述。
[0138] 参见图1所示，本发明实施例提供一种基于互乘法窗函数的三谱线插值FFT谐波 分析方法，包括以下步骤：
[0139] S1、构造互乘法窗函数：
[0140] 乘法窗函数的通用公式是由多个窗函数乘积产生的，乘法窗函数w(n)的通用公 式为：
[0142] 其中，n为采样点的序数，n为自然数；m为参加构造乘法窗的基本窗函数的个数，m 为正整数;Wi (n)为第i个基本窗函数表达式，i = 1，2…m力为参加第i个基本窗函数的 个数，称为wjn)的子阶次：
，c为互乘法窗函数的总阶次。根据Cl、c 2、(：3的 不同取值，形成不同的组合（c^cv-cj，即构造出不同形式的乘法窗函数。
[0143] 当^〇1)表达式相同时，即所采用的基本窗相同时，w(n)为自乘法窗函数，当wjn) 表达式不同时，即所采用的基本窗不同时，w(n)为互乘法窗函数。
[0144] 下面以三个常用窗函数为例，表达式分别*Wl(n)，w2(n)，w 3(n)。总阶次c = Cl+C2+C3。由于C越大，意味着参与乘法窗的基本窗函数越多，构成的乘法窗函数项数越多， 会增加计算的复杂度。一般给C 一个最大值，例如，令：c = 3。
[0145] 根据CpCyq的不同取值，形成不同的组合（c而…cm)，构造出以下9种乘法窗函 数：
[0146] 当Ci= 3、c 2= 0、c 3= 0时，构造出互乘300窗函数；
[0147] 当Ci = 2、c 2= 1、c 3= 0时，构造出互乘210窗函数；
[0148] 当Ci= 2、c 2= 0、c 3= 1时，构造出互乘201窗函数；
[0149] 当Ci= 1、c 2= 2、c 3= 0时，构造出互乘120窗函数；
[0150] 当Ci= 1、c 2= 0、c 3= 2时，构造出互乘102窗函数；
[0151] 当Ci= 0、c 2= 3、c 3= 0时，构造出互乘030窗函数；
[0152] 当(^= 1、c 2= 1、c 3= 1时，构造出互乘111窗函数；
[0153] 当Ci= 0、c 2= 2、c 3= 1时，构造出互乘021窗函数；
[0154] 当Ci= 0、c 2= 0、c 3= 3时，构造出互乘003窗函数；
[0155] 从上述9种乘法窗函数可以看出：互乘300窗函数、互乘030窗函数、互乘003窗 函数，这三种窗函数实际上属于自乘法窗函数（本发明不予讨论），其余6种为构造出的互 乘法窗函数。
[0156] S2、信号预处理：
[0157] 互感器采集电网信号，电网信号包括电流信号、电压信号，将互感器采集到的电网 信号x(n)，传输到上位机；对电网信号x(n)进行加互乘法窗函数w(n)进行截断，得到加窗 信号x w(n):
[0158] xw (n) = x (n) w (n) (2)
[0159] 对公式⑵的加窗信号进行FFT变换后，得到加窗FFT频谱：
[0161] 其中 ，W(?)为窗函数的频谱，k为正整数，X(k)表示第k次谐波的频谱，Ak为第 k次谐波的幅值，j表示虚数单位，e是自然对数的底数，％t为第k次谐波的初始相位，第一 次谐波为基波，fs为采样频率，f〇为基波频率，A f为离散频率间隔，且A f = f s/N。
[0162] 令：1^=h/A f，&为真实频谱的谱线位置，由于非同步采样或非整周期截断，h 一般为小数，不为整数。
[0163] 若忽略负频率点处旁瓣的影响，公式（3)变为：
[0165] S3、确定三根谱线：
[0166] 在S2得到的加窗FFT频谱峰值附近区域，心处频率点较大的三根谱线分别为：第 根，k 为正整数，k广k3的关系为：k2=k汴1，k3= k2+1，这三根谱线对应 的幅值分别为y2、y3。
[0167] 为方便计算，记变量a=k_k2,则-0. 5彡a彡〇.5;
[0168] 另记变量
[0169] S4、计算三谱线插值算法的修正公式：
[0170] 根据公式（4)和（5)可以得到：
[0172] 当N值比较大，例如：N>1000时，公式（6)可以化简为一个函数0 = g(a)，其 反函数为a =g4(f3)。由于所采用的余弦窗系数均为实系数，其频率响应是偶对称的，因 而g( ?)和g<( ?)均为奇函数。
[0173] 采用多项式逼近方法计算奇函数0 =g4(a)，表达式为：
[0174] a~pnX 0+p13X 03+…+plp0p(7)
[0175] pn，p13; 多项式逼近的奇次项系数，p是奇数。
[0176]根据公式（4)，求得电网信号第i次谐波的幅值Ai:
[0178] 其中，i为正整数，yi为加窗FFT后第i次谐波的幅值；
[0179] 以基波为例，考虑到y2是离真实谱线点最近谱线，给予较大权重，可以得到：
[0181] 类似公式（7)的逼近方法，当N比较大，例如N> 1000时，窗函数系数为实系数， 公式（9)可表示为：
[0182] Aj= N_1 (y3+2y2+y!)u(a )
[0183] 其中，u(a)为修正公式，且为偶函数，逼近多项式不含奇次项。
[0184] 三谱线修正逼近多项式如下：
[0185] u ( a ) = (p2〇+p22 a 2+…+p2d a d) (10)
[0186] 公式（10)中，p2(l，p22…p2d为多项式逼近的偶次项系数，d为拟合的最高阶次，且d 为偶数。
[0187] S5、计算基波参数：
[0188] 考虑到y2是离真实谱线点最近的两根谱线，给予较大权重，可以计算基波频率fV 基波幅值A ::
[0189] f0= k? Af = (k2+a)Af (11)
[0190] A! = N-1 (p2CI+p22 a 2+... +p2d a d) (12)
[0191] 根据公式（4)，还可以计算基波的相位：
[0193] 仿照基波参数的求取，根据公式（6)、（7)、（9)、（11)、（12)、（13)即可进行各次谐 波参数的分析。
[0194] S6、确定谐波参数：确定基波频率&后，在范围（kf ^5, kf^+5)内，重复步骤S3~ S5,直到所有谐波参数计算完毕；
[0195] S7、进行误差分析：分析基于互乘法窗函数的三谱线插值FFT算法的误差，并与常 规窗函数插值算法精度进行比较。
[0196] 下面通过一个具体实施例来进行详细说明。选定一个电网信号的谐波，谐波的波 形参见图1所示。
[0197] 本发明实施例提供一种基于互乘法窗函数的三谱线插值FFT谐波分析方法，具体 包括以下步骤：
[0198] 步骤1、互乘法窗函数的构造：
[0199] 以Hanning窗、Hamming窗、Blackman窗为例，给出互乘法窗函数的构造模式，构 造的6种乘法窗函数性能如表2所示。考虑到计算量的问题，乘法窗函数的阶次不适宜太 高，因此限定阶次c = 3。每种窗函数的子阶次可能取0、1、2三个值。为了书写的方面，下 面简化表达形式：Hanning简写为Hn，Hamming简写为Hm，Blackman简写为Bm。
[0200] 表1、基于常规函数构造互乘法窗函数的性能
[0202] 步骤2、信号的预处理：
[0203]将互感器采集到的离散电网信号x(n)，信号模型参见表2所示。为了验证所提算 法的精度，进行10次谐波仿真分析。信号模型为：
[0205] 其中：基波频率fQ为50. 5Hz，采样频率匕为5120Hz，采样点数N为1024。
[0206] 表2、电网信号的谐波参数
[0208] 加窗FFT变换后得到信号频谱：
[0209] 步骤3、确定三根谱线：
[0210] 在45~55Hz中寻求模最大三根谱线，分别为yp y2, y3。
[0211] 那么
[0212] 步骤4、计算三谱线插值算法的修正公式：
[0213] 根据【具体实施方式】中步骤S4的推导，表2中的4种窗函数的修正公式a = g4(f3)，u(a)分别如下：
[0214] (1)互乘210窗函数：
[0215] a = 1. 11458715 0 -0? 0908720 0 3+0. 01413650 0 5
[0216] u ( a ) = 1. 80783068+0. 41014621 a 2+〇. 05139583 a 4
[0217] 互乘210窗函数幅频特性参见图3所示。
[0218] (2)互乘201窗函数：
[0219] a = 1. 28076537 0 -0? 09843530 0 3+0. 01499840 0 5
[0220] u ( a ) = 1. 95369960+0. 38423865 a 2+〇. 04102745 a 4
[0221] 互乘201窗函数幅频特性参见图4所示。
[0222] (3)互乘120窗函数：
[0223] a = 1. 08516300 0 -0? 08842400 0 3+0. 0140181 0 5
[0224] u ( a ) = 1. 78667431+0. 41627357 a 2+〇. 05339022 a 4
[0225] 互乘120窗函数的幅频特性参见图5所示。
[0226] (4)互乘102窗函数：
[0227] a = 1. 42354957 0 -0? 10399918 0 3+0. 01605665 0 5
[0228] u ( a ) = 2. 07275558+0. 36590106 a 2+〇. 03458435 a 4
[0229] 互乘102窗函数的幅频特性参见图6所示。
[0230](5)互乘111窗函数：
[0231] a = 1. 2530640 0 -0. 09572246 0 3+0. 01451709 0 5
[0232]u(a) =1.93442536+0. 38883257a 2+〇.04234880a4
[0233] 互乘111窗函数的幅频特性参见图7所示。
[0234] (6)互乘021窗函数：
[0235] a = 1. 22447178 0 -0? 09322158 0 3+0. 01445149 0 5
[0236] u(a) = 1. 91516084+0. 39387119a 2+〇. 04376915a 4
[0237] 互乘021窗函数的幅频特性参见图8所示。
[0238] 步骤5、计算基波参数：
[0239]f0=k?Af=(k2+a)Af
[0240] Af N-1 (p2CI+p22 a 2+... +p2d a d)
[0242] 步骤6、确定谐波参数：确定基波频率&后，在范围（kfV5, kL+5)内，重复步骤 3~5,直到所有谐波参数计算完毕。
[0243] 步骤7、进行误差分析：分析基于互乘法窗函数的三谱线插值FFT算法的误差，并 与常规窗函数插值算法精度进行比较。误差分析比较结果参见表3~表5所示，其中，D Ai 表示基波和各次谐波幅值测量值的相对误差；Df(l表示基波频率测量值的相对误差；?表 示基波和各次谐波初始相位测量值的相对误差，均用百分比表示。
[0244] 表3、互乘法窗幅值、频率测量相对误差
[0246] 表4、互乘法窗相位测量相对误差
[0249] 而单独采用Hanning、Hamming和Blackman窗函数幅值、频率和相位测量误差如表 5所示。
[0250] 表5、Hn、Hm、Bm窗的幅值、频率和相位相对误差
[0252] 由此可见，本发明实施例中的互乘法窗函数在基本频率上，要比基本窗函数高 1~3个数量级，在幅值计算上高2~3个数量级，在相位计算上优势更为突出，高3~6个 数量级。由以上的仿真结果可以看出，计算结果普遍好于普通窗函数谱线插值算法。
[0253] 综上所述，本发明实施例中的基于互乘法窗函数谐波分析方法，首次 实现构造互 乘法窗函数，将互乘法窗函数应用到三谱线插值FFT算法中，计算谐波参数，能够获得比常 规窗函数更高的精度。
[0254] 本发明实施例还提供一种基于互乘法窗函数的三谱线插值FFT谐波分析系统，包 括互乘法窗函数构造单元、信号预处理单元、谱线确定单元、修正公式计算单元、基波参数 计算单元、谐波参数确定单元、误差分析单元，其中：
[0255] 互乘法函数构造单元，用于构造互乘法窗函数：
[0256] 乘法窗函数的通用公式是由多个窗函数乘积产生的，乘法窗函数w(n)的通用公 式为：
[0258] 其中，n为采样点的序数，n为自然数；m为参加构造乘法窗的基本窗函数的个数， m为正整数为第i个基本窗函数表达式，i = 1，2…m七为参加第i个基本窗函数 的个数，称为K (n)的子阶次
，c为互乘法窗函数的总阶次；根据Cl、c 2、(：3的 不同取值，形成不同的组合（ClcvCni)，即构造出不同形式的乘法窗函数；
[0259] 当力⑷表达式相同时，即所采用的基本窗相同时，w(n)为自乘法窗函数，当wjn) 表达式不同时，即所采用的基本窗不同时，w(n)为互乘法窗函数；
[0260] 下面以三个常用窗函数为例，表达式分别*Wl(n)，w2(n)，w 3(n)。总阶次c = Cl+C2+C3。由于C越大，意味着参与乘法窗的基本窗函数越多，构成的乘法窗函数项数越多， 会增加计算的复杂度。一般给C 一个最大值，例如，令：c = 3。
[0261] 根据Cl、c2、c 3的不同取值，形成不同的组合（c &???〇,构造出以下9种乘法窗函 数：
[0262] 当Ci= 3、c 2= 0、c 3= 0时，构造出互乘300窗函数；
[0263] 当Ci= 2、c 2= 1、c 3= 0时，构造出互乘210窗函数；
[0264] 当Ci= 2、c 2= 0、c 3= 1时，构造出互乘201窗函数；
[0265] 当Ci= 1、c 2= 2、c 3= 0时，构造出互乘120窗函数；
[0266] 当Ci= 1、c 2= 0、c 3= 2时，构造出互乘102窗函数；
[0267] 当Ci= 0、c 2= 3、c 3= 0时，构造出互乘030窗函数；
[0268] 当(^= l、c2= l、c3= 1时，构造出互乘111窗函数；
[0269] 当Ci= 0、c 2= 2、c 3= 1时，构造出互乘021窗函数；
[0270] 当Ci= 0、c 2= 0、c 3= 3时，构造出互乘003窗函数；
[0271] 从上述9种乘法窗函数可以看出：互乘300窗函数、互乘030窗函数、互乘003窗 函数，这三种窗函数实际上属于自乘法窗函数（本发明不予讨论），其余6种为构造出的互 乘法窗函数。
[0272] 信号预处理单元，用于对信号进行预处理：
[0273] 互感器采集电网信号，电网信号包括电流信号、电压信号；将互感器采集到的电网 信号x(n)，传输到上位机；对电网信号x(n)进行加互乘法窗函数w(n)进行截断，得到加窗 信号x w(n):
[0274] xw(n) = x(n)w(n) (2)
[0275] 对公式⑵的加窗信号进行FFT变换后，得到加窗FFT频谱：
[0277] 其中，W( ?)为窗函数的频谱，k为正整数，X(k)表示第k次谐波的频谱，Ak为第 k次谐波的幅值，j表示虚数单位，e是自然对数的底数，％t为第k次谐波的初始相位，第一 次谐波为基波，f s为采样频率，f 〇为基波频率，A f为离散频率间隔，且A f = f s/N ;
[0278] 令：kQ= fQ/Af，kQ为真实频谱的谱线位置，由于非同步采样或非整周期截断，kQ 一般为小数，不为整数。
[0279] 忽略负频率点处旁瓣的影响，公式（3)变为：
[0281] 谱线确定单元，用于确定三根谱线：
[0282] 在S2得到的加窗FFT频谱峰值附近区域，匕处频率点较大的三根谱线分别为：第 根，k 为正整数，k k 3的关系为：k 2= k i+1，k3= k 2+1，这三根谱线对应 的幅值分别为y2、
[0283] 记变量 a = k_k2,则-0? 5 彡 a 彡 〇? 5 ;
[0284] 另记变量，
[0285] 修正公式计算单元，用于计算三谱线插值算法的修正公式：
[0286] 根据公式⑷和（5)得到：
[0288] 当N值比较大，例如：N> 1000时，公式（6)可以化简为一个函数0 = g( a)，其 反函数为a =g4(f3)。由于所采用的余弦窗系数均为实系数，其频率响应是偶对称的，因 而g( ?)和g<( ?)均为奇函数。
[0289] 采用多项式逼近方法计算奇函数0 =g4(a)，表达式为：
[0290] a~pnX 0+p13X 03+…+plp0p(7)
[0291] pn，p13;…p lp为多项式逼近的奇次项系数，p是奇数；
[0292]根据公式（4)，求得电网信号第i次谐波的幅值Ai:
[0294] 其中，i为正整数，yi为加窗FFT后第i次谐波的幅值；
[0295] 以基波为例，考虑到y2是离真实谱线点最近谱线，给予较大权重，可以得到：
[0297] N> 1000时，窗函数系数为实系数，公式（9)表示为：
[0298] Aj= N_1 (y3+2y2+y1)u(a )
[0299] 其中，u(a)为修正公式，且为偶函数，逼近多项式不含奇次项；
[0300] 三谱线修正逼近多项式如下：
[0301] u ( a ) = (p2(l+p22 a 2+…+p2d a d) (10)
[0302] 公式（10)中，p2Q，p22~p2$多项式逼近的偶次项系数，d为拟合的最高阶次，且d 为偶数；
[0303] 基波参数计算单元，用于计算基波参数：
[0304]考虑到y2是离真实谱线点最近的两根谱线，给予较大权重，可以计算基波频率f。、 基波幅值A::
[0305] f0= k? Af = (k2+a)Af (11)
[0306] A! = N-1 (p2CI+p22 a 2+... +p2d a d) (12)
[0307]根据公式（4)，计算基波的相位：
[0309] 仿照基波参数的求取，根据公式（6)、（7)、（9)、（11)、（12)、（13)进行各次谐波参 数的分析；
[0310] 谐波参数确定单元，用于确定谐波参数：确定基波频率&后，在范围 (kfd-5, kf^+5)内，所述谱线确定单元、修正公式计算单元、基波参数计算单元重复进行计 算，直到所有谐波参数计算完毕；
[0311] 误差分析单元，用于进行误差分析：分析基于互乘法窗函数的三谱线插值FFT算 法的误差，并与常规窗函数插值算法精度进行比较。
[0312] 本领域的技术人员可以对本发明实施例进行各种修改和变型，倘若这些修改和变 型在本发明权利要求及其等同技术的范围之内，则这些修改和变型也在本发明的保护范围 之内。
[0313] 说明书中未详细描述的内容为本领域技术人员公知的现有技术。
【主权项】
1. 一种基于互乘法窗函数的三谱线插值FFT谐波分析方法，其特征在于，包括以下步 骤： 51、 构造互乘法窗函数： 乘法窗函数的通用公式是由多个窗函数乘积产生的，乘法窗函数w(n)的通用公式为：(1) 其中，η为采样点的序数，η为自然数;m为参加构造乘法窗的基本窗函数的个数，m为 正整数％ (η)为第i个基本窗函数表达式，i = 1，2···πι ;(^为参加第i个基本窗函数的个 数，称为Wi (η)的子阶次，c为互乘法窗函数的总阶次；根据Cl、c 2、C3的不同 取值，形成不同的组合（cl(V·· cm)，即构造出不同形式的乘法窗函数； 当Wi(Ii)表达式相同时，即所采用的基本窗相同时，w(n)为自乘法窗函数，当Wi(Ii)表 达式不同时，即所采用的基本窗不同时，w(n)为互乘法窗函数； 三个常用窗函数的表达式分别为W1(Ii)，W2 (n)，W3 (η)，总阶次C = CfCfC3，给C 一个最 大值，令：c = 3 ; 根据Cl、c2、C3的不同取值，形成不同的组合（c而…cm)，构造出以下9种乘法窗函数： 当C1= 3、C2= 0、c3= O时，构造出互乘300窗函数； 当C1= 2、C2= Uc3= 0时，构造出互乘210窗函数； 当C1= 2、C2= 0、c3= 1时，构造出互乘201窗函数； 当C1= Uc2= 2、C3= 0时，构造出互乘120窗函数； 当C1= Uc2= (Kc3= 2时，构造出互乘102窗函数； 当C1= (Kc2= 3、C3= 0时，构造出互乘030窗函数； 当C1= Uc2= Uc3= 1时，构造出互乘111窗函数； 当C1= (Kc2= 2、C3= 1时，构造出互乘021窗函数； 当C1= 0、c2= (Kc3= 3时，构造出互乘003窗函数； 其中，互乘300窗函数、互乘030窗函数、互乘003窗函数，这三种窗函数属于自乘法窗 函数，其余6种为构造出的互乘法窗函数； 52、 信号预处理： 互感器采集电网信号，将互感器采集到的电网信号x( n)，传输到上位机；对电网信号 X(H)进行加互乘法窗函数w (η)进行截断，得到加窗信号Xw (η): xw (η) = X (n) w (η) (2) 对公式（2)的加窗信号进行FFT变换后，得到加窗FFT频谱：其中，W( ·)为窗函数的频谱，k为正整数，X(k)表示第k次谐波的频谱，Ak为第k次 谐波的幅值，j表示虚数单位，e是自然对数的底数，妁为第k次谐波的初始相位，第一次谐 波为基波，fs为采样频率，f 0为基波频率，Δ f为离散频率间隔，且Λ f = f s/N ; 令：1?= fd/Δ f，1?为真实频谱的谱线位置， 忽略负频率点处旁瓣的影响，公式（3)变为：(4) 53、 确定三根谱线： 在S2得到的加窗FFT频谱峰值附近区域，1?处频率点较大的三根谱线分别为：第k i、 k2、k3根，k η k2、k3均为正整数，k广k 3的关系为：k 2= k i+1，k3= k 2+1，这三根谱线对应的 幅值分别为y2、y3; 记变量 α = k_k2，则-〇· 5 < α < 〇· 5 ; 另记变量（5) 54、 计算三谱线插值算法的修正公式： 根据公式（4)和（5)得到：(6) 采用多项式逼近方法计算奇函数β = Ρ(α)，表达式为： α ~ ρηΧ β+ρ13Χ β 3+…+ρ1ρβρ (7) Pn，P13;…P 1Ρ为多项式逼近的奇次项系数，P是奇数； 根据公式（4)，求得电网信号第i次谐波的幅值Ai:C8) 其中，i为正整数，Yi为加窗FFT后第i次谐波的幅值； 考虑到y2是离真实谱线点最近谱线，得到：(9) N > 1000时，窗函数系数为实系数，公式（9)表示为： A1= N (y3+2y2+y1)u( α ) 其中，u(a)为修正公式，且为偶函数，逼近多项式不含奇次项； 三谱线修正逼近多项式如下： u ( a ) = (ρ2(ι+ρ22 α 2+…+p2d a d) (10) 公式（10)中，P2(I，P22··· P2A多项式逼近的偶次项系数，d为拟合的最高阶次，且d为偶 数； 55、 计算基波参数： 计算基波频率fci、基波幅值A1: f0=k· Af = (k2+a)Af (11) A1 = N-1 (y^Sy^y!) (p2CI+p22 a 2+... +p2da d) (12) 根据公式（4)，计算基波的相位：(13) 仿照基波参数的求取，根据公式（6)、（7)、（9)、（11)、（12)、（13)进行各次谐波参数的 分析； 56、 确定谐波参数：确定基波频率&后，在范围内，重复步骤S3~S5,直 到所有谐波参数计算完毕； 57、 进行误差分析：分析基于互乘法窗函数的三谱线插值FFT算法的误差，并与常规窗 函数插值算法精度进行比较。2. 如权利要求1所述的基于互乘法窗函数的三谱线插值FFT谐波分析方法，其特征在 于：所述电网信号包括电流信号、电压信号。3. 如权利要求1所述的基于互乘法窗函数的三谱线插值FFT谐波分析方法，其特征在 于：所述真实频谱的谱线位置1?为小数。4. 一种基于互乘法窗函数的三谱线插值FFT谐波分析系统，其特征在于：包括互乘法 窗函数构造单元、信号预处理单元、谱线确定单元、修正公式计算单元、基波参数计算单元、 谐波参数确定单元、误差分析单元，其中： 所述互乘法函数构造单元，用于构造互乘法窗函数： 乘法窗函数的通用公式是由多个窗函数乘积产生的，乘法窗函数w(n)的通用公式为：(1) 其中，η为采样点的序数，η为自然数;m为参加构造乘法窗的基本窗函数的个数，m为 正整数％ (η)为第i个基本窗函数表达式，i = 1，2···πι ;(^为参加第i个基本窗函数的个 数，称为Wi(Il)的子阶次：c为互乘法窗函数的总阶次；根据Cl、c 2、C3的不同 取值，形成不同的组合（cl(V·· cm)，即构造出不同形式的乘法窗函数； 当Wi(Ii)表达式相同时，即所采用的基本窗相同时，w(n)为自乘法窗函数，当Wi(Ii)表 达式不同时，即所采用的基本窗不同时，w(n)为互乘法窗函数； 三个常用窗函数的表达式分别为W1(Ii)，W2 (n)，W3 (η)，总阶次C = CfCfC3，给C 一个最 大值，令：c = 3 ; 根据Cl、c2、C3的不同取值，形成不同的组合（c而…cm)，构造出以下9种乘法窗函数： 当C1= 3、C2= 0、c3= O时，构造出互乘300窗函数； 当C1= 2、C2= Uc3= 0时，构造出互乘210窗函数； 当C1= 2、C2= 0、c3= 1时，构造出互乘201窗函数； 当C1= Uc2= 2、C3= 0时，构造出互乘120窗函数； 当C1= Uc2= (Kc3= 2时，构造出互乘102窗函数； 当C1= (Kc2= 3、C3= 0时，构造出互乘030窗函数； 当C1= Uc2= Uc3= 1时，构造出互乘111窗函数； 当C1= (Kc2= 2、C3= 1时，构造出互乘021窗函数； 当C1= 0、c2= (Kc3= 3时，构造出互乘003窗函数； 其中，互乘300窗函数、互乘030窗函数、互乘003窗函数，这三种窗函数属于自乘法窗 函数，其余6种为构造出的互乘法窗函数； 所述信号预处理单元，用于对信号进行预处理： 互感器采集电网信号，将互感器采集到的电网信号x(n)，传输到上位机；对电网信号 X(H)进行加互乘法窗函数w (η)进行截断，得到加窗信号Xw (η): xw (η) = X (n) w (η) (2) 对公式（2)的加窗信号进行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), 所述谱线确定单元，用于确定三根谱线： 在S2得到的加窗FFT频谱峰值附近区域，1?处频率点较大的三根谱线分别为：第k i、 k2、k3根，k η k2、k3均为正整数，k广k 3的关系为：k 2= k i+1，k3= k 2+1，这三根谱线对应的 幅值分别为y2、y3; 记变量 α = k_k2，则-〇· 5 < α < 〇· 5 ; 另记变量（5); 所述修正公式计算单元，用于计算三谱线插值算法的修正公式： 根据公式（4)和（5)得到：(6) 采用多项式逼近方法计算奇函数β = Ρ(α)，表达式为： α ~ ρηΧ β+ρ13Χ β 3+…+ρ1ρβρ (7) Pn，P13;…P 1Ρ为多项式逼近的奇次项系数，P是奇数； 根据公式（4)，求得电网信号第i次谐波的幅值Ai:C8) 其中，i为正整数，Yi为加窗FFT后第i次谐波的幅值； 考虑到y2是离真实谱线点最近谱线，得到：(9) N > 1000时，窗函数系数为实系数，公式（9)表示为： A1= N (y3+2y2+y1)u( α ) 其中，u(a)为修正公式，且为偶函数，逼近多项式不含奇次项； 三谱线修正逼近多项式如下： u ( a ) = (ρ2(ι+ρ22 α 2+…+p2d a d) (10) 公式（10)中，P2(I，P22··· P2A多项式逼近的偶次项系数，d为拟合的最高阶次，且d为偶 数； 所述基波参数计算单元，用于计算基波参数： 计算基波频率fci、基波幅值A1: f0=k* Af= (k2+a)Af (11) A1 = N-1 (p2CI+p22 a 2+... +p2da d) (12) 根据公式（4)，计算基波的相位：(13) 仿照基波参数的求取，根据公式（6)、（7)、（9)、（11)、（12)、（13)进行各次谐波参数的 分析； 所述谐波参数确定单元，用于确定谐波参数：确定基波频率后，在范围 (kfd-5, kf^+5)内，所述谱线确定单元、修正公式计算单元、基波参数计算单元重复进行计 算，直到所有谐波参数计算完毕； 所述误差分析单元，用于进行误差分析：分析基于互乘法窗函数的三谱线插值FFT算 法的误差，并与常规窗函数插值算法精度进行比较。5. 如权利要求4所述的基于互乘法窗函数的三谱线插值FFT谐波分析系统，其特征在 于：所述电网信号包括电流信号、电压信号。6. 如权利要求4所述的基于互乘法窗函数的三谱线插值FFT谐波分析系统，其特征在 于：所述真实频谱的谱线位置1?为小数。
【专利摘要】本发明公开了一种基于互乘法窗函数的三谱线插值FFT谐波分析方法及系统，涉及谐波分析领域。该方法包括以下步骤：构造互乘法窗函数；信号预处理；确定三根谱线；计算三谱线插值算法的修正公式；计算基波参数；确定谐波参数；进行误差分析。本发明提出一种新的窗函数的构造方式，利用常规窗函数进行不同的乘法组合，首次实现构造互乘法窗函数，将互乘法窗函数应用到三谱线插值FFT算法中，进行谐波分析，计算谐波参数。多种常用的余弦窗函数计算实例表明，本发明构造出的互乘法窗函数相对于常规窗函数插值算法，有更高的准确度，能实现谐波的高精度测量。
【IPC分类】G01R23/16
【公开号】CN104897961
【申请号】CN201510335000
【发明人】张俊敏, 刘开培, 汪立, 王黎, 田微, 陈文娟
【申请人】中南民族大学
【公开日】2015年9月9日
【申请日】2015年6月17日
【技术领域】
[0001] 本发明涉及谐波分析领域，具体是涉及一种基于互乘法窗函数的三谱线插值FFT 谐波分析方法及系统。
【背景技术】
[0002] 随着大量非线性电力电子器件的使用，使得电网中谐波污染问题日益严重。谐波 问题不仅恶化电能质量，对电网的安全稳定和经济运行也造成较大影响。因此，对系统中谐 波参数进行高精度测量，对于减少谐波危害，维护电网安全稳定、高效运行是十分必要的。
[0003] 谐波检测的关键问题是：非同步采样下，如何解决频谱泄露和栅栏效应。目前，常 见的检测谐波和间谐波的方法主要有：傅里叶变换，机器学习法，小波变换，功率谱估计，希 尔伯特-黄变换。
[0004] 在傅里叶变换分析谐波中，所用方法是：运用各种特殊窗函数，对信号进行截断， 然后结合谱线插值FFT (Fast Fourier Transform，快速傅里叶变换）进行谐波分析。
[0005] 常用的窗函数有：汉宁（Hanning)窗函数、布莱克曼（Blackman)窗函数、布莱克曼 汉斯（Blackman-Harris)窗函数、纳托尔（Nuttall)窗函数、莱夫文森特（Rife-Vincent) 窗函数以及各种组合窗函数。
[0006] 在加窗基础上，D. Agrez和庞浩等人各自提出了双谱线的修正算法，Wu Jing、牛胜 锁和黄冬梅等人提出了三谱线修正算法。这些改进降低了频谱泄漏和栅栏效应的影响，提 高了谐波分析的准确性。
[0007] 然而，在工程实际使用中，常用的窗函数插值算法仍然无法满足高精度的谐波分 析要求。
【发明内容】
[0008] 本发明的目的是为了克服上述【背景技术】的不足，提供一种基于互乘法窗函数的三 谱线插值FFT谐波分析方法及系统，与常规窗函数比，具有更高的精度。
[0009] 本发明提供一种基于互乘法窗函数的三谱线插值FFT谐波分析方法，包括以下步 骤：
[0010] S1、构造互乘法窗函数：
[0011] 乘法窗函数的通用公式是由多个窗函数乘积产生的，乘法窗函数W(n)的通用公 式为：
[0013] 其中，n为采样点的序数，n为自然数；m为参加构造乘法窗的基本窗函数的个数， m为正整数为第i个基本窗函数表达式，i = 1，2…m七为参加第i个基本窗函数 的个数，称为Wi (n)的子阶次
，C为互乘法窗函数的总阶次；根据Cl、c 2、(33的 不同取值，形成不同的组合（c^cv-cj，即构造出不同形式的乘法窗函数；
[0014] 当^〇1)表达式相同时，即所采用的基本窗相同时，w(n)为自乘法窗函数，当wjn) 表达式不同时，即所采用的基本窗不同时，w(n)为互乘法窗函数；
[0015] 三个常用窗函数的表达式分别为％ (n)，w2 (n)，w3(n)，总阶次c = Wq，给c 一 个最大值，令：c = 3 ;
[0016] 根据Ci、c2、c3的不同取值，形成不同的组合（cffcj，构造出以下9种乘法窗函 数：
[0017] 当Ci= 3、c 2= 0、c 3= 0时，构造出互乘300窗函数；
[0018] 当(^=2、(32= l、c3=0时，构造出互乘210窗函数；
[0019] 当(^=2、(32=0、(33= 1时，构造出互乘201窗函数；
[0020] 当Ci= 1、c 2= 2、c 3= 0时，构造出互乘120窗函数；
[0021] 当~二l、c2=0、c3=2时，构造出互乘102窗函数；
[0022] 当Ci= 0、c 2= 3、c 3= 0时，构造出互乘030窗函数；
[0023] 当(^= l、c2= l、c3= 1时，构造出互乘111窗函数；
[0024] 当Ci= 0、c 2= 2、c 3= 1时，构造出互乘021窗函数；
[0025] 当Ci= 0、c 2= 0、c 3= 3时，构造出互乘003窗函数；
[0026] 其中，互乘300窗函数、互乘030窗函数、互乘003窗函数，这三种窗函数属于自乘 法窗函数，其余6种为构造出的互乘法窗函数；
[0027] S2、信号预处理：
[0028] 互感器采集电网信号，将互感器采集到的电网信号x(n)，传输到上位机；对电网 信号x(n)进行加互乘法窗函数w(n)进行截断，得到加窗信号\〇1):
[0029] xw(n) = x(n)w(n) (2)
[0030] 对公式⑵的加窗信号进行FFT变换后，得到加窗FFT频谱：
[0032] 其中，W( ?)为窗函数的频谱，k为正整数，X(k)表示第k次谐波的频谱，Ak为第 k次谐波的幅值，j表示虚数单位，e是自然对数的底数，外;为第k次谐波的初始相位，第一 次谐波为基波，fs为采样频率，f〇为基波频率，A f为离散频率间隔，且A f = f s/N ;
[0033] 令：k0= f。/Af，k。为真实频谱的谱线位置，
[0034]忽略负频率点处旁瓣的影响，公式（3)变为：
[0036] S3、确定三根谱线：
[0037] 在S2得到的加窗FFT频谱峰值附近区域，&处频率点较大的三根谱线分别为：第 V k2、k3根，k i、k2、k3均为正整数，k k 3的关系为：k2 = k i+1，k3= k 2+1，这三根谱线对 应的幅值分别为yp y2、
[0038] 记变量 a = k_k2,则-0? 5 彡 a 彡 〇? 5 ;
[0039] 另记变量
[0040] S4、计算三谱线插值算法的修正公式：
[0041] 根据公式⑷和（5)得到：
[0043] 采用多项式逼近方法计算奇函数0 =g4(a)，表达式为：
[0044] a ~ pnX 0+p13X 03+…+plp0p (7)
[0045] pn，p13; -p 1￡)为多项式逼近的奇次项系数，p是奇数；
[0046] 根据公式⑷，求得电网信号第i次谐波的幅值Ai:
[0048] 其中，i为正整数，yi为加窗FFT后第i次谐波的幅值；
[0049] 考虑到y2是离真实谱线点最近谱线，得到：
[0051] N > 1000时，窗函数系数为实系数，公式（9)表示为：
[0052] Aj= N_1 (y3+2y2+y1)u( a )
[0053] 其中，u(a)为修正公式，且为偶函数，逼近多项式不含奇次项；
[0054] 三谱线修正逼近多项式如下：
[0055] u ( a ) = (p2〇+p22 a 2+…+p2d a d) (10)
[0056] 公式（10)中，p2(l，py p2$多项式逼近的偶次项系数，d为拟合的最高阶次，且d 为偶数；
[0057] S5、计算基波参数：
[0058] 计算基波频率fQ、基波幅值Ay
[0059] f0= k? Af = (k2+a)Af (11)
[0060] A! = N-1 (p2CI+p22 a 2+... +p2d a d) (12)
[0061]根据公式（4)，计算基波的相位：
[0063] 仿照基波参数的求取，根据公式（6)、（7)、（9)、（11)、（12)、（13)进行各次谐波参 数的分析；
[0064] S6、确定谐波参数：确定基波频率fQ后，在范围（kf。_5, kfQ+5)内，重复步骤S3~ S5,直到所有谐波参数计算完毕；
[0065] S7、进行误差分析：分析基于互乘法窗函数的三谱线插值FFT算法的误差，并与常 规窗函数插值算法精度进行比较。
[0066] 在上述技术方案的基础上，所述电网信号包括电流信号、电压信号。
[0067] 在上述技术方案的基础上，所述真实频谱的谱线位置&为小数。
[0068] 本发明还提供一种基于互乘法窗函数的三谱线插值FFT谐波分析系统，包括互乘 法窗函数构造单元、信号预处理单元、谱线确定单元、修正公式计算单元、基波参数计算单 元、谐波参数确定单元、误差分析单元，其中：
[0069] 所述互乘法函数构造单元，用于构造互乘法窗函数：
[0070] 乘法窗函数的通用公式是由多个窗函数乘积产生的，乘法窗函数w(n)的通用公 式为：
[0072] 其中，n为采样点的序数，n为自然数；m为参加构造乘法窗的基本窗函数的个数， m为正整数为第i个基本窗函数表达式，i = 1，2…m七为参加第i个基本窗函数 的个数，称为K (n)的子阶次；
，c为互乘法窗函数的总阶次；根据Cl、c 2、(：3的 不同取值，形成不同的组合（c^cv-cj，即构造出不同形式的乘法窗函数；
[0073] 当^〇1)表达式相同时，即所采用的基本窗相同时，w(n)为自乘法窗函数，当wjn) 表达式不同时，即所采用的基本窗不同时，w(n)为互乘法窗函数；
[0074] 三个常用窗函数的表达式分别为％ (n)，w2 (n)，w3(n)，总阶次c = Ci+c2+c3，给c 一 个最大值，令：c = 3 ;
[0075] 根据Cl、c2、c3的不同取值，形成不同的组合（ Clcvcm)，构造出以下9种乘法窗函 数 ：
[0076] 当Ci= 3、c 2= 0、c 3= 0时，构造出互乘300窗函数；
[0077] 当Ci= 2、c 2= 1、c 3= 0时，构造出互乘210窗函数；
[0078] 当Ci= 2、c 2= 0、c 3= 1时，构造出互乘201窗函数；
[0079] 当Ci= 1、c 2= 2、c 3= 0时，构造出互乘120窗函数；
[0080] 当Ci= 1、c 2= 0、c 3= 2时，构造出互乘102窗函数；
[0081] 当Ci= 0、c 2= 3、c 3= 0时，构造出互乘030窗函数；
[0082] 当(^= l、c2= l、c3= 1时，构造出互乘111窗函数；
[0083] 当Ci= 0、c 2= 2、c 3= 1时，构造出互乘021窗函数；
[0084] 当Ci= 0、c 2= 0、c 3= 3时，构造出互乘003窗函数；
[0085] 其中，互乘300窗函数、互乘030窗函数、互乘003窗函数，这三种窗函数属于自乘 法窗函数，其余6种为构造出的互乘法窗函数；
[0086] 所述信号预处理单元，用于对信号进行预处理：
[0087] 互感器采集电网信号，将互感器采集到的电网信号x(n)，传输到上位机；对电网 信号x(n)进行加互乘法窗函数w(n)进行截断，得到加窗信号\〇1):
[0088] xw(n) = x(n)w(n) (2)
[0089] 对公式⑵的加窗信号进行FFT变换后，得到加窗FFT频谱：
[0091] 其中，W( ?)为窗函数的频谱，k为正整数，X(k)表示第k次谐波的频谱，Ak为第 k次谐波的幅值，j表示虚数单位，e是自然对数的底数，？^为第k次谐波的初始相位，第一 次谐波为基波，f s为采样频率，f〇为基波频率，A f为离散频率间隔，且A f = f s/N ;
[0092] 令：k0= f。/ A f，k。为真实频谱的谱线位置，
[0093] 忽略负频率点处旁瓣的影响，公式（3)变为：
[0095] 所述谱线确定单元，用于确定三根谱线：
[0096] 在S2得到的加窗FFT频谱峰值附近区域，h处频率点较大的三根谱线分别为：第 根，k 为正整数，k广k 3的关系为：k 2= k汴1，k3= k 2+1，这三根谱线对应 的幅值分别为y2、y3;
[0097] 记变量 a = k_k2，则 _0? 5 < a < 〇? 5 ;
[0098] 另记变量
[0099] 所述修正公式计算单元，用于计算三谱线插值算法的修正公式：
[0100] 根据公式⑷和（5)得到：
[0102] 采用多项式逼近方法计算奇函数0 =g4(a)，表达式为：
[0103] a ~ pnX 0+p13X 03+…+plp0p (7)
[0104] pn，p13;…p lp为多项式逼近的奇次项系数，p是奇数；
[0105] 根据公式（4)，求得电网信号第i次谐波的幅值Ai:
[0107] 其中，i为正整数，yi为加窗FFT后第i次谐波的幅值；
[0108] 考虑到y2是离真实谱线点最近谱线，得到：
[0110] N > 1000时，窗函数系数为实系数，公式（9)表示为：
[0111] Aj= N_1 (y3+2y2+y1)u( a )
[0112] 其中，u(a)为修正公式，且为偶函数，逼近多项式不含奇次项；
[0113] 三谱线修正逼近多项式如下：
[0114] u ( a ) = (p2Q+p22 a 2+…+p2d a d) (10)
[0115] 公式（10)中，p2Q，p22…p 2d为多项式逼近的偶次项系数，d为拟合的最高阶次，且d 为偶数；
[0116] 所述基波参数计算单元，用于计算基波参数：
[0117] 计算基波频率fQ、基波幅值Ay
[0118] f0= k? Af = (k2+a)Af (11)
[0119] Ai= N 1 (y3+2y2+y) (p20+p22 a 2+…+p2d a d) (12)
[0120]根据公式（4)，计算基波的相位：
[0122] 仿照基波参数的求取，根据公式（6)、（7)、（9)、（11)、（12)、（13)进行各次谐波参 数的分析；
[0123] 所述谐波参数确定单元，用于确定谐波参数：确定基波频率4后，在范围 (kfd-5, kf^+5)内，所述谱线确定单元、修正公式计算单元、基波参数计算单元重复进行计 算，直到所有谐波参数计算完毕；
[0124] 所述误差分析单元，用于进行误差分析：分析基于互乘法窗函数的三谱线插值 FFT算法的误差，并与常规窗函数插值算法精度进行比较。
[0125] 在上述技术方案的基础上，所述电网信号包括电流信号、电压信号。
[0126] 在上述技术方案的基础上，所述真实频谱的谱线位置k0为小数。
[0127] 与现有技术相比，本发明的优点如下：
[0128] 本发明提出一种新的窗函数的构造方式，利用常规窗函数进行不同的乘法组合， 首次实现构造互乘法窗函数，将互乘法窗函数应用到三谱线插值FFT算法中，进行谐波分 析，计算谐波参数。多种常用的余弦窗函数计算实例表明，本发明构造出的互乘法窗函数相 对于常规窗函数插值算法，有更高的准确度，实现了谐波的高精度测量，具有较高的实用价 值。
【附图说明】
[0129] 图1是本发明实施例中基于互乘法窗函数的三谱线插值FFT谐波分析方法的流程 图。
[0130] 图2是谐波信号的波形图。
[0131] 图3是互乘210窗函数的幅频特性图。
[0132] 图4是互乘201窗函数的幅频特性图。
[0133] 图5是互乘120窗函数的幅频特性图。
[0134] 图6是互乘102窗函数的幅频特性图。
[0135] 图7是互乘111窗函数的幅频特性图。
[0136] 图8是互乘021窗函数的幅频特性图。
【具体实施方式】
[0137] 下面结合附图及具体实施例对本发明作进一步的详细描述。
[0138] 参见图1所示，本发明实施例提供一种基于互乘法窗函数的三谱线插值FFT谐波 分析方法，包括以下步骤：
[0139] S1、构造互乘法窗函数：
[0140] 乘法窗函数的通用公式是由多个窗函数乘积产生的，乘法窗函数w(n)的通用公 式为：
[0142] 其中，n为采样点的序数，n为自然数；m为参加构造乘法窗的基本窗函数的个数，m 为正整数;Wi (n)为第i个基本窗函数表达式，i = 1，2…m力为参加第i个基本窗函数的 个数，称为wjn)的子阶次：
，c为互乘法窗函数的总阶次。根据Cl、c 2、(：3的 不同取值，形成不同的组合（c^cv-cj，即构造出不同形式的乘法窗函数。
[0143] 当^〇1)表达式相同时，即所采用的基本窗相同时，w(n)为自乘法窗函数，当wjn) 表达式不同时，即所采用的基本窗不同时，w(n)为互乘法窗函数。
[0144] 下面以三个常用窗函数为例，表达式分别*Wl(n)，w2(n)，w 3(n)。总阶次c = Cl+C2+C3。由于C越大，意味着参与乘法窗的基本窗函数越多，构成的乘法窗函数项数越多， 会增加计算的复杂度。一般给C 一个最大值，例如，令：c = 3。
[0145] 根据CpCyq的不同取值，形成不同的组合（c而…cm)，构造出以下9种乘法窗函 数：
[0146] 当Ci= 3、c 2= 0、c 3= 0时，构造出互乘300窗函数；
[0147] 当Ci = 2、c 2= 1、c 3= 0时，构造出互乘210窗函数；
[0148] 当Ci= 2、c 2= 0、c 3= 1时，构造出互乘201窗函数；
[0149] 当Ci= 1、c 2= 2、c 3= 0时，构造出互乘120窗函数；
[0150] 当Ci= 1、c 2= 0、c 3= 2时，构造出互乘102窗函数；
[0151] 当Ci= 0、c 2= 3、c 3= 0时，构造出互乘030窗函数；
[0152] 当(^= 1、c 2= 1、c 3= 1时，构造出互乘111窗函数；
[0153] 当Ci= 0、c 2= 2、c 3= 1时，构造出互乘021窗函数；
[0154] 当Ci= 0、c 2= 0、c 3= 3时，构造出互乘003窗函数；
[0155] 从上述9种乘法窗函数可以看出：互乘300窗函数、互乘030窗函数、互乘003窗 函数，这三种窗函数实际上属于自乘法窗函数（本发明不予讨论），其余6种为构造出的互 乘法窗函数。
[0156] S2、信号预处理：
[0157] 互感器采集电网信号，电网信号包括电流信号、电压信号，将互感器采集到的电网 信号x(n)，传输到上位机；对电网信号x(n)进行加互乘法窗函数w(n)进行截断，得到加窗 信号x w(n):
[0158] xw (n) = x (n) w (n) (2)
[0159] 对公式⑵的加窗信号进行FFT变换后，得到加窗FFT频谱：
[0161] 其中 ，W(?)为窗函数的频谱，k为正整数，X(k)表示第k次谐波的频谱，Ak为第 k次谐波的幅值，j表示虚数单位，e是自然对数的底数，％t为第k次谐波的初始相位，第一 次谐波为基波，fs为采样频率，f〇为基波频率，A f为离散频率间隔，且A f = f s/N。
[0162] 令：1^=h/A f，&为真实频谱的谱线位置，由于非同步采样或非整周期截断，h 一般为小数，不为整数。
[0163] 若忽略负频率点处旁瓣的影响，公式（3)变为：
[0165] S3、确定三根谱线：
[0166] 在S2得到的加窗FFT频谱峰值附近区域，心处频率点较大的三根谱线分别为：第 根，k 为正整数，k广k3的关系为：k2=k汴1，k3= k2+1，这三根谱线对应 的幅值分别为y2、y3。
[0167] 为方便计算，记变量a=k_k2,则-0. 5彡a彡〇.5;
[0168] 另记变量
[0169] S4、计算三谱线插值算法的修正公式：
[0170] 根据公式（4)和（5)可以得到：
[0172] 当N值比较大，例如：N>1000时，公式（6)可以化简为一个函数0 = g(a)，其 反函数为a =g4(f3)。由于所采用的余弦窗系数均为实系数，其频率响应是偶对称的，因 而g( ?)和g<( ?)均为奇函数。
[0173] 采用多项式逼近方法计算奇函数0 =g4(a)，表达式为：
[0174] a~pnX 0+p13X 03+…+plp0p(7)
[0175] pn，p13; 多项式逼近的奇次项系数，p是奇数。
[0176]根据公式（4)，求得电网信号第i次谐波的幅值Ai:
[0178] 其中，i为正整数，yi为加窗FFT后第i次谐波的幅值；
[0179] 以基波为例，考虑到y2是离真实谱线点最近谱线，给予较大权重，可以得到：
[0181] 类似公式（7)的逼近方法，当N比较大，例如N> 1000时，窗函数系数为实系数， 公式（9)可表示为：
[0182] Aj= N_1 (y3+2y2+y!)u(a )
[0183] 其中，u(a)为修正公式，且为偶函数，逼近多项式不含奇次项。
[0184] 三谱线修正逼近多项式如下：
[0185] u ( a ) = (p2〇+p22 a 2+…+p2d a d) (10)
[0186] 公式（10)中，p2(l，p22…p2d为多项式逼近的偶次项系数，d为拟合的最高阶次，且d 为偶数。
[0187] S5、计算基波参数：
[0188] 考虑到y2是离真实谱线点最近的两根谱线，给予较大权重，可以计算基波频率fV 基波幅值A ::
[0189] f0= k? Af = (k2+a)Af (11)
[0190] A! = N-1 (p2CI+p22 a 2+... +p2d a d) (12)
[0191] 根据公式（4)，还可以计算基波的相位：
[0193] 仿照基波参数的求取，根据公式（6)、（7)、（9)、（11)、（12)、（13)即可进行各次谐 波参数的分析。
[0194] S6、确定谐波参数：确定基波频率&后，在范围（kf ^5, kf^+5)内，重复步骤S3~ S5,直到所有谐波参数计算完毕；
[0195] S7、进行误差分析：分析基于互乘法窗函数的三谱线插值FFT算法的误差，并与常 规窗函数插值算法精度进行比较。
[0196] 下面通过一个具体实施例来进行详细说明。选定一个电网信号的谐波，谐波的波 形参见图1所示。
[0197] 本发明实施例提供一种基于互乘法窗函数的三谱线插值FFT谐波分析方法，具体 包括以下步骤：
[0198] 步骤1、互乘法窗函数的构造：
[0199] 以Hanning窗、Hamming窗、Blackman窗为例，给出互乘法窗函数的构造模式，构 造的6种乘法窗函数性能如表2所示。考虑到计算量的问题，乘法窗函数的阶次不适宜太 高，因此限定阶次c = 3。每种窗函数的子阶次可能取0、1、2三个值。为了书写的方面，下 面简化表达形式：Hanning简写为Hn，Hamming简写为Hm，Blackman简写为Bm。
[0200] 表1、基于常规函数构造互乘法窗函数的性能
[0202] 步骤2、信号的预处理：
[0203]将互感器采集到的离散电网信号x(n)，信号模型参见表2所示。为了验证所提算 法的精度，进行10次谐波仿真分析。信号模型为：
[0205] 其中：基波频率fQ为50. 5Hz，采样频率匕为5120Hz，采样点数N为1024。
[0206] 表2、电网信号的谐波参数
[0208] 加窗FFT变换后得到信号频谱：
[0209] 步骤3、确定三根谱线：
[0210] 在45~55Hz中寻求模最大三根谱线，分别为yp y2, y3。
[0211] 那么
[0212] 步骤4、计算三谱线插值算法的修正公式：
[0213] 根据【具体实施方式】中步骤S4的推导，表2中的4种窗函数的修正公式a = g4(f3)，u(a)分别如下：
[0214] (1)互乘210窗函数：
[0215] a = 1. 11458715 0 -0? 0908720 0 3+0. 01413650 0 5
[0216] u ( a ) = 1. 80783068+0. 41014621 a 2+〇. 05139583 a 4
[0217] 互乘210窗函数幅频特性参见图3所示。
[0218] (2)互乘201窗函数：
[0219] a = 1. 28076537 0 -0? 09843530 0 3+0. 01499840 0 5
[0220] u ( a ) = 1. 95369960+0. 38423865 a 2+〇. 04102745 a 4
[0221] 互乘201窗函数幅频特性参见图4所示。
[0222] (3)互乘120窗函数：
[0223] a = 1. 08516300 0 -0? 08842400 0 3+0. 0140181 0 5
[0224] u ( a ) = 1. 78667431+0. 41627357 a 2+〇. 05339022 a 4
[0225] 互乘120窗函数的幅频特性参见图5所示。
[0226] (4)互乘102窗函数：
[0227] a = 1. 42354957 0 -0? 10399918 0 3+0. 01605665 0 5
[0228] u ( a ) = 2. 07275558+0. 36590106 a 2+〇. 03458435 a 4
[0229] 互乘102窗函数的幅频特性参见图6所示。
[0230](5)互乘111窗函数：
[0231] a = 1. 2530640 0 -0. 09572246 0 3+0. 01451709 0 5
[0232]u(a) =1.93442536+0. 38883257a 2+〇.04234880a4
[0233] 互乘111窗函数的幅频特性参见图7所示。
[0234] (6)互乘021窗函数：
[0235] a = 1. 22447178 0 -0? 09322158 0 3+0. 01445149 0 5
[0236] u(a) = 1. 91516084+0. 39387119a 2+〇. 04376915a 4
[0237] 互乘021窗函数的幅频特性参见图8所示。
[0238] 步骤5、计算基波参数：
[0239]f0=k?Af=(k2+a)Af
[0240] Af N-1 (p2CI+p22 a 2+... +p2d a d)
[0242] 步骤6、确定谐波参数：确定基波频率&后，在范围（kfV5, kL+5)内，重复步骤 3~5,直到所有谐波参数计算完毕。
[0243] 步骤7、进行误差分析：分析基于互乘法窗函数的三谱线插值FFT算法的误差，并 与常规窗函数插值算法精度进行比较。误差分析比较结果参见表3~表5所示，其中，D Ai 表示基波和各次谐波幅值测量值的相对误差；Df(l表示基波频率测量值的相对误差；?表 示基波和各次谐波初始相位测量值的相对误差，均用百分比表示。
[0244] 表3、互乘法窗幅值、频率测量相对误差
[0246] 表4、互乘法窗相位测量相对误差
[0249] 而单独采用Hanning、Hamming和Blackman窗函数幅值、频率和相位测量误差如表 5所示。
[0250] 表5、Hn、Hm、Bm窗的幅值、频率和相位相对误差
[0252] 由此可见，本发明实施例中的互乘法窗函数在基本频率上，要比基本窗函数高 1~3个数量级，在幅值计算上高2~3个数量级，在相位计算上优势更为突出，高3~6个 数量级。由以上的仿真结果可以看出，计算结果普遍好于普通窗函数谱线插值算法。
[0253] 综上所述，本发明实施例中的基于互乘法窗函数谐波分析方法，首次 实现构造互 乘法窗函数，将互乘法窗函数应用到三谱线插值FFT算法中，计算谐波参数，能够获得比常 规窗函数更高的精度。
[0254] 本发明实施例还提供一种基于互乘法窗函数的三谱线插值FFT谐波分析系统，包 括互乘法窗函数构造单元、信号预处理单元、谱线确定单元、修正公式计算单元、基波参数 计算单元、谐波参数确定单元、误差分析单元，其中：
[0255] 互乘法函数构造单元，用于构造互乘法窗函数：
[0256] 乘法窗函数的通用公式是由多个窗函数乘积产生的，乘法窗函数w(n)的通用公 式为：
[0258] 其中，n为采样点的序数，n为自然数；m为参加构造乘法窗的基本窗函数的个数， m为正整数为第i个基本窗函数表达式，i = 1，2…m七为参加第i个基本窗函数 的个数，称为K (n)的子阶次
，c为互乘法窗函数的总阶次；根据Cl、c 2、(：3的 不同取值，形成不同的组合（ClcvCni)，即构造出不同形式的乘法窗函数；
[0259] 当力⑷表达式相同时，即所采用的基本窗相同时，w(n)为自乘法窗函数，当wjn) 表达式不同时，即所采用的基本窗不同时，w(n)为互乘法窗函数；
[0260] 下面以三个常用窗函数为例，表达式分别*Wl(n)，w2(n)，w 3(n)。总阶次c = Cl+C2+C3。由于C越大，意味着参与乘法窗的基本窗函数越多，构成的乘法窗函数项数越多， 会增加计算的复杂度。一般给C 一个最大值，例如，令：c = 3。
[0261] 根据Cl、c2、c 3的不同取值，形成不同的组合（c &???〇,构造出以下9种乘法窗函 数：
[0262] 当Ci= 3、c 2= 0、c 3= 0时，构造出互乘300窗函数；
[0263] 当Ci= 2、c 2= 1、c 3= 0时，构造出互乘210窗函数；
[0264] 当Ci= 2、c 2= 0、c 3= 1时，构造出互乘201窗函数；
[0265] 当Ci= 1、c 2= 2、c 3= 0时，构造出互乘120窗函数；
[0266] 当Ci= 1、c 2= 0、c 3= 2时，构造出互乘102窗函数；
[0267] 当Ci= 0、c 2= 3、c 3= 0时，构造出互乘030窗函数；
[0268] 当(^= l、c2= l、c3= 1时，构造出互乘111窗函数；
[0269] 当Ci= 0、c 2= 2、c 3= 1时，构造出互乘021窗函数；
[0270] 当Ci= 0、c 2= 0、c 3= 3时，构造出互乘003窗函数；
[0271] 从上述9种乘法窗函数可以看出：互乘300窗函数、互乘030窗函数、互乘003窗 函数，这三种窗函数实际上属于自乘法窗函数（本发明不予讨论），其余6种为构造出的互 乘法窗函数。
[0272] 信号预处理单元，用于对信号进行预处理：
[0273] 互感器采集电网信号，电网信号包括电流信号、电压信号；将互感器采集到的电网 信号x(n)，传输到上位机；对电网信号x(n)进行加互乘法窗函数w(n)进行截断，得到加窗 信号x w(n):
[0274] xw(n) = x(n)w(n) (2)
[0275] 对公式⑵的加窗信号进行FFT变换后，得到加窗FFT频谱：
[0277] 其中，W( ?)为窗函数的频谱，k为正整数，X(k)表示第k次谐波的频谱，Ak为第 k次谐波的幅值，j表示虚数单位，e是自然对数的底数，％t为第k次谐波的初始相位，第一 次谐波为基波，f s为采样频率，f 〇为基波频率，A f为离散频率间隔，且A f = f s/N ;
[0278] 令：kQ= fQ/Af，kQ为真实频谱的谱线位置，由于非同步采样或非整周期截断，kQ 一般为小数，不为整数。
[0279] 忽略负频率点处旁瓣的影响，公式（3)变为：
[0281] 谱线确定单元，用于确定三根谱线：
[0282] 在S2得到的加窗FFT频谱峰值附近区域，匕处频率点较大的三根谱线分别为：第 根，k 为正整数，k k 3的关系为：k 2= k i+1，k3= k 2+1，这三根谱线对应 的幅值分别为y2、
[0283] 记变量 a = k_k2,则-0? 5 彡 a 彡 〇? 5 ;
[0284] 另记变量，
[0285] 修正公式计算单元，用于计算三谱线插值算法的修正公式：
[0286] 根据公式⑷和（5)得到：
[0288] 当N值比较大，例如：N> 1000时，公式（6)可以化简为一个函数0 = g( a)，其 反函数为a =g4(f3)。由于所采用的余弦窗系数均为实系数，其频率响应是偶对称的，因 而g( ?)和g<( ?)均为奇函数。
[0289] 采用多项式逼近方法计算奇函数0 =g4(a)，表达式为：
[0290] a~pnX 0+p13X 03+…+plp0p(7)
[0291] pn，p13;…p lp为多项式逼近的奇次项系数，p是奇数；
[0292]根据公式（4)，求得电网信号第i次谐波的幅值Ai:
[0294] 其中，i为正整数，yi为加窗FFT后第i次谐波的幅值；
[0295] 以基波为例，考虑到y2是离真实谱线点最近谱线，给予较大权重，可以得到：
[0297] N> 1000时，窗函数系数为实系数，公式（9)表示为：
[0298] Aj= N_1 (y3+2y2+y1)u(a )
[0299] 其中，u(a)为修正公式，且为偶函数，逼近多项式不含奇次项；
[0300] 三谱线修正逼近多项式如下：
[0301] u ( a ) = (p2(l+p22 a 2+…+p2d a d) (10)
[0302] 公式（10)中，p2Q，p22~p2$多项式逼近的偶次项系数，d为拟合的最高阶次，且d 为偶数；
[0303] 基波参数计算单元，用于计算基波参数：
[0304]考虑到y2是离真实谱线点最近的两根谱线，给予较大权重，可以计算基波频率f。、 基波幅值A::
[0305] f0= k? Af = (k2+a)Af (11)
[0306] A! = N-1 (p2CI+p22 a 2+... +p2d a d) (12)
[0307]根据公式（4)，计算基波的相位：
[0309] 仿照基波参数的求取，根据公式（6)、（7)、（9)、（11)、（12)、（13)进行各次谐波参 数的分析；
[0310] 谐波参数确定单元，用于确定谐波参数：确定基波频率&后，在范围 (kfd-5, kf^+5)内，所述谱线确定单元、修正公式计算单元、基波参数计算单元重复进行计 算，直到所有谐波参数计算完毕；
[0311] 误差分析单元，用于进行误差分析：分析基于互乘法窗函数的三谱线插值FFT算 法的误差，并与常规窗函数插值算法精度进行比较。
[0312] 本领域的技术人员可以对本发明实施例进行各种修改和变型，倘若这些修改和变 型在本发明权利要求及其等同技术的范围之内，则这些修改和变型也在本发明的保护范围 之内。
[0313] 说明书中未详细描述的内容为本领域技术人员公知的现有技术。
【主权项】
1. 一种基于互乘法窗函数的三谱线插值FFT谐波分析方法，其特征在于，包括以下步 骤： 51、 构造互乘法窗函数： 乘法窗函数的通用公式是由多个窗函数乘积产生的，乘法窗函数w(n)的通用公式为：(1) 其中，η为采样点的序数，η为自然数;m为参加构造乘法窗的基本窗函数的个数，m为 正整数％ (η)为第i个基本窗函数表达式，i = 1，2···πι ;(^为参加第i个基本窗函数的个 数，称为Wi (η)的子阶次，c为互乘法窗函数的总阶次；根据Cl、c 2、C3的不同 取值，形成不同的组合（cl(V·· cm)，即构造出不同形式的乘法窗函数； 当Wi(Ii)表达式相同时，即所采用的基本窗相同时，w(n)为自乘法窗函数，当Wi(Ii)表 达式不同时，即所采用的基本窗不同时，w(n)为互乘法窗函数； 三个常用窗函数的表达式分别为W1(Ii)，W2 (n)，W3 (η)，总阶次C = CfCfC3，给C 一个最 大值，令：c = 3 ; 根据Cl、c2、C3的不同取值，形成不同的组合（c而…cm)，构造出以下9种乘法窗函数： 当C1= 3、C2= 0、c3= O时，构造出互乘300窗函数； 当C1= 2、C2= Uc3= 0时，构造出互乘210窗函数； 当C1= 2、C2= 0、c3= 1时，构造出互乘201窗函数； 当C1= Uc2= 2、C3= 0时，构造出互乘120窗函数； 当C1= Uc2= (Kc3= 2时，构造出互乘102窗函数； 当C1= (Kc2= 3、C3= 0时，构造出互乘030窗函数； 当C1= Uc2= Uc3= 1时，构造出互乘111窗函数； 当C1= (Kc2= 2、C3= 1时，构造出互乘021窗函数； 当C1= 0、c2= (Kc3= 3时，构造出互乘003窗函数； 其中，互乘300窗函数、互乘030窗函数、互乘003窗函数，这三种窗函数属于自乘法窗 函数，其余6种为构造出的互乘法窗函数； 52、 信号预处理： 互感器采集电网信号，将互感器采集到的电网信号x( n)，传输到上位机；对电网信号 X(H)进行加互乘法窗函数w (η)进行截断，得到加窗信号Xw (η): xw (η) = X (n) w (η) (2) 对公式（2)的加窗信号进行FFT变换后，得到加窗FFT频谱：其中，W( ·)为窗函数的频谱，k为正整数，X(k)表示第k次谐波的频谱，Ak为第k次 谐波的幅值，j表示虚数单位，e是自然对数的底数，妁为第k次谐波的初始相位，第一次谐 波为基波，fs为采样频率，f 0为基波频率，Δ f为离散频率间隔，且Λ f = f s/N ; 令：1?= fd/Δ f，1?为真实频谱的谱线位置， 忽略负频率点处旁瓣的影响，公式（3)变为：(4) 53、 确定三根谱线： 在S2得到的加窗FFT频谱峰值附近区域，1?处频率点较大的三根谱线分别为：第k i、 k2、k3根，k η k2、k3均为正整数，k广k 3的关系为：k 2= k i+1，k3= k 2+1，这三根谱线对应的 幅值分别为y2、y3; 记变量 α = k_k2，则-〇· 5 < α < 〇· 5 ; 另记变量（5) 54、 计算三谱线插值算法的修正公式： 根据公式（4)和（5)得到：(6) 采用多项式逼近方法计算奇函数β = Ρ(α)，表达式为： α ~ ρηΧ β+ρ13Χ β 3+…+ρ1ρβρ (7) Pn，P13;…P 1Ρ为多项式逼近的奇次项系数，P是奇数； 根据公式（4)，求得电网信号第i次谐波的幅值Ai:C8) 其中，i为正整数，Yi为加窗FFT后第i次谐波的幅值； 考虑到y2是离真实谱线点最近谱线，得到：(9) N > 1000时，窗函数系数为实系数，公式（9)表示为： A1= N (y3+2y2+y1)u( α ) 其中，u(a)为修正公式，且为偶函数，逼近多项式不含奇次项； 三谱线修正逼近多项式如下： u ( a ) = (ρ2(ι+ρ22 α 2+…+p2d a d) (10) 公式（10)中，P2(I，P22··· P2A多项式逼近的偶次项系数，d为拟合的最高阶次，且d为偶 数； 55、 计算基波参数： 计算基波频率fci、基波幅值A1: f0=k· Af = (k2+a)Af (11) A1 = N-1 (y^Sy^y!) (p2CI+p22 a 2+... +p2da d) (12) 根据公式（4)，计算基波的相位：(13) 仿照基波参数的求取，根据公式（6)、（7)、（9)、（11)、（12)、（13)进行各次谐波参数的 分析； 56、 确定谐波参数：确定基波频率&后，在范围内，重复步骤S3~S5,直 到所有谐波参数计算完毕； 57、 进行误差分析：分析基于互乘法窗函数的三谱线插值FFT算法的误差，并与常规窗 函数插值算法精度进行比较。2. 如权利要求1所述的基于互乘法窗函数的三谱线插值FFT谐波分析方法，其特征在 于：所述电网信号包括电流信号、电压信号。3. 如权利要求1所述的基于互乘法窗函数的三谱线插值FFT谐波分析方法，其特征在 于：所述真实频谱的谱线位置1?为小数。4. 一种基于互乘法窗函数的三谱线插值FFT谐波分析系统，其特征在于：包括互乘法 窗函数构造单元、信号预处理单元、谱线确定单元、修正公式计算单元、基波参数计算单元、 谐波参数确定单元、误差分析单元，其中： 所述互乘法函数构造单元，用于构造互乘法窗函数： 乘法窗函数的通用公式是由多个窗函数乘积产生的，乘法窗函数w(n)的通用公式为：(1) 其中，η为采样点的序数，η为自然数;m为参加构造乘法窗的基本窗函数的个数，m为 正整数％ (η)为第i个基本窗函数表达式，i = 1，2···πι ;(^为参加第i个基本窗函数的个 数，称为Wi(Il)的子阶次：c为互乘法窗函数的总阶次；根据Cl、c 2、C3的不同 取值，形成不同的组合（cl(V·· cm)，即构造出不同形式的乘法窗函数； 当Wi(Ii)表达式相同时，即所采用的基本窗相同时，w(n)为自乘法窗函数，当Wi(Ii)表 达式不同时，即所采用的基本窗不同时，w(n)为互乘法窗函数； 三个常用窗函数的表达式分别为W1(Ii)，W2 (n)，W3 (η)，总阶次C = CfCfC3，给C 一个最 大值，令：c = 3 ; 根据Cl、c2、C3的不同取值，形成不同的组合（c而…cm)，构造出以下9种乘法窗函数： 当C1= 3、C2= 0、c3= O时，构造出互乘300窗函数； 当C1= 2、C2= Uc3= 0时，构造出互乘210窗函数； 当C1= 2、C2= 0、c3= 1时，构造出互乘201窗函数； 当C1= Uc2= 2、C3= 0时，构造出互乘120窗函数； 当C1= Uc2= (Kc3= 2时，构造出互乘102窗函数； 当C1= (Kc2= 3、C3= 0时，构造出互乘030窗函数； 当C1= Uc2= Uc3= 1时，构造出互乘111窗函数； 当C1= (Kc2= 2、C3= 1时，构造出互乘021窗函数； 当C1= 0、c2= (Kc3= 3时，构造出互乘003窗函数； 其中，互乘300窗函数、互乘030窗函数、互乘003窗函数，这三种窗函数属于自乘法窗 函数，其余6种为构造出的互乘法窗函数； 所述信号预处理单元，用于对信号进行预处理： 互感器采集电网信号，将互感器采集到的电网信号x(n)，传输到上位机；对电网信号 X(H)进行加互乘法窗函数w (η)进行截断，得到加窗信号Xw (η): xw (η) = X (n) w (η) (2) 对公式（2)的加窗信号进行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), 所述谱线确定单元，用于确定三根谱线： 在S2得到的加窗FFT频谱峰值附近区域，1?处频率点较大的三根谱线分别为：第k i、 k2、k3根，k η k2、k3均为正整数，k广k 3的关系为：k 2= k i+1，k3= k 2+1，这三根谱线对应的 幅值分别为y2、y3; 记变量 α = k_k2，则-〇· 5 < α < 〇· 5 ; 另记变量（5); 所述修正公式计算单元，用于计算三谱线插值算法的修正公式： 根据公式（4)和（5)得到：(6) 采用多项式逼近方法计算奇函数β = Ρ(α)，表达式为： α ~ ρηΧ β+ρ13Χ β 3+…+ρ1ρβρ (7) Pn，P13;…P 1Ρ为多项式逼近的奇次项系数，P是奇数； 根据公式（4)，求得电网信号第i次谐波的幅值Ai:C8) 其中，i为正整数，Yi为加窗FFT后第i次谐波的幅值； 考虑到y2是离真实谱线点最近谱线，得到：(9) N > 1000时，窗函数系数为实系数，公式（9)表示为： A1= N (y3+2y2+y1)u( α ) 其中，u(a)为修正公式，且为偶函数，逼近多项式不含奇次项； 三谱线修正逼近多项式如下： u ( a ) = (ρ2(ι+ρ22 α 2+…+p2d a d) (10) 公式（10)中，P2(I，P22··· P2A多项式逼近的偶次项系数，d为拟合的最高阶次，且d为偶 数； 所述基波参数计算单元，用于计算基波参数： 计算基波频率fci、基波幅值A1: f0=k* Af= (k2+a)Af (11) A1 = N-1 (p2CI+p22 a 2+... +p2da d) (12) 根据公式（4)，计算基波的相位：(13) 仿照基波参数的求取，根据公式（6)、（7)、（9)、（11)、（12)、（13)进行各次谐波参数的 分析； 所述谐波参数确定单元，用于确定谐波参数：确定基波频率后，在范围 (kfd-5, kf^+5)内，所述谱线确定单元、修正公式计算单元、基波参数计算单元重复进行计 算，直到所有谐波参数计算完毕； 所述误差分析单元，用于进行误差分析：分析基于互乘法窗函数的三谱线插值FFT算 法的误差，并与常规窗函数插值算法精度进行比较。5. 如权利要求4所述的基于互乘法窗函数的三谱线插值FFT谐波分析系统，其特征在 于：所述电网信号包括电流信号、电压信号。6. 如权利要求4所述的基于互乘法窗函数的三谱线插值FFT谐波分析系统，其特征在 于：所述真实频谱的谱线位置1?为小数。
【专利摘要】本发明公开了一种基于互乘法窗函数的三谱线插值FFT谐波分析方法及系统，涉及谐波分析领域。该方法包括以下步骤：构造互乘法窗函数；信号预处理；确定三根谱线；计算三谱线插值算法的修正公式；计算基波参数；确定谐波参数；进行误差分析。本发明提出一种新的窗函数的构造方式，利用常规窗函数进行不同的乘法组合，首次实现构造互乘法窗函数，将互乘法窗函数应用到三谱线插值FFT算法中，进行谐波分析，计算谐波参数。多种常用的余弦窗函数计算实例表明，本发明构造出的互乘法窗函数相对于常规窗函数插值算法，有更高的准确度，能实现谐波的高精度测量。
【IPC分类】G01R23/16
【公开号】CN104897961
【申请号】CN201510335000
【发明人】张俊敏, 刘开培, 汪立, 王黎, 田微, 陈文娟
【申请人】中南民族大学
【公开日】2015年9月9日
【申请日】2015年6月17日
最新回复(0)