混合自卷积窗谐波算法研究与FPGA方案

赖增强，王贵忠，肖智宏，张国庆，李洪波，刘颖，关铁成

(1.哈尔滨工业大学 电气工程及自动化学院，黑龙江 哈尔滨 150090；2.哈工大(张家口) 工业技术研究院，河北 张家口 075421；3.国网经济技术研究院有限公司，北京 102209；4.国网营口供电公司，辽宁 营口 115002)

随着电力系统电力电子化进程的推进，大量电力电子装置的并网使用使得电网中的谐波检测与治理问题得到人们的广泛关注，谐波分析在电力系统中显得尤为重要。电力谐波分析算法一般分为如下4种[1-2]：①时域仿真分析法；②频域仿真分析法；③基于数学变换理论算法；④基于人工智能算法。迄今为止文献数量最多、研究最深入、运用最广泛当属基于数学变换理论算法，它主要包括快速傅里叶变换(fast Fourier transform, FFT)、小波变换分析法、矢量分析法、Prony法、广义S变换法以及相关改进算法。其中，FFT是把时域与频域相关联的算法，具有完备、正交等众多优势。尽管FFT不如小波变换那样特别适合处理突变和不平稳信号的分析，但它的改进算法因其计算方便特别适合现场可编程门阵列(field-programmable gate array, FPGA)等硬件编程实现，是电能质量装置中谐波分析的首选算法。然而FFT不足之处是非同步采样和非整数周期截断时易出现“频谱泄漏”与“栅栏效应”，致使不能精确地测量出电网的谐波数据[3]；同时对于电网的频率波动，即使是锁相环也难以实现电网的同步采样[4]，因此需要采用一种性能优良的窗函数对FFT算法进行加窗处理来抑制频谱效应。常用的传统窗函数有矩形窗、Hanning窗、Hamming窗和Blackman窗等[5-9]。

电力谐波测量标准要求幅值误差不应大于5%，相角误差不大于5°[10]。为了提高谐波分析精确度，文献[11-12]提出了Nuttall窗及其自卷积窗算法，文献[13]提出了Hanning自卷积窗算法，自卷积窗有利于提高窗的本身优势，但是同时也放大了窗的一些劣势，如Nuttall窗自卷积后主瓣宽度增加。文献[14]提出了不同窗混合卷积方法，即将两种窗的优势互补，进而提高计算精度。文献[15-16]提出了五点加权法，该方法虽精度高但因对FFT输出序列进行加权变换导致计算量大，难以硬件实现。为了得到频谱特性更为优越的窗函数，本文提出一种新型Nuttall混合自卷积窗函数，即先将Nuttall窗与Hanning窗混合卷积后再L阶自卷积，并对其进行双谱线插值与数值拟合修正。该算法利用Nuttall窗旁瓣峰值衰减速度快和Hanning窗主瓣较窄的特点，经自卷积进一步提升优势，从而得到的混合自卷积窗十分适用于谐波分析。最后通过仿真对比，验证混合自卷积窗谐波算法的准确性。同时对算法进行了硬件结构方案设计与FPGA关键技术研究，为谐波分析仪和电能质量装置的研制提供参考。

1 混合自卷积窗设计

1.1 Hanning窗与Nuttall窗

Hanning窗为一种升余弦窗，其离散表达式为

(1)

式中：n为采样样点序列号;N为采样样点总数;RN(n)为矩形窗函数。

式(1)进行离散傅里叶变换(discrete Fourier transform, DFT)后的频域表达式为

(2)

式中WR(·)为矩形窗的DFT表达式。

Nuttall窗也为一种余弦窗，函数表达式为

(3)

式中al为系数。

式(3)进行DFT后的频域表达式为

(4)

由于Nuttall具有多项多阶，即不同的L值具有不同的al系数。Nuttall窗的多项多阶与系数对应关系见表1。

表1 Nuttall窗的多项多阶与系数对应关系Tab.1 Corresponding relation between multiple orders and coefficients of Nuttall window

本文将采用4项3阶Nuttall窗，即

(5)

给出Hanning窗与Nuttall窗的幅频响应曲线，且幅频图中以归一化后的频率做横坐标，以最大值为0分贝值做纵坐标，如图1所示。

图1中，不仅表明了窗函数的幅频关系，更包含了一些十分重要的指标，如主瓣宽度、第一旁瓣衰减以及旁瓣峰值衰减等。根据图1，不难发现Nuttall窗旁瓣峰值衰减速度快及Hanning窗主瓣比较窄。实际上，窗函数第一旁瓣衰减大和旁瓣峰值衰减速度快最能有效抑制频谱效应，但是主瓣宽度较大往往又给信号分析带来不良影响，所以这是一个综合评估并折中选择的过程[17-18]。

