CN105717359A - 一种谐波分析算法 - Google Patents

Abstract

本发明公开了一种谐波分析算法,为增加其旁瓣抑制能力，本算法将Hamming窗特征系数a修改为25/46≈0.543478。本算法主要实现方法如下：构建长度为M的改进Hamming窗序列,M为2的整数次幂。将改进Hamming窗序列做自卷积，自卷积后在所得序列末尾补0。得到长度为N＝2M的2阶改进Hamming自卷积窗序列。对输入信号进行采样，采用长度为N的2阶改进Hamming自卷积窗对采样所得的无限长度的离散信号加窗截断。对加窗截断后的信号进行FFT计算，获得离散频谱。对离散频谱进行插值处理，计算频率偏移量。根据频率偏移量计算谐波频率、幅值和相角。

Description

一种谐波分析算法

技术领域

本发明涉及信号处理领域,具体涉及一种谐波分析算法。

背景技术

目前随着工业负荷的不断发展及电力电子设备的大量应用，非线性负荷己经成为电力系统负荷的重要组成部分。据统计，目前已有20％的电力负荷通过各种形式的功率交换来实现。大量非线性负荷的投入将给电网带来严重的谐波污染、影响电能质量。基于电力电子装置的新型有源滤波器是治理谐波的重要手段。然而要使有源滤波器等电力电子控制设备良好工作，实现对电力系统谐波的快速精确测量是其重要的前提。

FFT算法，即快速离散傅里叶算法是目前计算电力系统谐波的主要方法。常规FFT算法在计算电力系统谐波时，要求实现整周期采样，否则测量结果将受到栅栏效应和频谱泄露的影响。然而工程实际中，由于电网频率波动等因素的影响，完全严格的整周期采样难以实现，此时普通FFT算法将难以实现电力系统谐波的精确测量。

现阶段减少频谱泄漏和栅栏效应的主要方法是使用加窗插值算法。加窗插值算法通过加窗来抑制由旁瓣造成的谱间干扰；通过插值算法来计算被测信号在FFT计算点X(K)、X(K+1)之间的频率分量，从而解决栅栏效应。

综上所述，提供一种更精准的谐波分析算法是亟需解决的问题。

发明内容

本发明提供了一种谐波分析算法，能够实现对谐波的精确和快速检测。

本发明公开了如下技术方案：

一种谐波分析算法，包括步骤如下：

采集离散电学信号构成原始输入序列X[n]；

构建长度为M的Hamming窗序列，其中M为2的整数次幂，所述Hamming窗的系数a＝0.543478；

将所述Hamming窗序列做一次自卷积，并在卷积序列的末尾补一个0，得到长度为N的2阶离散Hamming自卷积窗，其中N＝2M；

使用所述2阶离散Hamming自卷积窗对所述原始输入序列X[n]加窗截断，并进行FFT计算，得到离散序列y[n]，获得谐波频率、幅值和相角。

优选的，在上述谐波分析算法中，所述使用所述2阶离散Hamming自卷积窗对所述原始输入序列加窗截断，并进行FFT计算，得到离散序列y[n]，获得谐波频率、幅值和相角具体是指：

进行FFT计算，获得离散频谱，对所述离散频谱进行插值处理，计算频率偏移量；

根据频率偏移量计算基频频率，根据基频频率计算各次谐波分量所在范围；

计算与各次谐波分量相邻的FFT计算点X(K)、X(K+1)，得到各次谐波的相位和幅值、相角。

由上述技术方案可见，本发明对Hamming窗的系数a的大小进行了改进，这是因为，谐波测量仪器对采样所得离散电学信号进行截断处理会造成频谱泄露。时域截断即将原输入序列X(n)乘以矩形窗序列，在频域原信号频谱X(w)将和窗频谱R(w)发生卷积运算。由于窗函数不可能无限宽，即R(w)不是冲激函数，在卷积运算后，原信号的频谱将发生展宽和拖尾，出现原来不存在的频率分量，造成频谱泄露和谱间干扰。而当整周期采样时，输入信号各频率分量旁瓣刚好在其它频率分量处的值为零。故当整周期采样时，各次谐波测量不受频谱泄露影响。

所以，加窗插值算法选择窗函数时，为了抑制非整周期采样过程中的频谱泄露，往往希望窗函数主瓣尽可能的窄，旁瓣衰减尽可能高。两项余弦窗的表达式如下式所示。

$w\left(n\right)=a-(1-a)\×\cos \left(\frac{2\mathrm{\π}n}{N}\right)$

其频域表达式幅值部分如下所示：

