CN107643446A - 一种基于主瓣宽度的多谱线插值谐波分析方法及系统 - Google Patents
一种基于主瓣宽度的多谱线插值谐波分析方法及系统 Download PDFInfo
- Publication number
- CN107643446A CN107643446A CN201710685207.0A CN201710685207A CN107643446A CN 107643446 A CN107643446 A CN 107643446A CN 201710685207 A CN201710685207 A CN 201710685207A CN 107643446 A CN107643446 A CN 107643446A
- Authority
- CN
- China
- Prior art keywords
- mrow
- frequency
- formula
- main lobe
- msub
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
Abstract
一种基于主瓣宽度的多谱线插值谐波分析方法及系统，涉及谐波分析领域。该方法包括以下步骤：信号预处理；谱线数目确定；权值确定；计算多谱线插值算法的修正公式；计算基波参数；确定谐波参数；进行误差分析。本发明从电力系统的电网信号(电流信号或电压信号)加窗后的频域表达式入手，首次根据所加窗函数的主瓣宽度确定参与运算的谱线数目，并推导出各个谱线权值的大小。在此基础上进行谐波分析，计算谐波相关参数。通过算例表明，本发明利用主瓣内所有谱线进行插值分析的方法，精度高，实现了谐波参数的高精度测量，具有较好的使用价值。
Description
技术领域
本发明涉及谐波分析领域，具体来讲涉及一种基于主瓣宽度的多谱线插值谐波分析方法及系统。
背景技术
电网中谐波的问题日益严重，不仅恶化电能质量，对电网的安全稳定和经济运行也造成较大影响。因此，对电网中谐波参数进行高精度分析与测量，对于提高电能质量，减少谐波危害，维护电网安全稳定、高效运行是十分必要的。
电力系统谐波分析最常用的算法是采用FFT(Fast Fourier Transform，快速傅里叶变换)分析算法。但由于电网频率的不固定，FFT很难做到同步采样和整周期截断，所带来的频谱泄露和栅栏效应可能会导致谐波的频率、幅值和相位等参数测量不准，无法满足测量要求。加窗插值FFT算法被广泛用来解决这一问题，具体做法是：运用各种窗函数对信号进行截断，然后结合谱线插值FFT进行谐波分析，从而提高测量精度。
常用窗函数例如汉宁(Hanning)窗函数、布莱克曼(Blackman)窗函数、布莱克曼汉斯(Blackman-Harris)窗函数、纳托尔(Nuttall)窗函数、莱夫文森特(Rife-Vincent)窗函数以及各种组合窗函数。随着复杂组合窗函数的出现，主瓣宽度越来越宽，主瓣内的所含谱线数目越来越多，常规的双插值和三插值会出现丢失有用信息的情况，导致谐波分析精度降低。
发明内容
针对现有技术中存在的缺陷，本发明的目的在于提供一种基于主瓣宽度的多谱线插值谐波分析方法及系统，有效提高谐波分析的精度。
为达到以上目的，本发明提供一种基于主瓣宽度的多谱线插值谐波分析方法，包括步骤：
S1.上位机接收到互感器采集到的离散的电网信号x(n)，n为采样点的序数，n为自然数；采用离散余弦窗函数w(n)对电网信号x(n)进行加窗截断，得到加窗信号xw(n)：
xw(n)＝x(n)w(n) 式(1)
对式(1)的加窗信号xw(n)进行FFT变换后，忽略负频率点处旁瓣的影响，得到加窗FFT频谱X(k)：
其中W(·)为窗函数的频谱，k为正整数，X(k)表示第k次谐波的频谱，Ak为第k次谐波的幅值，j表示虚数单位，e是自然对数的底数，为第k次谐波的初始相位，k0为真实频谱的谱线位置；
S2.整数次谐波之间间隔为50Hz，非同步采样带来的泄露频谱线在谐波间隔中最多有H1根，H1＝int(50/Δf)，其中Δf为离散频率间隔；主瓣内的谱线根数H2为H2＝int(Δω/2πΔf)，多谱线的根数H确定为H＝min(H1,H2)；
设峰值频率点左右多条的谱线分别为k1＜…＜kh＜k＜kh+1＜…＜kH，对应的谱线幅值分别为y1，…，yh，yh+1，…，yH；记α＝k-kh-0.5，由于0≤k-kh≤1，则-0.5≤α≤0.5，另记
其中α和β均为系数，ci为加权系数，i＝1…h，ci为yi的权值，令ci＝c(H-i+1)，根据离散余弦窗函数w(n)和频域主瓣表达式，加余弦窗时得到
S4.根据式(4)和式(5)计算多谱线插值算法的修正式；
S5.计算基波参数，包括基波频率f0和基波幅值A1，根据式(4)得出基波的相位
S6.确定基波频率f0后，在范围(kf0-5,kf0+5)内重复步骤S2～S5，直到所有谐波参数计算完毕。
在上述技术方案的基础上，离散余弦窗函数w(n)的表达式为
其中，N为采样点数，N为正整数，n＝0,1,2…N-1；Σ表示求和；m为窗函数的累加次数，m＝0,1,2…M-1；M为窗函数项数，M为正整数；bm为窗函数系数。
在上述技术方案的基础上，所述窗函数系数bm满足约束条件
在上述技术方案的基础上，对式(1)的加窗信号进行FFT变换后，得到加窗FFT频谱
其中，第一次谐波为基波，f0为基波频率，且Δf＝fs/N。
在上述技术方案的基础上，所述S4中，根据式(4)和式(5)得到
ci为加权系数，i＝1…h，ci为yi和yM-i的权值，当N大于1000时，式(7)简化为β＝g(α)，其反函数为α＝g-1(β)；采用多项逼近方法计算奇函数β＝g-1(α)，α≈p11×β+p13×β3+…+p1pβp，其中，p11，p13，…，p1p为多项式逼近的奇次项系数，由式(4)可知信号幅值为得到
通过逼近方法，当N大于1000时，窗函数系数为实系数，式(10)可表示为u(·)为偶函数，逼近多项式不含奇次项，则多谱线修正逼近多项式为
u(α)＝(p20+p22α2+…+p2dαd) 式(11)
其中，p20，p22，…，p2d为多项式逼近的偶次项系数。
在上述技术方案的基础上，所述S5中，
f0＝k·Δf＝(k2+α+0.5)Δf，
在上述技术方案的基础上，所述电网信号包括电流信号、电压信号，所述真实频谱的谱线位置k0为小数。
本发明还提供一种基于主瓣宽度的多谱线插值谐波分析系统，设置于上位机，包括信号预处理单元、谱线数目确定单元、权值确定单元、修正公式计算单元、基波参数计算单元、谐波参数确定单元；
所述信号预处理单元，用于对上位机中来自互感器的离散的电网信号x(n)进行预处理；采用离散余弦窗函数w(n)对电网信号x(n)进行加窗截断，得到加窗信号xw(n)
xw(n)＝x(n)w(n) 式(1)
对式(1)的加窗信号xw(n)进行FFT变换后，忽略负频率点处旁瓣的影响，得到加窗FFT频谱X(k)：
其中W(·)为窗函数的频谱，k为正整数，X(k)表示第k次谐波的频谱，Ak为第k次谐波的幅值，j表示虚数单位，e是自然对数的底数，为第k次谐波的初始相位，k0为真实频谱的谱线位置；
所述谱线数目确定单元，用于确定多根谱线；整数次谐波之间间隔为50Hz，非同步采样带来的泄露频谱线在谐波间隔中最多有H1根，H1＝int(50/Δf)，其中Δf为离散频率间隔；主瓣内的谱线根数H2为H2＝int(Δω/2πΔf)，多谱线的根数H确定为H＝min(H1,H2)；
权值确定单元，用于计算多谱线所对应得权值大小；
设峰值频率点左右多条的谱线分别为k1＜…＜kh＜k＜kh+1＜…＜kH，对应的谱线幅值分别为y1，…，yh，yh+1，…，yH；记α＝k-kh-0.5，由于0≤k-kh≤1，则-0.5≤α≤0.5，另记
其中α和β均为系数，ci为加权系数，i＝1…h，ci为yi的权值，令ci＝c(H-i+1)，根据离散余弦窗函数w(n)和频域主瓣表达式，加余弦窗时得到
所述修正公式计算单元，用于根据式(4)和式(5)计算多谱线插值算法的修正公式；
所述基波参数计算单元，用于计算基波参数，包括基波频率f0和基波幅值A1，以及基波的相位
所述谐波参数确定单元，用于确定基波频率f0后，在范围(kf0-5,kf0+5)内重复步骤S2～S5，直到所有计算得到所有谐波参数。
在上述技术方案的基础上，系统还包括误差分析单元，在窗函数主瓣宽度不断增加，参与运算的谱线数目也相应增加的情况下，进行多谱线插值计算的误差分析。
在上述技术方案的基础上，所述电网信号包括电流信号、电压信号，所述真实频谱的谱线位置k0为小数。
本发明的有益效果在于：通过从电力系统的电网信号加窗后的频域表达式入手，根据所加窗函数的主瓣宽度确定参与运算的谱线数目，进而推导出各个谱线权值的大小。并在此基础上进行谐波分析，计算谐波相关参数。通过算例表明，本发明利用主瓣内所有谱线进行插值分析，提高了精度，实现了谐波参数的高精度测量，具有较好的使用价值。
附图说明
图1为本发明基于主瓣宽度的多谱线插值谐波分析方法流程图；
图2为本发明实施例中谐波信号的波形图；
图3为1阶Hanning窗函数的幅频特性图；
图4为2阶Hanning窗函数的幅频特性图；
图5为3阶Hanning窗函数的幅频特性图；
图6为4阶Hanning窗函数的幅频特性图；
图7为本发明基于主瓣宽度的多谱线插值谐波分析系统示意图。
具体实施方式
以下结合附图及实施例对本发明作进一步详细说明。
如图1所示，本发明基于主瓣宽度的多谱线插值谐波分析方法，包括如下步骤：
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为窗函数系数。窗函数系数bm满足约束条件如下：
对式(1)的加窗信号xw(n)进行FFT变换后，得到加窗FFT频谱：
其中，W(·)为窗函数的频谱，k为正整数，X(k)表示第k次谐波的频谱，Ak为第k次谐波的幅值，j表示虚数单位，e是自然对数的底数，为第k次谐波的初始相位，第一次谐波为基波，fs为采样频率，f0为基波频率，Δf为离散频率间隔，且Δf＝fs/N。
令：k0＝f0/Δf，k0为真实频谱的谱线位置，真实频谱的谱线位置k0为小数，忽略负频率点处旁瓣的影响，式(3)变为：
S2.根据主瓣宽度确定谱线数目。
具体的，对于整数次谐波之间间隔为50Hz，那么非同步采样带来的泄露频谱线在谐波间隔中最多有H1根，H1＝int(50/Δf)。
若所加窗函数的主瓣宽度为Δω rad/s，主瓣内的谱线根数H2为H2＝int(Δω/2πΔf)，多谱线的根数H确定为H＝min(H1,H2)。
S3.确定权值。
具体的，令
对信号的非同步采样或非整周期数据截断时，由于栅栏效应，信号的峰值频率f0＝k0·Δf很难刚好位于离散谱线的频点上，即k0一般不为整数。设峰值频率点左右多条的谱线分别为k1＜…＜kh＜k＜kh+1＜…＜kH，对应的谱线幅值分别为y1，…，yh，yh+1，…，yH。
记α＝k-kh-0.5，由于0≤k-kh≤1，则-0.5≤α≤0.5，另记：
其中α和β均为系数，ci为加权系数，i＝1…h，ci为yi的权值。考虑到计算的方便，令ci＝c(H-i+1)，余弦窗函数的表达式如式(2)所示时，根据频域主瓣表达式，加余弦窗时候，可以求得：
S4.计算多普线插值算法的修正公式。
具体的，根据式(4)和式(5)得到
其中，ci为加权系数，i＝1…h，ci为yi和yM-i的权值。
当N大于1000时，式(7)可以简化为β＝g(α)，其反函数为α＝g-1(β)；由于所采用的余弦窗系数均为实系数，其频率响应是偶对称的，因而g(·)和g-1(·)均为奇函数。采用多项逼近方法计算奇函数β＝g-1(α)，表达式为：
α≈p11×β+p13×β3+…+p1pβp 式(8)
其中，p11，p13，…，p1p为多项式逼近的奇次项系数，由式(4)可知信号幅值Ai为：
可以得到：
通过类似于式(10)的逼近方法，当N大于1000时，窗函数系数为实系数，式(10)可表示为：u(·)为偶函数，逼近多项式不含奇次项，则多谱线修正逼近多项式为：
u(α)＝(p20+p22α2+…+p2dαd) 式(11)
其中，p20，p22，…，p2d为多项式逼近的偶次项系数。
S5.计算基波参数。
具体的，包括计算基波频率f0、基波幅值A1：
f0＝k·Δf＝(k2+α+0.5)Δf 式(12)
根据式(4)得出基波的相位
S6.确定谐波参数。
具体的，确定基波频率f0后，在范围(kf0-5,kf0+5)内重复步骤S2～S5，直到所有谐波参数计算完毕。
上述步骤S1～S6完成后，可以增加一个误差分析的步骤，当窗函数主瓣宽度不断增加，参与运算的谱线数目也相应增加，进行多谱线插值算法的误差分析，确保高精度。
下面通过一个具体实施例来进行详细说明。
如图2所示，是一个选定电网信号谐波的波形，基于主瓣宽度的多谱线插值谐波分析方法包括：
a1.信号预处理。具体的，将互感器采集到的离散电网信号x(n)，信号模型参见表1所示。为了验证所提算法的精度，进行10次谐波仿真分析。信号模型为：
其中，基波频率f0为50.5Hz，采样频率fs为5120Hz，采样点数N为1024。
表1、电网信号的谐波参数
自乘法窗函数的定义为：若干个相同窗函数进行相乘所构成的新的窗函数。乘法窗函数的公式：
其中，w(n)为基本窗函数，p为参加乘法的窗函数个数，称为乘法窗函数的阶数，p＝1时即为原窗函数。乘法窗函数的主瓣特性和旁瓣特性同基本窗函数选取和乘法窗的阶数相关。根据原窗函数的不同和乘法窗阶数的不同可以构造多种乘法窗函数。本实施例选取经典的Hanning窗函数作为原窗函数构造新的乘法窗函数。加窗FFT变换后得到信号频谱：
a2.确定谱线数目。由于频域分辨率Δf＝(2π)/N，P阶Hanning乘法窗函数的主瓣宽度(MB)为(4(p+1)π)/N，主瓣内谱线根数为2(p+1)。Hanning自乘法窗函数的频谱特性如表2所示。
表2、Hanning自乘法窗函数频谱特性
a3.计算相关谱线系数，参见表3。
表3、谱线的加权系数
| 窗类型 | c1 | c2 | c3 | c4 | c5 |
| 1阶 | 0.6329 | 1 | 0 | 0 | 0 |
| 2阶 | 0.4140 | 0.8272 | 1 | 0 | 0 |
| 3阶 | 0.3006 | 0.6329 | 0.9004 | 1 | 0 |
| 4阶 | 0.2344 | 0.5050 | 0.7571 | 0.9356 | 1 |
由于4阶Hanning乘法窗函数主主瓣内有10根谱线，两次谐波间隔中最多也为10根。
a4.根据步骤S4的推导，表2中的4种窗函数的修正公式α＝g-1(β)，u(α)分别如下：
(1)1阶Hanning自乘法窗函数：
α＝1.1301β-0.1836β3+0.0811β5
u(α)＝0.5355+0.1758α2+0.1043α4
如图3所示，为1阶Hanning自乘法窗函数幅频特性。
(2)2阶Hanning自乘法窗函数：
α＝1.5438β-0.3134β3+0.1618β5
u(α)＝0.7543+0.2249α2+0.1545α4
如图4所示，为2阶Hanning自乘法窗函数幅频特性。
(3)3阶Hanning自乘法窗函数：
α＝1.8496β-0.4045β3+0.2181β5
u(α)＝0.5745+0.1716α2+0.1053α4
如图5所示，为3阶Hanning自乘法窗函数的幅频特性。
(4)4阶Hanning自乘法窗函数：
α＝2.1895β-0.5473β3+0.2937β5
u(α)＝0.4668+0.1854α2+0.1512α4
如图6所示，为4阶Hanning自乘法窗函数的幅频特性。
a5.计算基波参数，公式如下：
f0＝k·Δf＝(k2+α+0.5)Δf；
a6.确定基波频率f0后，在范围(kf0-5,kf0+5)内重复步骤a2～a5，直到所有谐波参数计算完毕。
a7.进行误差分析。仿真实验的误差分析比较结果参见表3～表6所示，其中，DAi表示基波和各次谐波幅值测量值的相对误差；Df0表示基波频率测量值的相对误差；表示基波和各次谐波初始相位测量值的相对误差，均用百分比表示。
表4、Hanning各阶乘法窗的基波频率相对误差
| 阶次 | 1 | 2 | 3 | 4 |
| Df0 | -1.97E-7 | 1.43E-8 | 2.06E-9 | -2.55E-10 |
表5、Hanning窗各阶幅值相对误差
| 比较项 | 1阶 | 2阶 | 3阶 | 4阶 |
| DA1 | -1.23E-5 | -2.37E-6 | -7.65E-7 | -7.35E-8 |
| DA2 | -7.96E-5 | -6.25E-5 | -1.96E-5 | -7.53E-6 |
| DA3 | -4.53E-6 | -2.12E-6 | -6.49E-7 | -3.66E-7 |
| DA4 | 6.59E-5 | 4.63E-6 | 2.51E-6 | -3.66E-7 |
| DA5 | 2.58E-5 | -3.69E-5 | 1.38E-6 | 9.87E-7 |
| DA6 | 7.75E-5 | 3.45E-6 | 2.42E-6 | 9.57E-7 |
| DA7 | -2.61E-5 | 3.80E-6 | 1.17E-6 | -7.85E-7 |
| DA8 | -7.37E-5 | -4.25E-6 | -2.75E-6 | 9.94E-7 |
| DA9 | -5.55E-5 | -3.52E-6 | -1.84E-6 | -1.00E-6 |
| DA10 | -1.10E-4 | -3.53E-5 | -1.13E-5 | -1.09E-5 |
表6、Hanning窗各阶相位相对误差
从表4～表6数据可以看出，本发明实施例所推导的多谱线插值算法具有良好的精度。
如图7所示，本发明基于主瓣宽度的多谱线插值谐波分析系统，包括信号预处理单元、谱线数目确定单元、权值确定单元、修正公式计算单元、基波参数计算单元和谐波参数确定单元。
所述信号预处理单元，用于对上位机中来自互感器的离散的电网信号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为窗函数系数。窗函数系数bm满足约束条件如下：
对式(1)的加窗信号xw(n)进行FFT变换后，得到加窗FFT频谱：
其中，W(·)为窗函数的频谱，k为正整数，X(k)表示第k次谐波的频谱，Ak为第k次谐波的幅值，j表示虚数单位，e是自然对数的底数，为第k次谐波的初始相位，第一次谐波为基波，fs为采样频率，f0为基波频率，Δf为离散频率间隔，且Δf＝fs/N。
令：k0＝f0/Δf，k0为真实频谱的谱线位置，真实频谱的谱线位置k0为小数，忽略负频率点处旁瓣的影响，式(3)变为：
所述谱线数目确定单元，用于确定多根谱线。
具体的，对于整数次谐波之间间隔为50Hz，那么非同步采样带来的泄露频谱线在谐波间隔中最多有H1根，H1＝int(50/Δf)。
若所加窗函数的主瓣宽度为Δω rad/s，主瓣内的谱线根数H2为H2＝int(Δω/2πΔf)，多谱线的根数H确定为H＝min(H1,H2)。
所述权值确定单元，用于计算多谱线所对应得权值大小；具体的，令
对信号的非同步采样或非整周期数据截断时，由于栅栏效应，信号的峰值频率f0＝k0·Δf很难刚好位于离散谱线的频点上，即k0一般不为整数。设峰值频率点左右多条的谱线分别为k1＜…＜kh＜k＜kh+1＜…＜kH，对应的谱线幅值分别为y1，…，yh，yh+1，…，yH。
记α＝k-kh-0.5，由于0≤k-kh≤1，则-0.5≤α≤0.5，另记：
其中α和β均为系数，ci为加权系数，i＝1…h，ci为yi的权值。考虑到计算的方便，令ci＝c(H-i+1)，余弦窗函数的表达式如式(2)所示时，根据频域主瓣表达式，加余弦窗时候，可以求得：
所述修正公式计算单元，用于计算多谱线插值算法的修正公式。具体的，根据式(4)和式(5)得到
其中，ci为加权系数，i＝1…h，ci为yi和yM-i的权值。
当N大于1000时，式(7)可以简化为β＝g(α)，其反函数为α＝g-1(β)；由于所采用的余弦窗系数均为实系数，其频率响应是偶对称的，因而g(·)和g-1(·)均为奇函数。采用多项逼近方法计算奇函数β＝g-1(α)，表达式为：
α≈p11×β+p13×β3+…+p1pβp 式(8)
其中，p11，p13，…，p1p为多项式逼近的奇次项系数，由式(4)可知信号幅值Ai为：
可以得到：
通过类似于式(10)的逼近方法，当N大于1000时，窗函数系数为实系数，式(10)可表示为：u(·)为偶函数，逼近多项式不含奇次项，则多谱线修正逼近多项式为：
u(α)＝(p20+p22α2+…+p2dαd) 式(11)
其中，p20，p22，…，p2d为多项式逼近的偶次项系数。
所述基波参数计算单元，用于计算基波参数。具体的，包括计算基波频率f0、基波幅值A1：
f0＝k·Δf＝(k2+α+0.5)Δf 式(12)
根据式(4)得出基波的相位
所述谐波参数确定单元，用于确定谐波参数。具体的，确定基波频率f0后，在范围(kf0-5,kf0+5)内重复步骤S2～S5，直到所有谐波参数计算完毕。
本发明还可以设置一个误差分析单元，在窗函数主瓣宽度不断增加，参与运算的谱线数目也相应增加的情况下，进行多谱线插值计算的误差分析，验证本发明分析的精确度。
本发明不局限于上述实施方式，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也视为本发明的保护范围之内。本说明书中未作详细描述的内容属于本领域专业技术人员公知的现有技术。
Claims (10)
1.一种基于主瓣宽度的多谱线插值谐波分析方法，其特征在于，包括步骤：
S1.上位机接收到互感器采集到的离散的电网信号x(n)，n为采样点的序数，n为自然数；采用离散余弦窗函数w(n)对电网信号x(n)进行加窗截断，得到加窗信号xw(n)：
xw(n)＝x(n)w(n) 式(1)
对式(1)的加窗信号xw(n)进行FFT变换后，忽略负频率点处旁瓣的影响，得到加窗FFT频谱X(k)：
其中W(·)为窗函数的频谱，k为正整数，X(k)表示第k次谐波的频谱，Ak为第k次谐波的幅值，j表示虚数单位，e是自然对数的底数，为第k次谐波的初始相位，k0为真实频谱的谱线位置；
S2.整数次谐波之间间隔为50Hz，非同步采样带来的泄露频谱线在谐波间隔中最多有H1根，H1＝int(50/Δf)，其中Δf为离散频率间隔；主瓣内的谱线根数H2为H2＝int(Δω/2πΔf)，多谱线的根数H确定为H＝min(H1,H2)；
S3.令
设峰值频率点左右多条的谱线分别为k1＜…＜kh＜k＜kh+1＜…＜kH，对应的谱线幅值分别为y1，…，yh，yh+1，…，yH；记α＝k-kh-0.5，由于0≤k-kh≤1，则-0.5≤α≤0.5，另记
其中α和β均为系数，ci为加权系数，i＝1…h，ci为yi的权值，令ci＝c(H-i+1)，根据离散余弦窗函数w(n)和频域主瓣表达式，加余弦窗时得到
S4.根据式(4)和式(5)计算多谱线插值算法的修正式；
S5.计算基波参数，包括基波频率f0和基波幅值A1，根据式(4)得出基波的相位
S6.确定基波频率f0后，在范围(kf0-5,kf0+5)内重复步骤S2～S5，直到所有谐波参数计算完毕。
2.如权利要求1所述基于主瓣宽度的多谱线插值谐波分析方法，其特征在于：离散余弦窗函数w(n)的表达式为
其中，N为采样点数，N为正整数，n＝0,1,2...N-1；Σ表示求和；m为窗函数的累加次数，m＝0,1,2...M-1；M为窗函数项数，M为正整数；bm为窗函数系数。
3.如权利要求2所述基于主瓣宽度的多谱线插值谐波分析方法，其特征在于：所述窗函数系数bm满足约束条件
<mfenced open = "" close = "">
<mtable>
<mtr>
<mtd>
<mrow>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>m</mi>
<mo>=</mo>
<mn>0</mn>
</mrow>
<mrow>
<mi>M</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</munderover>
<msub>
<mi>b</mi>
<mi>m</mi>
</msub>
<mo>=</mo>
<mn>1</mn>
<mo>,</mo>
</mrow>
</mtd>
<mtd>
<mrow>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>m</mi>
<mo>=</mo>
<mn>0</mn>
</mrow>
<mrow>
<mi>M</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</munderover>
<msup>
<mrow>
<mo>(</mo>
<mo>-</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
<mi>m</mi>
</msup>
<msub>
<mi>b</mi>
<mi>m</mi>
</msub>
<mo>=</mo>
<mn>0.</mn>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
4.如权利要求2所述基于主瓣宽度的多谱线插值谐波分析方法，其特征在于：对式(1)的加窗信号进行FFT变换后，得到加窗FFT频谱
其中，第一次谐波为基波，f0为基波频率，且Δf＝fs/N。
5.如权利要求1所述基于主瓣宽度的多谱线插值谐波分析方法，其特征在于：所述S4中，根据式(4)和式(5)得到
ci为加权系数，i＝1…h，ci为yi和yM-i的权值，当N大于1000时，式(7)简化为β＝g(α)，其反函数为α＝g-1(β)；采用多项逼近方法计算奇函数β＝g-1(α)，α≈p11×β+p13×β3+…+p1pβp，其中，p11,p13,…,p1p为多项式逼近的奇次项系数，由式(4)可知信号幅值为得到
通过逼近方法，当N大于1000时，窗函数系数为实系数，式(10)可表示为u(·)为偶函数，逼近多项式不含奇次项，则多谱线修正逼近多项式为
u(α)＝(p20+p22α2+…+p2dαd) 式(11)
其中，p20,p22,…,p2d为多项式逼近的偶次项系数。
6.如权利要求5所述基于主瓣宽度的多谱线插值谐波分析方法，其特征在于：所述S5中，
f0＝k·Δf＝(k2+α+0.5)Δf，
<mrow>
<msub>
<mi>A</mi>
<mn>1</mn>
</msub>
<mo>=</mo>
<msup>
<mi>N</mi>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msup>
<mrow>
<mo>(</mo>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>i</mi>
<mo>=</mo>
<mi>m</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
<mi>M</mi>
</munderover>
<msub>
<mi>c</mi>
<mrow>
<mo>(</mo>
<mi>M</mi>
<mo>-</mo>
<mi>i</mi>
<mo>+</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
</msub>
<msub>
<mi>y</mi>
<mi>i</mi>
</msub>
<mo>+</mo>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>i</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>m</mi>
</munderover>
<msub>
<mi>c</mi>
<mi>i</mi>
</msub>
<msub>
<mi>y</mi>
<mi>i</mi>
</msub>
<mo>)</mo>
</mrow>
<mrow>
<mo>(</mo>
<msub>
<mi>p</mi>
<mn>20</mn>
</msub>
<mo>+</mo>
<msub>
<mi>p</mi>
<mn>22</mn>
</msub>
<msup>
<mi>&alpha;</mi>
<mn>2</mn>
</msup>
<mo>+</mo>
<mo>...</mo>
<mo>+</mo>
<msub>
<mi>p</mi>
<mrow>
<mn>2</mn>
<mi>d</mi>
</mrow>
</msub>
<msup>
<mi>&alpha;</mi>
<mi>d</mi>
</msup>
<mo>)</mo>
</mrow>
<mo>,</mo>
</mrow>
7.如权利要求1-6任一所述基于主瓣宽度的多谱线插值谐波分析方法，其特征在于：所述电网信号包括电流信号、电压信号，所述真实频谱的谱线位置k0为小数。
8.一种基于主瓣宽度的多谱线插值谐波分析系统，设置于上位机，其特征在于，包括信号预处理单元、谱线数目确定单元、权值确定单元、修正公式计算单元、基波参数计算单元、谐波参数确定单元；
所述信号预处理单元，用于对上位机中来自互感器的离散的电网信号x(n)进行预处理；采用离散余弦窗函数w(n)对电网信号x(n)进行加窗截断，得到加窗信号xw(n)
xw(n)＝x(n)w(n) 式(1)
对式(1)的加窗信号xw(n)进行FFT变换后，忽略负频率点处旁瓣的影响，得到加窗FFT频谱X(k)：
其中W(·)为窗函数的频谱，k为正整数，X(k)表示第k次谐波的频谱，Ak为第k次谐波的幅值，j表示虚数单位，e是自然对数的底数，为第k次谐波的初始相位，k0为真实频谱的谱线位置；
所述谱线数目确定单元，用于确定多根谱线；整数次谐波之间间隔为50Hz，非同步采样带来的泄露频谱线在谐波间隔中最多有H1根，H1＝int(50/Δf)，其中Δf为离散频率间隔；主瓣内的谱线根数H2为H2＝int(Δω/2πΔf)，多谱线的根数H确定为H＝min(H1,H2)；
权值确定单元，用于计算多谱线所对应得权值大小；
令
设峰值频率点左右多条的谱线分别为k1＜…＜kh＜k＜kh+1＜…＜kH，对应的谱线幅值分别为y1，…，yh，yh+1，…，yH；记α＝k-kh-0.5，由于0≤k-kh≤1，则-0.5≤α≤0.5，另记
其中α和β均为系数，ci为加权系数，i＝1…h，ci为yi的权值，令ci＝c(H-i+1)，根据离散余弦窗函数w(n)和频域主瓣表达式，加余弦窗时得到
所述修正公式计算单元，用于根据式(4)和式(5)计算多谱线插值算法的修正公式；
所述基波参数计算单元，用于计算基波参数，包括基波频率f0和基波幅值A1，以及基波的相位
所述谐波参数确定单元，用于确定基波频率f0后，在范围(kf0-5,kf0+5)内重复步骤S2～S5，直到所有计算得到所有谐波参数。
9.如权利要求8所述基于主瓣宽度的多谱线插值谐波分析系统，其特征在于：系统还包括误差分析单元，在窗函数主瓣宽度不断增加，参与运算的谱线数目也相应增加的情况下，进行多谱线插值计算的误差分析。
10.如权利要求8或9所述基于主瓣宽度的多谱线插值谐波分析系统，其特征在于：所述电网信号包括电流信号、电压信号，所述真实频谱的谱线位置k0为小数。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710685207.0A CN107643446B (zh) | 2017-08-11 | 2017-08-11 | 一种基于主瓣宽度的多谱线插值谐波分析方法及系统 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710685207.0A CN107643446B (zh) | 2017-08-11 | 2017-08-11 | 一种基于主瓣宽度的多谱线插值谐波分析方法及系统 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN107643446A true CN107643446A (zh) | 2018-01-30 |
| CN107643446B CN107643446B (zh) | 2019-11-08 |
Family
ID=61110305
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710685207.0A Active CN107643446B (zh) | 2017-08-11 | 2017-08-11 | 一种基于主瓣宽度的多谱线插值谐波分析方法及系统 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107643446B (zh) |
Cited By (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108427031A (zh) * | 2018-05-02 | 2018-08-21 | 三峡大学 | 基于多项式拟合及非干扰区域划分的间谐波检测方法 |
| CN108519512A (zh) * | 2018-03-23 | 2018-09-11 | 深圳市计量质量检测研究院 | 用于高精度谐波测量的方法和系统 |
| CN109884571A (zh) * | 2019-02-20 | 2019-06-14 | 武汉电力职业技术学院 | 一种基于非标准器多传感器融合的直流互感器计量方法 |
| CN110954746A (zh) * | 2019-11-27 | 2020-04-03 | 云南电网有限责任公司电力科学研究院 | 一种基于四项Nuttall余弦窗的六插值FFT算法 |
| CN111308198A (zh) * | 2020-03-10 | 2020-06-19 | 国网江苏省电力有限公司扬州供电分公司 | 一种基于Hanning窗的加窗插值DFT的谐波测量装置及测量方法 |
| CN113640579A (zh) * | 2021-10-13 | 2021-11-12 | 四川大学 | 基于双重谱线变换的谐波测量方法、电子设备及存储介质 |
| CN115825557A (zh) * | 2022-11-25 | 2023-03-21 | 国网四川省电力公司映秀湾水力发电总厂 | 基于谐波分量置零的广义谐波分析方法、装置及介质 |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| DE19515712A1 (de) * | 1994-04-28 | 1995-11-02 | Advantest Corp | Vorrichtung zum Kompensieren von Dopplerverschiebungen |
| JP3170600B2 (ja) * | 1997-03-26 | 2001-05-28 | 防衛庁技術研究本部長 | 波形解析方式 |
| CN104597321A (zh) * | 2015-01-28 | 2015-05-06 | 常洪山 | 基于四条离散傅里叶复数谱线的信号频率测量方法及装置 |
| CN104897961A (zh) * | 2015-06-17 | 2015-09-09 | 中南民族大学 | 基于互乘法窗函数的三谱线插值fft谐波分析方法及系统 |
| CN104897960A (zh) * | 2015-06-15 | 2015-09-09 | 中南民族大学 | 基于加窗四谱线插值fft的谐波快速分析方法及系统 |
-
2017
- 2017-08-11 CN CN201710685207.0A patent/CN107643446B/zh active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| DE19515712A1 (de) * | 1994-04-28 | 1995-11-02 | Advantest Corp | Vorrichtung zum Kompensieren von Dopplerverschiebungen |
| JP3170600B2 (ja) * | 1997-03-26 | 2001-05-28 | 防衛庁技術研究本部長 | 波形解析方式 |
| CN104597321A (zh) * | 2015-01-28 | 2015-05-06 | 常洪山 | 基于四条离散傅里叶复数谱线的信号频率测量方法及装置 |
| CN104897960A (zh) * | 2015-06-15 | 2015-09-09 | 中南民族大学 | 基于加窗四谱线插值fft的谐波快速分析方法及系统 |
| CN104897961A (zh) * | 2015-06-17 | 2015-09-09 | 中南民族大学 | 基于互乘法窗函数的三谱线插值fft谐波分析方法及系统 |
Non-Patent Citations (2)
| Title |
|---|
| TIANYUAN TAN ET AL.: "Harmonic analysis based on time domain mutual-multiplication window", 《J. MOD. POWER SYST. CLEAN ENERGY》 * |
| 康维 等: "一种改进FFT多谱线插值谐波分析方法", 《电测与仪表》 * |
Cited By (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108519512A (zh) * | 2018-03-23 | 2018-09-11 | 深圳市计量质量检测研究院 | 用于高精度谐波测量的方法和系统 |
| CN108519512B (zh) * | 2018-03-23 | 2020-05-22 | 深圳市计量质量检测研究院 | 用于高精度谐波测量的方法和系统 |
| CN108427031A (zh) * | 2018-05-02 | 2018-08-21 | 三峡大学 | 基于多项式拟合及非干扰区域划分的间谐波检测方法 |
| CN108427031B (zh) * | 2018-05-02 | 2020-08-04 | 三峡大学 | 基于多项式拟合及非干扰区域划分的间谐波检测方法 |
| CN109884571A (zh) * | 2019-02-20 | 2019-06-14 | 武汉电力职业技术学院 | 一种基于非标准器多传感器融合的直流互感器计量方法 |
| CN109884571B (zh) * | 2019-02-20 | 2021-06-08 | 武汉电力职业技术学院 | 一种基于非标准器多传感器融合的直流互感器计量方法 |
| CN110954746A (zh) * | 2019-11-27 | 2020-04-03 | 云南电网有限责任公司电力科学研究院 | 一种基于四项Nuttall余弦窗的六插值FFT算法 |
| CN111308198A (zh) * | 2020-03-10 | 2020-06-19 | 国网江苏省电力有限公司扬州供电分公司 | 一种基于Hanning窗的加窗插值DFT的谐波测量装置及测量方法 |
| CN111308198B (zh) * | 2020-03-10 | 2021-09-24 | 国网江苏省电力有限公司扬州供电分公司 | 一种基于Hanning窗的加窗插值DFT的谐波测量方法 |
| CN113640579A (zh) * | 2021-10-13 | 2021-11-12 | 四川大学 | 基于双重谱线变换的谐波测量方法、电子设备及存储介质 |
| CN115825557A (zh) * | 2022-11-25 | 2023-03-21 | 国网四川省电力公司映秀湾水力发电总厂 | 基于谐波分量置零的广义谐波分析方法、装置及介质 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN107643446B (zh) | 2019-11-08 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107643446A (zh) | 一种基于主瓣宽度的多谱线插值谐波分析方法及系统 | |
| CN104897961B (zh) | 基于互乘法窗函数的三谱线插值fft谐波分析方法及系统 | |
| CN104897960A (zh) | 基于加窗四谱线插值fft的谐波快速分析方法及系统 | |
| Duda | DFT interpolation algorithm for Kaiser–Bessel and Dolph–Chebyshev windows | |
| Luo et al. | Interpolated DFT algorithms with zero padding for classic windows | |
| CN109030941A (zh) | Hanning自乘卷积窗FFT三谱线插值谐波分析方法 | |
| CN109521275B (zh) | 一种同步相量确定方法、系统、装置及可读存储介质 | |
| CN105548739B (zh) | 一种避雷器运行状态信号处理方法 | |
| CN104122443B (zh) | Iec框架下的邻近谐波间谐波分离测量方法 | |
| CN103207319A (zh) | 数字化变电站电力信号非同步采样条件下的谐波测量方法 | |
| CN110728331A (zh) | 一种改进最小二乘支持向量机的谐波发射水平评估方法 | |
| CN105137180A (zh) | 基于六项余弦窗四谱线插值的高精度谐波分析方法 | |
| CN102818930A (zh) | 一种高精度快速计算电力谐波参数的方法 | |
| CN102495285B (zh) | 对称窗函数功率重心估计电力谐波参数的方法 | |
| CN105486921A (zh) | 凯撒三阶互卷积窗三谱线插值的谐波与间谐波检测方法 | |
| CN115575707A (zh) | 基于改进fft算法与小波变换结合的谐波检测装置及方法 | |
| Li et al. | Frequency estimation based on modulation FFT and MUSIC algorithm | |
| CN111579867A (zh) | 一种电力系统中谐波和间谐波的测量方法及装置 | |
| CN110954746A (zh) | 一种基于四项Nuttall余弦窗的六插值FFT算法 | |
| CN103543331A (zh) | 一种计算电信号谐波和间谐波的方法 | |
| Jiao et al. | An approach for electrical harmonic analysis based on interpolation DFT | |
| CN105223906B (zh) | 一种数控系统伺服驱动信号谐波频率的自动校正方法 | |
| Wen et al. | Performance comparison of windowed interpolation FFT and quasisynchronous sampling algorithm for frequency estimation | |
| Xuming et al. | Harmonic analysis based on Blackman-Harris self-multiplication window | |
| CN111579868B (zh) | 一种高次谐波的测量方法及装置 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |