1 引言

随着电力电子装置及非线性负荷的广泛应用，越来越多的谐波成分被注入电网，电网被污染的程度急剧上升.因此，在研究减小或消除电力系统谐波的影响时，准确地分析出谐波信号的相关参数显得尤为重要^[1-3].

目前对于谐波信号分析的方法主要有傅里叶变换算法、小波变换、Prony算法等.傅里叶变换无法检测间谐波，仅适用于分析平稳信号^[4-7].小波变换无法计算偶次谐波的参数，没有考虑实际电网的基频波动等^[8]. Prony算法在低信噪比情况下的准确度较低^[9].

在非同步采样时，频谱泄漏和栅栏效应是影响谐波分析准确性的主要原因之一，而选择合适的窗函数，进行加窗插值，可以有效地解决这一问题^[10-13].据此，提出一种将新型窗函数和插值结合的分析方法，对谐波信号进行分析，从而减小频谱泄漏或栅栏效应产生的误差，提高谐波分析的准确性.

理想窗函数应具有窄主瓣，低旁瓣，高旁瓣衰减率的特性.但经典窗函数的特性不能同时满足上述要求，只能在一定程度上改善窗的主瓣或旁瓣性能.而时域卷积运算虽然能大幅提升旁瓣性能，但使窗函数的主瓣变宽.为此，要得到理想窗函数特性，在综合考虑经典窗函数与卷积方法的基础上，设计一种新型窄主瓣、低旁瓣的梯形自卷积窗函数，从而可以分析出对频率分辨率要求更高的谐波成分.

为获得主瓣较窄、旁瓣最优的梯形窗，首先选择合适的梯形窗系数，即梯形窗上底长度与下底长度的比值；为使梯形窗的旁瓣性能得到更加明显的提升，需对该梯形窗进行时域卷积，得到低旁瓣的梯形自卷积窗.

时域连续的梯形窗表达式为

(1)

式中，T为梯形窗的时域长度；a为梯形窗上底长度(0 < a < T).

将式(1) 进行离散化：

(2)

式中，A为上底长度，N为时域离散长度.

对式(2) 进行Z变换后，并将z=e^jω代入得：

(3)

即：

(4)

当N=128，A=20时的梯形窗时域及频域曲线如图 1所示.其频域特性与时域长度N和上底长度A有关.梯形窗的幅频响应函数为

(5)

图 1 N=128，A=20时梯形窗时域及频域曲线 Figure 1 Trapezoidal window curve of time domain and frequency domain when N=128, A=20

图选项

梯形窗的主瓣宽度为

(6)

若A与N两者之间满足比例关系γ_T=A/N，则：

(7)

当时域长度N=128时，比例系数γ_T与梯形窗旁瓣峰值电平的关系如图 2所示.

由表 1可以看出经典余弦窗函数的主瓣宽度B_w、旁瓣峰值电平A_sp(dB)、旁瓣衰减速率D_ea(dB/oct)与参数α有如下近似关系：

(8)

(9)

(10)

参见图 2和表 1，梯形窗的上底长度与窗长度的比值γ_T∈[0.12，0.2]时，梯形窗获得较低的旁瓣峰值电平.当γ_T=0.16时，旁瓣峰值电平达到最低值A_sp=-32.73 dB，此时主瓣宽度为6.89π/N，梯形窗的旁瓣峰值电平接近于Hanning窗，主瓣宽度比Nuttall窗窄，旁瓣衰减速率低于Retangular窗，但并能满足理想窗函数的性能.

图 2 梯形窗旁瓣峰值电平与比例系数对应关系 Figure 2 Corresponding relation between trapezoidal window peak sidelobe level and proportion coefficient

图选项

定义梯形自卷积窗的定义为若干个相同的梯形窗进行时域卷积运算的结果，即：

(11)

式中，p为参与卷积运算的梯形窗的个数.

根据卷积定理，p阶梯形自卷积窗的频谱函数等于参与卷积梯形窗频谱函数的p次幂，即：

(12)

根据式(6)，由p个长度为N的梯形窗构建的p阶梯形自卷积窗的主瓣宽度：

(13)

由式(13) 可见，p阶梯形自卷积窗的主瓣宽度等于参与卷积的梯形窗主瓣宽度的p倍.当p阶梯形自卷积窗的序列长度M为定值时，由M=pN可知，N的取值与卷积阶数成反比，因而固定序列长度M的梯形自卷积窗的主瓣宽度取决于卷积阶数，卷积阶数越高，主瓣将越宽.

p阶梯形自卷积窗函数的旁瓣峰值电平A_sp(dB)、旁瓣衰减速率D_ea(dB/oct)与参数p有如下近似关系：

