CN113032716A - 基于加窗插值和Prony算法的谐波与间谐波分析方法 - Google Patents

Abstract

本发明公开了一种基于加窗插值和Prony算法的谐波与间谐波分析方法。该方法为：首先比较常用窗函数的频域特性，选择合适的窗函数；然后对待测信号进行加窗得到加窗信号，并进行FFT变换，得到信号的频谱；接着在频谱中搜索谱峰，计算每个谱峰附近两根谱线之间的相位差；最后判断主瓣干扰是否存在，若不存在，则使用插值算法校正对应的谐波或间谐波成分的参数；若存在，则不处理当前主瓣，利用插值算法算出的参数构建一个新信号，用原信号减去新信号得到残差信号，用Prony算法分析残差信号，得到频率相近的谐波和间谐波的参数。本发明线路简单，硬件成本低，可靠性高，能够有效分析出信号中谐波与间谐波的参数。

Description

基于加窗插值和Prony算法的谐波与间谐波分析方法

技术领域

本发明属于电力系统信号分析领域，特别是一种基于加窗插值和Prony算法的谐波与间谐波分析方法。

背景技术

电能是由电力企业向用户提供的一种特殊商品，在理想情况下，公用电网是以恒定的频率、标准的正弦波形、规定的电压水平和对称的三相电压电流向用户供电，由电厂、电网和用电方共同保证电能质量。但在实际电网中，由于系统中的发电机、变压器、输电线路等设备的固有缺陷及非线性负荷的大量使用等原因，使得电力企业提供的电压、电流信号波形和理想波形相比存在着一定程度的偏差，当偏差较大时，会给电网带来严重危害。电压、频率和波形是电能质量的三个因素，如果它们任意一个不在限定的范围内，就会给电力系统和用户带来不同程度的危害。

近年来，随着电力电子设备和冲击性负荷的增多，使得电网中谐波及间谐波的含量大量增加。谐波的存在会降低变压器容量，增加电机损耗，干扰电气设备正常工作，严重时甚至会损毁设备。与传统谐波相比，间谐波可以是直流到高次谐波间的任意频率，因此间谐波不但具有谐波的危害，其特性还会导致电压闪变、继电器误动作及引起无源滤波器过载等多种危害。同时，随着分布式电源的兴起和智能电网的快速发展，间谐波的来源越来越多，对电力系统的安全运行造成严重的威胁，因此对间谐波的治理迫在眉睫。然而，治理间谐波的前提条件是能够准确检测出间谐波参数，只有分析出间谐波信号，才能对其进行抑制和治理，进而改善供电环境，确保供电设备正常运转。因此，准确快速的检测出间谐波，对于抑制间谐波进而改善提高电能质量具有重要的现实意义。

目前文献中一般使用FFT插值算法来分析谐波信号，当信号中不含频率相近的谐波和间谐波信号时，现有FFT方法及其改进算法的具有较高的精度和稳定性，但当待分析信号中含有频率相近的谐波和间谐波成分时，该方法的频率分辨率不能满足需求。而Prony算法可使用一组指数函数的线性组合来拟合待分析信号，在理论上具有无限小的频率分辨率，当电网信号为宽带信号且所含谐波和间谐波成分较多时，算法准确性会受到一定影响，因此将两种算法结合可以弥补两种算法的局限性，同时，能保证精确度的情况下，准确分析出频率相近的间谐波信号。

发明内容

本发明的目的在于提供一种结构简单、运算速度高、可靠性高的基于加窗插值和Prony算法的电力系统谐波与间谐波分析方法。

实现本发明目的的技术解决方案为：一种基于加窗插值和Prony算法的谐波与间谐波分析方法，包括以下步骤：

步骤1、分析比较基本窗函数的频域特性，选择性能最优的窗函数；

步骤2、对待测信号进行加窗得到加窗信号，对加窗信号进行FFT变换，得到信号的频谱；

步骤3、在频谱中搜索谱峰，仅搜索幅值大于设定阈值的谱线，计算每个谱峰两根相邻谱线之间的相位差；

步骤4、判断主瓣干扰是否存在：

若主瓣干扰不存在，则使用插值算法校正对应的谐波或间谐波成分的参数；

若主瓣干扰存在，则不处理当前主瓣，利用插值算法算出的参数构建一个新信号，用原信号减去新信号得到残差信号，用Prony算法分析残差信号，得到含谐波和间谐波的信号的参数值。

