CN106053937A - 一种基于fft+ft的基波测频方法 - Google Patents

Abstract

本发明公开了一种基于FFT+FT的基波测频方法，通过先对离散序列添加窗函数，再基于FFT+FT原理对离散序列进行傅里叶变换，最后通过0.618算法搜索基于FFT+FT函数的最大幅值来得到基波频率值；在具体配置中，本发明选用Blackman窗函数，具有精度很高，且是单峰函数，可通过一维搜索算法，减少了计算量。

Description

一种基于FFT+FT的基波测频方法

技术领域

本发明属于电网测频技术领域，更为具体地讲，涉及一种基于FFT+FT的基波测频方法。

背景技术

随着我国经济的快速发展，现代电能是具有经济适用、清洁等特性的二次能源，既影响着社会进步、国民经济发展、国家战略的实施。电力系统频率一方面作为衡量电能质量的指标，需加以动态监测；另一方面作为实施安全稳定控制的重要状态反馈量，要求能实时重构，因此频率的准确测量至关重要。

目前测量电网基波频率的算法很多，主要使用的有零交法，通过计算信号的相邻过零点个数，计算基波周期，这种方法实现简单，在没有噪声和谐波的影响下，零点发生偏移，精度有限；其次，对信号进行FFT(DFT)计算得到基波频率，这种方法由于不能保证是同步采样，故存在严重的频谱泄露和栅栏效应，精度较低；另外，数学模型法主要有最小二乘法，可对采样数据进行数据拟合，但对于谐波和尖峰信号会产生较大误差，若进行滤波，成本增加而且滤波会造成延时。针对目前测频算法的现状，本文基于FFT+FT原理提出了一种的新的测频方法，并对该算法进行了优化，使改进的算法具有高精度、运算快、实时性好特性。

发明内容

本发明的目的在于克服现有技术的不足，提供一种基于FFT+FT的基波测频方法，通过0.618算法搜索基于FFT+FT函数的最大幅值来得到基波频率值，具有高精度、运算快、实时性好等特点。

为实现上述发明目的，本发明为一种基于FFT+FT的基波测频方法，其特征在于，包括以下步骤：

(1)、对采样点为N的离散序列添加Blackman窗，并标记为序列x(n)；

(2)、基于FFT+FT原理对序列x(n)进行傅里叶变换，得：

$X\left(f\right)=\underset{n=0}{\overset{N-1}{\Σ}}x\left(n\right)e^{-j2\mathrm{\π}n\frac{f}{f_s}}$

其中，f为连续频率，f_s为采样率；

(3)、设置初始参数值

设置三相电能质量分析仪的基波频率测量区间为[f₁,f₂]Hz，测量精度在cHz以内；再根据该设置将0.618算法的边界值初值设置为a＝f₁和b＝f₂，精度设置为L，L∈c；设置试探点λ＝a+0.382(b-a)和μ＝a+0.618(b-a)；

(4)、通过0.618算法搜索X(f)的最大幅值来得到基波频率值

将试探点λ和μ分别代入X(f)，计算出X(λ)和X(μ)，再对X(λ)和X(μ)求幅值，得到f(λ)和f(μ)；

当b-a＞L时，若f(λ)＞f(μ)时，则令b＝μ，μ＝λ，λ＝a+0.382(b-a)，再代入到X(f)进行重复迭代；

当b-a＞L时，若f(λ)＜f(μ)时，则令a＝λ，λ＝μ，μ＝a+0.618(b-a)，再代入到X(f)进行重复迭代；

当b-a＜L时，则停止计迭代，得到基波频率值为(a+b)/2。

本发明的发明目的是这样实现的：

本发明一种基于FFT+FT的基波测频方法，通过先对离散序列添加窗函数，再基于FFT+FT原理对离散序列进行傅里叶变换，最后通过0.618算法搜索基于FFT+FT函数的最大幅值来得到基波频率值；在具体配置中，本发明选用Blackman窗函数，具有精度很高，且是单峰函数，可通过一维搜索算法，减少了计算量。

