CN104914308A - 一种基于两条dft复数谱线的信号相位测量方法 - Google Patents

Abstract

本发明涉及一种基于两条DFT复数谱线的信号相位测量方法，属于信号参数测量技术领域。本发明的特征在于其处理步骤包含：将采样信号经过加窗处理后进行DFT变换，查找对应待测信号频率附近的两条复数谱线，基于两条谱线的复数值通过直接推导公式、或逼近多项式公式计算出中间参数，最终的幅值测量结果等于中间参数的模。本发明直接基于谱线复数进行计算，无需对每条谱线取模，减少了计算量，而且计算过程能够抵消其他频率信号的旁瓣干扰，提高了测量精度。

Description

一种基于两条DFT复数谱线的信号相位测量方法

技术领域

本发明涉及一种基于两条DFT复数谱线的信号幅值和相位测量方法，属于信号参数测量技术领域。

背景技术

当前，基于离散傅里叶变换DFT或其快速算法FFT分析频率信号的方法已经广泛使用。但是，DFT具有栏栅效应，即实际信号频率未必落在离散谱线上，由此需要采用插值算法估计实际信号的频率、幅值和相位。2003年《中国电机工程学报》23卷6期上发表的“应用FFT进行电力系统谐波分析的改进算法”文章中提出了对输入离散信号加窗傅里叶变换后，通过选择幅值最高和次高两条谱线，插值测量信号频率、幅值和相位的方法。如果两条谱线的离散频率序号分别对应k₁和k₂＝k₁+1，则实际信号频率对应的位置k₀满足k₁≤k₀≤k₂。引入一个辅助参数α＝k₀-k₁-0.5，忽略其他信号干扰，则α的数值范围是[-0.5,0.5]。由此，基于两条谱线幅值|Y(k₁)|和|Y(k₂)|计算信号幅度A可以按照下面插值公式计算：

$A=\frac{2\left(\right|Y\left(k_1\right)|+|Y\left(k_2\right)\left|\right)}{\left|W(2\mathrm{\π}\·(-\mathrm{\α}-0.5)/N)\right|+\left|W(2\mathrm{\π}\·(-\mathrm{\α}+0.5)/N)\right|}.$

对于一般的实系数窗函数，当N较大时，上式可以进一步简化为A＝(1/N)·(|Y(k₁)|+|Y(k₂)|)·v(α)的形式，v(α)是频偏参数α的函数、且与N无关。如果采用最高M次的逼近多项式计算函数，则信号幅度A的计算公式可以进一步表示为：

$A=1/N\·\left(\right|Y\left(k_1\right)|+|Y\left(k_2\right)\left|\right)\·(\underset{m=0}{\overset{M}{\mathrm{\Σ}}}{b}_m\·{\mathrm{\α}}^m)$

已有方法给出的相位计算公式为：

θ＝arg(Y(k_i))+π/2-arg(W(2π·(k_i-k₀)/N))

其中，i取1或2。

已有方法的不足在于信号幅值和相位的计算是相互独立的，其幅值计算需要计算实部和虚部的平方和、然后进行开方，其相位计算需要计算Y(k_i)和W(2π·(k_i-k₀)/N)两个复数的角度，所以计算量大。同时已有方法还容易受到其他频率信号的旁瓣干扰。

发明内容

本发明的目的是提供一种基于两条DFT复数谱线的信号相位测量方法，用以解决现有方法运动量大和旁瓣干扰的问题。

为实现上述目的，本发明的方案包括：

一种基于两条DFT复数谱线的信号相位测量方法，步骤如下：

步骤(1)：将采样率为F_S、采样点为连续截取的N点的采样信号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)：依据所需测量幅值和相位的信号的频率f₀所对应的离散频率序号值k₀，查找到临近k₀的两条谱线，其离散频率序号分别为k₁和k₂，其中k₀＝N·f₀/F_S，k₁等于不大于k₀的最大整数，即k₁＝floor(k₀)，k₂＝k₁+1；

步骤(4)：依据k₁和k₂对应的两条复数谱线Y(k₁)和Y(k₂)计算中间参数Y：

$Y=\frac{2\·(Y\left(k_2\right)-Y\left(k_1\right))}{W(2\mathrm{\π}\·(k_2-k_0)/N)-W(2\mathrm{\π}\·(k_1-k_0)/N)};$

步骤(5)：相位测量结果θ等于Y的幅角加上π/2，即