进一步地，步骤1所述的分析比较基本窗函数的频域特性，选择性能最优的窗函数，具体如下：

步骤1.1、hanning窗也称升余弦窗，长度为N的离散hanning窗时域表达式w_Hn(n)为：

其中，n＝0,1,2,…,N-1，w_Hn(n)的DTFT为W_Hn(n)：

式中，W_R(ω)为矩形窗的频谱函数：

步骤1.2、hamming窗也称改进升余弦窗，长度为N的离散hamming窗时域表达式w_Hm(n)为：

其中，n＝0,1,2,…,N-1，w_Hm(n)的DTFT为W_Hm(n)：

式中，W_R(ω)为矩形窗的频谱函数；

步骤1.3、长度为N的离散Blackman窗时域表达式w_Bm(n)为：

其中，n＝0,1,2,…,N-1，w_Bm(n)的DTFT为W_Bm(n)：

式中，W_R(ω)为矩形窗的频谱函数；

步骤1.4、长度为N的离散4项3阶Nuttall窗的时域表达式w_Nu(n)为：

其中，n＝0,1,2,…,N-1，4项最低旁瓣Nuttall窗的DTFT为W_Nu(n)：

式中，W_R(ω)为矩形窗的频谱函数；

步骤1.5、根据以下原则选取窗函数：1)能量尽可能集中在主瓣内；2)旁瓣峰值电平尽可能低且衰减快；3)窗函数时域、频域表达式尽可能简洁。

进一步地，步骤2所述的对待测信号进行加窗得到加窗信号，对加窗信号进行FFT变换，得到信号的频谱，具体如下：

步骤2.1、设定待分析信号为：

式中，I为谐波个数，i表示为第i次谐波，i＝1时，f₁为基波频率，f_i、A_i、

为第i次谐波或间谐波的频率、幅值、相位；

步骤2.2、对x(t)以采样频率f_S、采样长度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＝f_s/N；W(·)为Nuttall窗w(n)的频谱函数，余弦窗函数的频域表达式为W(λ)：

其中：

式中，a_h为窗函数的各项系数，H为窗函数多项式的项数，h表示多项式第h项，w为频域上的角频率；

设定信号的频点k_i所对应幅值最大的谱线为k_m，k_m左边的谱线为k_m-1，k_m右边的谱线为k_m+1，其中，k_m-1＝k_m-1，k_m+1＝k_m+1，则三条谱线的幅值分别为y₁＝|X(k_mΔf)|、y₂＝|X(k_m-1Δf)|、y₃＝|X(k_m+1Δf)|，频率分别为v(1)＝arg[X(k_mΔf)]、v(2)＝arg|X(k_m-1Δf)|、v(3)＝|X(k_m+1Δf)|。

进一步地，步骤3所述的在频谱中搜索谱峰，仅搜索幅值大于设定阈值的谱线，计算每个谱峰两根相邻谱线之间的相位差，具体如下：

步骤3.1、由式(14)和式(15)得：

步骤3.2、得k_m-1，k_m-1，k_m+1处谱线的相位分别为：

arg[X(k_mΔf)]＝θ-π(k_m-k_i)-π/2 (19)

arg[X(k_m-1Δf)]＝θ-π(k_m-1-k_i)-π/2 (20)

arg[X(k_m+1Δf)]＝θ-π(k_m+1-k_i)-π/2 (21)

由式(19)至式(21)知，k_m-1，k_m-1，k_m+1每相邻两条谱线之间的相位差分别为π，即对于信号加窗后而言，每相邻两条谱线之间的相位差为π，若主瓣存在其他谱线的干扰，则不成立，因此经过加窗信号主瓣的相邻两谱线之间的相位差是否为π来判断该峰值谱线是否有主瓣干扰的情况。

进一步地，步骤4所述的判断主瓣干扰是否存在：若主瓣干扰不存在，则使用插值算法校正对应的谐波或间谐波成分的参数，具体如下：