同时，本发明一种基于FFT+FT的基波测频方法还具有以下有益效果：

(1)、精度高

对信号不加噪声进行最大幅值搜索得到的基波频率值，与真值的差值几乎可以忽略。

(2)、抗噪能力强

对信号加载40dB的高斯白噪声，在不需要滤波的情况下基波测量值与真值的误差在0.0006Hz以内，加载随机噪声误差在0.006Hz以内，测量精度可达到基本仪器的测量要求，因此抗噪能力强。

(3)、抗谐波干扰

对电网信号加载高含量的谐波信号与随机噪声，在不需要滤波的情况下基波测量值与真值的误差在0.008Hz以内，测量精度可达到基本仪器的测量要求，具有抗谐波干扰的特性。

附图说明

图1是本发明基于FFT+FT的基波测频方法流程图；

图2是加不同窗函数计算基波频率与真值的差值图；

图3是基波42.5Hz处加三种窗函数的FFT+FT算法仿真图；

图4是基波57.5Hz处加三种窗函数的FFT+FT算法仿真图；

图5是对输入信号加载40dB的高斯白噪声在不同频率的理论与测量误差图；

图6是对输入信号加随机噪声在不同频率的理论与测量误差图；

图7是对输入信号加高含量谐波在不同频率的理论与测量误差图。

具体实施方式

下面结合附图对本发明的具体实施方式进行描述，以便本领域的技术人员更好地理解本发明。需要特别提醒注意的是，在以下的描述中，当已知功能和设计的详细描述也许会淡化本发明的主要内容时，这些描述在这里将被忽略。

实施例

为了方便描述，先对具体实施方式中出现的相关专业术语进行说明：

Blackman：布拉克曼窗；

FFT+FT离散频谱校正法；

图1是本发明基于FFT+FT的基波测频方法流程图。

在本实施例中，本发明一种基于FFT+FT的基波测频方法是针对三相电能质量分析仪进行基波频率测量，因此，需要设置三相电能质量分析仪的初始参数值，即：测频范围为42.5～57.5Hz，频率精度在±0.01Hz以内，下面对基波测频的具体方法进行详细说明，包括以下步骤：

S1、对采样点为N的离散序列添加窗函数，并标记为序列x(n)；

离散的频谱分析在时域被截断为有限长度，必然会产生频谱的能量泄露现象，而能量泄露与加窗函数的旁瓣息息相关，若旁瓣的幅度趋向于零，便可使能量相对集中在主瓣，此时更接近于真实频谱值。因此，时域采用不同的窗函数进行截取会产生不同的影响，结果如下表1所示。

表1五种窗函数基本参数的比较；

窗函数	旁瓣峰值/dB	主瓣峰值/
			矩形窗	-13	2
三角窗	-25	4
			汉宁窗	-31	4
海明窗	-41	4
			布拉克曼窗	-57	6
高斯窗	-55

表1

加窗函数是一种提高频谱分析质量的重要措施，选用恰当的窗函数能够达到抑制假频和提高信噪比。加窗函数主要有两个目的：第一可提高离散频谱幅值的精度；第二可提高信噪比和抑制假频，会出现假频的现象是由于常用的窗函数存在旁瓣。由表1可知采用合适的窗函数，可减少旁瓣的影响，将能量集中于主瓣中，但减少旁瓣必然会造成主瓣加宽，会导致主瓣的频谱泄露严重和频谱干扰，因此在选择加窗函数的时候，必须综合考虑这两个方面，做出合理的选择。

S2、离散时间序列{x(n)}(n＝0,1,…,N-1)DFT分别用实部和虚部来表示，即

$X_R\left(k\right)=\underset{n=0}{\overset{N-1}{\Σ}}x\left(n\right)cos(2\mathrm{\π}kn/N),k=0,1,2,...N/2-1---\left(1\right)$

$X_I\left(k\right)=-\underset{n=0}{\overset{N-1}{\Σ}}x\left(n\right)sin(2\mathrm{\π}kn/N),k=0,1,2,...N/2-1---\left(2\right)$