所述的步骤(4)采用逼近多项式计算中间参数Y，其计算公式为：

$Y\≈\frac{1}{N}\·(Y\left(k_2\right)-Y\left(k_1\right))\·(\underset{p=0}{\overset{P}{\mathrm{\Σ}}}{b}_p\·{\mathrm{\α}}^p+j\·\underset{q=0}{\overset{Q}{\mathrm{\Σ}}}{c}_q\·{\mathrm{\α}}^q),$

其中，α＝k₀-k₁-0.5，P和Q分别是实部和虚部逼近多项式的最高次数，b_p(p＝0:P)和c_q(q＝0:Q)分别是实部逼近多项式第p次项α^p和虚部逼近多项式第q次项α^q的系数。

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

$x\left(n\right)=A\·\sin (2\mathrm{\π}\frac{f_0}{F_S}n+\mathrm{\θ})$

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

$\stackrel{\‾}{X}\left(f\right)=\frac{A}{2j}{e}^{\mathrm{j\θ}}W\left(\frac{2\mathrm{\π}(f-f_0)}{f_s}\right)$

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

$Y\left(k\right)=\stackrel{\‾}{X}(k\·\mathrm{\Δf})=\frac{A}{2j}{e}^{\mathrm{j\θ}}W\left(\frac{2\mathrm{\π}(k\·\mathrm{\Δf}-f_0)}{F_S}\right)$

其中，离散频率间隔为Δf＝F_S/N。于是，

$Y\left(k_2\right)-Y\left(k_1\right)=\frac{A}{2}{e}^{j(\mathrm{\θ}-\frac{\mathrm{\π}}{2})}(W\left(\frac{2\mathrm{\π}(k_2-k_0)}{N}\right)-W\left(\frac{2\mathrm{\π}(k_1-k_0)}{N}\right))$

其中，离散频率间隔为Δf＝F_S/N。由此，

$Y=A{e}^{j(\mathrm{\θ}-\frac{\mathrm{\π}}{2})}=\frac{2\·(Y\left(k_2\right)-Y\left(k_1\right))}{W(2\mathrm{\π}\·(k_2-k_0)/N)-W(2\mathrm{\π}\·(k_1-k_0)/N)}.$

所以，直接采用复数谱线进行计算所得的幅值测量结果A等于中间参数Y的模，相位测量结果θ等于Y的幅角。

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

$\begin{array}{cc}w\left(n\right)=\underset{h=0}{\overset{H}{\mathrm{\Σ}}}{g}_h\·\cos (\frac{2\mathrm{\πh}}{N}\·n)& n=\mathrm{0,1},...,N-1\end{array}$

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

$W\left(\mathrm{\ω}\right)=N\·e^{-j\frac{\mathrm{N\ω}}{2}}\·\sin \frac{\mathrm{N\ω}}{2}\·F\left(\mathrm{\ω}\right)$

其中：

$F\left(\mathrm{\ω}\right)=\underset{h=0}{\overset{H}{\mathrm{\Σ}}}\frac{g_h}{2N}\·\frac{\sin \mathrm{\ω}}{\sin (\frac{\mathrm{\ω}}{2}-\frac{\mathrm{\πh}}{N})\sin (\frac{\mathrm{\ω}}{2}+\frac{\mathrm{\πh}}{N})}+j\·\underset{h=0}{\overset{H}{\mathrm{\Σ}}}{g}_h$

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

$F\left(\frac{2\mathrm{\πk}}{N}\right)\≈\frac{1}{\mathrm{\π}}\underset{h=0}{\overset{H}{\mathrm{\Σ}}}\frac{g_h\·k}{(k^2-h^2)}$

当时，上式取等号。依据常用余弦窗函数系数，在主瓣-H<k<H内，其相邻两条谱线W(ω)和的相位相差近似为π；而对应H<k<N/2的旁瓣内W(ω)和接近同相位。由此，对多数余弦窗函数频域的处理所得到的新的窗函数，能够进一步抑制旁瓣，因此可以减小其他频率信号及其DFT的负频率信号对待测频率信号谱线的影响，从而提高测量精度。

附图说明

图1是本发明相位测量的过程图；

图2是本发明实施例同时测量幅值和相位的过程图。

具体实施方式

下面结合附图对本发明做进一步详细的说明。

以下两个实施例均用于对50Hz附近的频率信号进行测量，并且在测量相位信号的同时，也给出了幅值。

实施例1

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

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

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