步骤4.1、设置可动态调整的参量，通过增加一个临界信号来获得临界值；IEC61000-3-6规定当信号的幅值小于基波幅值的0.1％时，忽略此信号，故将残余信号的终止条件б₂的值设定为基波幅值的0.1％；在谱线k_m处构建单频信号x₁(t)，信号x₁(t)的频率为k₁Δf，幅值为|X(k₁)|*4/N，相位为angle[X(k₁)]；设定在最大谱线k_m处存在信号x₂(t)，信号x₂(t)的频率为k_mΔf+5，幅值为基波幅值大小的0.1％，相位大小为0，将信号x₁(t)与x₂(t)叠加，进行加窗插值，获得最大谱线以及次最大谱线相位差作为主瓣谱线是否存在干扰的临界条件б₁的取值；

步骤4.2、判断无其他谐波谱线干扰，则用插值法进行计算，具体如下：

设

α＝k_i-k_m，其中k_i为真实频率处的谱线，k_m为离k_i谱线最近的幅值最大谱线，k_m-1＝k_m-1，k_m+1＝k_m+1，k_m-1、k_m、k_m+1谱线处的幅值分别为y₁、y₂、y₃，α的取值范围为(-0.5,0.5)，得：

记式(22)为β＝h(α)，其反函数记为α＝h^-1(β)，利用多项式拟合逼近，得α＝h^-1(β)逼近式为α＝H(β)，由α求出参数β，则第i次谐波的频率f_i修正公式为：

f_i＝k_iΔf＝(α+k_m)Δf (23)

步骤4.3、第i次谐波的幅值修正A_i是对k_m-1、k_m、k_m+1三根谱线的幅值进行加权平均,其计算公式为：

当N值大于阈值时，式(24)简化为：

A_i＝N^-1(y₁+2y₂+y₃)h(α) (25)

利用多项式拟合逼近，得A_i的逼近式为：

A_i＝N^-1(y₁+2y₂+y₃)m(α) (26)

式中，m(α)为h(α)的拟合逼近式；

步骤4.4、由式(13)得相位修正公式为：

其中，Nuttall窗函数的三谱线插值7阶逼近式为：

α＝h(β)＝1.724339β-0.085167λ³+0.014833λ⁵-0.002882λ⁷ (28)

m(α)＝1.724339+0.430783λ²+0.048030λ⁴+0.004225λ⁶ (29)。

进一步地，步骤4所述的若主瓣干扰存在，则不处理当前主瓣，利用插值算法算出的参数构建一个新信号，用原信号减去新信号得到残差信号，用Prony算法分析残差信号，得到含谐波和间谐波的信号的参数值，具体如下：

步骤5.1、采用p个具有任意幅值、相位、频率和衰减因子的指数函数的线性组合对振荡数据进行等间距采样，设N个原始信号为x(0),x(1),…,x(n),有数学模型：

式中，n＝0,1,…,N-1，A_i是振幅，f_i是频率，θ_i是初相，α_i是衰减因子，Δt是采样间隔，p为系统的估计阶数，

为原始信号的x(n)的近似估计值；

步骤5.2、定义样本函数r(i,j)为：

式中，x(·)为原始信号，x*(·)为其共轭值；

利用样本函数r(i,j)构建矩阵为：

步骤5.3、线性参数估计看做求解方程：

R_c[1 a₁ … a_p]^T＝[ε_P 0 … 0] (33)

式中，a_i为求解方程而定义的特征多项式，其中i＝1,2,…,p，

当ε_P趋于无穷小时，

就能近似表示为x(n)；

对R_c进行奇异值分解，得到样本矩阵的奇异值分布：

0＝…σ_pe＝σ_P+1≤σ_P≤…≤σ₂≤σ₁ (34)

式中，σ₁、…、σ_P为R_c的q个非0奇异值，σ_pe为R_c为0的奇异值；

由于实际中存在噪声，且M>p，因此用噪声空间取代样本矩阵p-M维零空间；p-M为估计阶数与实际阶数的差值；

定义系统阶数估计w_i为：

式中，p为系统的估计阶数，由于w_i单调递增，当i值从1向p递增，w_i的值向1趋近，当w_i大于限值0.995时，此时的i认定为系统的实际阶数M；

步骤5.4、用零空间代替系统的噪声空间，得到样本矩阵R_c的最优矩阵，预测参数的解的表达式为：

A_i＝S^-(M)(i+1,1)/S^-(M)(1,1) (36)

式中：

其中，

是酉矩阵第j列向量的加窗函数，S^-(M)是S^(M)的逆矩阵；

步骤5.5、由式(36)求解预测参数组成的多项式得到z_i，通过递推公式推导求出

