CN109782063A - 一种基于纳托尔自卷积窗三谱线插值fft的动态间谐波分析方法 - Google Patents

Abstract

本发明揭示了一种基于纳托尔自卷积窗三谱线插值FFT的动态间谐波分析方法：步骤1、对动态谐波信号进行采样，得到所需要的时域序列；步骤2、用纳托尔自卷积窗对采样的离散序列进行加窗截断，得到所需要的加窗信号；步骤3、对加窗截断信号做快速傅里叶变换，得到信号的频谱信息；步骤4、搜索最大谱线k_i1以及其左右谱线k_i0和k_i2，采集三条所需要的谱线，并确定所对应的幅值；步骤5、引入变量，通过算法求取插值修正公式；步骤6、用三谱线插值公式对所测量的频率、幅值以及相位进行修正。由于自卷积窗抑制频谱泄漏和栅栏效应的能力强，只需要三谱线插值就可以达到很高的精度，并且大大节约了计算时间，减少了计算量。

Description

一种基于纳托尔自卷积窗三谱线插值FFT的动态间谐波分析 方法

技术领域

本发明涉及动态间谐波分析技术领域。

背景技术

随着各种非线性负荷特别是电力电子设备在电力***中的广泛应用，电网 谐波污染日益严重。现代电力***中，电能质量的五个评价标准中不可或缺的 一个指标--谐波，对其进行无误差检测、滤除或者再利用对电力***潮流计算、 电网用电质量监管、电网安全运行等有十分重要价值，对电网谐波能量再利用、 保持清洁电网环境、保证***运行稳定具有重大意义。

电网对于谐波的主要研究方法包括，快速傅里叶变换(fast Fourier transform，FFT)、神经网络算法、瞬时无功功率理论算法、卡尔曼滤波算法等。 近年来，FFT算法计算量的节省越来越显著,并且其易于嵌入式实现而得到广泛 应用。但是直接采用FFT对信号进行频谱分析很难做到同步采样和整数周期截 断，由于电网运行情况复杂，电网频率出现小范围波动时，离散锁相环技术也 很难实现同步采样，由此产生的频谱泄漏和栅栏效应会导致谐波参数测量出现 误差。

发明内容

本发明所要解决的技术问题是实现一种提高动态微弱间谐波信号检测的精 确性的谐波获取处理方法。

为了实现上述目的，本发明采用的技术方案为：一种基于纳托尔自卷积窗 三谱线插值FFT的动态间谐波分析方法：

步骤1、对动态谐波信号进行采样，得到所需要的时域序列；

步骤2、用纳托尔自卷积窗对采样的离散序列进行加窗截断，得到所需要的 加窗信号；

步骤3、对加窗截断信号做快速傅里叶变换，得到信号的频谱信息；

步骤4、搜索最大谱线k_i1以及其左右谱线k_i0和k_i2，采集三条所需要的谱线， 并确定所对应的幅值；

步骤5、引入变量，通过算法求取插值修正公式；

步骤6、用三谱线插值公式对所测量的频率、幅值以及相位进行修正。

所述步骤1中，时域信号x(t)的表达式为

式中，H为谐波的个数；当h＝1时，f₁、A₁、分别为基波的频率、幅值和 初相位；当h≠1时，r_h表示第h项谐波次数，A_h、分别为第h项谐波的幅值和相 位。

所述步骤2含间谐波分量的信号x(t)以采样频率f_s均匀采样得到的离散时间 信号为

对信号x(n)加纳托尔自卷积窗得到频谱信息

所述步骤3中，对x_w(n)做离散快速傅里叶变换得到

简化后的加自卷积窗FFT变换表达式为

式中，Δf为抽样间隔；k为抽样频率点；N为数据截断长度。

所述步骤4中，只考虑一个抽样频率点，忽略其余抽样频率点的影响，则 有

设峰值频点附近抽样得到的最大谱线为第k_i1条，则它左右侧的谱线也含有 丰富的谱线信息，记为第k_i0条和第k_i2条，其大小关系为k_i0(＝k_i1-1)≤k_i1≤k_i2(＝k_i1+1)， 上述3条谱线的幅值分别记为y₀＝|y_w(k_i0Δf)|，y₁＝|y_w(k_i1Δf)|和y₂＝|y_w(k_i2Δf)|。

