CN108710029B - 一种信号谐波分量初相位的精确估计方法 - Google Patents

Abstract

一种基于计算机程序对信号谐波分量的初相位进行精确估计的方法，其特征在于，包括四个步骤，(1)读取待估计信号的采样序列x(n)n＝0，1，...，N‑1，设采样率为f_sHz，其感兴趣信号谐波分量的频率f₀已知或已测得；(2)计算序列x(n)幅频谱|X(f)|和相频谱φ(f)在全频域的γN个等间隔的离散值|X(k)|和φ(k)k＝0，1，...，γN‑1(其中γ＞1)；(3)计算幅频谱序列|X(k)|k＝0，1，...，γN‑1中幅度值大于一阈值th的序列频率序号集合，记为K_th；(4)求解一优化问题确定感兴趣以f₀为频率的谐波分量初相位的精确估值。简言之，搜索出一最优初相点

，使得以

为初相的以f₀为频率的构造正弦序列和待测序列的相频谱在K_th频率点集合中的差值的范数最小，则将之估计为信号感兴趣的以f₀为频率的谐波分量的初相

。

Description

一种信号谐波分量初相位的精确估计方法

技术领域

本申请涉及一种基于计算机程序对信号谐波分量的初相位进行精确估计的方法。实际中，很多场合都涉及基于计算机程序对信号谐波分量的初相位进行精确估计的问题，如在电力系统、激光测距、卫星导航等领域。

在电力系统中，谐波初相位的测量具有很多重要的作用。例如，通过测得电流与电压之间的初相位差

从而利用

求出功率因数，进而求出有用功率，实现对电能的计费功能。当电网并网合闸时，需要两电网的电力信号之间的初相位相同，这时需要精确测量两网工频信号的初相位和初相位差。

在激光测距和卫星导航等测距定位系统中，接收机同时接受两个基本点发射机传来的信号，两个信号的初相位差与它们的发射机与接收机的距离差成正比，故求得初相位差即可确定被测物体如飞机和船舰所处的位置。为了提高精度，需要利用高精度的谐波初相位估计算法实现高精度的距离测量。

两个同频率的信号谐波分量，它们的相位差和初相位差相同，故以下文中不做区分。以下文中初相位也简称初相，在不产生混淆的情况下，初相有时也简称相位。

背景技术

最常见的信号谐波分量初相位估计方法自然是在频域进行的。其一般原理是：设有一个待求谐波分量初相位的实时间信号x(t)，它的采样序列为x(n)n＝0，1，...，N-1，采样率为 f_sHz，则其傅里叶变换为

其中，|X(f)|称为幅频谱，

称为相频谱；如果信号x(t)中有一个以f₀为频率以

为初相的谐波分量，则称

为该 f₀谐波分量真正初相位

的估值。

以上估计信号谐波分量初相位的方法是利用离散时间傅里叶变换(Discrete-time Fourier Transform，DTFT)直接估计频率为f₀的谐波分量的初相位

然而，事实上，按照数字信号处理的理论，由于不可避免的加窗，频谱泄漏和频谱混叠是客观存在的，这会造成相频谱φ(f)在频率点f₀处的值

与真实初相位

有一点偏差；另一方面，信号x(t)或者x(n)中的随机噪声和其他干扰分量，也会造成相频谱φ(f)在频率点f₀处的值

与真实初相位

出现一点偏差。因此，以上方法在理论上就是一种不完美的情况。

衡量信号谐波分量初相位估计算法精确度的指标，可采用估计的初相位与实际初相位的均方误差(mean-square error，MSE)来度量。在一定信噪比水平下，一个信号谐波分量初相位估计方法能达到的均方误差有一个理论上的极限，称为克拉美罗界(Cramer-RaoLower Bound， CRB)。直接利用DTFT计算估计谐波分量初相位的方法，其均方误差MSE远没有达到CRB，因此，信号谐波分量初相位的准确估计还有很大的可改进的空间。

参考文献：

[1]张佑峰.谐波对相位测量的影响[J].计量技术，1998(2)：26-29.

[2]马胜前.受干扰波相位测量的一种新方法[J].仪器仪表学报，2002，23(2)：218-220.

[3]刘远社，张振亚.同频率正弦信号动态相位差的测量[J].仪器仪表学报，2005，26(z1)：92-93.

[4]Lapuh R.Accurate phase measurement with two sampling voltmeters[C].Instrumentation and Measurement Technology Conference，2001.Imtc2001.Proceedings of the，IEEE.IEEE，2001：645-647 vol.1.