FFT是当N＝2^M(M为正整数)的DFT可进行的递推一种快速算法。

以上公式变换中，频谱分辨率为Δf＝f_s/N，与采样率f_s成正比，与采样点数N成反比，因此当这两个参数一定时，频率分辨率无法再得到提高。

离散序列x(n)包含的频域信息为0～f_s/2，若用连续的傅里叶变换计算，把频谱曲线看做连续的，即公式(1)和(2)中的k看做区间[0,N/2]内的连续实数，公式可变为

$X_R\left(f\right)=\underset{n=0}{\overset{N-1}{\&Sigma;}}x\left(n\right)cos(2\mathrm{\&pi;}nf/f_s),0\&le;f<f_s/2---\left(3\right)$

$X_I\left(f\right)=-\underset{n=0}{\overset{N-1}{\&Sigma;}}x\left(n\right)sin(2\mathrm{\&pi;}nf/f_s),0\&le;f<f_s/2---\left(4\right)$

此时频率分辨率不受采样点数的影响，f是一个连续的频率变量。

然而计算机只能计算离散数据，故实际应用中得到的是离散频谱，因此可选择合适并且很窄的频率区间[f₁,f₂]，指定较小分辨率，从而大大提高分析精度。

此方法的实质是由连续傅里叶变换经过四个步骤到离散傅里叶变换DFT的过程。四个步骤包括将积分变成求和、时域离散化、频域离散化和时域截断。如果不进行频域离散化这一步骤，可将原DFT表达式进行转化

$2\mathrm{\&pi;k}/N\&DoubleLeftRightArrow;2\mathrm{\&pi;f}/f_s---\left(5\right)$

由上式可知k与N比值等价与连续频率f与采样率f_s的比值。由于采样序列x(n)有限，可以基于FFT+FT原理对序列x(n)进行傅里叶变换，得：

$X\left(f\right)=\underset{n=0}{\overset{N-1}{\&Sigma;}}x\left(n\right)e^{-j2\mathrm{\&pi;}n\frac{f}{f_s}}---\left(6\right)$

对指定的频率区间[f₁,f₂]，此区间在[0,f_s/2]内，对其进行n点的等间隔细化处理步骤如下：

(1)可得细化分辨率为

Δf＝(f₂-f₁)/L (7)

(2)从而可得计算细化的频率序列为

{f₁,f₁+Δf,f₁+2Δf,…,f₂} (8)

三相电能质量分析仪的测频范围为42.5～57.5Hz，频率精度为±0.01Hz，即在 有限范围内进行基波周期的测量，基波频率采取较小的步长0.015Hz进行等间隔取点计算1000个频率点值，即FFT+FT对此段测频范围进行细化，步长为0.001Hz，进行DFT运算。若测量每一个频率进行全部点计算，计算量大，很多点的结算并没有实际效果，工程实时性一般，因此需要采取选点计算进行比较，通过少量点的比较，得到测频值。但由于加不同的窗函数对频谱曲线影响不同，在选择搜索算法也会有所不同。首先研究电力信号加不同的窗函数进行时域截取后，进行FFT+FT运算后的精度影响；然后再研究信号加不同的窗函数后的频谱在测频范围内的特征，根据其不同的特性来决定搜索算法和窗函数的选定。下面进行具体说明：在本实施例中，通过加窗函数的FFT+FT测量值与真值误差分析

设采样频率f_s＝25600Hz，采样点数N＝5120，则频谱分辨率为f_s/N＝5Hz。设输入信号：x＝100*cos(2πf)，加高斯白噪声40dB，频率f＝42.5:0.015:47.5，总共统计1000个频率点的计算与真值误差，频谱[42.5,57.5]细化步长为0.001Hz，总共遍历15000点的计算。