所述步骤5中，由于0≤k_i-k_i1≤1，则引入2个变量α和β，令 α＝k_i-k_i1-0.5(α∈[-0.5，0.5])，设β＝(y₂-y₀)/y₁，则有：

则β＝f(α)，式中α为自变量，β为因变量，其反函数为α＝f^-1(β)；

可得p阶纳托尔自卷积窗的DTFT为

在α∈[-0.5,0.5]内取一组α值，得到对应的β，调用函数polyfit(β,α,m)进行多项式反拟合，可得函数α＝F(β)的逼近式，其中，m为拟合逼近多项式的阶数；

由此可得第i项频率的插值修正公式为

f_i＝k_iΔf＝(α+k_i1+0.5)Δf。

所述步骤6中，采用三谱线插值修正公式，三根谱线加权平均计算，得到 幅值修正公式为

由于N一般较大则可化简为

A_i＝N^-p(y₂+2y₁+y₀)g(α)

第i项谐波相位的修正公式为

式中h＝0,1,2。

本发明的优点在于采用具有最小旁瓣和最快衰减速率的窗函数纳托尔自卷 积窗。频谱发生泄漏时，峰值频点附近的旁瓣中也含有丰富的信息，为了进一 步减少误差，以及为了防止间谐波信息被淹没，采用三谱线插值算法。仿真实 验验证了纳托尔自卷积窗三谱线插值算法检测动态微弱间谐波信号的准确性和 及时性。

附图说明

下面对本发明说明书中每幅附图表达的内容作简要说明：

图1为4项5阶纳托尔窗的归一化对数频谱；

图2为4项3阶纳托尔窗的归一化对数频谱；

图3为4项最小旁瓣纳托尔窗的归一化对数频谱；

图4为3项最小旁瓣纳托尔窗的归一化对数频谱；

图5为纳托尔自卷积窗的幅频响应曲线。

具体实施方式

基于纳托尔自卷积窗三谱线插值FFT的动态间谐波分析方法步骤为：

对动态谐波信号进行采样，得到所需要的时域序列。

用纳托尔自卷积窗对采样的离散序列进行加窗截断，得到所需要的加窗信 号。

对加窗截断信号做快速傅里叶变换，得到信号的频谱信息。

搜索最大谱线k_i1以及其左右谱线k_i0和k_i2，采集三条所需要的谱线，并确定 所对应的幅值。

引入变量α和β，通过算法求取插值修正公式。

用三谱线插值公式对所测量的频率、幅值以及相位进行修正。

1、母窗的选择：

自卷积窗可以在母窗的基础上，增加近两倍的序列长度，因为可以检测更 多的有效数据个数，提高频谱分辨率，更便于频谱分析。因此，自卷积母窗的 选择至关重要，因为它的性能直接影响自卷积在动态间谐波中的灵敏性和有效 性。

母窗函数的选择，需要考虑诸多因素，最主要的是主瓣的宽度、旁瓣峰值 电平以及旁瓣衰减率。主瓣宽，则频率分辨率低；旁瓣大，泄漏就越多；旁瓣 衰减越快，对泄漏的抑制能力越强。

纳托尔(Nuttall)窗是一组典型的余弦组合窗,其旁瓣性能良好，其时域表 达形式如式(1)。

式中,n＝0,1,...,N'-1；M为窗函数的相数；N'为窗函数的长度；b_m满足约束条 件见式(2)。

纳托尔窗的频谱函数为式(3)。

式中，W_R(ω)为矩形窗的频谱函数，表达式为式(4)。

将式(4)带入式(3)得纳托尔窗具体运算频谱函数为

表1给出了典型的纳托尔窗函数及其他余弦窗函数。由于***采样很少能 做到同步采样，在非同步采样时，为了增加窗函数对频谱泄漏的抑制能力，所 选窗函数峰值电平越小越好，以及旁瓣衰减速率越大越好。表2给出了纳托尔 窗与经典余弦窗旁瓣特性的比较。当N＝50时，典型纳托尔窗的幅频相应曲线如 图1-4所示。其中，M_mag表示幅值，ω表示角频率。