进而根据式(30)求出b_i，将z_i和b_i带人式(38)计算幅值A_i、频率f_i和初相角θ_i，完成谐波分量特征参数的提取：

步骤5.6、结合步骤4求出的不含主瓣干扰的谐波信号的频率、幅值、相角，计算含谐波和间谐波的信号的频率、幅值、相角。

本发明与现有技术相比，其显著优点在于：(1)分析比较基本常用窗函数的频域特性，选择性能较好的窗函数，对待测信号进行加窗得到加窗信号，对加窗信号进行FFT变换，得到信号的频谱，结构简单；(2)在频谱中搜索谱峰，计算每个谱峰附近两根谱线之间的相位差；若主瓣干扰不存在，则使用插值算法校正对应的谐波或间谐波成分的参数；若主瓣干扰存在，则暂不处理当前主瓣，运算速度快、可靠性高；(3)利用插值算法算出的参数构建一个新信号，用原信号减去新信号得到残差信号，用Prony算法分析残差信号，准确快速的得到频率相近的谐波和间谐波的参数，从而有效的分析信号中谐波与间谐波的参数，对于抑制间谐波进而改善提高电能质量具有重要的现实意义。

附图说明

图1是本发明加窗插值和Prony算法的谐波与间谐波分析方法的流程示意图。

图2是本发明实施例中部分常用窗函数的归一化频率特性图。

图3是本发明实施例中含频率相近的谐波与间谐波信号的幅频谱图。

具体实施方式

下面结合附图及具体实施例对本发明作进一步详细说明。

结合图1，本发明基于加窗插值和Prony算法的谐波与间谐波分析方法，包括以下步骤：

步骤1、分析比较基本常用窗函数的频域特性，选择性能较好的窗函数，结合图2，具体如下：

步骤1.1、hanning窗也称升余弦窗，长度为N的离散hanning窗时域表达式w_Hn(n)为：

其中，n＝0,1,2,…,N-1，w_Hn(n)的DTFT为W_Hn(n)：

式中，W_R(ω)为矩形窗的频谱函数：

步骤1.2、hamming窗也称改进升余弦窗，长度为N的离散hamming窗时域表达式w_Hm(n)为：

其中，n＝0,1,2,…,N-1，w_Hm(n)的DTFT为W_Hm(n)：

式中，W_R(ω)为矩形窗的频谱函数；

步骤1.3、长度为N的离散Blackman窗时域表达式w_Bm(n)为：

其中，n＝0,1,2,…,N-1，w_Bm(n)的DTFT为W_Bm(n)：

式中，W_R(ω)为矩形窗的频谱函数；

步骤1.4、长度为N的离散4项3阶Nuttall窗的时域表达式w_Nu(n)为：

其中，n＝0,1,2,…,N-1，4项最低旁瓣Nuttall窗的DTFT为W_Nu(n)：

式中，W_R(ω)为矩形窗的频谱函数；

步骤1.5、根据以下原则选取窗函数：1)能量尽可能集中在主瓣内；2)旁瓣峰值电平尽可能低且衰减快；3)窗函数时域、频域表达式尽可能简洁。

但窗函数无法同时具备主瓣窄、旁瓣峰值电平低且旁瓣衰减快的条件。因此，在采用FFT谐波检测算法时要综合考虑实时性和准确度两方面的要求。如果仅要求精确计算频率，而不考虑幅值、初相角等参数的准确度，则可选用主瓣宽度比较窄而便于分辨的矩形窗，例如测量物体的自振频率等；如果分析窄带信号，且有较强的干扰噪声，则应选用旁瓣幅度小的窗函数，如Hanning窗、三角窗等。本发明分析信号中谐波与间谐波的参数，综合考虑选择Nuttall窗。

步骤2、对待测信号进行加窗得到加窗信号，对加窗信号进行FFT变换，得到信号的频谱，具体如下：

步骤2.1、设定待分析信号为：

式中，I为谐波个数，i表示为第i次谐波，i＝1时，f₁为基波频率，f_i、A_i、

为第i次谐波或间谐波的频率、幅值、相位；

步骤2.2、对x(t)以采样频率f_S、采样长度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＝f_s/N；W(·)为Nuttall窗w(n)的频谱函数，余弦窗函数的频域表达式为W(λ)：

其中：