[5]Kramer G，Klische W.Multi-channel synchronous digital phaserecorder[C].Frequency Control Symposium and PDA Exhibition，2001.Proceedingsof the 2001 IEEE International.IEEE，2001：144-151.

[6]Ciglaric S，Fefer D，Jeglic A.Evaluation of an alternativelydesigned digital phase angle standard[J].IEEE Transactions onInstrumentation&Measurement，2002，51(4)：845-848.

[7]Marcin M R.Digital receiver phase meter for LISA[J].IEEETransactions on Instrumentation&Measurement， 2005，54(6)：2446-2453.

发明目的

提出一种信号谐波分量初相位的精确估计方法，使初相位估计误差能够小于利用DTFT直接估计初相位方法的误差。

技术方案

提出一种基于计算机程序对信号谐波分量的初相位进行精确估计的方法，其特征在于，包括四个步骤，(1)读取待估计信号的采样序列x(n)n＝0，1，...，N-1，设采样率为f_sHz，其感兴趣信号谐波分量的频率f₀已知或已测得；(2)计算序列x(n)幅频谱|X(f)|和相频谱φ(f)在全频域的γN个等间隔的离散值|X(k)|和φ(k)k＝0，1，...，γN-1(其中γ＞1)；(3)计算幅频谱序列|X(k)|k＝0，1，...，γN-1中幅度值大于一阈值th的序列频率序号集合，记为K_th；(4) 求解一优化问题，以确定感兴趣的以f₀为频率的谐波分量初相位的精确估值，不妨记初相位为

其中，优化问题定义为

其中，

为自变量，

为构造序列

的相频谱序列

的子集，即

f₀已知，φ_th(k)为待分析原序列相频谱φ(k)k＝0，1，...，γN-1的子集，即φ_th(k)＝φ(k)k∈K_th，

为两序列

和φ_th(k)差异大小的范数，构造序列初相

的搜索范围可根据先验知识确定。简言之，搜索出一最优初相点

使得构造序列和待测序列的相频谱在K_th频率点集合中的差值的范数最小，则将之估计为信号感兴趣f₀谐波分量的初相

如图1所示。

以上方案的原理是：(a)对于谐波分量频率f₀已知的情况，其构造序列的相频谱序列

只与构造序列

的假定初相

有关；(b)构造序列的假定初相

作为自变量变化时，对应的构造序列相频谱序列

与原待分析序列x(n)的相频谱序列φ(k)的差异大小也发生变化，可看成它是假定初相

的函数；(c)当构造序列的假定相位

刚好等于信号谐波分量真实初相位

时，其相频谱序列

与原序列的相频谱序列φ(k)k＝0，1，..，γN-1的差异出现极小值，据此求一个最小化问题确定信号谐波分量初相位

的精确估值；(d)以上最小化问题定义中，用子集

和φ_th(k)分别取代

和φ(k)k＝0，1，..，γN-1，主要是为了改善数值性，原理一样。

根据以上提出的一种基于计算机程序对信号谐波分量初相位进行精确估计的方法，其第 (4)步最小化问题中

的搜索范围，根据先验知识确定，其特征在于，对|X(f)|的γN个等间隔离散值|X(k)|k＝0，1，..，γN-1，找出已知频率f₀邻域两个最高的相邻谱线，其频率点从左到右分别记为f_-1和f₁，将f_-1和f₁对应的两谐波分量用DTFT方法估计的初相分别记为

和

则将初相位搜索区间确定为

的确定方法如附图2流程框图所示。

根据以上提出的一种基于计算机程序对信号谐波分量的初相位进行精确估计的方法，第 (4)步求解一个最小化问题，其特征在于，在初相位搜索区间

内，按照二分法进行搜索，具体步骤为：(I)在

和

两个初相位点上，分别构造序列

计算幅值谱阈值th以上对应的相频谱序列

和

并计算它们与原序列x(n)的阈值th以上对应的相频谱序列φ_th(k)之差的范数

和

(II)令

构造序列

计算其阈值th上相频谱序列

并计算它与原序列x(n)的阈值上相频谱序列φ_th(k)之差的范数

(III)比较

与

若

则令

否则令

以缩小初相位范围

(IV)若

已经缩小到需要的相位精度，则结束搜索、输出相位估值

否则回到步骤(II)继续搜索。具体操作步骤如附图3所示。

