CN105137181A - 基于Nuttall-Kaiser组合窗双谱线插值的谐波分析方法 - Google Patents

Abstract

本发明涉及电力系统电能质量技术领域，特别是一种基于Nuttall-Kaiser组合窗双谱线插值的谐波分析算法。基于快速傅里叶变换（FFT）算法，将余弦组合函数和第1类变形零阶贝塞尔函数的频谱表达式进行组合，并将其代入加窗后信号的离散傅里叶表达式，然后采用双谱线插值方法，比较找出频谱图中每个峰附近的两根谱线，将其代入插值公式，利用MATLAB软件中的曲线拟合函数进行多项式拟合逼近，最后求出频率、幅值和相角修正公式，并求得频率、幅值和相角的值和相对误差。该算法能够对非整数周期截断造成的频谱泄露有着很好的抑制作用；并且在电力系统含有复杂谐波信号的条件下，仍然拥有较高的测量精度和抗噪性；相比于Nuttall四项三阶窗和Kaiser窗单一进行插值更具可靠性。

Description

基于Nuttall-Kaiser组合窗双谱线插值的谐波分析方法

技术领域

本发明涉及电力系统电能质量技术领域，特别是一种基于Nuttall-Kaiser组合窗双谱线插值的谐波分析算法。

背景技术

21世纪以来，随着智能电网技术的不断成熟，人们对电网参数的测量要求也在不断的提高。但是各种非线性电力负载和电力电子装置的广泛应用，正在不断的影响着电网的稳定运行。此外，广域测量系统(WAMS)作为当前智能电网的重要组成部分，它的基础是同步向量测量技术，通过同步向量测量单元(phasormeasurementunit，PMU)对电网参数进行实时测量，为智能电网的谐波分析提供数据支持，但由于装置本身以及授时装置的影响，PMU大都采用非同步采样的方式进行数据采集。在非同步采样的情况下，对信号进行快速傅里叶变换(fastFouriertransform，FFT)，时域截断引起的频谱泄露和频域离散化引起的栅栏效应会使谐波分析产生误差。

针对以上问题，国内外的研究学者采用了不同的方法来减少误差。硬件上可以采用锁相环电路来实现同步采样，而软件上则可以采用加窗插值FFT算法、小波分析算法、神经网络算法等进行误差的矫正。其中，加窗插值是当应用较为广泛的一种算法。1979年由V.K.Jain等提出的矩形窗插值算法，可以有效的提高计算精度。随着近年来窗函数的不断发展，大大小小的窗函数已经有20余种，Hanning窗、Blackman-Harris窗、Nuttall窗、Kaiser窗以及文献所提出的卷积窗等已经被应用到FFT谐波分析中。其中Nuttall窗是一种余弦组合窗，旁瓣峰值电平小且旁瓣渐进衰减速率大，可以很好地抑制临近泄露和远离泄露；而Kaiser窗可以定义一组可调的窗函数，其主瓣能量和旁瓣能量的比例近乎最大，且可自由选择主瓣宽度和旁瓣高度之间的比重。二者配合插值算法，对栅栏效应和频谱泄露所带来的误差都有很好的抑制作用。

发明内容

本发明的目的在于提供一种基于Nuttall-Kaiser组合窗双谱线插值的谐波分析算法。该算法综合了Nuttall四项三阶窗和Kaiser窗的优点，能对电力系统含有复杂谐波的信号进行分析，相对Nuttall四项三阶窗和Kaiser窗单一进行插值更具可靠性和抗噪性，提高了电力系统发生故障时的测量精度。

为实现上述目的，本发明的技术方案是：将余弦组合函数和第1类变形零阶贝塞尔函数的频谱表达式进行组合，并将组合后的函数代入加窗后信号的离散傅里叶表达式，然后采用双谱线插值方法，比较找出频谱图中每个峰附近的两根谱线，将两根谱线代入插值公式，并利用MATLAB软件中的曲线拟合函数polyfit(·)进行多项式拟合逼近，最后求出频率、幅值和相角修正公式，并求得频率、幅值和相角的值和相对误差。其具体步骤如下：

步骤1：假设电力系统电力信号含有谐波分量，其基波的频率、幅值和相角分别为f₁、A₁和θ₁，通过采样频率均匀采样得到的离散时间信号为：式中：i是谐波次数，当i≠1时，A_i、θ_i分别表示第i次谐波的幅值和相角；

步骤2：利用Nuttall-Kaiser组合窗各部分的表示：

Nuttall窗的时域表达式为： $w_N\left(n\right)=\underset{m=0}{\overset{M-1}{\Σ}}{(-1)}^m{b}_{m}cos\left(\frac{2\mathrm{\π}n\·m}{N}\right),$ 式中：M为窗函数的项数；n＝0,1,2,…,N-1；N为自然数；b_m满足约束条件：经过傅里叶变换后Nuttall窗的频谱函数表达式为： $W_N\left(w\right)=\underset{m=0}{\overset{M-1}{\Σ}}{(-1)}^m\frac{b_m}{2}\[W_R(w-\frac{2\mathrm{\π}}{N}m)+W_R(w+\frac{2\mathrm{\π}}{N}m)\],$ 式中 $W_R\left(w\right)=(sin\frac{Nw}{2}/sin\frac{w}{2})e^{-j\frac{N-1}{2}w};$

Kaiser窗的时域表达式为： $w_K\left(n\right)=\frac{I_0\[\mathrm{\β}\sqrt{1-{\left(\frac{n}{N/2}\right)}^2}\]}{I_0\left(\mathrm{\β}\right)},$ 0≤|n|≤N/2式中：I₀(β)是第1类变形零阶贝塞尔函数，β是窗函数的形状参数，经过傅里叶变换后Kaiser窗的频谱函数表达式为：

$W_K\left(w\right)=\frac{N}{I_0\left(\mathrm{\β}\right)}\·\frac{\sinh \left(\sqrt{{\mathrm{\β}}^2-{(Nw/2)}^2}\right)}{\sqrt{{\mathrm{\β}}^2-{(Nw/2)}^2}},$

将上式平移(N-1)/2，使其范围能与Nuttall窗一样满足[0,N-1]，得

$W_K\left(w\right)=\frac{N-1}{I_0\left(\mathrm{\β}\right)}\·\frac{\sinh \left(\sqrt{{\mathrm{\β}}^2-\[(N-1)w/2\]}\right)}{\sqrt{{\mathrm{\β}}^2-\[(N-1)w/2\]}}\·e^{-j\frac{N-1}{N}\mathrm{\π}k};$

Nuttall-Kaiser组合窗的频谱函数表达式为：W(w)＝0.5×[W_N(w)+W_K(w)]；

步骤3：用组合窗的窗函数w(n)对步骤1中的信号x(n)进行处理，可得加窗后信号的离散傅里叶表达式： $X\left(k\mathrm{\Δ}f\right)=\frac{A_i}{2j}{e}^{{\mathrm{j\θ}}_i}W\[\frac{2\mathrm{\π}(k\mathrm{\Δ}f-f_i)}{f_s}\],$ 令 $w=\frac{2\mathrm{\π}}{N}k,$ 可得：

$W\left(\frac{2\mathrm{\π}}{N}k\right)=0.5\×\{sin\mathrm{\π}k\·e^{-j\mathrm{\π}k}\[\underset{m=0}{\overset{M-1}{\Σ}}{(-1)}^m\frac{b_m}{2}\·\frac{sin\frac{2\mathrm{\π}k}{N}}{sin\frac{\mathrm{\π}(k-m)}{N}\·sin\frac{\mathrm{\π}(k+m)}{N}}\]+\frac{N-1}{I_0\left(\mathrm{\β}\right)}\·\frac{\sinh \left(\sqrt{{\mathrm{\β}}^2-{\[\frac{(N-1)}{N}\mathrm{\π}k\]}^2}\right)}{\sqrt{{\mathrm{\β}}^2-{\[\frac{(N-1)}{N}\mathrm{\π}k\]}^2}}\·e^{-j\frac{N-1}{N}\mathrm{\π}k}\},$

步骤4：非同步采样时，峰值频率f_i≠k_i△f，其中k_i为峰值点对应的谱线，假设频谱图中最大谱线为k_i1，次最大谱线为k_i2(＝k_i1+1)，则k_i1≤k_i≤k_i2，这两条谱线对应幅值记为y₁和y₂，设α＝k_i-k_i1-0.5，其中α的范围为(-0.5,0.5)，可得：

$\mathrm{\β}=\frac{\left|W\left(\frac{2\mathrm{\π}(-\mathrm{\α}+0.5)}{N}\right)\right|-\left|W\left(\frac{2\mathrm{\π}(-\mathrm{\α}-0.5)}{N}\right)\right|}{\left|W\left(\frac{2\mathrm{\π}(-\mathrm{\α}+0.5)}{N}\right)\right|+\left|W\left(\frac{2\mathrm{\π}(-\mathrm{\α}-0.5)}{N}\right)\right|};$

步骤5：令k＝-α±0.5代入

$W\left(\frac{2\mathrm{\π}}{N}k\right)=0.5\×\{sin\mathrm{\π}k\·e^{-j\mathrm{\π}k}\[\underset{m=0}{\overset{M-1}{\Σ}}{(-1)}^m\frac{b_m}{2}\·\frac{sin\frac{2\mathrm{\π}k}{N}}{sin\frac{\mathrm{\π}(k-m)}{N}\·sin\frac{\mathrm{\π}(k+m)}{N}}\]+\frac{N-1}{I_0\left(\mathrm{\β}\right)}\·\frac{\sinh \left(\sqrt{{\mathrm{\β}}^2-{\[\frac{(N-1)}{N}\mathrm{\π}k\]}^2}\right)}{\sqrt{{\mathrm{\β}}^2-{\[\frac{(N-1)}{N}\mathrm{\π}k\]}^2}}\·e^{-j\frac{N-1}{N}\mathrm{\π}k}\}$ 中，化简得到：

$\left|W\left(\frac{2\mathrm{\π}\left(-\mathrm{\α}\±0.5\right)}{N}\right)\right|\≈0.5\×\{\left|\sin \mathrm{\π}\left(-\mathrm{\α}\±0.5\right)\·\[\underset{m=0}{\overset{M-1}{\Σ}}{\left(-1\right)}^m\frac{b_m}{\mathrm{\π}}\frac{N\left(-\mathrm{\α}\±0.5\right)}{{\left(-\mathrm{\α}\±0.5\right)}^2-m^2}\]\right|+\left|\frac{N-1}{I_0\left(\mathrm{\β}\right)}\·\frac{\sinh \left(\sqrt{{\mathrm{\β}}^2-{\left(\mathrm{\π}k\right)}^2}\right)}{\sqrt{{\mathrm{\β}}^2-{\left(\mathrm{\π}k\right)}^2}}\right|\}$

步骤6：把步骤5的结果代入步骤4，再利用MATLAB软件中的曲线拟合函数polyfit(·)进行多项式拟合逼近，可得α＝h^-1(β)的逼近式为：α＝H(β)；

步骤7：由β可求出参数α，则频率的修正公式为：f_i＝k_i△f＝(α+k_i1+0.5)△f；

步骤8：幅值修正是对k₁和k₂谱线进行加权平均，其修正公式为：

$A_i=\frac{A_{i1}\left|W\left(\frac{2\mathrm{\π}(k_{i1}-k_i)}{N}\right)\right|+A_{i2}\left|W\left(\frac{2\mathrm{\π}(k_{i2}-k_i)}{N}\right)\right|}{\left|W\left(\frac{2\mathrm{\π}(k_{i1}-k_i)}{N}\right)\right|+\left|W\left(\frac{2\mathrm{\π}(k_{i2}-k_i)}{N}\right)\right|}=\frac{2(y_1+y_2)}{\left|W\left(\frac{2\mathrm{\π}(-\mathrm{\α}+0.5)}{N}\right)\right|+\left|W\left(\frac{2\mathrm{\π}(-\mathrm{\α}-0.5)}{N}\right)\right|}$

当N值较大时，上式可以化简为：其中， $v\left(a\right)=\frac{2N}{\left|W\left(\frac{2\mathrm{\π}(-\mathrm{\α}+0.5)}{N}\right)\right|+\left|W\left(\frac{2\mathrm{\π}(-\mathrm{\α}-0.5)}{N}\right)\right|},$ 式中的W()采用步骤5中未化简的公式；

同理，采用拟合函数进行多项式逼近，可得逼近式为：

其中， $g\left(a\right)=\frac{2N}{\left|W\left(\frac{2\mathrm{\π}(-\mathrm{\α}+0.5)}{N}\right)\right|+\left|W\left(\frac{2\mathrm{\π}(-\mathrm{\α}-0.5)}{N}\right)\right|},$ 式中的W()采用步骤5中化简的公式；

步骤9：相角的修正公式为：

${\mathrm{\θ}}_i=\arg \[X\left(k_{i1}\mathrm{\Δ}f\right)\]+\frac{\mathrm{\π}}{2}-\arg \[W\left(\frac{2\mathrm{\π}\left(k_{i1}\mathrm{\Δ}f-f_i\right)}{f_s}\right)\];$

步骤10：计算出各次谐波的频率、幅值和相角，并求取相对误差。

在本发明一实施例中，还包括以下步骤：在[-0.49,0.49]范围内每隔0.001取一个数组成一组α值，分别代入步骤4和步骤8得到对应的一组β和v(α)，然后在MATLAB中分别调用polyfit(β,α,m)和polyfit(α,v(α),m)，polyfit(β,α,m)m为5，polyfit(α,v(α),m)m为4，得到H(β)和g(α)的系数：

α＝H(β)＝1.2618269186β⁵+1.0368877011β³+3.6499846828β

g(α)＝0.1068160200α⁴+0.8628506225α²+3.6418169592。

相较于现有技术，本发明有以下有益效果：

1、能够对非整数周期截断造成的频谱泄露有着很好的抑制作用。

2、在电力系统含有复杂谐波信号的条件下，仍然拥有较高的测量精度和抗噪性。

3、相比于Nuttall四项三阶窗和Kaiser窗单一进行插值更具可靠性。

附图说明

图1是本发明实施例的工作流程图。

图2是Nuttall-Kaiser组合窗与Nuttall、Kaiser窗的归一化对数谱比较。

图3是Nuttall-Kaiser组合窗与Nuttall、Kaiser窗幅值相对误差曲线图。

图4是Nuttall-Kaiser组合窗与Nuttall、Kaiser窗相角相对误差曲线图。

具体实施方式

下面结合附图和具体实施方式对本发明做进一步说明。

本发明提供一种基于Nuttall-Kaiser组合窗双谱线插值的谐波分析算法，具体流程结合图1进行说明，并对电力系统电力信号进行分析。将余弦组合函数和第1类变形零阶贝塞尔函数的频谱表达式进行组合，并将其代入加窗后信号的离散傅里叶表达式，然后采用双谱线插值方法，比较找出频谱图中每个峰附近的两根谱线，将其代入插值公式，再利用MATLAB软件中的曲线拟合函数polyfit(·)进行多项式拟合逼近，最后求出频率、幅值和相角修正公式，并求得频率、幅值和相角的值和相对误差。具体步骤如下：

步骤1：假设信号含有谐波分量，其基波的频率、幅值和相角分别为、和，通过采样频率均匀采样得到的离散时间信号为：

$x\left(n\right)=\underset{i=1}{\overset{I}{\Σ}}{A}_{i}sin(2\mathrm{\π}\frac{{\mathrm{if}}_1}{f_s}n+{\mathrm{\θ}}_i)$

式中：i是谐波次数，当i≠1时，A_i、θ_i分别表示第i次谐波的幅值和相角，采样频率f_s取为1000Hz，截断信号的数据长度N为1024

步骤2：Nuttall-Kaiser组合窗各部分的表示：

Nuttall窗是一种余弦组合窗，其时域表达式为：

$w_N\left(n\right)=\underset{m=0}{\overset{M-1}{\Σ}}{(-1)}^m{b}_{m}cos\left(\frac{2\mathrm{\π}n\·m}{N}\right)$

式中：M为窗函数的项数；n＝0,1,2,…,N-1，N为自然数；b_m应满足约束条件， $\underset{m=0}{\overset{M-1}{\Σ}}{(-1)}^m{b}_m=0$

经过傅里叶变换，Nuttall窗的频谱函数表达式为：

$W_N\left(w\right)=\underset{m=0}{\overset{M-1}{\Σ}}{(-1)}^m\frac{b_m}{2}\[W_R(w-\frac{2\mathrm{\π}}{N}m)+W_R(w+\frac{2\mathrm{\π}}{N}m)\],$

式中 $W_R\left(w\right)=(sin\frac{Nw}{2}/sin\frac{w}{2})e^{-j\frac{N-1}{2}w};$

Kaiser窗为可自定义一组可调的窗函数，其时域表达式为：

$w_K\left(n\right)=\frac{I_0\[\mathrm{\β}\sqrt{1-{\left(\frac{n}{N/2}\right)}^2}\]}{I_0\left(\mathrm{\β}\right)},0\≤\left|n\right|\≤N/2,$

式中：I₀(β)是第1类变形零阶贝塞尔函数，β是窗函数的形状参数。经过傅里叶变换，Kaiser窗的频谱函数表达式为：

$W_K\left(w\right)=\frac{N}{I_0\left(\mathrm{\β}\right)}\·\frac{\sinh \left(\sqrt{{\mathrm{\β}}^2-{\left(Nw/2\right)}^2}\right)}{\sqrt{{\mathrm{\β}}^2-{\left(Nw/2\right)}^2}},$

将上式平移(N-1)/2，使其范围能与Nuttall窗一样满足[0,N-1]，可得：

$W_K\left(w\right)=\frac{N-1}{I_0\left(\mathrm{\β}\right)}\·\frac{\sinh \left(\sqrt{{\mathrm{\β}}^2-\[(N-1)w/2\]}\right)}{\sqrt{{\mathrm{\β}}^2-\[(N-1)w/2\]}}\·e^{-j\frac{N-1}{N}\mathrm{\π}k}$

步骤3：用组合窗的窗函数w(n)对步骤1形式的信号x(n)进行处理，可得加窗后信号的离散傅里叶表达式：

$X\left(k\mathrm{\Δ}f\right)=\frac{A_i}{2j}{e}^{{\mathrm{j\θ}}_i}W\[\frac{2\mathrm{\π}(k\mathrm{\Δ}f-f_i)}{f_s}\]$

式中W(·)是Nuttall-Kaiser组合窗的频谱函数表达式为：

W(w)＝0.5×[W_N(w)+W_K(w)]

令 $w=\frac{2\mathrm{\π}}{N}k,$ 可得：

$W\left(\frac{2\mathrm{\π}}{N}k\right)=0.5\×\{sin\mathrm{\π}k\·e^{-j\mathrm{\π}k}\[\underset{m=0}{\overset{M-1}{\Σ}}{(-1)}^m\frac{b_m}{2}\·\frac{sin\frac{2\mathrm{\π}k}{N}}{sin\frac{\mathrm{\π}(k-m)}{N}\·sin\frac{\mathrm{\π}(k+m)}{N}}\]+\frac{N-1}{I_0\left(\mathrm{\β}\right)}\·\frac{\sinh \left(\sqrt{{\mathrm{\β}}^2-{\[\frac{(N-1)}{N}\mathrm{\π}k\]}^2}\right)}{\sqrt{{\mathrm{\β}}^2-{\[\frac{(N-1)}{N}\mathrm{\π}k\]}^2}}\·e^{-j\frac{N-1}{N}\mathrm{\π}k}\}$

步骤4：非同步采样时，峰值频率f_i≠k_i△f，其中k_i为峰值点对应的谱线。假设频谱图中最大谱线为k_i1，次最大谱线为k_i2(＝k_i1+1)，则k_i1≤k_i≤k_i2，这两条谱线对应幅值记为y₁和y₂，设α＝k_i-k_i1-0.5，其中α的范围为(-0.5,0.5)，可得：

$\mathrm{\β}=\frac{\left|W\left(\frac{2\mathrm{\π}(-\mathrm{\α}+0.5)}{N}\right)\right|-\left|W\left(\frac{2\mathrm{\π}(-\mathrm{\α}-0.5)}{N}\right)\right|}{\left|W\left(\frac{2\mathrm{\π}(-\mathrm{\α}+0.5)}{N}\right)\right|+\left|W\left(\frac{2\mathrm{\π}(-\mathrm{\α}-0.5)}{N}\right)\right|}$

步骤5：令k＝-α±0.5代入步骤3，因为N一般较大，得

$\left|W\left(\frac{2\mathrm{\π}\left(-\mathrm{\α}\±0.5\right)}{N}\right)\right|\≈0.5\×\{\left|\sin \mathrm{\π}\left(-\mathrm{\α}\±0.5\right)\·\[\underset{m=0}{\overset{M-1}{\Σ}}{\left(-1\right)}^m\frac{b_m}{\mathrm{\π}}\frac{N\left(-\mathrm{\α}\±0.5\right)}{{\left(-\mathrm{\α}\±0.5\right)}^2-m^2}\]\right|+\left|\frac{N-1}{I_0\left(\mathrm{\β}\right)}\·\frac{\sinh \left(\sqrt{{\mathrm{\β}}^2-{\left(\mathrm{\π}k\right)}^2}\right)}{\sqrt{{\mathrm{\β}}^2-{\left(\mathrm{\π}k\right)}^2}}\right|\}$

步骤6：把步骤5式代入步骤4，利用MATLAB软件中的曲线拟合函数polyfit(·)进行多项式拟合逼近，可得α＝h^-1(β)的逼近式为：α＝H(β)

步骤7：由β可求出参数α，则频率的修正公式为：

f_i＝k_i△f＝(α+k_i1+0.5)△f

步骤8：幅值修正是对k₁和k₂谱线进行加权平均，其修正公式为：

$A_i=\frac{A_{i1}\left|W\left(\frac{2\mathrm{\π}(k_{i1}-k_i)}{N}\right)\right|+A_{i2}\left|W\left(\frac{2\mathrm{\π}(k_{i2}-k_i)}{N}\right)\right|}{\left|W\left(\frac{2\mathrm{\π}(k_{i1}-k_i)}{N}\right)\right|+\left|W\left(\frac{2\mathrm{\π}(k_{i2}-k_i)}{N}\right)\right|}=\frac{2(y_1+y_2)}{\left|W\left(\frac{2\mathrm{\π}(-\mathrm{\α}+0.5)}{N}\right)\right|+\left|W\left(\frac{2\mathrm{\π}(-\mathrm{\α}-0.5)}{N}\right)\right|}$

当N值较大时，上式可以化简为：

$A_i=\frac{(y_1+y_2)}{N}v\left(\mathrm{\α}\right);$ 其中， $v\left(a\right)=\frac{2N}{\left|W\left(\frac{2\mathrm{\π}(-\mathrm{\α}+0.5)}{N}\right)\right|+\left|W\left(\frac{2\mathrm{\π}(-\mathrm{\α}-0.5)}{N}\right)\right|},$ 式中的W()采用步骤5中未化简的公式；

同理，采用拟合函数进行多项式逼近，可得逼近式为：

其中， $g\left(a\right)=\frac{2N}{\left|W\left(\frac{2\mathrm{\π}(-\mathrm{\α}+0.5)}{N}\right)\right|+\left|W\left(\frac{2\mathrm{\π}(-\mathrm{\α}-0.5)}{N}\right)\right|},$ 式中的W()采用步骤5中化简的公式；

在[-0.49,0.49]范围内每隔0.001取一个数组成一组α值，分别代入步骤4和步骤8得到对应的一组β和v(α)，然后在MATLAB中分别调用polyfit(β,α,m)和polyfit(α,v(α),m)，取前者m为5，后者m为4，可推出相应的H(β)和g(α)的系数：

α＝H(β)＝1.2618269186β⁵+1.0368877011β³+3.6499846828β

g(α)＝0.1068160200α⁴+0.8628506225α²+3.6418169592；

步骤9：相角的修正公式为：

${\mathrm{\θ}}_i=\arg \[X\left(k_{i1}\mathrm{\Δ}f\right)\]+\frac{\mathrm{\π}}{2}-\arg \[W\left(\frac{2\mathrm{\π}(k_{i1}\mathrm{\Δ}f-f_i)}{f_s}\right)\];$

步骤10：计算出各次谐波的频率、幅值和相角，并求取相对误差。

Nuttall-Kaiser组合窗与Nuttall、Kaiser窗的归一化对数谱比较参见图2，其中幅值和相角的相对误差效果图分别示于附图3、4。可以看出采用本发明的技术方案相比于Nuttall四项三阶窗和Kaiser窗单一进行插值更具可靠性。

以上是本发明的较佳实施例，凡依本发明技术方案所作的改变，所产生的功能作用未超出本发明技术方案的范围时，均属于本发明的保护范围。

Claims (3)

1.一种基于Nuttall-Kaiser组合窗双谱线插值的谐波分析方法，其特征在于：将余弦组合函数和第1类变形零阶贝塞尔函数的频谱表达式进行组合，并将组合后的函数代入加窗后信号的离散傅里叶表达式，然后采用双谱线插值方法，比较找出频谱图中每个峰附近的两根谱线，将两根谱线代入插值公式，并利用MATLAB软件中的曲线拟合函数polyfit(·)进行多项式拟合逼近，最后求出频率、幅值和相角修正公式，并求得频率、幅值和相角的值和相对误差。

2.根据权利要求1所述的基于Nuttall-Kaiser组合窗双谱线插值的谐波分析方法，其特征在于：包括以下具体步骤：

步骤1：假设电力系统电力信号含有谐波分量，其基波的频率、幅值和相角分别为f₁、A₁和θ₁，通过采样频率均匀采样得到的离散时间信号为：式中：i是谐波次数，当i≠1时，A_i、θ_i分别表示第i次谐波的幅值和相角；

步骤2：利用Nuttall-Kaiser组合窗各部分的表示：

