CN113032716A - 基于加窗插值和Prony算法的谐波与间谐波分析方法 - Google Patents
基于加窗插值和Prony算法的谐波与间谐波分析方法 Download PDFInfo
- Publication number
- CN113032716A CN113032716A CN201911351006.2A CN201911351006A CN113032716A CN 113032716 A CN113032716 A CN 113032716A CN 201911351006 A CN201911351006 A CN 201911351006A CN 113032716 A CN113032716 A CN 113032716A
- Authority
- CN
- China
- Prior art keywords
- signal
- harmonic
- frequency
- window
- formula
- 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
- 238000004458 analytical method Methods 0.000 title claims abstract description 14
- 230000003595 spectral effect Effects 0.000 claims abstract description 98
- 238000001228 spectrum Methods 0.000 claims abstract description 21
- 238000000034 method Methods 0.000 claims abstract description 16
- 230000014509 gene expression Effects 0.000 claims description 21
- 238000005070 sampling Methods 0.000 claims description 20
- 239000011159 matrix material Substances 0.000 claims description 18
- 238000004364 calculation method Methods 0.000 claims description 4
- 238000013459 approach Methods 0.000 claims description 3
- 238000000354 decomposition reaction Methods 0.000 claims description 3
- 238000000605 extraction Methods 0.000 claims description 3
- 238000013178 mathematical model Methods 0.000 claims description 3
- 230000010355 oscillation Effects 0.000 claims description 3
- 239000004576 sand Substances 0.000 claims description 3
- 239000000126 substance Substances 0.000 claims description 3
- 238000006243 chemical reaction Methods 0.000 abstract 1
- 230000006378 damage Effects 0.000 description 3
- 230000002401 inhibitory effect Effects 0.000 description 2
- 230000005540 biological transmission Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/14—Fourier, Walsh or analogous domain transformations, e.g. Laplace, Hilbert, Karhunen-Loeve, transforms
- G06F17/141—Discrete Fourier transforms
- G06F17/142—Fast Fourier transforms, e.g. using a Cooley-Tukey type algorithm
Landscapes
- Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Mathematical Physics (AREA)
- General Physics & Mathematics (AREA)
- Computational Mathematics (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Pure & Applied Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Theoretical Computer Science (AREA)
- Algebra (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Discrete Mathematics (AREA)
- Operations Research (AREA)
- Complex Calculations (AREA)
Abstract
本发明公开了一种基于加窗插值和Prony算法的谐波与间谐波分析方法。该方法为:首先比较常用窗函数的频域特性,选择合适的窗函数;然后对待测信号进行加窗得到加窗信号,并进行FFT变换,得到信号的频谱;接着在频谱中搜索谱峰,计算每个谱峰附近两根谱线之间的相位差;最后判断主瓣干扰是否存在,若不存在,则使用插值算法校正对应的谐波或间谐波成分的参数;若存在,则不处理当前主瓣,利用插值算法算出的参数构建一个新信号,用原信号减去新信号得到残差信号,用Prony算法分析残差信号,得到频率相近的谐波和间谐波的参数。本发明线路简单,硬件成本低,可靠性高,能够有效分析出信号中谐波与间谐波的参数。
Description
技术领域
本发明属于电力***信号分析领域,特别是一种基于加窗插值和Prony算法的谐波与间谐波分析方法。
背景技术
电能是由电力企业向用户提供的一种特殊商品,在理想情况下,公用电网是以恒定的频率、标准的正弦波形、规定的电压水平和对称的三相电压电流向用户供电,由电厂、电网和用电方共同保证电能质量。但在实际电网中,由于***中的发电机、变压器、输电线路等设备的固有缺陷及非线性负荷的大量使用等原因,使得电力企业提供的电压、电流信号波形和理想波形相比存在着一定程度的偏差,当偏差较大时,会给电网带来严重危害。电压、频率和波形是电能质量的三个因素,如果它们任意一个不在限定的范围内,就会给电力***和用户带来不同程度的危害。
近年来,随着电力电子设备和冲击性负荷的增多,使得电网中谐波及间谐波的含量大量增加。谐波的存在会降低变压器容量,增加电机损耗,干扰电气设备正常工作,严重时甚至会损毁设备。与传统谐波相比,间谐波可以是直流到高次谐波间的任意频率,因此间谐波不但具有谐波的危害,其特性还会导致电压闪变、继电器误动作及引起无源滤波器过载等多种危害。同时,随着分布式电源的兴起和智能电网的快速发展,间谐波的来源越来越多,对电力***的安全运行造成严重的威胁,因此对间谐波的治理迫在眉睫。然而,治理间谐波的前提条件是能够准确检测出间谐波参数,只有分析出间谐波信号,才能对其进行抑制和治理,进而改善供电环境,确保供电设备正常运转。因此,准确快速的检测出间谐波,对于抑制间谐波进而改善提高电能质量具有重要的现实意义。
目前文献中一般使用FFT插值算法来分析谐波信号,当信号中不含频率相近的谐波和间谐波信号时,现有FFT方法及其改进算法的具有较高的精度和稳定性,但当待分析信号中含有频率相近的谐波和间谐波成分时,该方法的频率分辨率不能满足需求。而Prony算法可使用一组指数函数的线性组合来拟合待分析信号,在理论上具有无限小的频率分辨率,当电网信号为宽带信号且所含谐波和间谐波成分较多时,算法准确性会受到一定影响,因此将两种算法结合可以弥补两种算法的局限性,同时,能保证精确度的情况下,准确分析出频率相近的间谐波信号。
发明内容
本发明的目的在于提供一种结构简单、运算速度高、可靠性高的基于加窗插值和Prony算法的电力***谐波与间谐波分析方法。
实现本发明目的的技术解决方案为:一种基于加窗插值和Prony算法的谐波与间谐波分析方法,包括以下步骤:
步骤1、分析比较基本窗函数的频域特性,选择性能最优的窗函数;
步骤2、对待测信号进行加窗得到加窗信号,对加窗信号进行FFT变换,得到信号的频谱;
步骤3、在频谱中搜索谱峰,仅搜索幅值大于设定阈值的谱线,计算每个谱峰两根相邻谱线之间的相位差;
步骤4、判断主瓣干扰是否存在:
若主瓣干扰不存在,则使用插值算法校正对应的谐波或间谐波成分的参数;
若主瓣干扰存在,则不处理当前主瓣,利用插值算法算出的参数构建一个新信号,用原信号减去新信号得到残差信号,用Prony算法分析残差信号,得到含谐波和间谐波的信号的参数值。
进一步地,步骤1所述的分析比较基本窗函数的频域特性,选择性能最优的窗函数,具体如下:
步骤1.1、hanning窗也称升余弦窗,长度为N的离散hanning窗时域表达式wHn(n)为:
其中,n=0,1,2,…,N-1,wHn(n)的DTFT为WHn(n):
式中,WR(ω)为矩形窗的频谱函数:
步骤1.2、hamming窗也称改进升余弦窗,长度为N的离散hamming窗时域表达式wHm(n)为:
其中,n=0,1,2,…,N-1,wHm(n)的DTFT为WHm(n):
式中,WR(ω)为矩形窗的频谱函数;
步骤1.3、长度为N的离散Blackman窗时域表达式wBm(n)为:
其中,n=0,1,2,…,N-1,wBm(n)的DTFT为WBm(n):
式中,WR(ω)为矩形窗的频谱函数;
步骤1.4、长度为N的离散4项3阶Nuttall窗的时域表达式wNu(n)为:
其中,n=0,1,2,…,N-1,4项最低旁瓣Nuttall窗的DTFT为WNu(n):
式中,WR(ω)为矩形窗的频谱函数;
步骤1.5、根据以下原则选取窗函数:1)能量尽可能集中在主瓣内;2)旁瓣峰值电平尽可能低且衰减快;3)窗函数时域、频域表达式尽可能简洁。
进一步地,步骤2所述的对待测信号进行加窗得到加窗信号,对加窗信号进行FFT变换,得到信号的频谱,具体如下:
步骤2.1、设定待分析信号为:
步骤2.2、对x(t)以采样频率fS、采样长度N进行均匀采样,得到的离散信号为:
式中,n=0,1,…,N-1,N为采样点数;
步骤2.3、用Nuttall窗w(n)对式(11)的信号x(n)进行加窗处理,得到加窗后信号的离散傅里叶变换为X(kΔf):
式中,k为一谐波处的采样点值,Δf为频率分辨率,那么kΔf则为该谐波处的频率;用j表示虚数,i表示谐波次数;
步骤2.4、忽略负频点处频峰的旁瓣影响,则在正频点周边的离散谱线函数为:
式中,离散谱线频率分辨率Δf=fs/N;W(·)为Nuttall窗w(n)的频谱函数,余弦窗函数的频域表达式为W(λ):
其中:
式中,ah为窗函数的各项系数,H为窗函数多项式的项数,h表示多项式第h项,w为频域上的角频率;
设定信号的频点ki所对应幅值最大的谱线为km,km左边的谱线为km-1,km右边的谱线为km+1,其中,km-1=km-1,km+1=km+1,则三条谱线的幅值分别为y1=|X(kmΔf)|、y2=|X(km-1Δf)|、y3=|X(km+1Δf)|,频率分别为v(1)=arg[X(kmΔf)]、v(2)=arg|X(km-1Δf)|、v(3)=|X(km+1Δf)|。
进一步地,步骤3所述的在频谱中搜索谱峰,仅搜索幅值大于设定阈值的谱线,计算每个谱峰两根相邻谱线之间的相位差,具体如下:
步骤3.1、由式(14)和式(15)得:
步骤3.2、得km-1,km-1,km+1处谱线的相位分别为:
arg[X(kmΔf)]=θ-π(km-ki)-π/2 (19)
arg[X(km-1Δf)]=θ-π(km-1-ki)-π/2 (20)
arg[X(km+1Δf)]=θ-π(km+1-ki)-π/2 (21)
由式(19)至式(21)知,km-1,km-1,km+1每相邻两条谱线之间的相位差分别为π,即对于信号加窗后而言,每相邻两条谱线之间的相位差为π,若主瓣存在其他谱线的干扰,则不成立,因此经过加窗信号主瓣的相邻两谱线之间的相位差是否为π来判断该峰值谱线是否有主瓣干扰的情况。
进一步地,步骤4所述的判断主瓣干扰是否存在:若主瓣干扰不存在,则使用插值算法校正对应的谐波或间谐波成分的参数,具体如下:
步骤4.1、设置可动态调整的参量,通过增加一个临界信号来获得临界值;IEC61000-3-6规定当信号的幅值小于基波幅值的0.1%时,忽略此信号,故将残余信号的终止条件б2的值设定为基波幅值的0.1%;在谱线km处构建单频信号x1(t),信号x1(t)的频率为k1Δf,幅值为|X(k1)|*4/N,相位为angle[X(k1)];设定在最大谱线km处存在信号x2(t),信号x2(t)的频率为kmΔf+5,幅值为基波幅值大小的0.1%,相位大小为0,将信号x1(t)与x2(t)叠加,进行加窗插值,获得最大谱线以及次最大谱线相位差作为主瓣谱线是否存在干扰的临界条件б1的取值;
步骤4.2、判断无其他谐波谱线干扰,则用插值法进行计算,具体如下:
设α=ki-km,其中ki为真实频率处的谱线,km为离ki谱线最近的幅值最大谱线,km-1=km-1,km+1=km+1,km-1、km、km+1谱线处的幅值分别为y1、y2、y3,α的取值范围为(-0.5,0.5),得:
记式(22)为β=h(α),其反函数记为α=h-1(β),利用多项式拟合逼近,得α=h-1(β)逼近式为α=H(β),由α求出参数β,则第i次谐波的频率fi修正公式为:
fi=kiΔf=(α+km)Δf (23)
步骤4.3、第i次谐波的幅值修正Ai是对km-1、km、km+1三根谱线的幅值进行加权平均,其计算公式为:
当N值大于阈值时,式(24)简化为:
Ai=N-1(y1+2y2+y3)h(α) (25)
利用多项式拟合逼近,得Ai的逼近式为:
Ai=N-1(y1+2y2+y3)m(α) (26)
式中,m(α)为h(α)的拟合逼近式;
步骤4.4、由式(13)得相位修正公式为:
其中,Nuttall窗函数的三谱线插值7阶逼近式为:
α=h(β)=1.724339β-0.085167λ3+0.014833λ5-0.002882λ7 (28)
m(α)=1.724339+0.430783λ2+0.048030λ4+0.004225λ6 (29)。
进一步地,步骤4所述的若主瓣干扰存在,则不处理当前主瓣,利用插值算法算出的参数构建一个新信号,用原信号减去新信号得到残差信号,用Prony算法分析残差信号,得到含谐波和间谐波的信号的参数值,具体如下:
步骤5.1、采用p个具有任意幅值、相位、频率和衰减因子的指数函数的线性组合对振荡数据进行等间距采样,设N个原始信号为x(0),x(1),…,x(n),有数学模型:
步骤5.2、定义样本函数r(i,j)为:
式中,x(·)为原始信号,x*(·)为其共轭值;
利用样本函数r(i,j)构建矩阵为:
步骤5.3、线性参数估计看做求解方程:
Rc[1 a1 … ap]T=[εP 0 … 0] (33)
对Rc进行奇异值分解,得到样本矩阵的奇异值分布:
0=…σpe=σP+1≤σP≤…≤σ2≤σ1 (34)
式中,σ1、…、σP为Rc的q个非0奇异值,σpe为Rc为0的奇异值;
由于实际中存在噪声,且M>p,因此用噪声空间取代样本矩阵p-M维零空间;p-M为估计阶数与实际阶数的差值;
定义***阶数估计wi为:
式中,p为***的估计阶数,由于wi单调递增,当i值从1向p递增,wi的值向1趋近,当wi大于限值0.995时,此时的i认定为***的实际阶数M;
步骤5.4、用零空间代替***的噪声空间,得到样本矩阵Rc的最优矩阵,预测参数的解的表达式为:
Ai=S-(M)(i+1,1)/S-(M)(1,1) (36)
式中:
步骤5.6、结合步骤4求出的不含主瓣干扰的谐波信号的频率、幅值、相角,计算含谐波和间谐波的信号的频率、幅值、相角。
本发明与现有技术相比,其显著优点在于:(1)分析比较基本常用窗函数的频域特性,选择性能较好的窗函数,对待测信号进行加窗得到加窗信号,对加窗信号进行FFT变换,得到信号的频谱,结构简单;(2)在频谱中搜索谱峰,计算每个谱峰附近两根谱线之间的相位差;若主瓣干扰不存在,则使用插值算法校正对应的谐波或间谐波成分的参数;若主瓣干扰存在,则暂不处理当前主瓣,运算速度快、可靠性高;(3)利用插值算法算出的参数构建一个新信号,用原信号减去新信号得到残差信号,用Prony算法分析残差信号,准确快速的得到频率相近的谐波和间谐波的参数,从而有效的分析信号中谐波与间谐波的参数,对于抑制间谐波进而改善提高电能质量具有重要的现实意义。
附图说明
图1是本发明加窗插值和Prony算法的谐波与间谐波分析方法的流程示意图。
图2是本发明实施例中部分常用窗函数的归一化频率特性图。
图3是本发明实施例中含频率相近的谐波与间谐波信号的幅频谱图。
具体实施方式
下面结合附图及具体实施例对本发明作进一步详细说明。
结合图1,本发明基于加窗插值和Prony算法的谐波与间谐波分析方法,包括以下步骤:
步骤1、分析比较基本常用窗函数的频域特性,选择性能较好的窗函数,结合图2,具体如下:
步骤1.1、hanning窗也称升余弦窗,长度为N的离散hanning窗时域表达式wHn(n)为:
其中,n=0,1,2,…,N-1,wHn(n)的DTFT为WHn(n):
式中,WR(ω)为矩形窗的频谱函数:
步骤1.2、hamming窗也称改进升余弦窗,长度为N的离散hamming窗时域表达式wHm(n)为:
其中,n=0,1,2,…,N-1,wHm(n)的DTFT为WHm(n):
式中,WR(ω)为矩形窗的频谱函数;
步骤1.3、长度为N的离散Blackman窗时域表达式wBm(n)为:
其中,n=0,1,2,…,N-1,wBm(n)的DTFT为WBm(n):
式中,WR(ω)为矩形窗的频谱函数;
步骤1.4、长度为N的离散4项3阶Nuttall窗的时域表达式wNu(n)为:
其中,n=0,1,2,…,N-1,4项最低旁瓣Nuttall窗的DTFT为WNu(n):
式中,WR(ω)为矩形窗的频谱函数;
步骤1.5、根据以下原则选取窗函数:1)能量尽可能集中在主瓣内;2)旁瓣峰值电平尽可能低且衰减快;3)窗函数时域、频域表达式尽可能简洁。
但窗函数无法同时具备主瓣窄、旁瓣峰值电平低且旁瓣衰减快的条件。因此,在采用FFT谐波检测算法时要综合考虑实时性和准确度两方面的要求。如果仅要求精确计算频率,而不考虑幅值、初相角等参数的准确度,则可选用主瓣宽度比较窄而便于分辨的矩形窗,例如测量物体的自振频率等;如果分析窄带信号,且有较强的干扰噪声,则应选用旁瓣幅度小的窗函数,如Hanning窗、三角窗等。本发明分析信号中谐波与间谐波的参数,综合考虑选择Nuttall窗。
步骤2、对待测信号进行加窗得到加窗信号,对加窗信号进行FFT变换,得到信号的频谱,具体如下:
步骤2.1、设定待分析信号为:
步骤2.2、对x(t)以采样频率fS、采样长度N进行均匀采样,得到的离散信号为:
式中,n=0,1,…,N-1,N为采样点数;
步骤2.3、用Nuttall窗w(n)对式(11)的信号x(n)进行加窗处理,得到加窗后信号的离散傅里叶变换为X(kΔf):
式中,k为某一谐波处的采样点值,Δf为频率分辨率,那么kΔf则为该谐波处的频率。用j表示虚数,i表示谐波次数。
步骤2.4、忽略负频点处频峰的旁瓣影响,则在正频点周边的离散谱线函数为:
式中,离散谱线频率分辨率Δf=fs/N;W(·)为Nuttall窗w(n)的频谱函数,余弦窗函数的频域表达式为W(λ):
其中:
式中,ah为窗函数的各项系数,H为窗函数多项式的项数,h表示多项式第h项,w为频域上的角频率。
设定信号的频点ki所对应幅值最大的谱线为km,km左边的谱线为km-1,km右边的谱线为km+1,其中,km-1=km-1,km+1=km+1,则三条谱线的幅值分别为y1=|X(kmΔf)|、y2=|X(km-1Δf)|、y3=|X(km+1Δf)|,频率分别为v(1)=arg[X(kmΔf)]、v(2)=arg|X(km-1Δf)|、v(3)=|X(km+1Δf)|。
步骤3、在频谱中搜索谱峰,仅搜索幅值大于设定阈值的谱线,计算每个谱峰两根相邻谱线之间的相位差,具体如下:
步骤3.1、由式(14)和式(15)得:
步骤3.2、得km-1,km-1,km+1处谱线的相位分别为:
arg[X(kmΔf)]=θ-π(km-ki)-π/2 (19)
arg[X(km-1Δf)]=θ-π(km-1-ki)-π/2 (20)
arg[X(km+1Δf)]=θ-π(km+1-ki)-π/2 (21)
由式(19)至式(21)知,km-1,km-1,km+1每相邻两条谱线之间的相位差分别为π,即对于信号加窗后而言,每相邻两条谱线之间的相位差为π,若主瓣存在其他谱线的干扰,则不成立,因此经过加窗信号主瓣的相邻两谱线之间的相位差是否为π来判断该峰值谱线是否有主瓣干扰的情况。
步骤4、判断主瓣干扰是否存在:
(1)若主瓣干扰不存在,则使用插值算法校正对应的谐波或间谐波成分的参数,具体如下:
步骤4.1、设置可动态调整的参量,通过增加一个临界信号来获得临界值;IEC61000-3-6规定当信号的幅值小于基波幅值的0.1%时,忽略此信号,故将残余信号的终止条件б2的值设定为基波幅值的0.1%;在谱线km处构建单频信号x1(t),信号x1(t)的频率为k1Δf,幅值为|X(k1)|*4/N,相位为angle[X(k1)];设定在最大谱线km处存在信号x2(t),信号x2(t)的频率为kmΔf+5,幅值为基波幅值大小的0.1%,相位大小为0,将信号x1(t)与x2(t)叠加,进行加窗插值,获得最大谱线以及次最大谱线相位差作为主瓣谱线是否存在干扰的临界条件б1的取值;
步骤4.2、判断无其他谐波谱线干扰,则用插值法进行计算,具体如下:
设α=ki-km,其中ki为真实频率处的谱线,km为离ki谱线最近的幅值最大谱线,km-1=km-1,km+1=km+1,km-1、km、km+1谱线处的幅值分别为y1、y2、y3,α的取值范围为(-0.5,0.5),得:
记式(22)为β=h(α),其反函数记为α=h-1(β),利用多项式拟合逼近,可得α=h-1(β)逼近式为α=H(β),由α可求出参数β,则第i次谐波的频率fi修正公式为:
fi=kiΔf=(α+km)Δf (23)
步骤4.3、第i次谐波的幅值修正Ai是对km-1、km、km+1三根谱线的幅值进行加权平均,其计算公式为:
当N值较大时,式(24)可简化为:
Ai=N-1(y1+2y2+y3)h(α) (25)
利用多项式拟合逼近,可得Ai的逼近式为:
Ai=N-1(y1+2y2+y3)m(α) (26)
式中,m(α)为h(α)的拟合逼近式。
步骤4.4、由式(13)得相位修正公式为:
其中,Nuttall窗函数的三谱线插值7阶逼近式为:
α=h(β)=1.724339β-0.085167λ3+0.014833λ5-0.002882λ7 (28)
m(α)=1.724339+0.430783λ2+0.048030λ4+0.004225λ6 (29)
(2)若主瓣干扰存在,则不处理当前主瓣,利用插值算法算出的参数构建一个新信号,用原信号减去新信号得到残差信号,用Prony算法分析残差信号,得到含谐波和间谐波的信号的参数值,如图3所示,具体如下:
步骤5.1、采用p个具有任意幅值、相位、频率和衰减因子的指数函数的线性组合对振荡数据进行等间距采样,设N个原始信号为x(0),x(1),…,x(n),有数学模型:
步骤5.2、定义样本函数r(i,j)为:
式中,x(·)为原始信号,x*(·)为其共轭值。
利用样本函数r(i,j)构建矩阵为:
步骤5.3、线性参数估计看做求解方程:
Rc[1 a1 … ap]T=[εP 0 … 0] (33)
对Rc进行奇异值分解,得到样本矩阵的奇异值分布:
0=…σpe=σP+1≤σP≤…≤σ2≤σ1 (34)
式中,σ1、…、σP为Rc的q个非0奇异值,σpe为Rc为0的奇异值。
由于实际中存在噪声,且M>p,因此用噪声空间取代样本矩阵p-M维零空间;p-M为估计阶数与实际阶数的差值。
定义***阶数估计wi为:
式中,p为***的估计阶数,由于wi单调递增,当i值从1向p递增,wi的值向1趋近,当wi大于限值0.995时,此时的i认定为***的实际阶数M;
步骤5.4、用零空间代替***的噪声空间,得到样本矩阵Rc的最优矩阵,预测参数的解的表达式为:
Ai=S-(M)(i+1,1)/S-(M)(1,1) (36)
式中:
步骤5.6、结合步骤4求出的不含主瓣干扰的谐波信号的频率、幅值、相角,至此可以全部计算含谐波和间谐波的信号的频率、幅值、相角。
综上所述,本发明线路简单,硬件成本低,可靠性高,能够有效分析出信号中谐波与间谐波的参数。
Claims (6)
1.一种基于加窗插值和Prony算法的谐波与间谐波分析方法,其特征在于,包括以下步骤:
步骤1、分析比较基本窗函数的频域特性,选择性能最优的窗函数;
步骤2、对待测信号进行加窗得到加窗信号,对加窗信号进行FFT变换,得到信号的频谱;
步骤3、在频谱中搜索谱峰,仅搜索幅值大于设定阈值的谱线,计算每个谱峰两根相邻谱线之间的相位差;
步骤4、判断主瓣干扰是否存在:
若主瓣干扰不存在,则使用插值算法校正对应的谐波或间谐波成分的参数;
若主瓣干扰存在,则不处理当前主瓣,利用插值算法算出的参数构建一个新信号,用原信号减去新信号得到残差信号,用Prony算法分析残差信号,得到含谐波和间谐波的信号的参数值。
2.根据权利要求1所述的基于加窗插值和Prony算法的谐波与间谐波分析方法,其特征在于,步骤1所述的分析比较基本窗函数的频域特性,选择性能最优的窗函数,具体如下:
步骤1.1、hanning窗也称升余弦窗,长度为N的离散hanning窗时域表达式wHn(n)为:
其中,n=0,1,2,…,N-1,wHn(n)的DTFT为WHn(n):
式中,WR(ω)为矩形窗的频谱函数:
步骤1.2、hamming窗也称改进升余弦窗,长度为N的离散hamming窗时域表达式wHm(n)为:
其中,n=0,1,2,…,N-1,wHm(n)的DTFT为WHm(n):
式中,WR(ω)为矩形窗的频谱函数;
步骤1.3、长度为N的离散Blackman窗时域表达式wBm(n)为:
其中,n=0,1,2,…,N-1,wBm(n)的DTFT为WBm(n):
式中,WR(ω)为矩形窗的频谱函数;
步骤1.4、长度为N的离散4项3阶Nuttall窗的时域表达式wNu(n)为:
其中,n=0,1,2,…,N-1,4项最低旁瓣Nuttall窗的DTFT为WNu(n):
式中,WR(ω)为矩形窗的频谱函数;
步骤1.5、根据以下原则选取窗函数:1)能量尽可能集中在主瓣内;2)旁瓣峰值电平尽可能低且衰减快;3)窗函数时域、频域表达式尽可能简洁。
3.根据权利要求1所述的基于加窗插值和Prony算法的谐波与间谐波分析方法,其特征在于,步骤2所述的对待测信号进行加窗得到加窗信号,对加窗信号进行FFT变换,得到信号的频谱,具体如下:
步骤2.1、设定待分析信号为:
步骤2.2、对x(t)以采样频率fS、采样长度N进行均匀采样,得到的离散信号为:
式中,n=0,1,…,N-1,N为采样点数;
步骤2.3、用Nuttall窗w(n)对式(11)的信号x(n)进行加窗处理,得到加窗后信号的离散傅里叶变换为X(kΔf):
式中,k为一谐波处的采样点值,Δf为频率分辨率,那么kΔf则为该谐波处的频率;用j表示虚数,i表示谐波次数;
步骤2.4、忽略负频点处频峰的旁瓣影响,则在正频点周边的离散谱线函数为:
式中,离散谱线频率分辨率Δf=fs/N;W(·)为Nuttall窗w(n)的频谱函数,余弦窗函数的频域表达式为W(λ):
其中:
式中,ah为窗函数的各项系数,H为窗函数多项式的项数,h表示多项式第h项,w为频域上的角频率;
设定信号的频点ki所对应幅值最大的谱线为km,km左边的谱线为km-1,km右边的谱线为km+1,其中,km-1=km-1,km+1=km+1,则三条谱线的幅值分别为y1=|X(kmΔf)|、y2=|X(km-1Δf)|、y3=|X(km+1Δf)|,频率分别为v(1)=arg[X(kmΔf)]、v(2)=arg|X(km-1Δf)|、v(3)=|X(km+1Δf)|。
4.根据权利要求3所述的基于加窗插值和Prony算法的谐波与间谐波分析方法,其特征在于,步骤3所述的在频谱中搜索谱峰,仅搜索幅值大于设定阈值的谱线,计算每个谱峰两根相邻谱线之间的相位差,具体如下:
步骤3.1、由式(14)和式(15)得:
步骤3.2、得km-1,km-1,km+1处谱线的相位分别为:
arg[X(kmΔf)]=θ-π(km-ki)-π/2 (19)
arg[X(km-1Δf)]=θ-π(km-1-ki)-π/2 (20)
arg[X(km+1Δf)]=θ-π(km+1-ki)-π/2 (21)
由式(19)至式(21)知,km-1,km-1,km+1每相邻两条谱线之间的相位差分别为π,即对于信号加窗后而言,每相邻两条谱线之间的相位差为π,若主瓣存在其他谱线的干扰,则不成立,因此经过加窗信号主瓣的相邻两谱线之间的相位差是否为π来判断该峰值谱线是否有主瓣干扰的情况。
5.根据权利要求4所述的基于加窗插值和Prony算法的谐波与间谐波分析方法,其特征在于,步骤4所述的判断主瓣干扰是否存在:若主瓣干扰不存在,则使用插值算法校正对应的谐波或间谐波成分的参数,具体如下:
步骤4.1、设置可动态调整的参量,通过增加一个临界信号来获得临界值;IEC61000-3-6规定当信号的幅值小于基波幅值的0.1%时,忽略此信号,故将残余信号的终止条件б2的值设定为基波幅值的0.1%;在谱线km处构建单频信号x1(t),信号x1(t)的频率为k1Δf,幅值为|X(k1)|*4/N,相位为angle[X(k1)];设定在最大谱线km处存在信号x2(t),信号x2(t)的频率为kmΔf+5,幅值为基波幅值大小的0.1%,相位大小为0,将信号x1(t)与x2(t)叠加,进行加窗插值,获得最大谱线以及次最大谱线相位差作为主瓣谱线是否存在干扰的临界条件б1的取值;
步骤4.2、判断无其他谐波谱线干扰,则用插值法进行计算,具体如下:
设α=ki-km,其中ki为真实频率处的谱线,km为离ki谱线最近的幅值最大谱线,km-1=km-1,km+1=km+1,km-1、km、km+1谱线处的幅值分别为y1、y2、y3,α的取值范围为(-0.5,0.5),得:
记式(22)为β=h(α),其反函数记为α=h-1(β),利用多项式拟合逼近,得α=h-1(β)逼近式为α=H(β),由α求出参数β,则第i次谐波的频率fi修正公式为:
fi=kiΔf=(α+km)Δf (23)
步骤4.3、第i次谐波的幅值修正Ai是对km-1、km、km+1三根谱线的幅值进行加权平均,其计算公式为:
当N值大于阈值时,式(24)简化为:
Ai=N-1(y1+2y2+y3)h(α) (25)
利用多项式拟合逼近,得Ai的逼近式为:
Ai=N-1(y1+2y2+y3)m(α) (26)
式中,m(α)为h(α)的拟合逼近式;
步骤4.4、由式(13)得相位修正公式为:
其中,Nuttall窗函数的三谱线插值7阶逼近式为:
α=h(β)=1.724339β-0.085167λ3+0.014833λ5-0.002882λ7 (28)
m(α)=1.724339+0.430783λ2+0.048030λ4+0.004225λ6 (29)。
6.根据权利要求4所述的基于加窗插值和Prony算法的谐波与间谐波分析方法,其特征在于,步骤4所述的若主瓣干扰存在,则不处理当前主瓣,利用插值算法算出的参数构建一个新信号,用原信号减去新信号得到残差信号,用Prony算法分析残差信号,得到含谐波和间谐波的信号的参数值,具体如下:
步骤5.1、采用p个具有任意幅值、相位、频率和衰减因子的指数函数的线性组合对振荡数据进行等间距采样,设N个原始信号为x(0),x(1),…,x(n),有数学模型:
步骤5.2、定义样本函数r(i,j)为:
式中,x(·)为原始信号,x*(·)为其共轭值;
利用样本函数r(i,j)构建矩阵为:
步骤5.3、线性参数估计看做求解方程:
Rc[1 a1…ap]T=[εP 0…0] (33)
对Rc进行奇异值分解,得到样本矩阵的奇异值分布:
0=…σpe=σP+1≤σP≤…≤σ2≤σ1 (34)
式中,σ1、…、σP为Rc的q个非0奇异值,σpe为Rc为0的奇异值;
由于实际中存在噪声,且M>p,因此用噪声空间取代样本矩阵p-M维零空间;p-M为估计阶数与实际阶数的差值;
定义***阶数估计wi为:
式中,p为***的估计阶数,由于wi单调递增,当i值从1向p递增,wi的值向1趋近,当wi大于限值0.995时,此时的i认定为***的实际阶数M;
步骤5.4、用零空间代替***的噪声空间,得到样本矩阵Rc的最优矩阵,预测参数的解的表达式为:
Ai=S-(M)(i+1,1)/S-(M)(1,1) (36)
式中:
步骤5.6、结合步骤4求出的不含主瓣干扰的谐波信号的频率、幅值、相角,计算含谐波和间谐波的信号的频率、幅值、相角。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201911351006.2A CN113032716B (zh) | 2019-12-24 | 2019-12-24 | 基于加窗插值和Prony算法的谐波与间谐波分析方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201911351006.2A CN113032716B (zh) | 2019-12-24 | 2019-12-24 | 基于加窗插值和Prony算法的谐波与间谐波分析方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN113032716A true CN113032716A (zh) | 2021-06-25 |
CN113032716B CN113032716B (zh) | 2024-06-18 |
Family
ID=76452146
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201911351006.2A Active CN113032716B (zh) | 2019-12-24 | 2019-12-24 | 基于加窗插值和Prony算法的谐波与间谐波分析方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN113032716B (zh) |
Cited By (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN114659791A (zh) * | 2022-02-28 | 2022-06-24 | 广东机电职业技术学院 | 汽轮机故障检测方法、***、装置和存储介质 |
CN115508618A (zh) * | 2022-10-13 | 2022-12-23 | 四川大学 | 一种基于时域Hermite插值的准同步谐波分析装置及方法 |
CN115684719A (zh) * | 2023-01-03 | 2023-02-03 | 国网江西省电力有限公司电力科学研究院 | 一种新能源并网***的宽频带耦合谐波分析方法及*** |
CN115825557A (zh) * | 2022-11-25 | 2023-03-21 | 国网四川省电力公司映秀湾水力发电总厂 | 基于谐波分量置零的广义谐波分析方法、装置及介质 |
CN116735957A (zh) * | 2023-06-07 | 2023-09-12 | 四川大学 | 计及主瓣重叠干扰的近频谐波与间谐波测量方法及*** |
WO2024087237A1 (zh) * | 2022-10-27 | 2024-05-02 | 苏州大学 | 一种电网谐波与间谐波的检测方法 |
Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20160274155A1 (en) * | 2013-12-16 | 2016-09-22 | State Grid Corporation Of China (Sgcc) | Method for acquiring parameters of dynamic signal |
CN108037361A (zh) * | 2017-12-05 | 2018-05-15 | 南京福致通电气自动化有限公司 | 一种基于滑动窗dft的高精度谐波参数估计方法 |
CN109541312A (zh) * | 2019-01-09 | 2019-03-29 | 华北电力大学 | 一种新能源汇集地区次同步谐波检测方法 |
CN109557367A (zh) * | 2018-10-23 | 2019-04-02 | 中国农业大学 | 一种高频分辨率谐波和间谐波Prony方法及装置 |
-
2019
- 2019-12-24 CN CN201911351006.2A patent/CN113032716B/zh active Active
Patent Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20160274155A1 (en) * | 2013-12-16 | 2016-09-22 | State Grid Corporation Of China (Sgcc) | Method for acquiring parameters of dynamic signal |
CN108037361A (zh) * | 2017-12-05 | 2018-05-15 | 南京福致通电气自动化有限公司 | 一种基于滑动窗dft的高精度谐波参数估计方法 |
CN109557367A (zh) * | 2018-10-23 | 2019-04-02 | 中国农业大学 | 一种高频分辨率谐波和间谐波Prony方法及装置 |
CN109541312A (zh) * | 2019-01-09 | 2019-03-29 | 华北电力大学 | 一种新能源汇集地区次同步谐波检测方法 |
Non-Patent Citations (2)
Title |
---|
张绍勇;余志顽;余冬荣;: "改进的间谐波分析法――加窗插值MUSIC法", 南京师范大学学报(工程技术版), no. 04, 20 December 2011 (2011-12-20) * |
裴源;宁薇薇;: "基于Prony改进算法的电力***谐波参数估计", 电力职业技术学刊, no. 02, 15 June 2010 (2010-06-15) * |
Cited By (9)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN114659791A (zh) * | 2022-02-28 | 2022-06-24 | 广东机电职业技术学院 | 汽轮机故障检测方法、***、装置和存储介质 |
CN114659791B (zh) * | 2022-02-28 | 2023-07-04 | 广东机电职业技术学院 | 汽轮机故障检测方法、***、装置和存储介质 |
CN115508618A (zh) * | 2022-10-13 | 2022-12-23 | 四川大学 | 一种基于时域Hermite插值的准同步谐波分析装置及方法 |
CN115508618B (zh) * | 2022-10-13 | 2023-10-03 | 四川大学 | 一种基于时域Hermite插值的准同步谐波分析装置及方法 |
WO2024087237A1 (zh) * | 2022-10-27 | 2024-05-02 | 苏州大学 | 一种电网谐波与间谐波的检测方法 |
CN115825557A (zh) * | 2022-11-25 | 2023-03-21 | 国网四川省电力公司映秀湾水力发电总厂 | 基于谐波分量置零的广义谐波分析方法、装置及介质 |
CN115684719A (zh) * | 2023-01-03 | 2023-02-03 | 国网江西省电力有限公司电力科学研究院 | 一种新能源并网***的宽频带耦合谐波分析方法及*** |
CN116735957A (zh) * | 2023-06-07 | 2023-09-12 | 四川大学 | 计及主瓣重叠干扰的近频谐波与间谐波测量方法及*** |
CN116735957B (zh) * | 2023-06-07 | 2024-02-27 | 四川大学 | 计及主瓣重叠干扰的近频谐波与间谐波测量方法及*** |
Also Published As
Publication number | Publication date |
---|---|
CN113032716B (zh) | 2024-06-18 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN113032716B (zh) | 基于加窗插值和Prony算法的谐波与间谐波分析方法 | |
Zygarlicki et al. | A reduced Prony's method in power-quality analysis—parameters selection | |
Su et al. | Power harmonic and interharmonic detection method in renewable power based on Nuttall double‐window all‐phase FFT algorithm | |
CN107102255B (zh) | 单一adc采集通道动态特性测试方法 | |
CN104897960B (zh) | 基于加窗四谱线插值fft的谐波快速分析方法及*** | |
Wen et al. | Hanning self-convolution window and its application to harmonic analysis | |
Wen et al. | Novel three-point interpolation DFT method for frequency measurement of sine-wave | |
CN106771586B (zh) | 一种直流控制保护板卡的回路信号分析方法及装置 | |
CN113361331B (zh) | 基于加窗插值fft的工频干扰消除方法、***和介质 | |
CN110579684A (zh) | 一种基于融合算法的小电流接地***选线方法 | |
CN114002475B (zh) | 一种避雷器阻性电流在线监测方法 | |
de Almeida Coelho et al. | Power measurement using stockwell transform | |
CN112730982A (zh) | 一种混合直流输电***的谐波检测方法 | |
CN107315714B (zh) | 一种去卷积功率谱估计方法 | |
CN109557367B (zh) | 一种高频分辨率谐波和间谐波Prony方法及装置 | |
CN117388574A (zh) | 基于msd混合卷积窗的高频谐波分析方法、***、设备及存储介质 | |
Rodrigues et al. | Low-cost embedded measurement system for power quality frequency monitoring | |
Tarasiuk | A few remarks about assessment methods of electric power quality on ships–present state and further development | |
Ferrero et al. | Employment of interpolated DFT-based PMU algorithms in three-phase systems | |
CN112946374B (zh) | 基于卷积窗函数的三相不平衡度检测方法及装置 | |
Firouzjah et al. | A predictive current control method for shunt active filter with windowing based wavelet transform in harmonic detection | |
CN113189398A (zh) | 一种置零点频域加窗的高阶谐波分析方法及装置 | |
CN112836390B (zh) | 一种变流器故障检测方法、***及存储介质 | |
CN112034232A (zh) | 一种供电***电压暂降检测方法 | |
Singh et al. | Improvement in taylor weighted least square based dynamic synchrophasor estimation algorithm for real time application |
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 |