$\left\{\begin{array}{c}W\left(w\right)=a\×R_N\left(w\right)-0.5\×(1-a)\×\left(R_N\right(w-\frac{2\mathrm{\π}}{N})+R_N(w+\frac{2\mathrm{\π}}{N}\left)\right)\\ R_N\left(w\right)=\frac{\sin (wN/2)}{\sin (w/2)}\end{array}\right.$

其中R_N(w)为矩形窗，通常Hanning窗的系数a＝0.5，Hamming窗的系数a＝0.54。现为了增加窗函数对最大旁瓣的抑制能力，令原本最大傍瓣，即w＝5π/N处的W(w)为零，由此通过解方程可以得到傍瓣抑制能力最强的两项窗，此时a＝25/46≈0.543478。

为了进一步增强窗函数的傍瓣抑制能力，本方法将改进Hamming在时域做自卷积处理。根据傅里叶变换的有关知识，在时域进行卷积计算，则在频域将进行相乘运算。若改进Hamming在时域做n次自卷积运算，则其傍瓣衰减速率(dB/oct)将提升n倍。然后同时其主瓣宽度也将增加n倍。当主瓣宽度增加时，为了防止主瓣造成谱间干扰，采样周期数将增加，这会影响算法的实时性。故综合考虑主瓣、旁瓣两方面影响，本发明对改进Hamming窗做一次自卷积运算。设改进Hamming窗的长度为M点，且在自卷积过后补零一点，最终得到长度为N(N＝2M)的改进Hamming自卷积窗。2阶卷积在减小傍瓣峰值、增加傍瓣衰减速率的同时，其主瓣宽度不至于过宽，保持了算法的实时性。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，对于本领域普通技术人员而言，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明提供的谐波分析算法的一个实施例。

具体实施方式

本发明提供的一种谐波分析算法，能够实现对谐波的精确和快速检测。

为了使本技术领域的人员更好地理解本说明中的技术方案，下面将结合本说明实施例中的附图，对本说明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本说明一部分实施例，而不是全部的实施例。基于本说明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都应当属于本说明保护的范围。

请参考图1，图一为本发明提供的谐波分析算法的一个实施例，如图1所示，

在步骤S01中，设置采样频率Fs＝3200Hz,对离散电学信号进行采样，得到原始输入序列x[n]。根据采样定理该采样频率最高可以测量31次谐波。

在步骤S02中，构建长度为256点的Hamming窗序列，改进Hamming窗函数的系数a＝0.543478。

在步骤S03中，对改进后的Hamming窗进行卷积运算，并补零，得到长度为512点的改进二阶Hamming自卷积窗W(n)。

在步骤S04中，使用改进后的Hamming自卷积窗函数W(n)对原始输入序列x[n]进行加窗截断:Wp(n)＝X(n)×W(n),得到截断后序列X0(n)。

在步骤S05中，对截断后序列X0(n)进行FFT计算，得到序列y[n],y[n]＝FFT[X(n)×W(n)]。

在步骤S06中，判断FFT计算点y(8)和y(10)的绝对值幅值大小，得到基波频率。这么做的原因是，因为电网基频频率一般在49Hz到51Hz之间波动，故电网基波在FFT计算点y[8]至y[10]之间。比较y[8]、y[10]幅值。若两种相等，则认为此时电网基波频率为50Hz；若y[8]幅值大于y[10]幅值，则认为基波频率在y[8]、y[9]之间；若y[8]幅值小于y[10]幅值，则认为基波频率在y[10]、y[9]之间。基波频率的计算如下：

当abs(y(10))＞abs(y(8))时，按照如下方式计算β，β＝(abs(y(10))-abs(y(9))/(abs(y(10))+abs(y(9)),当abs(y(10))＜abs(y(8))时，按照如下方式计算β，β＝(abs(y(9))-abs(y(8))/(abs(y(9))+abs(y(8))，然后均按照如下计算频率偏移量α＝α₇β⁷+α₅β⁵+α₃β³+α₁β¹，然后再分别按照如下方式计算基频频率f：前者f＝(8+0.5+a)×6.25，后者为f＝(7+0.5+a)×6.25。

在步骤S07中，根据基波频率，计算各次谐波分量所在范围。

在使用窗函数抑制傍瓣的干扰后，本算法使用插值算法消除栅栏效应的影响，并计算其实电网波动后的频率。FFT算法是在频域上对信号频谱进行采样，因此使用FFT计算信号频谱只能得到间隔为F的离散频谱。由于FFT计算所得频谱是离散的，若输入信号频率在FFT计算频率间隔处，则FFT算法无法直接计算该输入信号频谱。通常文献把这种效应称为栅栏效应。

假设输入周期信号的一个频率分量为：x(t)＝Ae^j2πft。对x(t)以采样频率fs进行采样，并用改进Hamming自卷积窗截断。FFT频率间隔为F，F＝fs/N。设该输入信号分量频率f＝(L+r)F，其中L为正整数，0≤r<1，即f在FFT计算频率X(L)，X(L+1)之间。则:

$\left\{\begin{array}{c}\left|X\left(l\right)\right|=A(a\frac{\sin \left(\frac{r\mathrm{\&pi;}}{2}\right)}{\sin \left(\frac{r\mathrm{\&pi;}}{N}\right)}+\frac{1-a}{2}{(\frac{\sin \left(\frac{(r+2)\mathrm{\&pi;}}{2}\right)}{\sin \left(\frac{(r+2)\mathrm{\&pi;}}{N}\right)}+\frac{\sin \left(\frac{r(r-2)\mathrm{\&pi;}}{2}\right)}{\sin \left(\frac{(r-2)\mathrm{\&pi;}}{N}\right)})}^2\\ \left|X(l+1)\right|=A(a\frac{\sin \left(\frac{(r-1)\mathrm{\&pi;}}{2}\right)}{\sin \left(\frac{(r-1)\mathrm{\&pi;}}{N}\right)}+\frac{1-a}{2}{(\frac{\sin \left(\frac{(r+1)\mathrm{\&pi;}}{2}\right)}{\sin \left(\frac{(r+1)\mathrm{\&pi;}}{N}\right)}+\frac{\sin \left(\frac{r(r-3)\mathrm{\&pi;}}{2}\right)}{\sin \left(\frac{(r-3)\mathrm{\&pi;}}{N}\right)})}^2\end{array}\right.$

设则β可通过FFT计算得到，且由频率偏移量r唯一决定，故可以使用β计算频率偏移量r。为了便于嵌入式设备计算，本发明使用多项式拟合的方法计算r。为了减少计算量，本发明引入中间变量α＝r-0.5。此时α是β的奇函数，故多项式拟合只有奇数项。多项式拟合结果如下：

α＝2.46354451884363*+0.30298811703459*β^3+0.13630447415588*β^5+0.08610324419670*β^7。此时r＝α+0.5，f＝(l+α+0.5)×fs/N。

在步骤S08中，计算与各次谐波分量范围相邻的FFT计算点y(l)和y(l+1)，得到各次谐波相位、幅值。

在计算得到频率偏移量r后，将r带回|x(l)|、|x(l+1)|的表达式，就可以分别根据|x(l)|、|x(l+1)|计算得到频率分量x(t)＝Ae^j2πft的幅值A1、A2。本发明取A1、A2的加权平均值A为最终的计算结果。因为FFT把各次谐波分解成正负两个频率分量，故谐波幅值A谐波＝2A。

本发明各次谐波相位可如下计算：Φ＝arg(X(l))+π/2-πr。

以上所述仅是本发明的具体实施方式，使本领域技术人员能够理解或实现本发明。对这些实施例的多种修改对本领域的技术人员来说将是显而易见的，本文中所定义的一般原理可以在不脱离本发明的精神或范围的情况下，在其它实施例中实现。因此，本发明将不会被限制于本文所示的这些实施例，而是要符合与本文所公开的原理和新颖特点相一致的最宽的范围。

以上所述仅是本发明的具体实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (2)

1.一种谐波分析算法，其特征在于，包括步骤如下：

采集离散电学信号构成原始输入序列X[n]；

构建长度为M的Hamming窗序列，其中M为2的整数次幂，所述Hamming窗的系数a＝0.543478；

将所述Hamming窗序列做一次自卷积，并在卷积序列的末尾补一个0，得到长度为N的2阶离散Hamming自卷积窗，其中N＝2M；

使用所述2阶离散Hamming自卷积窗对所述原始输入序列X[n]加窗截断，并进行FFT计算，得到离散序列y[n]，获得谐波频率、幅值和相角。

2.根据权利要求1所述的谐波分析算法，其特征在于，所述使用所述2阶离散Hamming自卷积窗对所述原始输入序列加窗截断，并进行FFT计算，得到离散序列y[n]，获得谐波频率、幅值和相角具体是指：

进行FFT计算，获得离散频谱，对所述离散频谱进行插值处理，计算频率偏移量；

根据频率偏移量计算基频频率，根据基频频率计算各次谐波分量所在范围；

计算与各次谐波分量相邻的FFT计算点X(K)、X(K+1)，得到各次谐波的相位和幅值、相角。