(14)

(15)

由此可见，梯形自卷积窗的旁瓣衰减速率和旁瓣峰值电平与卷积阶数成正比关系，即梯形卷积窗的旁瓣性能随着卷积阶数的增加而得到提高.

表 2列出1~4阶新型梯形自卷积窗(N=128、γ_T=0.16) 主瓣宽度、旁瓣峰值电平A_sp(dB)和旁瓣衰减速率D_ea(dB/oct).

由于卷积阶数越高，主瓣将越宽.因此，为保证计算速度，在保证采样长度的前提下，选择卷积阶数较低的窗函数.比较式(8)~(10) 和式(13)~(15)，结合表 1和2可以看出，新型2阶梯形自卷积窗比经典余弦窗函数具有较窄主瓣、较低旁瓣旁瓣峰值电平，可以更好地减小频谱泄漏和栅栏效应的影响，因此选择新型2阶梯形自卷积窗作为窗函数对信号进行加窗处理.

表 1 经典窗函数性能参数比较 Table 1 Comparison of performance parameters of cosine combined window

表选项

表 2 新型自卷积窗函数性能参数比较 Table 2 Comparison of performance parameters of new convolution window function

表选项

根据式(12)，新型2阶梯形自卷积窗的频域表达式可表示为

(16)

设x(t)是由基波、谐波组合而成的复杂电网谐波信号，以采样频率f_s对信号进行等间距采样，得到的离散时间信号x(n)为

(17)

式中：M为信号所含整数次频率成分的最高次数；m为谐波次数；A_m和φ_m为第m次谐波的幅值和相位；f₀为基波频率；f_s为采样频率.

首先将采样信号经复调制细化法(zoomed fast Fourier transformation，ZFFT)对基波附近频率间隔较近的谐波成分进行分析，复调制细化法原理如图 3所示.然后将原信号减去该谐波信号，得到剩余信号的表达式为

(18)

图 3 复调制细化法原理流程 Figure 3 ZFFT principle process

图选项

对剩余信号加窗处理，得x_w(n)=x^*(n)w_Tra-2(n)，对加窗后信号x_w(n)进行离散傅里叶变换，忽略负频点-f_m处频谱的旁瓣影响，那么在正频点f_m附近的频谱函数可以表示为

(19)

式中：f_m为第m次谐波频率，W(·)为窗函数频谱的表达式.

对式(19) 进行离散化，得到离散傅里叶变换的表达式为

(20)

式中：△f=f_s/N，f_m=k_m△f，N为数据截断长度.

为了便于分析，本文以第m次谐波进行分析，根据式(20)，令

(21)

根据对称性，在双谱线插值法的基础上，引入两条对称谱线，构成四谱线插值法.在不增加计算量的条件下，从而提高谐波信号分析的准确性.设峰值频率k_m△f附近4条谱线的频率分别为k₁△f、k₂△f、k₃△f、k₄△f，则有k₁ < k₂≤k_m≤k₃ < k₄，其中各值之间的关系为k₂=k₁+1，k₃=k₂+1，k₄=k₃+1.令这4条谱线的频点所对应的幅值分别为y₁、y₂、y₃、y₄，则y₁=|X_w(k₁△f)|，y₂=|X_w(k₂△f)|，y₃=|X_w(k₃△f)|、y₄=|X_w(k₄△f)|.引入参数δ，令δ=k_m-k₂，δ∈(-0.5，0.5)，求取δ的值即可求得频率、幅值和相位的相关数值.

由式(21) 得：

假设：

(22)

令r=2y₂+y₁，s=2y₃+y₄，则：

(23)

其中，y₁、y₂、y₃、y₄的值可以通过传统的傅里叶变换算法(fast Fourier transformation，FFT)获得，即幅值比β也就可以确定了.当N较大时，式(22) 一般可以简化为β=g(δ)，其反函数记为δ=g^-1(β)，当窗函数w(n)为实系数时，其幅频响应W(2πf)为偶函数，所以函数g(·)及其反函数g^-1(·)都是奇函数，因而δ=g^-1(β)只含有奇次项，其表达式为

(24)

式中：c₁，c₃，…，c_2α+1为2α+1次逼近多项式的奇次项系数.

求得δ后，可以求出实际信号的频率为

(25)

为了更加准确地计算出实际信号的幅值和相位，需要将y₁、y₂、y₃、y₄进行幅值修正，通过对这四根谱线的幅值进行加权平均来计算出实际的峰值点的幅值，经过对比试验，这4根谱线加权值分别为1，2，2，1.