由图1-4和表2可以看出，4项5阶纳托尔窗最适合作为母窗，对频谱泄漏 抑制能力很强。它的旁瓣衰减率为42dB/oct，旁瓣峰值电平为-60.95dB。

综合以上考虑，自卷积窗的母窗选择4项5阶纳托尔窗。

表1典型纳托尔窗函数：

表2几种余弦组合窗的旁瓣特性：

2、纳托尔自卷积窗

在时域运算时，纳托尔自卷积窗的定义为若干个相同的Nuttall窗进行时 域卷积运算得到的窗函数，即

式中，p为参与运算的母窗个数；W_N(m)为参与自卷积的Nuttall窗；为 纳托尔自卷积窗。

在频域运算时，纳托尔自卷积窗的定义为Nuttall窗的频谱函数的p次幂即 为p阶纳托尔自卷积窗的频谱函数，即

式中，W_N(ω)为Nuttall窗的频谱函数；p为参与运算的母窗个数。

将式(3)代入式(7)，得p阶离散纳托尔自卷积窗的频谱响应为

将式(4)代入式(8)，得具体频谱运算式为

在进行自卷积仿真时需要注意的问题有2个：

(1)补零操作。窗函数进行自卷积时，序列长度增加至近两倍，但不是整 两倍，因此需要进行补零操作。原来的两个长度为m的Nuttall窗序列进行自卷 积运算时，得到的自卷积窗函数长度为2m-1。在仿真时，需要将在该序列首端 或者末尾补零，以得到长度为2m的序列。同理，将p个长度为m的Nuttall窗序 列进行p-1自卷积运算时，需要补p-1个0，以得到长度为pm的自卷积窗函数。

(2)仿真时区分时域运算和频域运算。Nuttall窗进行自卷积时，时域运算 是对p个母窗函数进行自卷积，如式(6)。在频域运算则是对Nuttall窗的频谱 函数求p次幂，如式(7)。

窗函数的主瓣宽度定义为，在W_N(ω)频谱图中，当|W_N(ω)|＝0时，距离原点最 近的2个零值点距离即为窗函数的主瓣宽度。

|W_N(ω)|＝0，即式(3)等于0。式(3)需要满足的条件为

主瓣宽度B_N,即为d＝±1时的宽度

式中，M为窗函数的项数；N'为未卷积时Nuttall窗的长度。

同理，p阶纳托尔自卷积窗的主瓣宽度为

式中，M为窗函数的项数；N为纳托尔自卷积窗的长度。为主瓣宽度满 足约束条件

式中，f_i为谐波频率；f_i与f_i+1为相邻谐波；f_s为采样频率。

由式(12)可知：

(1)当N一定时，纳托尔自卷积窗的主瓣宽度等于p个Nuttall窗的主瓣宽 度。当p＝1时，卷积窗就是母窗Nuttall窗。

(2)当N变化时，自卷积窗的主瓣宽度与母窗的长度N'成反比。

(3)自卷积窗并不是无条件可以任意卷积，需要满足约束条件，这是因为 当N一定时，卷积阶数过大，则主瓣越宽，此时旁瓣对主瓣的干扰很大，无法 准确测量，一般取p＝2,3。

以p＝2时的4项5阶纳托尔自卷积窗为例来说明。取长度N＝50时，纳托尔 自卷积窗的幅频响应曲线如图5所示。可见，其旁瓣电平峰值为-121.9dB，旁 瓣衰减速率为84dB/oct。其旁瓣衰减速率优于Rife-Vincent窗和Nuttall窗。

3、三谱线插值FFT算法

输入时域信号x(t)的表达式为

式中，H为谐波的个数；当h＝1时，f₁、A₁、分别为基波的频率(Hz)、幅 值(V)和初相位(°)；当h≠1时，r_h表示第h项谐波次数，A_h、分别为第h项谐波 的幅值(V)和相位(°)。

含间谐波分量的信号x(t)以采样频率f_s均匀采样得到的离散时间信号为

对信号x(n)加纳托尔自卷积窗得到对x_w(n)做离散快速傅里叶变换(DTFT)得到

忽略其他次数谐波以及负频率点旁瓣的影响，得到简化后的加自卷积窗FFT 变换表达式为

式中，Δf为抽样间隔；k为抽样频率点；N为数据截断长度。

