CN105372492B

CN105372492B - 基于三条dft复数谱线的信号频率测量方法

Info

Publication number: CN105372492B
Application number: CN201410435495.0A
Authority: CN
Inventors: 庞浩
Original assignee: Changzhou Hao Yun Industrial Control Technology Co Ltd
Current assignee: Changzhou Haoyun Industrial Control Technology Co ltd
Priority date: 2014-08-31
Filing date: 2014-08-31
Publication date: 2018-06-22
Anticipated expiration: 2034-08-31
Also published as: CN105372492A

Abstract

本发明涉及一种基于三条DFT复数谱线的信号频率测量方法，属于信号参数测量技术领域。本发明的特征在于其处理步骤包含：采样信号经过加窗处理后进行DFT变换，在对应频率范围内通过比较谱线幅值或者实部虚部绝对值之和搜索最高谱线及其相邻的两条谱线，基于三条谱线对应复数的实部虚部的运算求取中间参数，再通过求解方程、或反函数公式、或逼近多项式获得频偏参数，最终测量出信号频率。本发明在搜索最高谱线、以及计算中间参数时减少了乘法和开方计算，简化了信号频率测量算法的实现；同时采用相邻谱线取差的形式，抑制了其他信号旁瓣的影响，提高了测量准确度。

Description

基于三条DFT复数谱线的信号频率测量方法

技术领域

本发明涉及一种基于三条DFT复数谱线的信号频率测量方法，属于信号参数测量技术领域。

背景技术

当前，基于离散傅里叶变换DFT或其快速算法FFT分析频率信号的方法已经广泛使用。但是，DFT具有栏栅效应，即实际信号频率未必落在离散谱线上，由此需要采用插值算法估计实际信号的频率。2012年《中国电机工程学报》32卷16期上发表的“基于三谱线插值FFT的电力谐波分析算法”文章中提出了对输入离散信号加窗傅里叶变换后，通过采用幅值最高谱线及其两边相邻的谱线，三条谱线共同插值测量信号频率的方法。如果三条谱线的离散频率序号分别对应k ₁-1、k ₁和k ₁+1，三条谱线的幅值|Y ₁ |=|Y(k ₁-1)|、|Y ₂|=|Y(k ₁)|和|Y ₃|=|Y(k ₁+1)|满足：

，

其中，γ= k ₀- k ₁，k ₀对应实际信号频率。当N较大时，上式可以简化表示为λ=g(γ)，其反函数记为。该文进一步提出采用多项式逼近方法计算：

。

已有方法的不足在于基于谱线幅值进行最高谱线搜索。而且，获得谱线幅值需要计算实部和虚部的平方和、然后进行开方，该计算量较大。已有方法对中间参数γ的计算也基于谱线幅值实现，同时虽然采用了多条谱线，并未有效抑制其他频点信号的旁瓣干扰，测量准确度低。

发明内容

本发明的目的是在搜索最高谱线、以及计算中间参数λ时减少乘法和开方计算，减少运算量；同时，设计更合理的中间参数γ的计算公式，实现旁瓣干扰的抑制，提高测量准确度。

本发明为解决上述技术问题而提供一种基于三条DFT复数谱线的信号频率测量方法，该测量方法的步骤如下：

步骤（1）：将采样率为F _S、采样点为连续截取的采样信号x(n)，进行加窗处理得到加窗信号y(n)，加窗处理公式为：

，

其中为N点的窗函数序列，n = 0:(N-1)；

步骤（2）：对加窗信号y(n)进行离散傅里叶DFT变换，得到离散频谱Y(k)，其中离散频率序号k = 0:(N-1)；

步骤（3）：在设定频率范围所对应的离散频率序号范围[k _s, k _e]内，搜索最大的谱线作为最高谱线，其离散频率序号为k ₁，同时选择该谱线两侧、离散频率序号为k ₁-1和k ₁+1的相邻两条谱线，其中Re(Z)表示取复数Z的实部值，Im(Z)表示取复数Z的虚部值；

步骤（4）：依据k ₁-1、k ₁和k ₁+1对应的三条谱线Y ₁ =Y(k ₁-1)、Y ₂ =Y(k ₁)和Y ₃ =Y(k ₁+1)计算中间参数λ：

步骤（5）：求解如下方程中的频偏参数γ：

，

其中，W(ω)是窗函数w(n)的离散时间傅里叶变换DTFT的结果，归一角频率ω=2π f/F _S=2πk/N；

步骤（6）：依据频偏参数γ计算被测信号频率f _γ，计算公式为：

。

