CN104391178A - 一种基于Nuttall窗的时移相位差稳态谐波信号校正方法 - Google Patents

Abstract

本发明涉及一种基于Nuttall窗的时移相位差稳态谐波信号校正方法，属于谐波分析领域；该方法采用Nuttall窗对谐波信号进行加权，同时结合相位差校正原理求得频率校正量对谐波进行校正和分析。本发明提出的一种基于Nuttall窗的时移相位差稳态谐波信号校正方法能很好地抑制由于频谱泄漏引起的基波及整数次谐波之间的相互干扰，与谐波分析的比值校正法和加其它窗函数相比，信号截断时间相同时，基于Nuttall窗的时移相位差校正法可以获得更高的计算精度。

Description

一种基于Nuttall窗的时移相位差稳态谐波信号校正方法

技术领域

本发明属于谐波分析领域，涉及一种基于Nuttall窗的时移相位差稳态谐波信号校正方法。

背景技术

快速傅里叶变换(FFT)易于在嵌入式数字信号系统中实现，是谐波分析的主要方法。但运用FFT进行电力谐波分析时很难做到同步采样和整周期截断，造成的频谱泄漏将影响谐波分析精度。为减小谐波分析误差，国内外学者提出了加窗插值的FFT算法,例如矩形窗、Hanning窗、Blackman窗、Rife-Vincent(I)窗和Nuttall窗等，在一定程度抑制了频谱泄漏，改善了谐波分析精度。

在谐波信号离散频谱分析中，对窗函数的要求是主瓣窄、旁瓣低、旁瓣衰减速度快，具体对某一种窗函数，这几个要求相互矛盾，很难同时满足。主瓣与频率分辨率有关，主瓣宽、频率分辨率低；旁瓣与泄漏直接有关，旁瓣峰值电平高，频谱泄漏多；旁瓣衰减速度与长范围泄漏有关，旁瓣衰减速度快，可以有效的抑制长范围泄漏。

在这三项指标中，窗函数最主要的指标是旁瓣峰值电平和旁瓣衰减速度，4项3阶的Nuttall窗的旁瓣峰值电平为-82.6dB，而其旁瓣衰减速度达到30dB/oct，与其它窗函数相比，Nuttall窗是旁瓣性能优良的窗函数，它们的旁瓣特性特别适合于周期谐波信号的频谱分析，因此提出一种基于Nuttall窗的时移相位差稳态谐波信号校正方法，将该方法能够显著地提高计算精度。

发明内容

有鉴于此，本发明的目的在于提供一种基于Nuttall窗的时移相位差稳态谐波信号校正方法，该方法结合相位差校正原理，采用Nuttall窗对谐波信号进行加权，该算法能很好地抑制由于频谱泄漏引起的基波及整数次谐波之间的相互干扰。

为达到上述目的，本发明提供如下技术方案：

步骤一：谐波信号x(t)为以采样频率f_s将x(t)离散化得序列x(n)，其中f_m为频率，A_m为幅值，为初相位，P为最高谐波次数，m＝1,2…p，Ω为模拟角频率，ω为数字角频率，ω_m＝Ω_mT_s，T_s为采样周期，T_s＝1/f_s；

步骤二：x₁(n)、x₂(n)为序列x(n)中的两段序列；用长度为N的Nuttall窗序列w_Nuttall(n)对x₁(n)、x₂(n)加权截断，得离散加窗信号x_1w(n)、x_2w(n)，x_1w(n)＝x₁(n)·w_1Nuttall(n)，x_2w(n)＝x₂(n)·w_2Nuttall(n)，n＝0,1,2,…,N-1；

步骤三：FFT计算出各频点处对应谱线的幅值和相角；

步骤四：搜索出x₁(n)和x₂(n)各次谐波的峰值谱线；

步骤五：通过以下公式，求出x₁(n)和x₂(n)两段序列峰值谱线处的相位差，