在实际测量间谐波时，为了得到某一项的谐波，只考虑一个点，忽略其余 点的影响，则有

非同步采样时，频谱发生泄漏，由于频谱泄露和栅栏效应，此时采样值f_i＝k_iΔf很少正好落在峰值谱线频率点上。设峰值频点附近抽样得到的最大谱线 为第k_i1条，此时它左右侧的谱线也含有丰富的谱线信息，记为第k_i0条和第k_i2条， 其大小关系为k_i0(＝k_i1-1)≤k_i1≤k_i2(＝k_i1+1)，这3条谱线的幅值分别记为y₀＝|y_w(k_i0Δf)|， y₁＝|y_w(k_i1Δf)|和y₂＝|y_w(k_i2Δf)|。

由于0≤k_i-k_i1≤1，为了便于计算，引入2个变量α和β，令 α＝k_i-k_i1-0.5(α∈[-0.5，0.5])，设β＝(y₂-y₀)/y₁，，将式(18)带入β可得：

记式(19)为β＝f(α)。此时α为自变量，β为因变量，其反函数为α＝f^-1(β)。

由式(17)可得p阶纳托尔自卷积窗的DTFT为

将式(20)带入式(19)，在α∈[-0.5,0.5]内取一组α值，得到对应的β，调用函 数polyfit(β,α,m)进行多项式反拟合，可得函数α＝F(β)的逼近式。其中，m为拟合 逼近多项式的阶数。

由此，由式(17)可得第i项频率的插值修正公式为

f_i＝k_iΔf＝(α+k_i1+0.5)Δf (21)

为了使测量误差最小，需要对幅值频率进行加窗插值修正，最大谱线k_i1中 含有丰富的信息，但是其左右两根谱线k_i0和k_i2依然含有信息，为了尽量减少误 差，采用三谱线插值修正公式，三根谱线加权平均计算，得到幅值修正公式为

由于N一般较大，式(22)可化简为

A_i＝N^-p(y₂+2y₁+y₀)g(α) (23)

第i项谐波相位的修正公式为

式中h＝0,1,2。

4、谐波参数插值修正

在α∈[-0.5,0.5]内取一组值，分别由式(19)、(23)得到对应的β和g(α)值，调 用函数polyfit(β,α,m)进行多项式反拟合，可得函数α＝F(β)的逼近式，再调用函数 polyfit(α,g(α),m)，求出g(α)＝ψ(α)。

由于式(13)对于自卷积窗的限制，自卷积窗并不是无条件可以任意卷积， 需要满足约束条件，如果无限制卷积，旁瓣对主瓣的干扰很大，此时的频谱泄 露和栅栏效应将会更加严重，不仅无法提高测量精度，测量结果误差也会比普 通加窗测量误差大很多倍，因此一般取p＝2,3。

考虑到计算速度等因素，本论文取p＝2，来说明纳托尔自卷积窗对于动态 间谐波以及加白噪声时，测量的优越性。

当p＝2时，基于4项5阶纳托尔自卷积窗的三谱线插值修正公式为

α＝2.017780726618934β-0.14416912519588β³+

0.020680914246738β⁵-0.003614378944554β⁷ (25)

g(α)＝21.93553832925077+2.726188198628158α²

+0.176288053666237α⁴+0.008047452065558α⁶ (26)

将式(25)、(26)带入式(23)、(24)就可以得到第i项谐波修正后的幅值、 相位。

5、仿真分析

为验证本发明算法的正确性，分别对含非整数次谐波信号和21次复杂谐波 信号进行仿真分析。

A、含非整数次谐波的FFT分析

在幅值较小的非整数次谐波信号附近存在幅值较大的间谐波时，在非同步 采样情况下，频谱泄漏和栅栏效应将导致大量的弱信号分量被淹没，致使弱间 谐波信号未被检测出，或者检测误差较大。由于间谐波对电网危害极大，这将 造成严重的后果。

采用本发明的纳托尔自卷积窗三谱线插值算法，对表3中的信号进行分析 检测(表3中的数据与文献[14]中的含间谐波信号相同)，采样频率f_s＝1500H_Z,采 样点数N＝1024。表4为频率、幅值和相位的测量结果以及相对误差，其中 ae-b＝a×10^-b。