加各种窗函数进行时域信号的截取，然后进行FFT+FT分析，仿真误差结果如图2所示；其中，加矩形窗、加三角窗、加汉明窗在一些频率处测量频率值与真值误差较大，不能满足仪器对频率的测量要求；加汉宁窗、加Blackman窗、加高斯窗所得到的频率测量值与真值误差较小，可满足测量要求，并且抗噪能力强，这三种窗函数的旁瓣的衰减较大，可减少频谱泄露以及正负频率的频谱相互影响。但由于对此频率段的细化倍数为1000倍，计算点数很多，运算速度以及搜索出最大值有待优化。因此需考虑加窗之后数据的趋势，研究取点计算的比较算法，找到最大值所对应的频率值，即测量频率。

从图2可知，汉宁窗、Blackman窗、高斯窗的测量精度较好。其中汉宁窗的主瓣包括4根谱线，Blackman窗主瓣包含6根谱线，高斯窗主瓣很宽，其中谱线的宽度为实际分辨率。考虑边界情况，即最坏的条件，当基波频率为42.5Hz、57.5Hz时，分别进行FFT+FT仿真计算，其仿真计算结果如图3、4所示；

其中，加载汉宁窗的频谱响应可能为多峰函数，加载Blackman窗与高斯窗的频谱相应都是单峰函数，由于三种窗函数都是对称窗，图3、图4为频率取边 界值时，即为最差的一种情况；由此可知，在测频范围内，加汉宁窗函数可会出现多峰情况，加Blackman窗与高斯窗由于边界处这两种窗函数包含的谱线为其主瓣的一半，其他频率值处包含的谱线多于边界处，故在其他频率处也为单峰函数。对于搜索算法，单峰要比多峰的算法具有实现简单、速度快、实时性好、精度高等特点，因此选择Blackman窗与高斯窗最佳；但由于高斯窗主瓣太宽，主瓣谱线太容易受到干扰，且高斯窗更适用于非周期函数的时域截取，Blackman窗主瓣宽度有限，故选Blackman窗更合适。

综上，Blackman窗具有进行FFT+FT计算精度高，旁瓣衰减快，具有对称性且在测频范围内可保证为单峰函数的性质，方便后续进行一维搜索择点运算，减少运算量。

S3、设置初始参数值

根据三相电能质量分析仪的基波频率测量区间42.5～57.5Hz，测量精度在0.01Hz以内，设置将0.618算法的边界值初值设置为a＝42.5Hz和b＝57.5Hz，精度设置为L＝0.001Hz；设置试探点λ＝a+0.382(b-a)和μ＝a+0.618(b-a)；

S4、通过0.618算法搜索X(f)的最大幅值来得到基波频率值

0.618算法是一维搜索算法中很成熟的方法，适用于单峰函数。单峰函数函数具有的重要性质：通过计算区间[a,b]内两个不同点处的函数值，就能确定一个包含极小点的子区间。

将试探点λ和μ分别代入X(f)，计算出X(λ)和X(μ)，再对X(λ)和X(μ)求幅值，得到f(λ)和f(μ)；

当b-a＞L时，若f(λ)＞f(μ)时，则令b＝μ，μ＝λ，λ＝a+0.382(b-a)，再代入到X(f)进行重复迭代；

当b-a＞L时，若f(λ)＜f(μ)时，则令a＝λ，λ＝μ，μ＝a+0.618(b-a)，再代入到X(f)进行重复迭代；

当b-a＜L时，则停止计迭代，得到基波频率值为(a+b)/2。

实例仿真

图5是对信号加载40dB的高斯白噪声的不同频率的理论与测量误差图。

图6是对信号随机噪声的不同频率的理论与测量误差图。

图7是对信号加载随机噪声、高含量谐波的不同频率的理论与测量误差图。

一、对输入信号加高斯白噪声

在本实施例中，基于加Blackman窗函数的电力信号经FFT+FT变换的频谱曲线可确定为单峰函数，利用0.618算法搜索，可以通过两个试探点计算少量点数进行比较，达到快速收敛频率区间，准确定位在精度允许范围内最大值的位置。

本实施例中通过此方法对输入信号测量值和真值的误差仿真，其频谱范围为[42.5,57.5]Hz,取频率点1000个进行计算，即步长为0.015Hz。设采样频率f_s＝25600Hz，采样点数为N＝5120，则频谱分辨率为5Hz。