步骤六：通过以下公式，求出各次谐波的频率校正量δ，

其中，m＝1,2…p；L_m＝m·L，L为最接近NT_s/T₁的正整数；N为时间窗t_p内采样点数；

步骤七：通过以下公式，分别求出各次谐波的频率、幅值、相位，

$f_m=\frac{{\mathrm{\ω}}_m}{2\mathrm{\π}{T}_s}=\frac{(L_m+{\mathrm{\δ}}_m)\mathrm{\Δ\ω}}{2\mathrm{\π}{T}_s}=\frac{(L_m+{\mathrm{\δ}}_m)}{2\mathrm{\π}{T}_{s}N}=(L_m+{\mathrm{\δ}}_m)\mathrm{\Δf};$

其中，△ω＝2π/N，△f＝1/NT_s

$A_m=\frac{2\left|X_w\left(L_m\right)\right|}{\underset{\mathrm{\η}=0}{\overset{k}{\mathrm{\Σ}}}\frac{\left|{\mathrm{\α}}_{\mathrm{\η}}\right|}{2}\left\{W_R\right[\frac{2\mathrm{\π}}{N}(\mathrm{\η}-{\mathrm{\λ}}_m)]+W_R[\frac{2\mathrm{\π}}{N}(-\mathrm{\η}-{\mathrm{\λ}}_m)\left]\right\}};$

其中，a_η为Nutall四项三阶窗函数的系数，λ_m为各次谐波的频谱偏差量；

其中，

进一步，所述步骤二中的Nuttall窗为一种余弦组合窗，其时域表达式为其中n＝0,1,2,…,N-1，M为窗函数的项数；b_m为窗函数中系数且满足条件

本发明的有益效果在于：本发明提出的一种基于Nuttall窗的时移相位差稳态谐波信号校正方法，该算法能很好地抑制由于频谱泄漏引起的基波及整数次谐波之间的相互干扰。与谐波分析的比值校正法和加其它窗函数相比，信号截断时间相同时，基于Nuttall窗的时移相位差校正法可以获得更高的计算精度。

附图说明

为了使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明作进一步的详细描述，其中：

图1为本发明所述流程图；

图2为谐波幅值相对误差比较；

图3为谐波相位相对误差比较。

具体实施方式

下面将结合附图，对本发明的优选实施例进行详细的描述。

本发明提供的一种基于Nuttall窗的时移相位差稳态谐波信号校正方法，该方法包括以下步骤：

步骤一：谐波信号以采样频率f_s将x(t)离散化得序列x(n)，

步骤二：x₁(n)、x₂(n)为序列x(n)中的两段序列；用长度为N的Nuttall窗序列w_Nuttall(n)对x₁(n)、x₂(n)加权截断，得离散加窗信号x_1w(n)、x_2w(n)，x_1w(n)＝x₁(n)·w_1Nuttall(n)，x_2w(n)＝x₂(n)·w_2Nuttall(n)，n＝0,1,2,…,N-1；

步骤三：FFT计算出各频点处对应谱线的幅值和相角；

步骤四：搜索出x₁(n)和x₂(n)各次谐波的峰值谱线；

步骤五：求出x₁(n)和x₂(n)两段序列峰值谱线处的相位差；

步骤六：求出各次谐波的频率校正量δ；

步骤七：分别求出各次谐波的频率、幅值、相位。

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

$w_H\left(n\right)=\underset{m=0}{\overset{M-1}{\mathrm{\Σ}}}{(-1)}^m{b}_m\cos \left(\frac{2\mathrm{\π}}{N}\mathrm{mn}\right),n=\mathrm{0,1,2},...,N-1---\left(1\right)$

式中，M为窗函数的项数；窗函数中系数b_m应满足约束条件

典型的Nuttall窗函数的系数如表1所示，表2给出了Nuttall窗的旁瓣特性。

表1 Nuttall窗的系数