-32 dB—Hanning窗第一旁瓣衰减；-83 dB—Nuttall窗第一旁瓣衰减；-18 dB/oct—Hanning窗旁瓣峰值衰减速度；-30 dB/oct—Nuttall窗旁瓣峰值衰减速度，dB/oct是衰减斜率的单位分频斜率，用来反映分频点以下频响曲线的下降斜率,用分贝/倍频程表示，即dB/oct。

图1 Hanning窗与 Nuttall窗的幅频响应对比
Fig.1 Contrast of amplitude frequency response ofHanning and Nuttall window

1.2 混合自卷积窗函数设计

Nuttall混合自卷积窗函数，是先将Nuttall窗与Hanning窗混合卷积，得到混合窗函数

wNHC(n)=wN(n)*wH(n).

(6)

再将得到的wNHC(n)进行L个时域卷积计算，即L阶自卷积，推导出混合自卷积窗函数

(7)

图2给出了不同L阶混合自卷积窗函数的幅频响应曲线。从图2可以看出，随着L的增加，混合自卷积窗函数的主瓣宽度缓慢增宽，但旁瓣特性却迅速提升，所以对于电力谐波信号，选择合适的L阶混合自卷积窗函数对信号加窗，可有效提升谐波的FFT分析精度。结合窗函数选择原则，若无法同时满足主瓣宽度与旁瓣特性，常常以牺牲主瓣的宽度来重视旁瓣特性的要求[19-20]，故以下将选择L=4作为混合自卷积窗函数的阶数来做进一步的研究。

2 混合自卷积窗双谱线插值算法研究

电力信号实际上并非同步采样，而且很可能是非整周波采样，因此在离散频率点上不一定存在峰值频率为基波整数倍的谐波,但双谱线插值算法能很好地解决这一问题[11,21]。

图2 不同L阶混合自卷积窗函数的幅频响应对比Fig.2 Contrast of amplitude frequency ofNHSC window function with different L orders

对时域信号x(t)进行离散采样，则

(8)

式中：f0为频域离散信号x(n)的初始频率；A为x(t)幅值；φ0为x(t)初相角;fs为采样频率;t为时间变量。

对信号x(n)加窗及快速傅里叶变换后经离散时间傅里叶变换(discrete-time Fourier transform，DTFT)离散抽样，若抽样时间间隔为Δf=fs/N,则

(9)

式中：X(kΔf) 为第k条谱线的DTFT的表达式;W(·)为窗函数的DTFT表达式。

设k1、k2分别为峰值点k0左右两侧幅值为最大和次大的谱线，可分别表示为y1=|X(k1Δf)|和y2=|X(k2Δf)|，且满足k1≤k0≤k2=k1+1。同时引入参数α，并令α=k0-k1-0.5，则α的取值范围为[-0.5,0.5]。定义β=(y2-y1)/(y2+y1)，则

(10)

由式(10)定义α和β的函数关系，记做β=g(α)或α=g-1(β)，并得到了频率修正式f0=k0Δf=(α+k1+0.5)Δf，同时定义两谱线加权平均值为幅值修正式，即

(11)

式中A1、A2分别为第k1和k2条谱线所对应的幅值。

当数据阶段长度N很大时，式(11)还可做进一步简化，即

A=N-1(y1+y2)v(α).

(12)

式中v(α)为下文式(14)的偶函数形式。同时，由式(9)，可得初相位修正式为

(13)

综上所述，双谱线插值算法计算繁琐，特别是应用在嵌入式实时监测系统中，加之混合自卷积窗对系统的运算量提出了很高的要求。为此，采用一种多项式逼近的数值分析算法求取α=g-1(β)。设α系列α1,α3,…,α2m+1取值范围为[-0.5,0.5]，代入式(10)，得到与之对应的β值，通过数据拟合从而得到

(14)

同时初相角的修正式也可化简为

(15)

针对混合自卷积窗，为了保证精度的同时提升其计算速度，本文将采用双谱线插值五次多项式逼近的修正式，同时给出修正式α、β的拟合曲线，如图3所示；不仅如此，为了对比分析，也分别给出了性能较优的Hanning窗和Nuttall窗的修正式，即：

a)混合自卷积窗函数及修正式

(16)

b)Hanning窗及修正式

(17)

c)Nuttall窗及修正式

(18)

图3 混合自卷积窗的修正式拟合曲线Fig.3 Correction formula fitting curve of NHSC window

3 仿真验证

