CN103941093B - 一种双向dft对称补偿相位测量方法及其装置 - Google Patents

Abstract

本发明公开了一种双向DFT对称补偿相位测量方法及其装置，方法包括：分别对样本做前向和后向两次点数均为N的DFT变换；提取相位测量值和对这两个相位测量值做叠加平均即得{x(n)}中心样点处的相位测量结果。装置包括：模数转化器、DSP器件、驱动和显示模块，待测频率为f₀的信号经过模数转化器采样得到样本序列x(n)，以并行数字输入的形式进入DSP器件，经过DSP器件的内部算法处理，得到信号的参数估计；驱动和显示模块显示出相位测量结果、测相方差和克拉美罗理论下限。本方法对相位进行了无偏估计，填补了对信号相位参数做有效精确估计方法的空白，实现了对复指数信号的高精确相位测量。

Description

一种双向DFT对称补偿相位测量方法及其装置

技术领域

本发明涉及数字信号处理领域，尤其涉及一种双向DFT对称补偿相位测量方法及其装置。

背景技术

经典复指数信号a·exp[j(ω₀n+θ₀)]的相位估计是在通信^[1]、GPS导航^[2]、电力谐波分析^[3]、光学工程^[4]等领域广泛遇到的重要问题，引起国内外相关学者的普遍关注。

相位估计可分为模拟方法和数字方法两种。数字化相位估计因其准确度高、测量过程灵活且不需要引入复杂电路等优势而成为相位估计方法的主流，吸引着学术界和工程界不断致力于其算法的研究。对于数字化相位估计法，衡量其测量精度的指标，就是在给定样本长度M和信噪比(Signal to Noise Ratio,SNR)条件下，由多次测量结果统计算出的均方误差(Mean Squared Error,MSE)。在信号检测与估计理论中，克拉美罗下限(Cramer-RaoLower Bound,CRLB)是衡量测相方法的尺度^[5,6]。即在给定样本长度M和信噪比条件下，由任何测量方法得到的均方误差都不能低于CRLB，只能逼近CRLB，因而测量方法MSE逼近CRLB的程度反映了该方法的测量能力。

现有的数字测量方法有：希尔伯特变换法^[7]、正弦曲线拟合法^[8]，以及各种基于离散傅里叶变换(Discrete Fourier Transform,DFT)谱校正的相位测量法^[9-11]。其中，相比于希尔伯特变换法、正弦曲线拟合法，离散傅立叶变换更常被用于工程界测相应用中。这是因为：(1)离散傅立叶变换结果为复数，含有丰富相位信息；(2)DFT有快速算法(即快速傅里叶变换，Fast Fourier Transformation,FFT)，FFT大大加快了DFT的运算速度，非常易于处理实时信号；(3)可以借助一些内插、修正等措施对DFT的谱分析结果做进一步校正而提升DFT的测相精度。如文献[12]利用插值离散傅立叶变换对信号进行初相位估计；文献[13]对不同窗函数下插值离散傅立叶变换法估计信号初相位的性能作了分析；文献[14]利用离散傅立叶变换频谱的相位信息来提高初相位的估计准确度，等等。此外，近似细化法^[15]、比值法^[16]、相位差法^[17,18]、余弦窗校正法^[19]都是基于传统DFT的相位测量法。

但是文献[7-19]所提出的这些测相法，都需先估测信号频率值后再利用频率值去校正相位，因而频率估计误差会带到相位估计误差中去而导致估计性能变差。正如文献[19]所指出：“频率估计的偏差会带入到相位估计中，从而增加了相位估计的不确定性”毫无疑问，这种不确定性使得推导相位估计的理论性能变得非常困难，因而文献[7-19]的相位估 计法都没有推导出相位估计方差的理论表达式。以上基于DFT的各种测相方法，由于相位测量过程中无法避免频率估计，故其测相涉及估计三个参数[a,ω₀,θ₀]^T，因而相位估计方差也被Rife在文献[5]所推导的三参数CRLB(记为CRLB₃)所限定。

现有技术中提出的全相位FFT(all-phase FFT,apFFT)测相法^[20,21]，其测量精度高于文献[7-19]各种测量方法。因而如何设计出测量精度比apFFT测相法更高的测相方法，是一项具有挑战性而又急需完成的工作。

另外，给定样本长度M和信噪比条件下，测相方法的测相均方误差的理论解析式及其均方误差的测相理论下限也是急需推导揭示出来的。这是因为，在实际工程测量中，若把握了均方误差与样本长度M和信噪比的内在联系与规律，则可以有针对性地确定系统的样本存储空间和调整测量现场的噪声环境，从而有利于缩短基于测相的工程开发进度。