$w\left(n\right)=0.5-0.5\&CenterDot;\cos \left(\frac{2\mathrm{\&pi;n}}{N}\right),n=0:(N-1);$

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

步骤(3)：依据所需测量幅值和相位的信号的频率f₀所对应的离散频率序号值k₀，查找到临近k₀的两条谱线，其离散频率序号分别为k₁和k₂，其中k₀＝N·f₀/F_S，k₁等于不大于k₀的最大整数，即k₁＝floor(k₀)，k₂＝k₁+1；

步骤(4)：依据k₁和k₂对应的两条复数谱线Y(k₁)和Y(k₂)计算中间参数Y：

$\begin{array}{c}Y=\frac{2\&CenterDot;(Y\left(k_2\right)-Y\left(k_1\right))}{W(2\mathrm{\&pi;}\&CenterDot;(k_2-k_0)/N)-W(2\mathrm{\&pi;}\&CenterDot;(k_1-k_0)/N)}\\ =\frac{4\mathrm{\&pi;}}{3N}\&CenterDot;e^{-\mathrm{j\&pi;}(\mathrm{\&alpha;}-0.5)}\&CenterDot;\frac{({\mathrm{\&alpha;}}^2-0.25)({\mathrm{\&alpha;}}^2-2.25)}{\cos \left(\mathrm{\&pi;\&alpha;}\right)}\&CenterDot;(Y\left(k_2\right)-Y\left(k_1\right))\end{array},$

其中，α＝k₀-k₁-0.5；

步骤(5)：对应频率f₀的被测信号的幅值测量结果A等于中间参数Y的模，相位测量结果θ等于Y的幅角，即：

$A=\left|Y\right|,\mathrm{\&theta;}=\arg \left(Y\right)+\frac{\mathrm{\&pi;}}{2}.$

实施例2

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

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

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

$w\left(n\right)=0.42659-0.49656\&CenterDot;\cos \left(\frac{2\mathrm{\&pi;n}}{N}\right)+0.076849\&CenterDot;\cos \left(\frac{4\mathrm{\&pi;n}}{N}\right),n=0:(N-1);$

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

步骤(3)：依据所需测量幅值和相位的信号的频率f₀所对应的离散频率序号值k₀，查找到临近k₀的两条谱线，其离散频率序号分别为k₁和k₂，其中k₀＝N·f₀/F_S，k₁等于不大于k₀的最大整数，即k₁＝floor(k₀)，k₂＝k₁+1；

步骤(4)：采用逼近多项式计算中间参数Y，实部和虚部逼近多项式的最高次数分别为7次和6次，实际采用的计算公式为：

$\begin{array}{c}Y\&ap;\frac{1}{N}\&CenterDot;(Y\left(k_2\right)-Y\left(k_1\right))\\ \&times;((8.407008293\&CenterDot;\mathrm{\&alpha;}-10.414137183\&CenterDot;{\mathrm{\&alpha;}}^3+1.967257952\&CenterDot;{\mathrm{\&alpha;}}^5+0.040949919\&CenterDot;{\mathrm{\&alpha;}}^7)\\ +j\&CenterDot;(2.676033776-12.118716462\&CenterDot;{\mathrm{\&alpha;}}^2+5.739097723\&CenterDot;{\mathrm{\&alpha;}}^4-0.315597560\&CenterDot;{\mathrm{\&alpha;}}^6))\end{array},$

其中，α＝k₀-k₁-0.5；

步骤(5)：对应频率f₀的被测信号的幅值测量结果A等于中间参数Y的模，相位测量结果θ等于Y的幅角加上π/2，即：

$A=\left|Y\right|,\mathrm{\&theta;}=\arg \left(Y\right)+\frac{\mathrm{\&pi;}}{2}.$

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

$x\left(n\right)=\underset{i=1}{\overset{9}{\mathrm{\&Sigma;}}}{A}_i\sin (2\mathrm{\&pi;}\frac{i\&CenterDot;f_1}{F_S}n+{\mathrm{\&theta;}}_i)$

其中，基波和各次谐波的幅值分别是：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°。仿真测试中需要测量50.1Hz基波信号的幅值和相位。基波频率f₀所对应的离散频率序号值k₀＝17.1008，选择临近k₀的两条谱线的离散频率序号k₁＝17和k₂＝18。

采用哈宁窗的第一个实施方式中，离散频率序号范围17和18的两条谱线的复数值为：Y(k₁)＝-10.9858392-j126.688437，Y(k₂)＝6.36734959+j73.4295187。由此，