上述二分法求解最小化问题，基于的原理是：构造序列的假定相位

作为自变量变化时，对应的构造序列阈值th上相频谱序列

与原序列x(n)的阈值th上相频谱序列φ_th(k) 之差的范数作为函数是单峰对称的，如附图4所示，因此，通过二分法寻找范数的最小值，相对于通过枚举法寻找，速度提高很多。

有益效果

进一步对附图4说明如下：发明人在相位范围[-10°，10°]内按0.1°的初相位间隔构造一系列正弦信号(采样点数为1000、采样频率为1000Hz、信号频率为50Hz)，并补0做11000点的FFT。另有一个采样点数为1000、采样频率为1000Hz、信号频率为50Hz、初相为0°的带加性白噪声(高斯分布，信噪比为35dB)的正弦信号，同样补0做11000点的FFT相频谱序列。取它们幅频谱大于最大幅度10％的频率点构成K_th。K_th上系列构造信号相频谱与待分析信号相频谱的差值平方和得到了附图4。从图4中可以看出，当假定相位取为0°时，所构造正弦波的相频谱序列和带加性白噪声的50Hz正弦信号的相频谱序列，在K_th上具有最小的差值平方和(2范数)。此例验证了本发明方案所依据原理的正确性。

为进一步验证本发明的效果，产生不同信噪比水平的仿真信号进行了对比性实验。用 MATLAB产生形如

的序列模拟待分析信号。其中，

为仿真信号的真实初始相位，在范围[-10°，10°]内随机取值；f_s＝1kHz为采样频率；f₀为信号频率，A为信号幅值，不失一般性，取A＝1 f₀＝50Hz；ω(n)为均值为零、方差为σ²的白噪声；补零后信号点数为γN＝11N。

实验中，变化ω(n)的大小以产生不同噪声水平的信号，实验分别在[45dB，70dB]范围内每隔5dB在该信噪比水平下随机产生10000个信号，分别用本发明的方法和直接用DTFT计算相位法估计这些信号的相位，并比较它们的效果，两种方法的估计误差均方值MSE比较如表1和表2。其中，表1为信号点数N＝1000即整周期采样时的效果比较；表2为信号点数 N＝1024即非整周期采样时的效果比较。从仿真结果中可以看出，虽然在整周期采样时，本发明的估计精确度不如DTFT直接计算相位法，但是在非整周期采样时，本发明的估计精确度远好于DTFT直接计算相位法。实际中，非整周期采样是一般情况。

表1 本发明方法与DTFT直接计算法估计相位的均方误差比较(整周期采样)

信噪比(dB)	70	65	60	55	50	45
							本发明方法	0.0010	0.0016	0.0029	0.0051	0.0093	0.0168
DTFT直接计算法	0.0006	0.0010	0.0017	0.0032	0.0056	0.0100

表2 本发明方法与DTFT直接计算法估计相位的均方误差比较(非整周期采样)

信噪比(dB)	70	65	60	55	50	45
							本发明方法	0.0023	0.0038	0.0068	0.013	0.0227	0.0392
DTFT直接计算法	0.1391	0.1391	0.1391	0.1391	0.1392	0.1394

附图说明

图1本发明信号谐波分量初相位估计方法流程图

图2本发明初相位初始搜索范围确定的方法流程图

图3本发明在初相位搜索区间二分法求解最小化问题流程图

图4本发明构造信号相频谱和原信号相频谱在K_th上的差值平方和与构造信号初相关系示意图

图5激光相位测距原理

实施例

结合相位式激光测距给出本发明的一个实施例。如附图5所示，激光相位测距的原理，是用一频率比光波频率低得多的正弦波电信号加到激光器上使激光光强按电信号的规律变化，这一过程称为调制。把经调制了的激光发射出去，在被测目标处置一反射器，使光波按原路返回。回来的被调制的光波就有一个相位延迟φ，从图5中可以看出它等于整数n个2π周期再加上一个尾数Δφ，即：φ＝2nπ+Δφ，则光波在两点之间来回飞行的时间

(f为正弦波电信号的频率)，被测距离

(c为光在介质中传播的速度)。