表2 Nuttall窗的旁瓣特性

选用旁瓣峰值电平小且旁瓣渐进衰减速率大的窗函数，可很好地抑制邻近泄漏和远离泄漏，提高谐波分析的准确度。由表1可见，在上述几种的余弦组合窗中，4项3阶Nuttall窗具有理想的旁瓣特性，旁瓣峰值电平为-82.6dB，旁瓣衰减速率为30dB/oct。

Nuttall窗谐波信号

设一频率为f_m、幅值为A_m、初相位为最高谐波次数为P的谐波信号x(t)为

以采样频率f_s将上式离散化得序列x(n)

其中Ω为模拟角频率，ω为数字角频率，ω_m＝Ω_mT_s，T_s＝1/f_s为采样周期；x(n)的频谱为

用长度为N的Nuttall窗序列w_Nuttall(n)对x(n)加权截断，得离散加窗信号x_w(n)

x_w(n)＝x(n)·w_Nuttall(n)  n＝0,1,2,…,N-1   (5)

根据频域卷积定理，时域相乘对应于频域卷积，因此，加Nuttall窗信号的x_w(n)的DTFT为

x_w(n)的ω≥0部分的频谱分量为

对信号进行DFT变换求离散频谱X_w(k)，相当于在DTFT频域X_w(e^jω)中以△ω＝2π/N(对应的△f＝f_s/N＝1/NT_s)的间隔抽样：

考虑采样不同步，即时间窗t_p＝NT_s不为信号基波周期T₁(T₁＝1/f₁)的整数倍，设

$\mathrm{\λ}=\frac{t_p}{T_1}=\frac{{\mathrm{NT}}_s}{T_1}=L+\mathrm{\δ}---\left(9\right)$

其中，L为最接近NT_s/T₁的正整数，δ为由非同步采样造成的频率偏差。

又△ω＝2π/N则结合(9)式得

$\frac{{\mathrm{\ω}}_m}{\mathrm{\Δ\ω}}=\frac{{\mathrm{\Ω}}_m{T}_s}{\mathrm{\Δ\ω}}=\frac{2\mathrm{\π}{f}_m{T}_s}{2\mathrm{\π}/N}=\frac{{\mathrm{NT}}_s}{m{T}_1}=\mathrm{mL}+\mathrm{m\δ}=L_m+{\mathrm{\δ}}_m---\left(10\right)$

其中L_m＝m·L，δ_m＝m·δ。

对于m＝1,2…p，△f＝1/NT_s，可得第m次谐波的频率为

$f_m=\frac{{\mathrm{\ω}}_m}{2\mathrm{\π}{T}_s}=\frac{(L_m+{\mathrm{\δ}}_m)\mathrm{\Δ\ω}}{2\mathrm{\π}{T}_s}=\frac{(L_m+{\mathrm{\δ}}_m)}{2\mathrm{\π}{T}_{s}N}=(L_m+{\mathrm{\δ}}_m)\mathrm{\Δf}---\left(11\right)$

由上式得

${\mathrm{\ω}}_m=2\mathrm{\π}{f}_m{T}_s=2\mathrm{\π}(L_m+{\mathrm{\δ}}_m)\frac{T_s}{{\mathrm{NT}}_s}=\frac{2\mathrm{\π}}{N}(L_m+{\mathrm{\δ}}_m)---\left(12\right)$

又

$\mathrm{\ω}-{\mathrm{\ω}}_m=k\frac{2\mathrm{\π}}{N}-\frac{2\mathrm{\π}}{N}(L_m+{\mathrm{\δ}}_m)=-{\mathrm{\δ}}_m\frac{2\mathrm{\π}}{N}---\left(13\right)$

式中，k与L_m表示各次谐波的峰值谱线，所以它们的值应是相等的。

将(13)式代入(4)式得

$\begin{array}{c}{W}_{\mathrm{Nuttall}}\left(e^{j(\mathrm{\ω}-{\mathrm{\ω}}_m)}\right)=\\ \underset{\mathrm{\η}=0}{\overset{3}{\mathrm{\Σ}}}\frac{\left|a_{\mathrm{\η}}\right|}{2}\left\{W_R^0\right[\frac{2\mathrm{\π}}{N}(\mathrm{\η}-{\mathrm{\δ}}_m)]+W_R^0[\frac{2\mathrm{\π}}{N}(-\mathrm{\η}-{\mathrm{\δ}}_m)\left]\right\}{e}^{j\frac{N-1}{N}{\mathrm{\π\δ}}_m}\end{array}$

则由(8)式和(14)式，可得第m次谐波的幅值

$A_m=\frac{2\left|X_w\left(L_m\right)\right|}{\underset{\mathrm{\η}=0}{\overset{\mathrm{\κ}}{\mathrm{\Σ}}}\frac{\left|{\mathrm{\α}}_{\mathrm{\η}}\right|}{2}\left\{W_R\right[\frac{2\mathrm{\π}}{N}(\mathrm{\η}-{\mathrm{\λ}}_m)]+W_R[\frac{2\mathrm{\π}}{N}(-\mathrm{\η}-{\mathrm{\λ}}_m)\left]\right\}}---\left(15\right)$

其中，a_η为Nutall四项三阶窗函数的系数，λ_m为各次谐波的频谱偏差量；

将(14)式代入(8)式，得

取上式两边的相角，可得第m次谐波的相位

由(11)式、(15)式、(17)式，可得基波及各次谐波的频率、幅值和相位信息，而从上述各式的推导中可以看出，最重要的是频率偏差量的求取方法。

求取频率偏差的方法有：多点卷积幅值校正法、能量重心校正法、峰值搜寻法、比值公式法和相位差法等，前4种在估计某次谐波参数时，均要求在该次谐波主瓣宽度内的两根或两根以上的谱线，其它次谐波及该次谐波自身负频分量的频谱泄漏均为1或很小，这样的条件较为严格。对于普通的余弦窗，只能通过加长时间来满足，因此增加了计算量和响应时间。而相位差校正法只要求该次谐波幅值最大的谱线处其它谐波对其干扰为0或很小，此条件在加余弦窗时较易满足。

取时间窗t_p＝τT₀(τ为所取的工频周期数取为正整数,T₀为工频周期0.02S)，在时间窗t_p内采样点数为N，则频率分辨率△f＝1/t_p＝1/τT₀，采样间隔T_s＝τT₀/N。

以T_s对谐波信号等间隔采样点，取0至N-1点构成时域序列x₁(n)，取第至点构成时域序列x₂(n)，则序列x₂(n)比序列x₁(n)滞后的时间为x₂(n)、x₁(n)对应的频域初相角的关系为

将(11)式代入上式，得

分别对x₁(n)和x₂(n)加Nuttall窗后均作N点DFT，各次谐波对应的第L_m条谱线，由(17)式得

式(20)和式(21)相减，得

根据(19)式和(22)式，得频率校正量

将δ_m代入(11)式、(15)式、(17)式，即可得基波及各次谐波的频率、幅值和相位信息。

实施例

采用的谐波信号模型为

式中：基波频率f₁为50.5Hz，A_m和分别为第m次此而谐波的幅值和初相角，采样频率为5kHz,截断信号的第一段和第二段的数据长度N均为1024点，约10个周期的数据，基波和各次谐波的幅值和相角如表3。

表3 基波及谐波参数

选取Nuttall窗长度M＝1024，分别采用3项1阶、3项最小旁瓣、4项1阶、4项3阶Nuttall窗，数据长度为N＝2M＝2048，幅值和相位的误差仿真结果如表4、表5。

表4 幅值相对误差比较(百分比)

表5 相位相对误差比较(百分比)