进一步地，所述的步骤（3）搜索最高和次高谱线的处理为：在设定频率范围所对应的离散频率序号范围[k _s, k _e]内，搜索|Y(k)|最大的谱线作为最高谱线，其离散频率序号为k ₁，同时选择该谱线两侧、离散频率序号为k ₁-1和k ₁+1的相邻两条谱线。

进一步地，所述的步骤（5）采用中间参数实测值λ与频偏参数γ关系函数的反函数公式计算频偏参数γ，该反函数为：

。

进一步地，所述的步骤（5）采用中间参数实测值λ与频偏参数γ关系函数的反函数的逼近多项式公式计算频偏参数γ，该逼近多项式公式为：

，

其中，M是逼近多项式的最高次数，a _m（m=0:M）是多项式第m次项λ ^m的系数。

本发明信号频率测量方法的设计原理是：假设一个频率为f ₀、幅值为、初相位为的单一频率信号x(t)，在经过了采样率为Fs的模数变换后得到如下形式的离散信号：

。

如果所加窗函数的时域形式为w(n)，其离散时间傅里叶变换DTFT得到的连续频谱为W(ω)，则忽略负频点-f ₀处频峰的旁瓣影响，在正频点f ₀附近的连续频谱函数可以表达为：

，

上式进行离散抽样，即可得到离散傅立叶变换DFT的表达式为：

，

其中，离散频率间隔为Δf=F _S/N。

余弦窗函数是DFT最为常用的一类窗函数。对应余弦窗函数的统一时域形式为：

，

余弦窗w(n)的离散时间傅里叶变换DTFT结果为：

，

其中：

，

在信号DTFT频谱曲线的主瓣内，且当N较大时，近似有：

，

当时，上式取等号。依据常用窗函数系数，在主瓣-H<k<H内，其相邻两条谱线和的相位相差近似为π；而对应H<k<N/2的旁瓣内和接近同相位。由此，当谱线Y(k ₁)和Y(k ₂)的相位相差0或π时，有：

，

由此，不必计算谱线幅值，直接通过谱线实部和虚部的绝对值之和，即可实现最高谱线的搜索、以及中间参数λ的计算。

进一步地，通过对多数窗函数进一步进行频域处理所得到的新的窗函数，将进一步抑制其旁瓣，由此可以减小被测信号DFT之后负频点、谐波频点对待测信号强度的旁瓣干扰。由此，本发明以(Y(k ₁)-Y(k ₁-1))=Y ₂-Y ₁和(Y(k ₁)-Y(k ₁+1))= Y ₂-Y ₃的比值为基础，识别实际信号频点位置。同时，为了减少开方运算，采用谱线的实部和虚部分别进行计算，即：

。

本发明方法基于上述原理设计，减少了乘法和开方运算，简化了信号频率测量算法实现。计算公式中采用相邻谱线取差的形式，进一步抑制其他信号旁瓣对测量的影响，提高了测量准确度。

附图说明

图1是本发明的基于三条DFT复数谱线的信号频率测量方法的计算流程图。

具体实施方式

下面结合附图1所示的计算流程图对本发明的两个具体实施方式作进一步的说明。这两个实施例应用于对50Hz附近频率信号进行测量。第一个具体实施方式采用哈宁（Hanning）窗，其具体步骤如下：

步骤（1）：将采样率Fs=1500Hz、连续截取N=512点的信号x(n)，进行加窗处理得到加窗信号y(n)，加窗处理公式为：

，

其中w(n)选择N=512点的Hanning窗函数序列，即：

，n = 0:(N-1)；

步骤（3）：在离散频率序号范围[15, 23]，即对应频率[43.945, 67.383]Hz的范围内，搜索最大的谱线作为最高谱线，其离散频率序号为k ₁，同时选择该谱线两侧、离散频率序号为k ₁-1和k ₁+1的相邻两条谱线；

；

步骤（5）：求解如下方程中的频偏参数γ：

；

由此，计算频偏参数γ的函数公式为：

，

。

第二个具体实施方式采用布莱克曼（BlackMan）窗，其具体步骤如下：

步骤（1）：将采样率Fs=1500Hz、连续截取 N=512点的信号x(n)，进行加窗处理得到加窗信号y(n)，加窗处理公式为：

，

其中w(n)选择N=512点的布莱克曼（BlackMan）窗函数序列，即：

，n = 0:(N-1)；

；

步骤（5）：采用中间参数实测值λ与频偏参数γ关系函数的反函数的逼近多项式公式计算频偏参数γ，采用最高为7次的逼近多项式公式，形式如下：

，

。

依据第一个和第二个实施方式，分别输入相同的一组仿真测试数据，以验证两个实施例的计算结果。该输入信号x(n)是基波频率f ₁为50.1 Hz、包含2至9次谐波的信号，具体形式为：