输入信号：x＝100*cos(2πf)，加高斯白噪声40dB，0.618算法的精度L＝0.001。FFT+FT仿真结果如图5所示；

0.618算法的迭代次数为20次，收敛速度较快。由不同频率的理论与测量误差图可知，在40dB高斯白噪声影响下，用0.618算法进行测频算法，得到与频率的真值误差在0.0006Hz以内，符合频率的测量要求，且抗噪能力强。

二、对输入信号加随机噪声

对三相电能质量分析仪的噪声进行采集，经过数据处理得到噪声的平均值为0.1327和方差为10.1618，根据噪声的平均值和方差，在输入信号上加载随机噪声，进行仿真计算，得到下表2；

表2是实际噪声影响下的测频算法结果；

表2

如图6所示，通过实际检测的随机噪声的影响可知，0.618算法测量频率点150个，得到的测频误差基本上在0.006Hz以内。虽然0.618算法会比遍历的误差较大，但相差很小，可满足频率测量精度的要求，而且此方法实际计算点数即迭代次数大约有20次，足以体现一维搜索0.618算法收敛速度较快，且与每次遍历测量基波频率计算15000点计算量相比，很大程度的加快了运算速度，更加优化了测频算法。

三、对输入信号加谐波

由于在实际电网中一般会产生奇次谐波，常见的3次谐波，因此仿真信号加入较高含量的3次谐波进行检测谐波对算法的影响。仿真信号：x＝100*cos(2πf)+20*cos(2π3f)，噪声加载实际检测参数。

对于基波信号加20％的3次谐波以及噪声影响，仿真频率是42.5:0.15:57.5，测试100个点的频率误差，由图7可知，基本在0.008Hz范围以内，本算法对于含有谐波和噪声的信号具有很高的识别度，而且不需要进行滤波，只需要对多关心的频率范围进行频谱细化，即可得到精度较高的测量值。

综上可见，优化的基于FFT+FT的算法具有精度高、抗噪性强、不需滤波和不受谐波影响、可任意设定测频范围、实时性强等优点，可应用的范围广。

尽管上面对本发明说明性的具体实施方式进行了描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

Claims (3)

1.一种基于FFT+FT的基波测频方法，其特征在于，包括以下步骤：

(1)、对采样点为N的离散序列添加窗函数，并标记为序列x(n)；

(2)、基于FFT+FT原理对序列x(n)进行傅里叶变换，得：

$X\left(f\right)=\underset{n=0}{\overset{N-1}{\&Sigma;}}x\left(n\right)e^{-j2\mathrm{\&pi;}n\frac{f}{f_s}}$

其中，f为连续频率，f_s为采样率；

(3)、设置初始参数值

设置三相电能质量分析仪的基波频率测量区间为[f₁,f₂]Hz，测量精度在cHz以内；再根据该设置将0.618算法的边界边界值初值设置为a＝f₁和b＝f₂，精度设置为L，L∈c；设置试探点λ＝a+0.382(b-a)和μ＝a+0.618(b-a)；

(4)、通过0.618算法搜索X(f)的最大幅值来得到基波频率值

将试探点λ和μ分别代入X(f)，计算出X(λ)和X(μ)，再对X(λ)和X(μ)求幅值，得到f(λ)和f(μ)；

当b-a＞L时，若f(λ)＞f(μ)时，则令b＝μ，μ＝λ，λ＝a+0.382(b-a)，再代入到X(f)进行重复迭代；

当b-a＞L时，若f(λ)＜f(μ)时，则令a＝λ，λ＝μ，μ＝a+0.618(b-a)，再代入到X(f)进行重复迭代；

当b-a＜L时，则停止计迭代，得到基波频率值为(a+b)/2。

2.根据权利要求1所述的基于FFT+FT的基波测频方法，其特征在于，所述的窗函数选择Blackman窗。

3.根据权利要求1所述的基于FFT+FT的基波测频方法，其特征在于，所述的精度L＝0.001Hz。