发明内容

本发明提供了一种双向DFT对称补偿相位测量方法及其装置，本发明实现了对给定样本复指数信号样本的高精确相位测量，并且可定量计算出相位估计的理论均方误差，详见下文描述：

一种双向DFT对称补偿相位测量方法，所述方法包括以下步骤：

（1）分别对样本做前向和后向两次点数均为N的DFT变换；

（2）从前向DFT谱的峰值谱位置k=k^*和后向DFT谱的对称峰值谱位置k=N-k^*中分别提取出对应的相位测量值和

（3）对这两个相位测量值和做叠加平均即得{x(n)}中心样点x(0)处的相位测量结果。

所述相位测量值具体为：

所述相位测量值具体为：

其中，θ₀为信号相位，δ为频偏值，k^*为峰值谱位置值。

通过所述方法获得的测相方差为：

其中，ρ为信噪比值，sinc(δ)=sin(δπ)/(δπ)，δ为频偏值。

通过所述方法获得的克拉美罗理论下限CRLB₂为：

其中，ρ为信噪比值，M=2N-1，a为信号幅度，σ²为方差。

一种双向DFT对称补偿相位测量装置，包括：模数转化器、DSP器件、驱动和显示模块，将待测频率为f₀的信号x(t)经过所述模数转化器采样得到样本序列x(n)，以并行数字输入的形式进入所述DSP器件，经过所述DSP器件的内部算法处理，得到信号的参数估计；最后借助所述驱动和显示模块显示出相位测量结果、测相方差和克拉美罗理论下限。

本发明提供的技术方案的有益效果是：

第一、本专利对相位进行了无偏估计，填补了对信号相位参数做有效精确估计方法的空白，实现了对复指数信号的高精确相位测量。原有的数字信号处理的测相法，都需先估测信号频率值后再利用频率值去校正相位，因而频率估计的误差会带到相位估计误差中去而导致估计性能变差。本专利使用了与apFFT同样长度的样本数据，从前、后向DFT中综合了可相互补偿且精确的相位信息，经过理论推导，本方法可以完全抵消因未知频偏δ而造成的测相误差，有噪无噪时均为无偏估计，实现了高精度的相位测量。

第二、文献[7-19]的相位估计法都没有推导出相位估计方差的理论表达式，本发明可定量计算出相位估计的理论均方误差。由式(7)推导出的均方误差可知，该方差与频偏δ、信噪比ρ和样本量N有关，因此可进行定量计算，利于工程应用。

第三、本发明提供了双向DFT对称补偿测相法的均方误差理论下限。对于单频信号测相，当信噪比ρ足够大时，对比双向DFT对称补偿测相法和克拉美罗限可知：本发明对应的克拉美罗限仅为传统三参数模型克拉美罗限的1/4。

第四、本发明计算量较小，估计方法简单且效率高，资源耗费少，大大节省了硬件成本。对于M=2N-1的样本x(n)，截取分为仅包含唯一共同样点x(0)的前、后两段，分别对样本做前向和后向两次点数均为N的DFT变换得到前向、后向序列峰值谱位置的相位谱值，此过程可利用前、向后向DFT的性质和傅里叶变换的性质，这就减少了很多不必要的计算，大大加快了计算速度，效率得以提高。

附图说明

图1为一种双向DFT对称补偿相位测量方法的流程示意图；

图2分别为δ=0.1、δ=0.2、δ=0.3、δ=0.4、δ=0.5的两种测相方差曲线及克拉美罗 界曲线对照图；

图3为一种双向DFT对称补偿相位测量装置的结构示意图；

图4为DSP内部程序流图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面对本发明实施方式作进一步地详细描述。

本发明提出的双向DFT对称补偿法即符合此要求，该方法使用了与apFFT同样长度的样本数据，从前、后向DFT中综合了可相互补偿且精确的相位信息。本发明还推导出了新方法相位估计的理论均方误差和适用于本方法的测相方差的克拉美罗理论下限。由于本专利方法的相位估计过程避开了频率估计，故与传统DFT测相法不同，该理论下限符合两参数(记为CRLB₂)估计模型，而不是传统的三参数估计模型，其均方误差突破了传统三参数CRLB₃。本专利还通过实验证明了其均方误差介于apFFT测相方差和CRLB₂之间，在信噪比较高时测相精度高于apFFT测相法，因而有较高应用价值。

101：分别对样本做前向和后向两次点数均为N的DFT变换；

如图1所示，令样本长度为M，N=(M+1)/2，给定复指数序列{x(n)=a·exp[j(ω₀n+θ₀)]，-N+1≤n≤N-1}，其中，a为信号幅度，θ_o为信号相位。本发明旨在估计出中心样点x(0)的相位θ₀x。