表3含间谐波信号的频率、幅值和相位：

表4仿真参数结果比较：

表5含间谐波时不同算法幅值误差：

表6含间谐波时不同算法相位误差：

在含有相同的间谐波信号的输电线路中，采用本发明的纳托尔自卷积三谱 线插值算法与传统算法对比可得

频率：本算法所有频率相对误差检测精度达到10^-12，间谐波相对于基波和谐 波而言幅值小，信号微弱。由表4可知，对于1.6次和3.5次间谐波本算法让 然可以保持很高的精确度。

幅值：如表5所示，本论文对于间谐波的幅值检测精度为10^-9。比文献[14] 中，对于1.6次和3.5次间谐波改进的四项5阶Nuttall窗四谱线检测算法幅 值相对误差检测精度高6个精确度，其余谐波也均高出1-5个检测精度。

相位：如表6所示，本论文对于间谐波的相位检测精度为10^-9。比文献[14] 中，对于1.6次间谐波改进的四项5阶Nuttall窗四谱线检测算法幅值相对误 差检测精度高6个精确度，其余谐波也均高出1-5个检测精度。

B、动态微弱间谐波的FFT检测

频率的变化将会导致频谱泄漏和栅栏效应更加严重，尤其是对于1.6次和 3.5次微弱的间谐波信号，频率的波动很容易“淹没”此类微弱间谐波信号。本 发明算法对频谱波动的含有微弱间谐波的信号进行加窗插值检测，充分利用最 大谱线以及最大谱线左右两侧的谱线进行三谱线插值修正，获得动态信号的频 率、幅值和相位信息。设基波频率的变化范围为49.5-51.5Hz，变动步长为0.01 Hz，采用基于纳托尔自卷积窗的三谱线插值FFT算法对表3中含有间谐波的信 号进行检测，得到各次谐波的频率、幅值和相位误差结果如下：

即使在频率波动时，检测微弱的间谐波动态频率误差也不会超过10^-7；检测 微弱的间谐波动态幅值误差不会超过10^-5；检测微弱的间谐波动态幅值误差不会 超过10^-4。由此可见，本发明算法能有效抑制基波波动的情况下间谐波的频谱泄 漏和栅栏效应，具有很高的间谐波检测精度。

C、白噪声影响动态间谐波的FFT分析

在实际输电线路中，信号中往往含有一定的白噪声。在其他条件不变的情 况下，在表3间谐波信号中加入不同信噪比SNR(signal-noise ratio)的高斯 白噪声，验证本发明算法的抗噪性和实际可行性。

表7 SNR＝70dB的仿真参数结果比较：

表8 SNR＝50dB的仿真参数结果比较：

表9 SNR＝30dB的仿真参数结果比较：

设备的信噪比SNR越高表明它产生的杂音越少。表7中，当SNR＝70dB时， 间谐波的频率误差不低于10^-6，幅值误差不低于10^-5，相位误差不低于10^-3。表8 中，当SNR＝50dB时，间谐波的频率误差不低于10^-5，幅值误差不低于10^-4，相位 误差不低于10^-3。表9中，在不考虑相位时，当SNR＝30dB时，间谐波的频率误 差不低于10^-4，幅值误差不低于10^-3。

当逐步降低信噪比SNR时，即杂音越来越大时，本发明检测精度虽然随着 噪音的增大会有所下降，但是均在较高的检测精度之内，结果也进一步表明本 发明算法具有较好的抗噪性。在一些实际测量情况下，如不要求相位精确度很 高的测量情况，本发明可以应用于噪音较大的测量环境。

综上所述，本发明利用4项5阶纳托尔窗进行2阶自卷积运算构造纳托尔 自卷积窗。经仿真实验表明纳托尔自卷积窗的旁瓣峰值电平、旁瓣衰减速率均 正比于卷积阶数，在非同步采样下，较余弦组合窗具有更好的频谱泄漏抑制能 力。本发明采用三谱线插值修正算法，对比纳托尔四谱线插值算法，虽然四谱 线可以更多的利用旁瓣中的信息，但是本文中由于自卷积窗抑制频谱泄漏和栅 栏效应的能力强，只需要三谱线插值就可以达到很高的精度，并且大大节约了 计算时间，减少了计算量。仿真结果验证了自卷积窗能更精确的检测动态微弱 间谐波信号的频率、幅值和相位，抗噪性能显著提高。