式中，a_h为窗函数的各项系数，H为窗函数多项式的项数，h表示多项式第h项，w为频域上的角频率。

设定信号的频点k_i所对应幅值最大的谱线为k_m，k_m左边的谱线为k_m-1，k_m右边的谱线为k_m+1，其中，k_m-1＝k_m-1，k_m+1＝k_m+1，则三条谱线的幅值分别为y₁＝|X(k_mΔf)|、y₂＝|X(k_m-1Δf)|、y₃＝|X(k_m+1Δf)|，频率分别为v(1)＝arg[X(k_mΔf)]、v(2)＝arg|X(k_m-1Δf)|、v(3)＝|X(k_m+1Δf)|。

步骤3、在频谱中搜索谱峰，仅搜索幅值大于设定阈值的谱线，计算每个谱峰两根相邻谱线之间的相位差，具体如下：

步骤3.1、由式(14)和式(15)得：

步骤3.2、得k_m-1，k_m-1，k_m+1处谱线的相位分别为：

arg[X(k_mΔf)]＝θ-π(k_m-k_i)-π/2 (19)

arg[X(k_m-1Δf)]＝θ-π(k_m-1-k_i)-π/2 (20)

arg[X(k_m+1Δf)]＝θ-π(k_m+1-k_i)-π/2 (21)

由式(19)至式(21)知，k_m-1，k_m-1，k_m+1每相邻两条谱线之间的相位差分别为π，即对于信号加窗后而言，每相邻两条谱线之间的相位差为π，若主瓣存在其他谱线的干扰，则不成立，因此经过加窗信号主瓣的相邻两谱线之间的相位差是否为π来判断该峰值谱线是否有主瓣干扰的情况。

步骤4、判断主瓣干扰是否存在：

(1)若主瓣干扰不存在，则使用插值算法校正对应的谐波或间谐波成分的参数，具体如下：

步骤4.1、设置可动态调整的参量，通过增加一个临界信号来获得临界值；IEC61000-3-6规定当信号的幅值小于基波幅值的0.1％时，忽略此信号，故将残余信号的终止条件б₂的值设定为基波幅值的0.1％；在谱线k_m处构建单频信号x₁(t)，信号x₁(t)的频率为k₁Δf，幅值为|X(k₁)|*4/N，相位为angle[X(k₁)]；设定在最大谱线k_m处存在信号x₂(t)，信号x₂(t)的频率为k_mΔf+5，幅值为基波幅值大小的0.1％，相位大小为0，将信号x₁(t)与x₂(t)叠加，进行加窗插值，获得最大谱线以及次最大谱线相位差作为主瓣谱线是否存在干扰的临界条件б₁的取值；

步骤4.2、判断无其他谐波谱线干扰，则用插值法进行计算，具体如下：

设

α＝k_i-k_m，其中k_i为真实频率处的谱线，k_m为离k_i谱线最近的幅值最大谱线，k_m-1＝k_m-1，k_m+1＝k_m+1，k_m-1、k_m、k_m+1谱线处的幅值分别为y₁、y₂、y₃，α的取值范围为(-0.5,0.5)，得：

记式(22)为β＝h(α)，其反函数记为α＝h^-1(β)，利用多项式拟合逼近，可得α＝h^-1(β)逼近式为α＝H(β)，由α可求出参数β，则第i次谐波的频率f_i修正公式为：

f_i＝k_iΔf＝(α+k_m)Δf (23)

步骤4.3、第i次谐波的幅值修正A_i是对k_m-1、k_m、k_m+1三根谱线的幅值进行加权平均,其计算公式为：

当N值较大时，式(24)可简化为：

A_i＝N^-1(y₁+2y₂+y₃)h(α) (25)

利用多项式拟合逼近，可得A_i的逼近式为：

A_i＝N^-1(y₁+2y₂+y₃)m(α) (26)

式中，m(α)为h(α)的拟合逼近式。

步骤4.4、由式(13)得相位修正公式为：

其中，Nuttall窗函数的三谱线插值7阶逼近式为：

α＝h(β)＝1.724339β-0.085167λ³+0.014833λ⁵-0.002882λ⁷ (28)

m(α)＝1.724339+0.430783λ²+0.048030λ⁴+0.004225λ⁶ (29)