Nuttall窗的时域表达式为：式中：M为窗函数的项数；n＝0,1,2,…,N-1；N为自然数；b_m满足约束条件：经过傅里叶变换后Nuttall窗的频谱函数表达式为： $W_N\left(w\right)=\underset{m=0}{\overset{M-1}{\Σ}}{(-1)}^m\frac{b_m}{2}\[W_R(w-\frac{2\mathrm{\π}}{N}m)+W_R(w+\frac{2\mathrm{\π}}{N}m)\],$ 式中 $W_R\left(w\right)=(sin\frac{Nw}{2}/sin\frac{w}{2})e^{-j\frac{N-1}{2}w};$

Kaiser窗的时域表达式为：0≤|n|≤N/2，式中：I₀(β)是第1类变形零阶贝塞尔函数，β是窗函数的形状参数，经过傅里叶变换后Kaiser窗的频谱函数表达式为：

$W_K\left(w\right)=\frac{N}{I_0\left(\mathrm{\β}\right)}\·\frac{\sinh \left(\sqrt{{\mathrm{\β}}^2-{(Nw/2)}^2}\right)}{\sqrt{{\mathrm{\β}}^2-{(Nw/2)}^2}},$

将上式平移(N-1)/2，使其范围能与Nuttall窗一样满足[0,N-1]，得

$W_K\left(w\right)=\frac{N-1}{I_0\left(\mathrm{\β}\right)}\·\frac{\sinh \left(\sqrt{{\mathrm{\β}}^2-\[(N-1)w/2\]}\right)}{\sqrt{{\mathrm{\β}}^2-\[(N-1)w/2\]}}\·e^{-j\frac{N-1}{N}\mathrm{\π}k};$

Nuttall-Kaiser组合窗的频谱函数表达式为：W(w)＝0.5×[W_N(w)+W_K(w)]；

步骤3：用组合窗的窗函数w(n)对步骤1中的信号x(n)进行处理，可得加窗后信号的离散傅里叶表达式： $X\left(k\mathrm{\Δ}f\right)=\frac{A_i}{2j}{e}^{{\mathrm{j\θ}}_i}W\[\frac{2\mathrm{\π}(k\mathrm{\Δ}f-f_i)}{f_s}\],$

令 $w=\frac{2\mathrm{\π}}{N}k,$ 可得：

$W\left(\frac{2\mathrm{\π}}{N}k\right)=0.5\×\{sin\mathrm{\π}k\·e^{-j\mathrm{\π}k}\[\underset{m=0}{\overset{M-1}{\Σ}}{(-1)}^m\frac{b_m}{2}\·\frac{sin\frac{2\mathrm{\π}k}{N}}{sin\frac{\mathrm{\π}(k-m)}{N}\·\sin \frac{\mathrm{\π}(k+m)}{N}}\]+\frac{N-1}{I_0\left(\mathrm{\β}\right)}\·\frac{\sinh \left(\sqrt{{\mathrm{\β}}^2-{\[\frac{(N-1)}{N}\mathrm{\π}k\]}^2}\right)}{\sqrt{{\mathrm{\β}}^2-{\[\frac{(N-1)}{N}\mathrm{\π}k\]}^2}}\·e^{-j\frac{N-1}{N}\mathrm{\π}k}\},$

步骤4：非同步采样时，峰值频率f_i≠k_i△f，其中k_i为峰值点对应的谱线，假设频谱图中最大谱线为k_i1，次最大谱线为k_i2(＝k_i1+1)，则k_i1≤k_i≤k_i2，这两条谱线对应幅值记为y₁和y₂，设α＝k_i-k_i1-0.5，其中α的范围为(-0.5,0.5)，可得：

$\mathrm{\β}=\frac{\left|W\left(\frac{2\mathrm{\π}(-\mathrm{\α}+0.5)}{N}\right)\right|-\left|W\left(\frac{2\mathrm{\π}(-\mathrm{\α}-0.5)}{N}\right)\right|}{\left|W\left(\frac{2\mathrm{\π}(-\mathrm{\α}+0.5)}{N}\right)\right|+\left|W\left(\frac{2\mathrm{\π}(-\mathrm{\α}-0.5)}{N}\right)\right|};$

步骤5：令k＝-α±0.5代入

$W\left(\frac{2\mathrm{\π}}{N}k\right)=0.5\×\{sin\mathrm{\π}k\·e^{-j\mathrm{\π}k}\[\underset{m=0}{\overset{M-1}{\Σ}}{(-1)}^m\frac{b_m}{2}\·\frac{sin\frac{2\mathrm{\π}k}{N}}{sin\frac{\mathrm{\π}(k-m)}{N}\·\sin \frac{\mathrm{\π}(k+m)}{N}}\]+\frac{N-1}{I_0\left(\mathrm{\β}\right)}\·\frac{\sinh \left(\sqrt{{\mathrm{\β}}^2-{\[\frac{(N-1)}{N}\mathrm{\π}k\]}^2}\right)}{\sqrt{{\mathrm{\β}}^2-{\[\frac{(N-1)}{N}\mathrm{\π}k\]}^2}}\·e^{-j\frac{N-1}{N}\mathrm{\π}k}\}$