以样本长度M=7为例，N=(M+1)/2=4，分别对样本做前向和后向两次点数均为N的DFT变换，该步骤的详细操作为本领域技术人员所公知，本发明实施例对此不做赘述。

102：从前向DFT谱的峰值谱位置k=k^*和后向DFT谱的对称峰值谱位置k=N-k^*中分别提取出对应的相位测量值和

假定频率分辨率Δω=2π/N，对于频率ω₀不妨做如下设定：

ω₀＝βΔω＝(k^*+δ)Δω，k^*∈z⁺，|δ|≤0.5 (1)

其中，β为频率分辨率倍数，k^*为峰值谱位置值，δ为频偏值，z⁺为正整数集。把该序列截取分为仅包含唯一共同样点x(0)的前、后两段{x(n)，0≤n≤N-1}与{x(n)，-N+1≤n≤0}。分别对样本为N的序列进行DFT变换，得到前向序列峰值谱位置的相位谱值为：

类似地，后向序列峰值谱位置的相位谱值为

103：对这两个相位测量值和做叠加平均即得{x(n)}中心样点x(0)处的相位测量结果。

联立式(2)和(3)，将与进行叠加取平均，可得

式(4)表明，与叠加后可完全抵消因未知频偏δ而造成的测相误差。

无噪时，相位估计的每一步都是没做近似的恒等推导，因而根据式(4)所得到的相位估计值即为真实相位θ₀；有噪时，在信噪比SNR值并不特别低(例如下文实验中信噪比SNR值大于10dB)的情况下，本方法提出的测相估计为无偏估计。

一、均方误差

假设w(n)为均值为0，方差为σ²的复高斯白噪声，则信号模型为x(n)=a·exp[j(ω₀n+θ₀)]+w(n)，因而信噪比SNR值ρ=a²/σ²，当ρ不特别低且N足够大时(例如，信噪比SNR值ρ=20dB，N=32)，对噪声取方差，可得前向DFT的测相方差为：

其中，sinc(δ)=sin(δπ)/(δπ)。由于前后向DFT的对称性，可得后向DFT的测相方差为：

由于前后向数据段仅有1个重合样点x(0)，故与可认为相互独立，因而推出测相方差为：

二、均方误差的理论下限

对于本专利的复信号模型，利用联合概率密度函数和Fisher信息矩阵，根据经典参数估计理论，可以得到本发明测相方差克拉美罗理论下限CRLB₂为

式(8)中，ρ为信噪比，M=2N-1，即为本发明所耗费的样本个数。式(8)所推出的克拉美 罗限是衡量测相方法的一个参考指标。对于本发明提出的测相方法而言，若其测量方差远大于CRLB₂，说明其测相精度低；若其测量方法略大于CRLB₂，则说明其测量精度高。因此后面给出的实验中，会将本发明的测相曲线与CRLB₂理论曲线做对照。

此外，与本发明无需估计频率而直接测相的情况不同，传统的各类需先估计频率，再根据频率估计相位的方法也存在一个理论下限CRLB₃，由于M=2N-1，根据文献推出的理论式子，可推导出CRLB₃为

比较CRLB₂和CRLB₃，有

从式(10)可以看出，由于本发明的测相无需依据频率估计结果，其测相参数估计模型可简化为两参数模型，对应的克拉美罗限仅为传统三参数模型克拉美罗限的1/4（即对应差别10log₁₀4≈6dB）。

下面以具体的实验来验证本发明的可行性，详见下文描述：

令N=32，对{x(n)=a·exp[j(3+δ)2π/(Nn)+π/3)]，-N+1≤n≤N-1}，分别用本方法、apFFT测相法进行测相。设定其频偏值δ在0.1～0.5间变化，对于每种频偏情况，做1000次Monte-Carlo测相模拟，并统计测相方差。图2(a)～图2(e)给出了测相方差曲线，并给出了CRLB₂和CRLB₃曲线做对照。

从图2(a)～图2(e)，可总结如下规律：

1、本方法与全相位FFT测相法一样，其测相理论方差是被两参数模型的CRLB₂界定，而不是被CRLB₃界定。例如，图2(a)～图2(c)给出的δ=0.1～0.3的情况，本专利方法(‘o’标记)和apFFT测相法(‘□’标记)的测相方差都突破了CRLB₃，并且逼近CRLB₂。因而两种方法都是高精度测相法。