(2)若主瓣干扰存在，则不处理当前主瓣，利用插值算法算出的参数构建一个新信号，用原信号减去新信号得到残差信号，用Prony算法分析残差信号，得到含谐波和间谐波的信号的参数值，如图3所示，具体如下：

步骤5.1、采用p个具有任意幅值、相位、频率和衰减因子的指数函数的线性组合对振荡数据进行等间距采样，设N个原始信号为x(0),x(1),…,x(n)，有数学模型：

式中，n＝0,1,…,N-1，A_i是振幅，f_i是频率，θ_i是初相，α_i是衰减因子，Δt是采样间隔，p为系统的估计阶数，

为原始信号的x(n)的近似估计值；

步骤5.2、定义样本函数r(i,j)为：

式中，x(·)为原始信号，x*(·)为其共轭值。

利用样本函数r(i,j)构建矩阵为：

步骤5.3、线性参数估计看做求解方程：

R_c[1 a₁ … a_p]^T＝[ε_P 0 … 0] (33)

式中，a_i为求解方程而定义的特征多项式，其中i＝1,2,…,p，

当ε_P趋于无穷小时，

就能近似表示为x(n)。

对R_c进行奇异值分解，得到样本矩阵的奇异值分布：

0＝…σ_pe＝σ_P+1≤σ_P≤…≤σ₂≤σ₁ (34)

式中，σ₁、…、σ_P为R_c的q个非0奇异值，σ_pe为R_c为0的奇异值。

由于实际中存在噪声，且M>p，因此用噪声空间取代样本矩阵p-M维零空间；p-M为估计阶数与实际阶数的差值。

定义系统阶数估计w_i为：

式中，p为系统的估计阶数，由于w_i单调递增，当i值从1向p递增，w_i的值向1趋近，当w_i大于限值0.995时，此时的i认定为系统的实际阶数M；

步骤5.4、用零空间代替系统的噪声空间，得到样本矩阵R_c的最优矩阵，预测参数的解的表达式为：

A_i＝S^-(M)(i+1,1)/S^-(M)(1,1) (36)

式中：

其中，

是酉矩阵第j列向量的加窗函数，S^-(M)是S^(M)的逆矩阵；

步骤5.5、由式(36)求解预测参数组成的多项式得到z_i，通过递推公式推导求出

进而根据式(30)求出b_i，将z_i和b_i带人式(38)计算幅值A_i、频率f_i和初相角θ_i，完成谐波分量特征参数的提取：

步骤5.6、结合步骤4求出的不含主瓣干扰的谐波信号的频率、幅值、相角，至此可以全部计算含谐波和间谐波的信号的频率、幅值、相角。

综上所述，本发明线路简单，硬件成本低，可靠性高，能够有效分析出信号中谐波与间谐波的参数。

Claims (6)

1.一种基于加窗插值和Prony算法的谐波与间谐波分析方法，其特征在于，包括以下步骤：

步骤1、分析比较基本窗函数的频域特性，选择性能最优的窗函数；

步骤2、对待测信号进行加窗得到加窗信号，对加窗信号进行FFT变换，得到信号的频谱；

步骤3、在频谱中搜索谱峰，仅搜索幅值大于设定阈值的谱线，计算每个谱峰两根相邻谱线之间的相位差；

步骤4、判断主瓣干扰是否存在：

若主瓣干扰不存在，则使用插值算法校正对应的谐波或间谐波成分的参数；

若主瓣干扰存在，则不处理当前主瓣，利用插值算法算出的参数构建一个新信号，用原信号减去新信号得到残差信号，用Prony算法分析残差信号，得到含谐波和间谐波的信号的参数值。

2.根据权利要求1所述的基于加窗插值和Prony算法的谐波与间谐波分析方法，其特征在于，步骤1所述的分析比较基本窗函数的频域特性，选择性能最优的窗函数，具体如下：

步骤1.1、hanning窗也称升余弦窗，长度为N的离散hanning窗时域表达式w_Hn(n)为：

其中，n＝0,1,2,…,N-1，w_Hn(n)的DTFT为W_Hn(n)：

式中，W_R(ω)为矩形窗的频谱函数：

步骤1.2、hamming窗也称改进升余弦窗，长度为N的离散hamming窗时域表达式w_Hm(n)为：

其中，n＝0,1,2,…,N-1，w_Hm(n)的DTFT为W_Hm(n)：

式中，W_R(ω)为矩形窗的频谱函数；

