本发明公开了一种基于互乘法窗函数的三谱线插值FFT谐波分析方法及系统，涉及谐波分析领域。该方法包括以下步骤：构造互乘法窗函数；信号预处理；确定三根谱线；计算三谱线插值算法的修正公式；计算基波参数；确定谐波参数；进行误差分析。本发明提出一种新的窗函数的构造方式，利用常规窗函数进行不同的乘法组合，首次实现构造互乘法窗函数，将互乘法窗函数应用到三谱线插值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)的通用公式为：
[0012]
[0013] 其中，n为采样点的序数，n为自然数；m为参加构造乘法窗函数的基本窗函数的个数，m为正整数；wi(n)为第i个基本窗函数表达式，i＝1,2...m；ci为参加第i个基本窗函数的个数，称为wi(n)的子阶次； c为互乘法窗函数的总阶次；根据c1、c2...cm的不同取值，形成不同的组合(c1c2...cm)，即构造出不同形式的乘法窗函数；
[0014] 当wi(n)表达式相同时，即所采用的基本窗函数相同时，w(n)为自乘法窗函数，当wi(n)表达式不同时，即所采用的基本窗函数不同时， w(n)为互乘法窗函数；
[0015] 三个常用窗函数的表达式分别为w1(n)，w2(n)，w3(n)，总阶次 c＝c1+c2+c3，给c一个最大值，令：c＝3；
[0016] 根据c1、c2、c3的不同取值，形成不同的组合(c1c2c3)，构造出以下9种乘法窗函数：
[0017] 当c1＝3、c2＝0、c3＝0时，构造出互乘300窗函数；
[0018] 当c1＝2、c2＝1、c3＝0时，构造出互乘210窗函数；
[0019] 当c1＝2、c2＝0、c3＝1时，构造出互乘201窗函数；
[0020] 当c1＝1、c2＝2、c3＝0时，构造出互乘120窗函数；
[0021] 当c1＝1、c2＝0、c3＝2时，构造出互乘102窗函数；
[0022] 当c1＝0、c2＝3、c3＝0时，构造出互乘030窗函数；
[0023] 当c1＝1、c2＝1、c3＝1时，构造出互乘111窗函数；
[0024] 当c1＝0、c2＝2、c3＝1时，构造出互乘021窗函数；
[0025] 当c1＝0、c2＝0、c3＝3时，构造出互乘003窗函数；
[0026] 其中，互乘300窗函数、互乘030窗函数、互乘003窗函数，这三种窗函数属于自乘法窗函数，其余6种为构造出的互乘法窗函数；
[0027] S2、信号预处理：
[0028] 互感器采集电网信号，将互感器采集到的电网信号x(n)，传输到上位机；对电网信号x(n)进行加互乘法窗函数w(n)进行截断，得到加窗信号xw(n)：
[0029] xw(n)＝x(n)w(n) (2)
[0030] 对公式(2)的加窗信号进行FFT变换后，得到加窗FFT频谱：
[0031]
[0032] 其中，W(·)为窗函数的频谱，k为正整数，X(k)表示第k次谐波的频谱，Ak为第k次谐波的幅值，j表示虚数单位，e是自然对数的底数， 为第k次谐波的初始相位，第一次谐波为基波，fs为采样频率，f0为基波频率，Δf为离散频率间隔，且Δf＝fs/N；
[0033] 令：k0＝f0/Δf，k0为真实频谱的谱线位置，
[0034] 忽略负频率点处旁瓣的影响，公式(3)变为：
[0035]
[0036] S3、确定三根谱线：
[0037] 在S2得到的加窗FFT频谱峰值附近区域，k0处频率点较大的三根谱线分别为：第k1、k2、k3根，k1、k2、k3均为正整数，k1～k3的关系为：k2＝k1+1，k3＝k2+1，这三根谱线对应的幅值分别为y1、y2、y3；
[0038] 记变量α＝k-k2，则-0.5≤α≤0.5；
[0039] 另记变量
[0040] S4、计算三谱线插值算法的修正公式：
[0041] 根据公式(4)和(5)得到：
[0042]
[0043] 采用多项式逼近方法计算奇函数α＝g-1(β)，表达式为：
[0044] α≈p11×β+p13×β3+...+p1pβp (7)
[0045] p11，p13；…p1p为多项式逼近的奇次项系数，p是奇数；
[0046] 根据公式(4)，求得电网信号第i次谐波的幅值Ai：
[0047]
[0048] 其中，i为正整数，yi为加窗FFT后第i次谐波的幅值；
[0049] 考虑到y2是离真实谱线点最近谱线，得到：
[0050]
[0051] N＞1000时，窗函数系数为实系数，公式(9)表示为：
[0052] A1＝N-1(y3+2y2+y1)u(α)
[0053] 其中，u(α)为修正公式，且为偶函数，逼近多项式不含奇次项；
[0054] 三谱线修正逼近多项式如下：
[0055] u(α)＝(p20+p22α2+…+p2dαd) (10)
[0056] 公式(10)中，p20,p22…p2d为多项式逼近的偶次项系数，d为拟合的最高阶次，且d为偶数；
[0057] S5、计算基波参数：
[0058] 计算基波频率f0、基波幅值A1：
[0059] f0＝k·Δf＝(k2+α)Δf (11)
[0060] A1＝N-1(y3+2y2+y1)(p20+p22α2+…+p2dαd) (12)
[0061] 根据公式(4)，计算基波的相位：
[0062]
[0063] 仿照基波参数的求取，根据公式(6)、(7)、(9)、(11)、(12)、 (13)进行各次谐波参数的分析；
[0064] S6、确定谐波参数：确定基波频率f0后，在范围(kf0-5,kf0+5) 内，重复步骤S3～S5，直到所有谐波参数计算完毕；
[0065] S7、进行误差分析：分析基于互乘法窗函数的三谱线插值FFT 算法的误差，并与常规窗函数插值算法精度进行比较。
[0066] 在上述技术方案的基础上，所述电网信号包括电流信号、电压信号。
[0067] 在上述技术方案的基础上，所述真实频谱的谱线位置k0为小数。
[0068] 本发明还提供一种基于互乘法窗函数的三谱线插值FFT谐波分析系统，包括互乘法窗函数构造单元、信号预处理单元、谱线确定单元、修正公式计算单元、基波参数计算单元、谐波参数确定单元、误差分析单元，其中：
[0069] 所述互乘法窗函数构造单元，用于构造互乘法窗函数：
[0070] 乘法窗函数的通用公式是由多个窗函数乘积产生的，乘法窗函数 w(n)的通用公式为：
[0071]
[0072] 其中，n为采样点的序数，n为自然数；m为参加构造乘法窗函数的基本窗函数的个数，m为正整数；wi(n)为第i个基本窗函数表达式，i＝1,2...m；ci为参加第i个基本窗函数的个数，称为wi(n)的子阶次； c为互乘法窗函数的总阶次；根据c1、c2...cm的不同取值，形成不同的组合(c1c2...cm)，即构造出不同形式的乘法窗函数；
[0073] 当wi(n)表达式相同时，即所采用的基本窗函数相同时，w(n)为自乘法窗函数，当wi(n)表达式不同时，即所采用的基本窗函数不同时， w(n)为互乘法窗函数；
[0074] 三个常用窗函数的表达式分别为w1(n)，w2(n)，w3(n)，总阶次 c＝c1+c2+c3，给c一个最大值，令：c＝3；
[0075] 根据c1、c2、c3的不同取值，形成不同的组合(c1c2c3)，构造出以下9种乘法窗函数：
[0076] 当c1＝3、c2＝0、c3＝0时，构造出互乘300窗函数；
[0077] 当c1＝2、c2＝1、c3＝0时，构造出互乘210窗函数；
[0078] 当c1＝2、c2＝0、c3＝1时，构造出互乘201窗函数；
[0079] 当c1＝1、c2＝2、c3＝0时，构造出互乘120窗函数；
[0080] 当c1＝1、c2＝0、c3＝2时，构造出互乘102窗函数；
[0081] 当c1＝0、c2＝3、c3＝0时，构造出互乘030窗函数；
[0082] 当c1＝1、c2＝1、c3＝1时，构造出互乘111窗函数；
[0083] 当c1＝0、c2＝2、c3＝1时，构造出互乘021窗函数；
[0084] 当c1＝0、c2＝0、c3＝3时，构造出互乘003窗函数；
[0085] 其中，互乘300窗函数、互乘030窗函数、互乘003窗函数，这三种窗函数属于自乘法窗函数，其余6种为构造出的互乘法窗函数；
[0086] 所述信号预处理单元，用于对信号进行预处理：
[0087] 互感器采集电网信号，将互感器采集到的电网信号x(n)，传输到上位机；对电网信号x(n)进行加互乘法窗函数w(n)进行截断，得到加窗信号xw(n)：
[0088] xw(n)＝x(n)w(n) (2)
[0089] 对公式(2)的加窗信号进行FFT变换后，得到加窗FFT频谱：
[0090]
[0091] 其中，W(·)为窗函数的频谱，k为正整数，X(k)表示第k次谐波的频谱，Ak为第k次谐波的幅值，j表示虚数单位，e是自然对数的底数， 为第k次谐波的初始相位，第一次谐波为基波，fs为采样频率，f0为基波频率，Δf为离散频率间隔，且Δf＝fs/N；
[0092] 令：k0＝f0/Δf，k0为真实频谱的谱线位置，
[0093] 忽略负频率点处旁瓣的影响，公式(3)变为：
[0094]
[0095] 所述谱线确定单元，用于确定三根谱线：
[0096] 在所述信号预处理单元得到的加窗FFT频谱峰值附近区域，k0处频率点较大的三根谱线分别为：第k1、k2、k3根，k1、k2、k3均为正整数，k1～k3的关系为：k2＝k1+1，k3＝k2+1，这三根谱线对应的幅值分别为y1、y2、y3；
[0097] 记变量α＝k-k2，则-0.5≤α≤0.5；
[0098] 另记变量
[0099] 所述修正公式计算单元，用于计算三谱线插值算法的修正公式：
[0100] 根据公式(4)和(5)得到：
[0101]
[0102] 采用多项式逼近方法计算奇函数α＝g-1(β)，表达式为：
[0103] α≈p11×β+p13×β3+…+p1pβp (7)
[0104] p11，p13；…p1p为多项式逼近的奇次项系数，p是奇数；
[0105] 根据公式(4)，求得电网信号第i次谐波的幅值Ai：
[0106]
[0107] 其中，i为正整数，yi为加窗FFT后第i次谐波的幅值；
[0108] 考虑到y2是离真实谱线点最近谱线，得到：
[0109]
[0110] N＞1000时，窗函数系数为实系数，公式(9)表示为：
[0111] A1＝N-1(y3+2y2+y1)u(α)
[0112] 其中，u(α)为修正公式，且为偶函数，逼近多项式不含奇次项；
[0113] 三谱线修正逼近多项式如下：
[0114] u(α)＝(p20+p22α2+…+p2dαd) (10)
[0115] 公式(10)中，p20,p22…p2d为多项式逼近的偶次项系数，d为拟合的最高阶次，且d为偶数；
[0116] 所述基波参数计算单元，用于计算基波参数：
[0117] 计算基波频率f0、基波幅值A1：
[0118] f0＝k·Δf＝(k2+α)Δf (11)
[0119] A1＝N-1(y3+2y2+y1)(p20+p22α2+…+p2dαd) (12)
[0120] 根据公式(4)，计算基波的相位：
[0121]
[0122] 仿照基波参数的求取，根据公式(6)、(7)、(9)、(11)、(12)、 (13)进行各次谐波参数的分析；
[0123] 所述谐波参数确定单元，用于确定谐波参数：确定基波频率f0后，在范围(kf0-5,kf0+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)的通用公式为：
[0141]
[0142] 其中，n为采样点的序数，n为自然数；m为参加构造乘法窗函数的基本窗函数的个数，m为正整数；wi(n)为第i个基本窗函数表达式，i＝1,2...m；ci为参加第i个基本窗函数的个数，称为wi(n)的子阶次； c为互乘法窗函数的总阶次。根据c1、c2...cm的不同取值，形成不同的组合(c1c2...cm)，即构造出不同形式的乘法窗函数。
[0143] 当wi(n)表达式相同时，即所采用的基本窗函数相同时，w(n)为自乘法窗函数，当wi(n)表达式不同时，即所采用的基本窗函数不同时， w(n)为互乘法窗函数。
[0144] 下面以三个常用窗函数为例，表达式分别为w1(n)，w2(n)，w3(n)。总阶次c＝c1+c2+c3。由于c越大，意味着参与乘法窗的基本窗函数越多，构成的乘法窗函数项数越多，会增加计算的复杂度。一般给c一个最大值，例如，令：c＝3。
[0145] 根据c1、c2、c3的不同取值，形成不同的组合(c1c2c3)，构造出以下9种乘法窗函数：
[0146] 当c1＝3、c2＝0、c3＝0时，构造出互乘300窗函数；
[0147] 当c1＝2、c2＝1、c3＝0时，构造出互乘210窗函数；
[0148] 当c1＝2、c2＝0、c3＝1时，构造出互乘201窗函数；
[0149] 当c1＝1、c2＝2、c3＝0时，构造出互乘120窗函数；
[0150] 当c1＝1、c2＝0、c3＝2时，构造出互乘102窗函数；
[0151] 当c1＝0、c2＝3、c3＝0时，构造出互乘030窗函数；
[0152] 当c1＝1、c2＝1、c3＝1时，构造出互乘111窗函数；
[0153] 当c1＝0、c2＝2、c3＝1时，构造出互乘021窗函数；
[0154] 当c1＝0、c2＝0、c3＝3时，构造出互乘003窗函数；
[0155] 从上述9种乘法窗函数可以看出：互乘300窗函数、互乘030窗函数、互乘003窗函数，这三种窗函数实际上属于自乘法窗函数(本发明不予讨论)，其余6种为构造出的互乘法窗函数。
[0156] S2、信号预处理：
[0157] 互感器采集电网信号，电网信号包括电流信号、电压信号，将互感器采集到的电网信号x(n)，传输到上位机；对电网信号x(n)进行加互乘法窗函数w(n)进行截断，得到加窗信号xw(n)：
[0158] xw(n)＝x(n)w(n) (2)
[0159] 对公式(2)的加窗信号进行FFT变换后，得到加窗FFT频谱：
[0160]
[0161] 其中，W(·)为窗函数的频谱，k为正整数，X(k)表示第k次谐波的频谱，Ak为第k次谐波的幅值，j表示虚数单位，e是自然对数的底数， 为第k次谐波的初始相位，第一次谐波为基波，fs为采样频率，f0为基波频率，Δf为离散频率间隔，且Δf＝fs/N。
[0162] 令：k0＝f0/Δf，k0为真实频谱的谱线位置，由于非同步采样或非整周期截断，k0一般为小数，不为整数。
[0163] 若忽略负频率点处旁瓣的影响，公式(3)变为：
[0164]
[0165] S3、确定三根谱线：
[0166] 在S2得到的加窗FFT频谱峰值附近区域，k0处频率点较大的三根谱线分别为：第k1、k2、k3根，k1、k2、k3均为正整数，k1～k3的关系为：k2＝k1+1，k3＝k2+1，这三根谱线对应的幅值分别为y1、y2、y3。
[0167] 为方便计算，记变量α＝k-k2，则-0.5≤α≤0.5；
[0168] 另记变量
[0169] S4、计算三谱线插值算法的修正公式：
[0170] 根据公式(4)和(5)可以得到：
[0171]
[0172] 当N值比较大，例如：N＞1000时，公式(6)可以化简为一个函数β＝g(α)，其反函数为α＝g-1(β)。由于所采用的余弦窗系数均为实系数，其频率响应是偶对称的，因而g(·)和g-1(·)均为奇函数。
[0173] 采用多项式逼近方法计算奇函数α＝g-1(β)，表达式为：
[0174] α≈p11×β+p13×β3+…+p1pβp (7)
[0175] p11，p13；…p1p为多项式逼近的奇次项系数，p是奇数。
[0176] 根据公式(4)，求得电网信号第i次谐波的幅值Ai：
[0177]
[0178] 其中，i为正整数，yi为加窗FFT后第i次谐波的幅值；
[0179] 以基波为例，考虑到y2是离真实谱线点最近谱线，给予较大权重，可以得到：
[0180]
[0181] 类似公式(7)的逼近方法，当N比较大，例如N＞1000时，窗函数系数为实系数，公式(9)可表示为：
[0182] A1＝N-1(y3+2y2+y1)u(α)
[0183] 其中，u(α)为修正公式，且为偶函数，逼近多项式不含奇次项。
[0184] 三谱线修正逼近多项式如下：
[0185] u(α)＝(p20+p22α2+…+p2dαd) (10)
[0186] 公式(10)中，p20,p22…p2d为多项式逼近的偶次项系数，d为拟合的最高阶次，且d为偶数。
[0187] S5、计算基波参数：
[0188] 考虑到y2是离真实谱线点最近的两根谱线，给予较大权重，可以计算基波频率f0、基波幅值A1：
[0189] f0＝k·Δf＝(k2+α)Δf (11)
[0190] A1＝N-1(y3+2y2+y1)(p20+p22α2+…+p2dαd) (12)
[0191] 根据公式(4)，还可以计算基波的相位：
[0192]
[0193] 仿照基波参数的求取，根据公式(6)、(7)、(9)、(11)、(12)、 (13)即可进行各次谐波参数的分析。
[0194] S6、确定谐波参数：确定基波频率f0后，在范围(kf0-5,kf0+5) 内，重复步骤S3～S5，直到所有谐波参数计算完毕；
[0195] S7、进行误差分析：分析基于互乘法窗函数的三谱线插值FFT 算法的误差，并与常规窗函数插值算法精度进行比较。
[0196] 下面通过一个具体实施例来进行详细说明。选定一个电网信号的谐波，谐波的波形参见图2所示。
[0197] 本发明实施例提供一种基于互乘法窗函数的三谱线插值FFT谐波分析方法，具体包括以下步骤：
[0198] 步骤1、互乘法窗函数的构造：
[0199] 以Hanning窗、Hamming窗、Blackman窗为例，给出互乘法窗函数的构造模式，构造的6种乘法窗函数性能如表1所示。考虑到计算量的问题，乘法窗函数的阶次不适宜太高，因此限定阶次c＝3。每种窗函数的子阶次可能取0、1、2三个值。为了书写的方面，下面简化表达形式：Hanning简写为Hn，Hamming简写为Hm，Blackman 简写为Bm。
[0200] 表1、基于常规函数构造互乘法窗函数的性能
[0201]Hn(c1) Hm(c2) Bm(c3) 主瓣宽度(×π) 副瓣衰减(dB)
2 1 0 0.0054 -62.8
2 0 1 0.0059 -84.5
1 2 0 0.0051 -64.3
1 0 2 0.0061 -100.6
1 1 1 0.0056 -84.9
0 2 1 0.0056 -86.1
[0202] 步骤2、信号的预处理：
[0203] 将互感器采集到的离散电网信号x(n)，信号模型参见表2所示。为了验证所提算法的精度，进行10次谐波仿真分析。信号模型为：
[0204]
[0205] 其中：基波频率f0为50.5Hz，采样频率fs为5120Hz，采样点数 N为1024。
[0206] 表2、电网信号的谐波参数
[0207]
[0208] 加窗FFT变换后得到信号频谱：
[0209] 步骤3、确定三根谱线：
[0210] 在45～55Hz中寻求模最大三根谱线，分别为y1，y2，y3。
[0211] 那么
[0212] 步骤4、计算三谱线插值算法的修正公式：
[0213] 根据具体实施方式中步骤S4的推导，表1中的6种窗函数的修正公式α＝g-1(β)，u(α)分别如下：
[0214] (1)互乘210窗函数：
[0215] α＝1.11458715β-0.0908720β3+0.01413650β5
[0216] u(α)＝1.80783068+0.41014621α2+0.05139583α4
[0217] 互乘210窗函数幅频特性参见图3所示。
[0218] (2)互乘201窗函数：
[0219] α＝1.28076537β-0.09843530β3+0.01499840β5
[0220] u(α)＝1.95369960+0.38423865α2+0.04102745α4
[0221] 互乘201窗函数幅频特性参见图4所示。
[0222] (3)互乘120窗函数：
[0223] α＝1.08516300β-0.08842400β3+0.0140181β5
[0224] u(α)＝1.78667431+0.41627357α2+0.05339022α4
[0225] 互乘120窗函数的幅频特性参见图5所示。
[0226] (4)互乘102窗函数：
[0227] α＝1.42354957β-0.10399918β3+0.01605665β5
[0228] u(α)＝2.07275558+0.36590106α2+0.03458435α4
[0229] 互乘102窗函数的幅频特性参见图6所示。
[0230] (5)互乘111窗函数：
[0231] α＝1.2530640β-0.09572246β3+0.01451709β5
[0232] u(α)＝1.93442536+0.38883257α2+0.04234880α4
[0233] 互乘111窗函数的幅频特性参见图7所示。
[0234] (6)互乘021窗函数：
[0235] α＝1.22447178β-0.09322158β3+0.01445149β5
[0236] u(α)＝1.91516084+0.39387119α2+0.04376915α4
[0237] 互乘021窗函数的幅频特性参见图8所示。
[0238] 步骤5、计算基波参数：
[0239] f0＝k·Δf＝(k2+α)Δf
[0240] A1＝N-1(y3+2y2+y1)(p20+p22α2+…+p2dαd)
[0241]
[0242] 步骤6、确定谐波参数：确定基波频率f0后，在范围 (kf0-5,kf0+5)内，重复步骤3～5，直到所有谐波参数计算完毕。
[0243] 步骤7、进行误差分析：分析基于互乘法窗函数的三谱线插值FFT 算法的误差，并与常规窗函数插值算法精度进行比较。误差分析比较结果参见表3～表5所示，其中，DAi表示基波和各次谐波幅值测量值的相对误差；Df0表示基波频率测量值的相对误差； 表示基波和各次谐波初始相位测量值的相对误差，均用百分比表示。
[0244] 表3、互乘法窗幅值、频率测量相对误差
[0245]
[0246]
[0247] 表4、互乘法窗相位测量相对误差
[0248]
[0249] 而单独采用Hanning、Hamming和Blackman窗函数幅值、频率和相位测量误差如表5所示。
[0250] 表5、Hn、Hm、Bm窗的幅值、频率和相位相对误差
[0251]
[0252] 由此可见，本发明实施例中的互乘法窗函数在基本频率上，要比基本窗函数高1～3个数量级，在幅值计算上高2～3个数量级，在相位计算上优势更为突出，高3～6个数量级。
由以上的仿真结果可以看出，计算结果普遍好于普通窗函数谱线插值算法。
[0253] 综上所述，本发明实施例中的基于互乘法窗函数谐波分析方法，首次实现构造互乘法窗函数，将互乘法窗函数应用到三谱线插值FFT 算法中，计算谐波参数，能够获得比常规窗函数更高的精度。
[0254] 本发明实施例还提供一种基于互乘法窗函数的三谱线插值FFT 谐波分析系统，包括互乘法窗函数构造单元、信号预处理单元、谱线确定单元、修正公式计算单元、基波参数计算单元、谐波参数确定单元、误差分析单元，其中：
[0255] 互乘法窗函数构造单元，用于构造互乘法窗函数：
[0256] 乘法窗函数的通用公式是由多个窗函数乘积产生的，乘法窗函数 w(n)的通用公式为：
[0257]
[0258] 其中，n为采样点的序数，n为自然数；m为参加构造乘法窗函数的基本窗函数的个数，m为正整数；wi(n)为第i个基本窗函数表达式，i＝1,2...m；ci为参加第i个基本窗函数的个数，称为wi(n)的子阶次； c为互乘法窗函数的总阶次；根据c1、c2...cm的不同取值，形成不同的组合(c1c2…cm)，即构造出不同形式的乘法窗函数；
[0259] 当wi(n)表达式相同时，即所采用的基本窗函数相同时，w(n)为自乘法窗函数，当wi(n)表达式不同时，即所采用的基本窗函数不同时， w(n)为互乘法窗函数；
[0260] 下面以三个常用窗函数为例，表达式分别为w1(n)，w2(n)，w3(n)。总阶次c＝c1+c2+c3。由于c越大，意味着参与乘法窗的基本窗函数越多，构成的乘法窗函数项数越多，会增加计算的复杂度。一般给c一个最大值，例如，令：c＝3。
[0261] 根据c1、c2、c3的不同取值，形成不同的组合(c1c2c3)，构造出以下9种乘法窗函数：
[0262] 当c1＝3、c2＝0、c3＝0时，构造出互乘300窗函数；
[0263] 当c1＝2、c2＝1、c3＝0时，构造出互乘210窗函数；
[0264] 当c1＝2、c2＝0、c3＝1时，构造出互乘201窗函数；
[0265] 当c1＝1、c2＝2、c3＝0时，构造出互乘120窗函数；
[0266] 当c1＝1、c2＝0、c3＝2时，构造出互乘102窗函数；
[0267] 当c1＝0、c2＝3、c3＝0时，构造出互乘030窗函数；
[0268] 当c1＝1、c2＝1、c3＝1时，构造出互乘111窗函数；
[0269] 当c1＝0、c2＝2、c3＝1时，构造出互乘021窗函数；
[0270] 当c1＝0、c2＝0、c3＝3时，构造出互乘003窗函数；
[0271] 从上述9种乘法窗函数可以看出：互乘300窗函数、互乘030窗函数、互乘003窗函数，这三种窗函数实际上属于自乘法窗函数(本发明不予讨论)，其余6种为构造出的互乘法窗函数。
[0272] 信号预处理单元，用于对信号进行预处理：
[0273] 互感器采集电网信号，电网信号包括电流信号、电压信号；将互感器采集到的电网信号x(n)，传输到上位机；对电网信号x(n)进行加互乘法窗函数w(n)进行截断，得到加窗信号xw(n)：
[0274] xw(n)＝x(n)w(n) (2)
[0275] 对公式(2)的加窗信号进行FFT变换后，得到加窗FFT频谱：
[0276]
[0277] 其中，W(·)为窗函数的频谱，k为正整数，X(k)表示第k次谐波的频谱，Ak为第k次谐波的幅值，j表示虚数单位，e是自然对数的底数， 为第k次谐波的初始相位，第一次谐波为基波，fs为采样频率，f0为基波频率，Δf为离散频率间隔，且Δf＝fs/N；
[0278] 令：k0＝f0/Δf，k0为真实频谱的谱线位置，由于非同步采样或非整周期截断，k0一般为小数，不为整数。
[0279] 忽略负频率点处旁瓣的影响，公式(3)变为：
[0280]
[0281] 谱线确定单元，用于确定三根谱线：
[0282] 在信号预处理单元得到的加窗FFT频谱峰值附近区域，k0处频率点较大的三根谱线分别为：第k1、k2、k3根，k1、k2、k3均为正整数，k1～k3的关系为：k2＝k1+1，k3＝k2+1，这三根谱线对应的幅值分别为y1、y2、y3；
[0283] 记变量α＝k-k2，则-0.5≤α≤0.5；
[0284] 另记变量
[0285] 修正公式计算单元，用于计算三谱线插值算法的修正公式：
[0286] 根据公式(4)和(5)得到：
[0287]
[0288] 当N值比较大，例如：N＞1000时，公式(6)可以化简为一个函数β＝g(α)，其反函数为α＝g-1(β)。由于所采用的余弦窗系数均为实系数，其频率响应是偶对称的，因而g(·)和g-1(·)均为奇函数。
[0289] 采用多项式逼近方法计算奇函数α＝g-1(β)，表达式为：
[0290] α≈p11×β+p13×β3+…+p1pβp (7)
[0291] p11，p13；…p1p为多项式逼近的奇次项系数，p是奇数；
[0292] 根据公式(4)，求得电网信号第i次谐波的幅值Ai：
[0293]
[0294] 其中，i为正整数，yi为加窗FFT后第i次谐波的幅值；
[0295] 以基波为例，考虑到y2是离真实谱线点最近谱线，给予较大权重，可以得到：
[0296]
[0297] N＞1000时，窗函数系数为实系数，公式(9)表示为：
[0298] A1＝N-1(y3+2y2+y1)u(α)
[0299] 其中，u(α)为修正公式，且为偶函数，逼近多项式不含奇次项；
[0300] 三谱线修正逼近多项式如下：
[0301] u(α)＝(p20+p22α2+…+p2dαd) (10)
[0302] 公式(10)中，p20,p22…p2d为多项式逼近的偶次项系数，d为拟合的最高阶次，且d为偶数；
[0303] 基波参数计算单元，用于计算基波参数：
[0304] 考虑到y2是离真实谱线点最近的两根谱线，给予较大权重，可以计算基波频率f0、基波幅值A1：
[0305] f0＝k·Δf＝(k2+α)Δf (11)
[0306] A1＝N-1(y3+2y2+y1)(p20+p22α2+…+p2dαd) (12)
[0307] 根据公式(4)，计算基波的相位：
[0308]
[0309] 仿照基波参数的求取，根据公式(6)、(7)、(9)、(11)、(12)、 (13)进行各次谐波参数的分析；
[0310] 谐波参数确定单元，用于确定谐波参数：确定基波频率f0后，在范围(kf0-5,kf0+5)内，所述谱线确定单元、修正公式计算单元、基波参数计算单元重复进行计算，直到所有谐波参数计算完毕；
[0311] 误差分析单元，用于进行误差分析：分析基于互乘法窗函数的三谱线插值FFT算法的误差，并与常规窗函数插值算法精度进行比较。
[0312] 本领域的技术人员可以对本发明实施例进行各种修改和变型，倘若这些修改和变型在本发明权利要求及其等同技术的范围之内，则这些修改和变型也在本发明的保护范围之内。
[0313] 说明书中未详细描述的内容为本领域技术人员公知的现有技术。