$\begin{array}{c}Y=\frac{4\mathrm{\&pi;}}{3N}{\&CenterDot;e}^{-\mathrm{j\&pi;}(\mathrm{\&alpha;}-0.5)}\&CenterDot;\frac{({\mathrm{\&alpha;}}^2-0.25)({\mathrm{\&alpha;}}^2-2.25)}{\cos \left(\mathrm{\&pi;\&alpha;}\right)}\&CenterDot;(Y\left(k_2\right)-Y\left(k_1\right))\\ =-0.392337042712-j0.919821553707\end{array}$

最终，幅值的测量结果为A＝|Y|＝1.0000000229，相对误差0.00000229％；相位的测量结果为-0.403170967rad，即-23.09999483°，绝对误差0.00000517°。

采用布莱克曼窗的第二个实施方式中，离散频率序号范围17和18的两条谱线的复数值为：Y(k₁)＝-9.38422261-j108.21320679，Y(k₂)＝6.10560557+j70.41558734。由此，

$\begin{array}{c}Y\&ap;\frac{1}{N}\&CenterDot;(Y\left(k_2\right)-Y\left(k_1\right))\\ \&times;((8.407008293\&CenterDot;\mathrm{\&alpha;}-10.414137183\&CenterDot;{\mathrm{\&alpha;}}^3+1.967257952\&CenterDot;{\mathrm{\&alpha;}}^5+0.040949919\&CenterDot;{\mathrm{\&alpha;}}^7)\\ +j\&CenterDot;\left(2.676033776-12.118716462\&CenterDot;{\mathrm{\&alpha;}}^2+5.739097723\&CenterDot;{\mathrm{\&alpha;}}^4-0.315597560\&CenterDot;{\mathrm{\&alpha;}}^6)\right)\\ =-0.392343927737-j0.919820343965\end{array}$

最终，幅值的测量结果为A＝|Y|＝1.0000016114，相对误差0.00016114％；相位的测量结果为-0.4031777747rad，即-23.10038489°，绝对误差-0.00038489°。

以上给出了具体的实施方式，但本发明不局限于所描述的实施方式。本发明的基本思路在于上述基本方案，对本领域普通技术人员而言，根据本发明的教导，设计出各种变形的模型、公式、参数并不需要花费创造性劳动。在不脱离本发明的原理和精神的情况下对实施方式进行的变化、修改、替换和变型仍落入本发明的保护范围内。

Claims (2)

1.一种基于两条DFT复数谱线的信号相位测量方法，其特征在于，步骤如下：

步骤(1)：将采样率为F_S、采样点为连续截取的N点的采样信号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)：依据所需测量幅值和相位的信号的频率f₀所对应的离散频率序号值k₀，查找到临近k₀的两条谱线，其离散频率序号分别为k₁和k₂，其中k₀＝N·f₀/F_S，k₁等于不大于k₀的最大整数，即k₁＝floor(k₀)，k₂＝k₁+1；

步骤(4)：依据k₁和k₂对应的两条复数谱线Y(k₁)和Y(k₂)计算中间参数Y：

$Y=\frac{2\&CenterDot;(Y\left(k_2\right)-Y\left(k_1\right))}{W(2\mathrm{\&pi;}\&CenterDot;(k_2-k_0)/N)-W(2\mathrm{\&pi;}\&CenterDot;(k_1-k_0)/N)};$

步骤(5)：相位测量结果θ等于Y的幅角加上π/2，即

2.根据权利要求1所述的一种基于两条DFT复数谱线的信号相位测量方法，其特征在于：所述的步骤(4)采用逼近多项式计算中间参数Y，其计算公式为：

$Y\&ap;\frac{1}{N}\&CenterDot;(Y\left(k_2\right)-Y\left(k_1\right))\&CenterDot;(\underset{p=0}{\overset{P}{\mathrm{\&Sigma;}}}{b}_p\&CenterDot;{\mathrm{\&alpha;}}^p+j\&CenterDot;\underset{q=0}{\overset{Q}{\mathrm{\&Sigma;}}}{c}_q\&CenterDot;{\mathrm{\&alpha;}}^q),$

其中，α＝k₀-k₁-0.5，P和Q分别是实部和虚部逼近多项式的最高次数，b_p(p＝0:P)和c_q(q＝0:Q)分别是实部逼近多项式第p次项α^p和虚部逼近多项式第q次项α^q的系数。