但是目前任何测量相位的方法都不能测出整数个2π周期n的值，只能测定不到2π的尾数Δφ。解决问题的办法是用高、低频率不同的两种正弦波电信号分别去调制激光做两次测量。第一次测量用调制频率足够低(一个周期就足够长)的激光去测量，使激光往返一次回来的被调制的光波相位延迟φ＜2π，测得一个距离。第二次测量用调制频率较高的激光去测量距离的尾数，然后，将两次测量组合起来就可以达到高精度的测量结果。这实际上就是用精度相同、长度和刻度不同的两把尺子组合测量距离。由于相位测距具体实施细节对相位测量的精度并不是本发明的研究目的，为了便于说明本发明的效果，本实施例采用的测量距离和发射信号频率等参数能够使发射波往返仅用不到一个周期，仅比较一次测量的测量精度。

已知实际测量距离为30.0m，发射信号频率f＝2MHz，发射信号初相为0°，则理想情况下回波信号的相位

用如下仿真信号来模拟此场景下的回波信号：

其中，信号幅值A＝1，f₀＝2MHz，

采样点数N＝1024，采样频率f_s＝40MHz，ω(n)是均值为0、方差为σ²的白噪声，此处取方差σ²＝10^-5。利用本发明所述的信号相位估计方法，按照技术方案的步骤，估计此信号的相位。信号相位

估计精度误差界设定为 (10^-6)°。

首先，读取信号x(n)，对x(n)补0做11264点的FFT，得到幅频谱序列|X(k)|和相频谱序列φ(k)k＝0，1，..，11263。计算序列|X(k)|的幅度大于max(|X(k)|)的10％的频率范围K_th，得到φ_th(k)。这里依据先验知识确定初始构造相位为

和

构造两正弦序列

并求出相应的阈值上相频谱序列

和

由于

则确定新的相位搜索范围

因为

继续比较

以不断缩小相位范围直到

经过29次迭代后，达到估计精度，则结束搜索、输出相位估值

由此方法估计得到相位为144.0147°，进而推出距离为30.0031m。若采用DTFT直接计算法，得到的估计相位为144.0539°，由此推出距离为30.0112m。此例中，本方法将误差从0.04％降到了0.01％。

Claims (3)

1.一种基于计算机程序对信号谐波分量的初相位进行精确估计的方法，其特征在于，包括四个步骤，(1)读取待估计信号的采样序列x(n)n＝0，1，...，N-1，设采样率为f_sHz，其感兴趣信号谐波分量的频率f₀已知或已测得；(2)计算序列x(n)幅频谱|X(f)|和相频谱φ(f)在全频域的γN个等间隔的离散值|X(k)|和φ(k)k＝0，1，...，γN-1(其中γ＞1)；(3)计算幅频谱序列|X(k)|k＝0，1，...，γN-1中幅度值大于一阈值th的序列频率序号集合，记为K_th；(4) (4)求解一优化问题，以确定感兴趣的以f_o为频率的谐波分量初相位的精确估值，不妨记初相位为

其中，优化问题定义为

其中，

为自变量，

为构造序列

的相频谱序列

的子集，即

f₀已知，φ_th(k)为待分析原序列相频谱φ(k)k＝0，1，...，γN-1的子集，即φ_th(k)＝φ(k)k∈K_th，

为两序列

和φ_th(k)差异大小的范数，构造序列初相

的搜索范围可根据先验知识确定。

2.根据以上权利要求1所述的一种基于计算机程序对信号谐波分量初相位进行精确估计的方法，其第(4)步最小化问题中

的搜索范围，根据先验知识确定，其特征在于，对|X(f)|的γN个等间隔离散值|X(k)|k＝0，1，..，γN-1，找出已知频率f₀邻域两个最高的相邻谱线，其频率点从左到右分别记为f_-1和f₁，将f_-1和f₁对应的两谐波分量用DTFT方法估计的初相分别记为

和

则将初相位搜索区间确定为

3.根据以上权利要求1所述的一种基于计算机程序对信号谐波分量的初相位进行精确估计的方法，第(4)步求解一个最小化问题，其特征在于，在初相位搜索区间

内，按照二分法进行搜索，具体步骤为：(I)在

和

两个初相位点上，分别构造序列

计算幅值谱阈值th以上对应的相频谱序列

和

并计算它们与原序列x(n)的阈值th以上对应的相频谱序列φ_th(k)之差的范数

和

(II) 令

构造序列

计算其阈值th上相频谱序列

并计算它与原序列x(n)的阈值上相频谱序列φ_th(k)之差的范数

(III)比较

与

若

则令

否则令

以缩小初相位范围

(IV)若

已经缩小到需要的相位精度，则结束搜索、输出相位估值

否则回到步骤(II)继续搜索。

Images