2、本方法与全相位FFT测相法一样，测相精度与频偏值δ有关。即δ越小，方差越小，精度越高。从图2(a)～图2(e)，随着δ值从0.1增加到0.5，对称补偿DFT测相法与全相位FFT测相法的测相方差曲线偏离CRLB₂越来越远。

3、在信噪比不太低的区域，图2(a)～图2(e)中测相法的实测方差曲线与由式(7)算出的 理论方差曲线吻合。这证明了式(7)理论表达式的正确性。

4、在信噪比不太低的区域，本方法的测相方差曲线比apFFT方差曲线更靠近CRLB₂，因而具有更高的估计精度。

参见图3，该双向DFT对称补偿相位测量装置包括：A/D（模数转化器）、DSP器件、驱动和显示模块，将待测频率为f₀的信号x(t)经过A/D（模数转化器）采样得到样本序列x(n)，以并行数字输入的形式进入DSP器件，经过DSP器件的内部算法处理，得到信号的参数估计；最后借助驱动和显示模块显示出相位的估计值、相位估计的理论均方误差和均方误差的理论下限。

其中，图4的DSP（Digital Signal Processor，数字信号处理器）为核心器件，在信号参数估计过程中，完成如下主要功能：

(1)调用核心算法，完成接收信号的参数估计处理；

(2)根据实际需要调整采样率f_s，使之在该采样率条件下的频率分辨率Δf=f_s/N满足其实际频偏值f₀/Δf-[f₀/Δf]与事先设定的频偏值δ值一致(其中“[.]”表示四舍五入取整)。同时设定测量的重复次数i(如i=1000)进行均方误差的计算。

例如：若信号频率f₀为3.3Hz，δ设为0.3，令N=256，采集2N-1=511个样点，则通过DSP与A/D转换器的反馈通路，将采样速率f_s调整为256Hz，这时恰好满足f₀/Δf-[f₀/Δf]=3.3-3=0.3。

(3)将相位估计结果实时输出至驱动和显示模块。

需指出，由于采用了数字化的估计方法，因而决定了图3系统的复杂度、实时程度和稳定度的主要因素并不是图3中DSP器件的外围连接，而是DSP内部程序存储器所存储的核心估计算法。

本发明将所提出的“双向DFT对称补偿相位测量法”这一核心估计算法植入DSP器件内，基于此完成高精度、低复杂度、高效的相位估计。图4流程分为如下几个步骤：

(1)首先需根据具体应用要求（例如：医学和军事等的具体测量要求），设置信号的采样点数2N-1，频偏值δ和重复测量的次数i，并根据具体需要设定精度要求。该步骤是从工程方面提出具体需求，以使得后续流程有针对性地进行处理。

(2)然后，DSP器件从I/O端口读采样数据，进入内部RAM。

(3)后续的“去直流处理”，是为了消除待测信号中的直流成分的影响。否则，直流成分的存在，会降低测量精度。直流成分很容易测出，仅需计算样点的平均值即可得到。

(4)按图1本发明的处理过程进行相位测量是DSP算法最核心的部分，运行该算法后，即可得到相位测量值。

(5)判断本方法是否满足工程需求，若不满足，程序返回，重新根据要求设定采样频率进行下一轮相位测量。

(6)直至测量结果符合工程要求。重复进行以上测量过程i次。

(7)通过DSP的输出总线输出至外部显示驱动设备，将相位测量结果、相位估计的理论均方误差和均方误差的理论下限进行数码显示。

需指出，由于采用了DSP器件实现，使得整个参数估计操作变得更为灵活，可根据信号所包含的各种分量的具体情况，通过编程灵活改变算法的内部参数设置，如采样点数2N-1，采样率f_s等。

本发明实施例对各器件的型号除做特殊说明的以外，其他器件的型号不做限制，只要能完成上述功能的器件均可。

参考文献

[1]X.Yi,W.Shieh,Y.Tang.Phase estimation for coherent optical OFDM[J].Photonics Technology Letters,IEEE,2007,19(12):919-921.

[2]A.Carta,N.Locci,C.Muscas,et al..A flexible GPS-based system forsynchronized phasor measurement in electric distribution networks[J].Instrumentation and Measurement,IEEE Transactions on,2008,57(11):2450-2456.

[3]唐求,王耀南,郭斯羽.电力系统谐波及其检测方法研究[J].电子测量与仪器学报,2009,23(5):29-33.

[4]J.Kong,Z.Ou,W.Zhang.Phase-measurement sensitivity beyond thestandard quantum limit in an interferometer consisting of a parametricamplifier and a beam splitter[J].Physical Review A,2013,87(2):023825-1-023825-6.

[5]D.Rife,R.Boorstyn.Single tone parameter estimation from discrete-time observations[J].Information Theory,IEEE Transactions on,1974,20(5):591-598.