从而有，

(26)

对于一般的实系数窗函数，当N较大时，式(26) 可表示为

(27)

其中，u(·)是偶函数，其结果中不含奇次项，则A_m可以表示为

(28)

式中：a₀，a₂，…，a_2l为2l次逼近多项式的偶次幂系数.

信号的相位由式(21) 可得：

(29)

由式(16) 可得：

(30)

(31)

(32)

将式(30)~(32) 代入式(22) 和式(26)，得到δ=g^-1(β)的逼近式和A_m=N^-1(y₂+2y₁+2y₃+y₄)u(δ)的逼近式.为了使δ和u(δ)的值在数据拟合中保持小数点后八位精度，对δ和u(δ)用Matlab多项式逼近拟合函数ployfit进行拟合.文中采用对δ数值拟合时采用修正公式幂次数选为7次，u(δ)的数据拟合采用修正公式幂次数选为6次，其推倒流程如图 4所示，得到δ=g^-1(β)和u(δ)的系数修正公式为

(33)

(34)

图 4 公式推导流程图 Figure 4 Program flow diagram of formulas derivation

图选项

频率修正公式为

(35)

幅值修正公式为

(36)

相位修正公式为

(37)

仿真模型为采用满足电网复杂波形的前11次谐波信号进行仿真验证：

(38)

式中：基波频率f₀为50 Hz，采样频率f_s为5 120 Hz，采样点N取1 024.复解析带通滤波器中心频率f_e为50 Hz，选抽比D为100，滤波器半阶数M为100.本文设置仿真信号各次谐波成分的幅值及相位如表 3所示.

表 3 仿真信号的构成 Table 3 Components of the simulated signal

表选项

本文算法的主要过程如下，首先采用复调制细化法对采样信号频率间隔较近的谐波成份进行分析，然后将剩余信号加新型2阶梯形自卷积窗处理，并采用离散傅里叶变换运算得到离散频谱，找出相应的峰值谱线，利用四谱线插值推导出谐波参数修正公式，最后分析出谐波相关参数.

为验证本文提出算法的准确性，本文针对直接FFT、加Hanning窗、Nuttall窗和与本文提出的新型2阶梯形自卷积窗四谱线插值法进行对比分析.不同方法分析出各次谐波的幅值和相位的相对误差如表 4和5所示图表中采用的算数符号a·E-b表示为a×10^-b，其中a与b为常数.

由图 5~8可以看出，采样信号直接采用傅里叶变换进行分析时频谱泄漏现象较为严重，同时出现了原有信号不存在频率为224.9 Hz虚假谐波成份.通过加Nuttall窗插值法不能准确分析出与基波频率较近的谐波成份，并且基波成份的幅值受临近谐波能量泄露的影响明显增大.采用复调制细化法能较为准确地分析出基波频率附近的谐波成份，加新型2阶梯形自卷积窗使分析的结果受频谱泄漏和栅栏效应的影响较小.

图 5 FFT分析 Figure 5 FFT analysis

图选项

图 6 加Nuttall窗插值分析 Figure 6 Nuttall window interpolation analysis

图选项

图 7 复调制细化法分析 Figure 7 Zoomed-FFT analysis

图选项

图 8 加新型2阶梯形自卷积窗分析 Figure 8 New 2 step self-convolution window analysis

图选项

由表 4和5可以看出采用加新型2阶梯形自卷积窗分析出信号的各次谐波的幅值和相位的值都有了明显的改善，幅值和相位相对误差集中在10^-6左右，与其他方法相比提高了1~2个数量级，总体来说本文提出的方法在进行谐波信号分析时能够获得比较准确的分析结果.

表 4 不同方法谐波信号的幅值相对误差(%)比较 Table 4 Comparisons of harmonic amplitude relative errors in different algorithms

表选项

表 5 不同方法谐波信号的相位相对误差(%)比较 Table 5 Comparisons of harmonic phase relative errors in different algorithms

表选项

为了减小谐波分析的误差，提高谐波分析的准确性，首先采用复调制细化法可以有效地分析出基波附近频率间隔较近的谐波成份，消除受临近谐波的影响，使分析出基波的参数结果更加准确，避免了虚假成分的存在.结合自卷的特性，在插值算法的基础上，提出一种新型2阶梯形自卷积窗四谱线插值的电网谐波分析方法，新型2阶梯形自卷积窗具有主瓣窄、旁瓣峰值电平低、旁瓣衰减速率较快的特点，有效地减小由频谱泄漏和栅栏效应引起的误差，采用四谱线插值法推导出谐波参数的修正公式，使分析出的谐波参数更加准确.