表中E_A表示基波及各次谐波的幅值测量值相对于真值的误差百分比；表示基波及各次谐波的相位测量值相对于真值的误差，图2和图3为误差曲线，可以看出采用4项1阶、4项3阶Nuttall窗函数时的计算结果比采用3项1阶、3项最小旁瓣Nuttall窗函数时的计算结果具有更高的精度和稳定性。

本发明提供的一种基于Nuttall窗的时移相位差稳态谐波信号校正方法简单，精度较高，适合各种对称窗函数。

最后说明的是，以上优选实施例仅用以说明本发明的技术方案而非限制，尽管通过上述优选实施例已经对本发明进行了详细的描述，但本领域技术人员应当理解，可以在形式上和细节上对其作出各种各样的改变，而不偏离本发明权利要求书所限定的范围。

Claims (2)

1.一种基于Nuttall窗的时移相位差稳态谐波信号校正方法，其特征在于：该方法包括以下步骤：

步骤一：谐波信号x(t)，将x(t)离散化得序列x(n)，其中，f_m为频率，A_m为幅值，为初相位，P为最高谐波次数，m＝1,2…p，Ω为模拟角频率，ω为数字角频率，ω_m＝Ω_mT_s，T_s为采样周期，f_s为采样频率，T_s＝1/f_s；

步骤二：x₁(n)、x₂(n)为序列x(n)中的两段序列，用长度为N的Nuttall窗序列w_Nuttall(n)对x₁(n)、x₂(n)加权截断，得离散加窗信号x_1w(n)、x_2w(n)，x_1w(n)＝x₁(n)·w_1Nuttall(n)，x_2w(n)＝x₂(n)·w_2Nuttall(n)，n＝0,1,2,…,N-1；

步骤三：FFT计算出各频点处对应谱线的幅值和相角；

步骤四：搜索出x₁(n)和x₂(n)各次谐波的峰值谱线L_m；

步骤五：通过以下公式，求出x₁(n)和x₂(n)两段序列峰值谱线处的相位差，

步骤六：通过以下公式，求出各次谐波的频率校正量δ_m，

其中，m＝1,2…p；L_m＝m·L，L为最接近NT_s/T₁的正整数；N为时间窗t_p内采样点数；

步骤七：通过以下公式，分别求出各次谐波的频率、幅值、相位，

$f_m=\frac{{\mathrm{\ω}}_m}{2\mathrm{\π}{T}_s}=\frac{(L_m+{\mathrm{\δ}}_m)\mathrm{\Δ\ω}}{2\mathrm{\π}{T}_s}=\frac{(L_m+{\mathrm{\δ}}_m)2\mathrm{\π}}{2\mathrm{\π}{T}_{s}N}=(L_m+{\mathrm{\δ}}_m)\mathrm{\Δf},$

其中，△ω＝2π/N，△f＝1/NT_s；

$A_m=\frac{2\left|X_w\left(L_m\right)\right|}{\underset{\mathrm{\η}=0}{\overset{k}{\mathrm{\Σ}}}\frac{\left|a_{\mathrm{\η}}\right|}{2}\left\{W_R\right[\frac{2\mathrm{\π}}{N}(\mathrm{\η}-{\mathrm{\λ}}_m)]+W_R[\frac{2\mathrm{\π}}{N}(-\mathrm{\η}-{\mathrm{\λ}}_m)\left]\right\}},$

其中，a_η为Nutall四项三阶窗函数的系数，λ_m为各次谐波的频谱偏差量；

其中，

2.根据权利要求1所述的一种基于Nuttall窗的时移相位差稳态谐波信号校正方法，其特征在于：所述步骤二中的Nuttall窗为一种余弦组合窗，其时域表达式为其中n＝0,1,2,…,N-1，M为窗函数的项数；b_m为窗函数中系数且满足条件 $\underset{m=0}{\overset{M-1}{\mathrm{\Σ}}}{(-1)}^m{b}_m=0.$