步骤1.3、长度为N的离散Blackman窗时域表达式w_Bm(n)为：

其中，n＝0,1,2,…,N-1，w_Bm(n)的DTFT为W_Bm(n)：

式中，W_R(ω)为矩形窗的频谱函数；

步骤1.4、长度为N的离散4项3阶Nuttall窗的时域表达式w_Nu(n)为：

其中，n＝0,1,2,…,N-1，4项最低旁瓣Nuttall窗的DTFT为W_Nu(n)：

式中，W_R(ω)为矩形窗的频谱函数；

步骤1.5、根据以下原则选取窗函数：1)能量尽可能集中在主瓣内；2)旁瓣峰值电平尽可能低且衰减快；3)窗函数时域、频域表达式尽可能简洁。

3.根据权利要求1所述的基于加窗插值和Prony算法的谐波与间谐波分析方法，其特征在于，步骤2所述的对待测信号进行加窗得到加窗信号，对加窗信号进行FFT变换，得到信号的频谱，具体如下：

步骤2.1、设定待分析信号为：

式中，I为谐波个数，i表示为第i次谐波，i＝1时，f₁为基波频率，f_i、A_i、

为第i次谐波或间谐波的频率、幅值、相位；

步骤2.2、对x(t)以采样频率f_S、采样长度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＝f_s/N；W(·)为Nuttall窗w(n)的频谱函数，余弦窗函数的频域表达式为W(λ)：

其中：

式中，a_h为窗函数的各项系数，H为窗函数多项式的项数，h表示多项式第h项，w为频域上的角频率；

设定信号的频点k_i所对应幅值最大的谱线为k_m，k_m左边的谱线为k_m-1，k_m右边的谱线为k_m+1，其中，k_m-1＝k_m-1，k_m+1＝k_m+1，则三条谱线的幅值分别为y₁＝|X(k_mΔf)|、y₂＝|X(k_m-1Δf)|、y₃＝|X(k_m+1Δf)|，频率分别为v(1)＝arg[X(k_mΔf)]、v(2)＝arg|X(k_m-1Δf)|、v(3)＝|X(k_m+1Δf)|。

4.根据权利要求3所述的基于加窗插值和Prony算法的谐波与间谐波分析方法，其特征在于，步骤3所述的在频谱中搜索谱峰，仅搜索幅值大于设定阈值的谱线，计算每个谱峰两根相邻谱线之间的相位差，具体如下：

步骤3.1、由式(14)和式(15)得：

步骤3.2、得k_m-1，k_m-1，k_m+1处谱线的相位分别为：

arg[X(k_mΔf)]＝θ-π(k_m-k_i)-π/2 (19)

arg[X(k_m-1Δf)]＝θ-π(k_m-1-k_i)-π/2 (20)

arg[X(k_m+1Δf)]＝θ-π(k_m+1-k_i)-π/2 (21)

由式(19)至式(21)知，k_m-1，k_m-1，k_m+1每相邻两条谱线之间的相位差分别为π，即对于信号加窗后而言，每相邻两条谱线之间的相位差为π，若主瓣存在其他谱线的干扰，则不成立，因此经过加窗信号主瓣的相邻两谱线之间的相位差是否为π来判断该峰值谱线是否有主瓣干扰的情况。

5.根据权利要求4所述的基于加窗插值和Prony算法的谐波与间谐波分析方法，其特征在于，步骤4所述的判断主瓣干扰是否存在：若主瓣干扰不存在，则使用插值算法校正对应的谐波或间谐波成分的参数，具体如下：

步骤4.1、设置可动态调整的参量，通过增加一个临界信号来获得临界值；IEC61000-3-6规定当信号的幅值小于基波幅值的0.1％时，忽略此信号，故将残余信号的终止条件б₂的值设定为基波幅值的0.1％；在谱线k_m处构建单频信号x₁(t)，信号x₁(t)的频率为k₁Δf，幅值为|X(k₁)|*4/N，相位为angle[X(k₁)]；设定在最大谱线k_m处存在信号x₂(t)，信号x₂(t)的频率为k_mΔf+5，幅值为基波幅值大小的0.1％，相位大小为0，将信号x₁(t)与x₂(t)叠加，进行加窗插值，获得最大谱线以及次最大谱线相位差作为主瓣谱线是否存在干扰的临界条件б₁的取值；