其中，基波和各次谐波的幅值分别是：1, 0.02, 0.1, 0.01, 0.05, 0.0, 0.02,0.0, 0.01；初始相位分别是-23.1°, 115.6°, 59.3°, 52.4°, 123.8°, 161.8°, -31.8°,119.9°, -63.7°。

第一个实施方式中，离散频率序号范围15至23的9条谱线计算结果依次是：0.152783905+j1.76229599, 4.70210906+j54.2254920, -10.9858392-j126.688437,6.36734959+j73.4295187, -0.221526634-j2.55345712, -0.0512202109-j0.589197696,-0.0199885230-j0.228698046, -0.00996303987-j0.112669252。不论幅值比较还是采用实部虚部绝对值之和，均可得到最高谱线对应k ₁=17，同时选择该谱线两侧、离散频率序号为k ₁-1=16和k ₁+1=18的相邻两条谱线。于是，中间参数计算结果为λ=0.05039996845，频偏参数γ=2∙λ= 0.1007999369，最终测量得到的信号频率为50.099999815Hz，测量误差-0.000000369%。

第二个实施方式中，离散频率序号范围15至23的9条谱线计算结果依次是：-0.63151373-j7.27542732, 4.87561219+j56.23242377, -9.38422261-j108.21320679,6.10560557+j70.41558734, -1.11464837-j12.84931264, 0.00734999+j0.0889938, -0.00022265+j0.00101658, 0.00048733+j0.00853355。不论幅值比较还是采用实部虚部绝对值之和，均可得到最高谱线对应k ₁=17，同时选择该谱线两侧、离散频率序号为k ₁-1=16和k ₁+1=18的相邻两条谱线。于是，中间参数计算结果为λ=0.0413416339。带入逼近多项式后，得到频偏参数γ=0.1007992629，最终测量得到的信号频率为50.099997840Hz，测量误差-0.000004311%。

Claims

1.一种基于三条DFT复数谱线的信号频率测量方法，其特征在于该信号频率测量方法包含如下步骤：

步骤(1)：将采样率为F_S、采样点为连续截取的采样信号x(n)，进行加窗处理得到加窗信号y(n)，加窗处理公式为：

y(n)＝x(n)·w(n)，

其中w(n)为N点的窗函数序列，n＝0:(N-1)；

步骤(2)：对加窗信号y(n)进行离散傅里叶DFT变换，得到离散频谱Y(k)，其中离散频率序号k＝0:(N-1)；

步骤(3)：在设定频率范围所对应的离散频率序号范围[k_s,k_e]内，搜索(|Re(Y(k))|+|Im(Y(k))|)最大的谱线作为最高谱线，其离散频率序号为k₁，同时选择该谱线两侧、离散频率序号为k₁-1和k₁+1的相邻两条谱线，其中Re(Z)表示取复数Z的实部值，Im(Z)表示取复数Z的虚部值；

步骤(4)：依据k₁-1、k₁和k₁+1对应的三条谱线Y₁＝Y(k₁-1)、Y₂＝Y(k₁)和Y₃＝Y(k₁+1)计算中间参数λ：

步骤(5)：求解如下方程中的频偏参数γ：

其中，W(ω)是窗函数w(n)的离散时间傅里叶变换DTFT的结果，归一角频率ω＝2πf/F_S＝2πk/N，f：信号频率；

步骤(6)：依据频偏参数γ计算被测信号频率f_γ，计算公式为：

f_γ＝(k₁+γ)·Fs/N。

2.根据权利要求1所述的基于三条DFT复数谱线的信号频率测量方法，其特征在于：所述的步骤(3)搜索最高和次高谱线的处理为：在设定频率范围所对应的离散频率序号范围[k_s,k_e]内，搜索|Y(k)|最大的谱线作为最高谱线，其离散频率序号为k₁，同时选择该谱线两侧、离散频率序号为k₁-1和k₁+1的相邻两条谱线。

3.根据权利要求1所述的基于三条DFT复数谱线的信号频率测量方法，其特征在于：所述的步骤(5)采用中间参数实测值λ与频偏参数γ关系函数λ＝g(γ)的反函数公式计算频偏参数γ，该反函数为：

γ＝g^-1(λ)。

4.根据权利要求1所述的基于三条DFT复数谱线的信号频率测量方法，其特征在于：所述的步骤(5)采用中间参数实测值λ与频偏参数γ关系函数λ＝g(γ)的反函数g^-1(λ)的逼近多项式公式计算频偏参数γ，该逼近多项式公式为：

其中，M是逼近多项式的最高次数，a_m，m＝0:M，是多项式第m次项λ^m的系数。