上面结合附图对本发明进行了示例性描述，显然本发明具体实现并不受上 述方式的限制，只要采用了本发明的方法构思和技术方案进行的各种非实质性 的改进，或未经改进将本发明的构思和技术方案直接应用于其它场合的，均在 本发明的保护范围之内。

Claims (7)

1.一种基于纳托尔自卷积窗三谱线插值FFT的动态间谐波分析方法，其特征在于：

步骤1、对动态谐波信号进行采样，得到所需要的时域序列；

步骤2、用纳托尔自卷积窗对采样的离散序列进行加窗截断，得到所需要的加窗信号；

步骤3、对加窗截断信号做快速傅里叶变换，得到信号的频谱信息；

步骤4、搜索最大谱线k_i1以及其左右谱线k_i0和k_i2，采集三条所需要的谱线，并确定所对应的幅值；

步骤5、引入变量，通过算法求取插值修正公式；

步骤6、用三谱线插值公式对所测量的频率、幅值以及相位进行修正。

2.根据权利要求1所述的基于纳托尔自卷积窗三谱线插值FFT的动态间谐波分析方法，其特征在于：

所述步骤1中，时域信号x(t)的表达式为

式中，H为谐波的个数；当h＝1时，f₁、A₁、分别为基波的频率、幅值和初相位；当h≠1时，r_h表示第h项谐波次数，A_h、分别为第h项谐波的幅值和相位。

3.根据权利要求2所述的基于纳托尔自卷积窗三谱线插值FFT的动态间谐波分析方法，其特征在于：

所述步骤2中，含间谐波分量的信号x(t)以采样频率f_s均匀采样得到的离散时间信号为

对信号x(n)加纳托尔自卷积窗得到频谱信息

4.根据权利要求3所述的基于纳托尔自卷积窗三谱线插值FFT的动态间谐波分析方法，其特征在于：

所述步骤3中，对x_w(n)做离散快速傅里叶变换得到

简化后的加自卷积窗FFT变换表达式为

式中，Δf为抽样间隔；k为抽样频率点；N为数据截断长度。

5.根据权利要求4所述的基于纳托尔自卷积窗三谱线插值FFT的动态间谐波分析方法，其特征在于：

所述步骤4中，只考虑一个抽样频率点，忽略其余抽样频率点的影响，则有

设峰值频点附近抽样得到的最大谱线为第k_i1条，则它左右侧的谱线也含有丰富的谱线信息，记为第k_i0条和第k_i2条，其大小关系为k_i0(＝k_i1-1)≤k_i1≤k_i2(＝k_i1+1)，上述3条谱线的幅值分别记为y₀＝|y_w(k_i0Δf)|，y₁＝|y_w(k_i1Δf)|和y₂＝|y_w(k_i2Δf)|。

6.根据权利要求5所述的基于纳托尔自卷积窗三谱线插值FFT的动态间谐波分析方法，其特征在于：

所述步骤5中，由于0≤k_i-k_i1≤1，则引入2个变量α和β，令α＝k_i-k_i1-0.5(α∈[-0.5,0.5])，设β＝(y₂-y₀)/y₁，则有：

则β＝f(α)，式中α为自变量，β为因变量，其反函数为α＝f^-1(β)；

可得p阶纳托尔自卷积窗的DTFT为

在α∈[-0.5,0.5]内取一组α值，得到对应的β，调用函数polyfit(β,α,m)进行多项式反拟合，可得函数α＝F(β)的逼近式，其中，m为拟合逼近多项式的阶数；

由此可得第i项频率的插值修正公式为

f_i＝k_iΔf＝(α+k_i1+0.5)Δf。

7.根据权利要求6所述的基于纳托尔自卷积窗三谱线插值FFT的动态间谐波分析方法，其特征在于：

所述步骤6中，采用三谱线插值修正公式，三根谱线加权平均计算，得到幅值修正公式为

由于N一般较大则可化简为

A_i＝N^-p(y₂+2y₁+y₀)g(α)

第i项谐波相位的修正公式为

式中h＝0,1,2。