步骤4.2、判断无其他谐波谱线干扰，则用插值法进行计算，具体如下：

设

α＝k_i-k_m，其中k_i为真实频率处的谱线，k_m为离k_i谱线最近的幅值最大谱线，k_m-1＝k_m-1，k_m+1＝k_m+1，k_m-1、k_m、k_m+1谱线处的幅值分别为y₁、y₂、y₃，α的取值范围为(-0.5,0.5)，得：

记式(22)为β＝h(α)，其反函数记为α＝h^-1(β)，利用多项式拟合逼近，得α＝h^-1(β)逼近式为α＝H(β)，由α求出参数β，则第i次谐波的频率f_i修正公式为：

f_i＝k_iΔf＝(α+k_m)Δf (23)

步骤4.3、第i次谐波的幅值修正A_i是对k_m-1、k_m、k_m+1三根谱线的幅值进行加权平均,其计算公式为：

当N值大于阈值时，式(24)简化为：

A_i＝N^-1(y₁+2y₂+y₃)h(α) (25)

利用多项式拟合逼近，得A_i的逼近式为：

A_i＝N^-1(y₁+2y₂+y₃)m(α) (26)

式中，m(α)为h(α)的拟合逼近式；

步骤4.4、由式(13)得相位修正公式为：

其中，Nuttall窗函数的三谱线插值7阶逼近式为：

α＝h(β)＝1.724339β-0.085167λ³+0.014833λ⁵-0.002882λ⁷ (28)

m(α)＝1.724339+0.430783λ²+0.048030λ⁴+0.004225λ⁶ (29)。

6.根据权利要求4所述的基于加窗插值和Prony算法的谐波与间谐波分析方法，其特征在于，步骤4所述的若主瓣干扰存在，则不处理当前主瓣，利用插值算法算出的参数构建一个新信号，用原信号减去新信号得到残差信号，用Prony算法分析残差信号，得到含谐波和间谐波的信号的参数值，具体如下：

步骤5.1、采用p个具有任意幅值、相位、频率和衰减因子的指数函数的线性组合对振荡数据进行等间距采样，设N个原始信号为x(0),x(1),…,x(n),有数学模型：

式中，n＝0,1,…,N-1，A_i是振幅，f_i是频率，θ_i是初相，α_i是衰减因子，Δt是采样间隔，p为系统的估计阶数，

为原始信号的x(n)的近似估计值；

步骤5.2、定义样本函数r(i,j)为：

式中，x(·)为原始信号，x*(·)为其共轭值；

利用样本函数r(i,j)构建矩阵为：

步骤5.3、线性参数估计看做求解方程：

R_c[1 a₁…a_p]^T＝[ε_P 0…0] (33)

式中，a_i为求解方程而定义的特征多项式，其中i＝1,2,…,p，

当ε_P趋于无穷小时，

就能近似表示为x(n)；

对R_c进行奇异值分解，得到样本矩阵的奇异值分布：

0＝…σ_pe＝σ_P+1≤σ_P≤…≤σ₂≤σ₁ (34)

式中，σ₁、…、σ_P为R_c的q个非0奇异值，σ_pe为R_c为0的奇异值；

由于实际中存在噪声，且M>p，因此用噪声空间取代样本矩阵p-M维零空间；p-M为估计阶数与实际阶数的差值；

定义系统阶数估计w_i为：

式中，p为系统的估计阶数，由于w_i单调递增，当i值从1向p递增，w_i的值向1趋近，当w_i大于限值0.995时，此时的i认定为系统的实际阶数M；

步骤5.4、用零空间代替系统的噪声空间，得到样本矩阵R_c的最优矩阵，预测参数的解的表达式为：

A_i＝S^-(M)(i+1,1)/S^-(M)(1,1) (36)

式中：

其中，

是酉矩阵第j列向量的加窗函数，S^-(M)是S^(M)的逆矩阵；

步骤5.5、由式(36)求解预测参数组成的多项式得到z_i，通过递推公式推导求出

进而根据式(30)求出b_i，将z_i和b_i带人式(38)计算幅值A_i、频率f_i和初相角θ_i，完成谐波分量特征参数的提取：

步骤5.6、结合步骤4求出的不含主瓣干扰的谐波信号的频率、幅值、相角，计算含谐波和间谐波的信号的频率、幅值、相角。

Images