中，化简得到：

$\left|W\left(\frac{2\mathrm{\π}\left(-\mathrm{\α}\±0.5\right)}{N}\right)\right|\≈0.5\×\left\{\right|sin\mathrm{\π}\left(-\mathrm{\α}\±0.5\right)\·\[\underset{m=0}{\overset{M-1}{\Σ}}{(-1)}^m\frac{b_m}{\mathrm{\π}}\frac{N\left(-\mathrm{\α}\±0.5\right)}{{\left(-\mathrm{\α}\±0.5\right)}^2-m^2}\]|+|\frac{N-1}{I_0\left(\mathrm{\β}\right)}\·\frac{\sinh \left(\sqrt{{\mathrm{\β}}^2-{\left(\mathrm{\π}k\right)}^2}\right)}{\sqrt{{\mathrm{\β}}^2-{\left(\mathrm{\π}k\right)}^2}}\left|\right\}$

步骤6：把步骤5的结果代入步骤4，再利用MATLAB软件中的曲线拟合函数polyfit(·)进行多项式拟合逼近，可得α＝h^-1(β)的逼近式为：α＝H(β)；

步骤7：由β可求出参数α，则频率的修正公式为：f_i＝k_i△f＝(α+k_i1+0.5)△f；

步骤8：幅值修正是对k₁和k₂谱线进行加权平均，其修正公式为：

$A_i=\frac{A_{i1}\left|W\left(\frac{2\mathrm{\π}\left(k_{i1}-k_i\right)}{N}\right)\right|+A_{i2}\left|W\left(\frac{2\mathrm{\π}\left(k_{i2}-k_i\right)}{N}\right)\right|}{\left|W\left(\frac{2\mathrm{\π}\left(k_{i1}-k_i\right)}{N}\right)\right|+\left|W\left(\frac{2\mathrm{\π}\left(k_{i2}-k_i\right)}{N}\right)\right|}=\frac{2\left(y_1+y_2\right)}{\left|W\left(\frac{2\mathrm{\π}\left(-\mathrm{\α}+0.5\right)}{N}\right)\right|+\left|W\left(\frac{2\mathrm{\π}\left(-\mathrm{\α}-0.5\right)}{N}\right)\right|},$

当N值较大时，上式可以化简为：其中：

$v\left(a\right)=\frac{2N}{\left|W\left(\frac{2\mathrm{\π}(-\mathrm{\α}+0.5)}{N}\right)\right|+\left|W\left(\frac{2\mathrm{\π}(-\mathrm{\α}-0.5)}{N}\right)\right|};$

同理，采用拟合函数进行多项式逼近，可得逼近式为：其中： $g\left(a\right)=\frac{2N}{\left|W\left(\frac{2\mathrm{\π}(-\mathrm{\α}+0.5)}{N}\right)\right|+\left|W\left(\frac{2\mathrm{\π}(-\mathrm{\α}-0.5)}{N}\right)\right|};$

步骤9：相角的修正公式为：

${\mathrm{\θ}}_i=\arg \[X\left(k_{i1}\mathrm{\Δ}f\right)\]+\frac{\mathrm{\π}}{2}-\arg \[W\left(\frac{2\mathrm{\π}(k_{i1}\mathrm{\Δ}f-f_i)}{f_s}\right)\];$

步骤10：计算出各次谐波的频率、幅值和相角，并求取相对误差。

3.根据权利要求2所述的基于Nuttall-Kaiser组合窗双谱线插值的谐波分析方法，其特征在于：还包括以下步骤：在[-0.49,0.49]范围内每隔0.001取一个数组成一组α值，分别代入步骤4和步骤8得到对应的一组β和v(α)，然后在MATLAB中分别调用polyfit(β,α,m)和polyfit(α,v(α),m)，polyfit(β,α,m)m为5，polyfit(α,v(α),m)m为4，得到H(β)和g(α)的系数：

α＝H(β)＝1.2618269186β⁵+1.0368877011β³+3.6499846828β

g(α)＝0.1068160200α⁴+0.8628506225α²+3.6418169592。