[6]齐国清.离散实正弦信号参数估计的Cramer Rao方差下限[J].数据采集与处理,2003,8(2):151-155.

[7]王凤鹏,邹万芳,尹真,等.希尔伯特变换实时全息干涉条纹相位提取[J].光电工程,2009,36(4):92-96.

[8]梁志国,张大治,孙璟宇,等.四参数正弦波曲线拟合的快速算法[J].计测技术,2006,26(1):4-7.

[9]丁康,钟舜聪.通用的离散频谱相位差校正方法[J].电子学报,2003,31(1):142-145.

[10]段虎明,秦树人,李宁.离散频谱的校正方法综述[J].振动与冲击,2008,26(11):138-145.

[11]齐国清.利用FFT相位差校正信号频率和初相估计的误差分析[J].数据采集与处理,2003,18(1):7-11.

[12]Jain V.K.,Collins W.L.,Davis D.C.High-accuracy analogmeasurements via interpolated FFT[J].IEEE Trans on IM,1979,28(2):113-122.

[13]Dusan Agrez.Improving phase estimation with leakage minimization[J].IEEE Trans on IM,2005,54(4):1347-1353.

[14]齐国清，贾欣乐.基于DFT相位的正弦波频率和初相的高精度估计方法[J].电子学报,2001,29(9):1164-1167.

[15]霍兵勇,易伟建.密集频率数字信号的判定和校正方法[J].振动与冲击,2013,32(2):171-174.

[16]黄翔东,朱晴晴,杜宇彬,等.有限带宽信号的双谱线比值压缩感知重构算法[J].信号处理,2012,28(6):793-798.

[17]安国臣,张秀清,王晓君,等.离散频谱的加权相位差校正方法[J].信号处理,2012,28(1):99-104.

[18]胡文彪,夏立,向东阳,吴正国.一种改进的基于相位差法的频谱校正方法[J].振动与冲击,2012,31(1):162-166.

[19]M.Novotny,D.Slepicka,M.Sedlacek.Uncertainty analysis of the RMSvalue and phase in the frequency domain by noncoherent sampling[J].Instrumentation and Measurement,IEEE Transactions on,2007,56(3):983-989.

[20]王兆华,黄翔东,杨尉.全相位FFT相位测量法[J].世界科技研究与发展,2007,29(4):28-32.

[21]Huang X,Wang Z,Ren L,et al.A novel high-accuracy digitalizedmeasuring phase method[C]//Signal Processing,2008.ICSP2008.9th InternationalConference on.IEEE,2008:120-123.

本领域技术人员可以理解附图只是一个优选实施例的示意图，上述本发明实施例序号仅仅为了描述，不代表实施例的优劣。

以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (4)

1.一种双向DFT对称补偿相位测量方法，其特征在于，所述方法包括以下步骤：

(1)分别对样本做前向和后向两次点数均为N的DFT变换；

(2)从前向DFT谱的峰值谱位置k＝k^*和后向DFT谱的对称峰值谱位置k＝N-k^*中分别提取出对应的相位测量值和

(3)对这两个相位测量值和做叠加平均即得{x(n)}中心样点x(0)处的相位测量结果；

所述相位测量值具体为：

所述相位测量值具体为：

其中，θ₀为信号相位，δ为频偏值，k^*为峰值谱位置值；

将与进行叠加取平均，可得

2.根据权利要求1所述的一种双向DFT对称补偿相位测量方法，其特征在于，通过所述方法获得的测相方差为：

其中，ρ为信号功率与噪声功率的比值，sinc(δ)＝sin(δπ)/(δπ)，δ为频偏值。

3.根据权利要求1所述的一种双向DFT对称补偿相位测量方法，其特征在于，通过所述方法获得的克拉美罗理论下限CRLB₂为：

${\mathrm{CRLB}}_2=\frac{{\mathrm{\σ}}^2}{2M\·a^2}=\frac{1}{2M(a^2/{\mathrm{\σ}}^2)}=\frac{1}{2M\mathrm{\ρ}}$

其中，ρ为信号功率与噪声功率的比值，M＝2N-1，a为信号幅度，σ²为方差。

4.一种用于权利要求1所述的双向DFT对称补偿相位测量方法的装置，包括：模数转化器、DSP器件、驱动和显示模块，其特征在于，

将待测频率为f₀的信号x(t)经过所述模数转化器采样得到样本序列x(n)，以并行数字输入的形式进入所述DSP器件，经过所述DSP器件的内部算法处理，得到信号的参数估计；最后借助所述驱动和显示模块显示出相位测量结果、测相方差和克拉美罗理论下限。