为了验证理论分析，采用MATLAB编程仿真，对表2中一组谐波数据信号分别用混合自卷积窗、Nuttall窗和Hanning窗进行处理。通过观察各类窗双谱线插值后的数据结果，对比分析各类窗的优劣。

设信号基波频率为49.5 Hz，采样频率为3 200 Hz。

MATLAB仿真流程如图4所示，其中输入信号为x(n)，分别为混合自卷积窗、Nuttall窗、Hanning窗,同时令基波频率、幅值及相位角分别为f0,A0,φ0，其他各次谐波频率、幅值及相位角分别为fi,Ai,φi。

表2 谐波数据表Tab.2 Harmonic data table

图4 仿真程序流程Fig.4 Simulation program flow chart

通过程序仿真，得到的具体数据见表3—表5，其中相位及幅值误差如图5和图6所示。

表3 混合自卷积窗谐波数据Tab.3 Harmonic data of NHSC window

表4 Hanning窗谐波数据Tab.4 Harmonic data of Hanning window

表5 Nuttall窗谐波数据Tab.5 Harmonic data of Nuttall window

图5 3类窗谐波分析相位误差对比Fig.5 Phase errors of harmonic analysis of three kinds of window functions

从表3—表5数据及图5和图6误差中可以分析得出：无论在相位还是幅值检测上，较之性能较好的Nuttall窗与Hanning窗，混合自卷积窗更能体现出误差小、精度高等强大优势。其中，相位误差稳定在0.000 01°，幅值误差几乎为0，完全满足谐波测量标准要求，该算法对弱幅值谐波也有较强的检测精度。

图6 3类窗谐波分析幅值误差对比Fig.6 Amplitude errors of harmonic analysis of three kinds of window functions

4 硬件结构与FPGA设计方案

FPGA作为集微处理器、存储设备和输入输出接口等模块于一身的系统级集成电路，因其功能强大、保密性能好以及良好的可重复编程等特性，深受工业界的青睐。因此本算法的硬件设计方案采用FPGA作为主控芯片。本算法的硬件全局结构如图7所示。

图7中，标准100 V/5 V电压互感器和5 A/2.5 mA电流互感器分别把三相电压信号和三相电流信号转换为6路小电压及小电流信号，并通过信号调制电路得到适合数模转换的量程；其次FPGA控制数模转换后的数据信号缓存到输入双口数据存储器中，然后运用数据处理模块对数据进行类型转换、混合自卷积窗处理、FFT运算以及有效值处理等，并且整个处理过程是由有限状态机控制，并缓存到输出双口数据存储器中，作进一步处理运算；最后进行存储，并利用串口将处理好的数据输出到上位机液晶屏中显示，从而实现电能质量谐波分析与监测的功能。

基于混合自卷积窗函数双谱线插值算法的硬件移植，包含混合自卷积窗和FFT运算两大部分。主要思路是运用摩尔状态机的控制将256个点的采样数据缓存到双口数据存储器中进行加窗、双峰谱线插值修正和FFT运算，其具体FFT程序控制流程如图8所示。

对于混合自卷积窗及插值修正，可以事先利用MATLAB获得修正后的混合自卷积窗的.mif数据文件，并将其存储到硬件电路的芯片中，可利用Quartus II 13.0设置为数据存储，具体程序架构如图9所示。

在实际Verilog编程中，流水线乘法器均被混合自卷积窗与采集数据所采用，因此可采用乘法单元的知识产权核来设计，便能将此算法在FPGA上实现电力谐波分析。

图7 硬件装置结构Fig.7 Hardware device structure

图8 FFT程序流程Fig.8 FFT program flow chart

图9 混合自卷积加窗插值程序架构Fig.9 NHSC Window interpolation program architecture

5 结束语

本文提出了一种新型的混合自卷积窗，分析了Hanning窗与Nuttall窗各自的优势，且利用二者混合卷积后再L阶自卷积的方式，使得新型窗—混合自卷积窗旁瓣特性大幅度提升，且性能不断优化。在此基础上，因考虑实际检测谐波是非周期且非同步采样，故对混合自卷积窗进行了双谱线插值，进一步提高了FFT算法的准确度，且抑制了一些频谱效应。对所采用的算法进行对比仿真验证，结果表明本文提出的基于混合自卷积窗的FFT插值修正算法误差小、精度高，适合电力系统谐波分析。最后给出了FPGA的硬件设计方案，运用MATLAB获取一系列窗的数值表，与采集后离散化数据做乘法加窗处理，在保证算法计算速度的同时，比起传统FFT算法，提高了测量精度，为相关谐波分析仪与电能质量装置的研制提供了一定的理论与技术支撑。