CN107305223A - 一种改进的相位差频率估计方法 - Google Patents

Abstract

本发明公开了一种改进的相位差频率估计方法，包括以下步骤：输入2N‑1个样本,对样本作无窗全相位快速傅立叶变换得频谱Y(k)，样本中选N个点作无窗快速傅立叶变换得频谱X(k)；对X(k)和Y(k)搜索，得谱峰位置，获得谱峰相位和得差值相位调整得除τ＝(N‑1)/2得频偏估计 构造序列样本时域调制，获得频偏估计与相加得频偏估计 频偏估计 代入公式求频率估计。本发明引入了频移补偿措施，提出了相位差整周模糊的解决措施，提高了其实用性。

Description

一种改进的相位差频率估计方法

技术领域

本发明涉及数字信号处理领域，更具体的说，是涉及一种改进的相位差频率估计方法。

背景技术

频率估计是数字信号处理的重要内容，特别是对含有高斯白噪声的信号进行频率估计一直是信号处理的经典课题，频率估计不仅在理论上，在实际应用中也有着重要的研究价值(如多普勒效应检测^[1]、阵列波达方向估计^[2]、振动分析中的转速测量^[3]等都可转化为频率估计问题)。由于FFT(Fast Fourier Transform,快速傅立叶变换)具有高运算效率的优势，因而基于FFT的频率估计一直是学术界长久不衰的研究课题。

目前，频率估计主要有内插法和相位差估计法。为消除FFT的栅栏效应，人们常常通过对FFT峰值谱附近的谱线做内插^[4-7]来提高估计精度，然而文献^[5]指出，FFT内插估计器通常是有偏的。为消除该固有偏差，通常需要采取频移、迭代等措施，对FFT谱峰附近的谱值做进一步修正和细化，如文献^[8]提出修正Rife算法，用该算法得到的初始值做一次牛顿迭代后即可得到高精度频率估计结果，其频率估计方差不但不随信号本身频率波动，而且还接近克拉美-罗限(Cramer-Rao lower Bound,CRB)^[9]；Tsui在文献^[10]指出，获得FFT谱峰位置后，若借助比值内插措施进一步算出距离谱峰左、右各0.5个频率分辨率位置的DTFT(Discrete Time FourierTransform)谱值，就可估算出更精确的谱峰位置，并将以上过程迭代多次，则不仅可消除固有偏差，而且在任意频偏情况下其频率估计方差可完全到达克拉美-罗限。

与内插法不同，相位差估计法^[11-14]仅需单根峰值谱线上的相位信息，即可得到频偏估计结果。且相位差估计器无需引入频移、修正措施，其频率估计就是无偏的^[13]。在对相位参数敏感的场合，如电力系统谐波分析^[15]、光学工程^[16]等，相位差法得到了广泛应用。

全相位FFT(all-phase FFT,apFFT)因兼备优良的抑制谱泄漏性能和高精度测相性能(即“相位不变性”^[17])，故非常适合做频率估计。文献^[11]从apFFT与FFT的相位差中提取了精确的频率信息，提出一种相位差频率估计法(为突出其特征，称之为经典apFFT/FFT相位差法)，该方法的价值已在低频实信号频率估计^[12]、超声波时间测量^[18]、水声波束生成^[19]等领域得到证实。但是，经典apFFT/FFT相位差法在估计精度上仍有提升的空间，其RMSE曲线距离克拉美-罗限仍具有一定的距离；此外，为了提升实际工程应用的价值，该估计器也应该具备对多频信号频率估计的适应性。

发明内容

本发明的目的是为了克服现有技术中的不足，提供一种改进的相位差频率估计方法，引入了频移补偿措施对原估计器做根本改进，提出了相位差整周模糊的解决措施，提高了其实用性。

本发明的目的是通过以下技术方案实现的。

本发明的一种改进的相位差频率估计方法，包括以下步骤：

(1)输入M＝2N-1个样本x(-N+1)～x(N-1),N为2的正整数次幂，对样本x(-N+1)～x(N-1)作无窗全相位快速傅立叶变换得到频谱Y(k)，同时从样本x(-N+1)～x(N-1)中选取N个点x(0)～x(N-1)作无窗快速傅立叶变换得到频谱X(k)；

(2)分别对频谱X(k)和频谱Y(k)进行搜索，得到其对应的谱峰位置k^*，进而获得谱峰对应的相位和求得两者的差值进行相位调整得到再除以τ＝(N-1)/2得到频偏估计

(3)若(ζ为给定小阈值)，令同时构造序列-N+1≤n≤N-1，并用序列s(n)对样本x(-N+1)～x(N-1)进行时域调制，将调制后的样本重新执行步骤(1)和步骤(2)，重新获得频偏估计再与相加得到校正后的频偏估计若则频偏估计

(4)将步骤(3)得到的频偏估计代入公式求出频率估计

所述步骤(2)中对频谱X(k)和频谱Y(k)搜索得到m个谱峰位置k^*时(m＝2,3……)，重复步骤(2)至步骤(4)算出各个谱峰位置对应的频率估计(i＝0,1,2....)。

所述步骤(2)中相位调整后得到为

与现有技术相比，本发明的技术方案所带来的有益效果是：

(1)本发明仅利用单根谱线信息便可以精确地估计频率，通过改变加窗方式和引入频移补偿两种措施，提升了估计器精度，使之在单频情况下RMSE逼近克拉美-罗界，与经典apFFT/FFT相位差频率估计器相比，其性能有了较大的提升，更能满足实际工程的需要；

(2)本发明仅需要做一次迭代便可以实现频率的精确估计，计算复杂度低，省去了每次迭代需要做FFT的大量计算，更加节省大量的计算成本，符合工程需要；

(3)本发明通过相位调整，解决了观测相位差可能存在的整周模糊问题，增强了实用性；

(4)本发明不仅适用单频信号频率估计，还实现了估计器对多频信号频率估计的适应性，估计单频信号以及多频率复合信号时都获得了很高的精度，对于单频信号，RMSE曲线已经逼近克拉美-罗界，对于多频复合信号，在低信噪比区域，频率估计精度高，具有很强的抗干扰能力，因此具有广泛的应用前景。

附图说明

图1是本发明的流程示意图；

图2是本发明的硬件实施图；

图3是全相位FFT谱分析流程示意图(N＝4)；

图4是双窗apFFT和加窗FFT的振幅谱；

图5是无窗apFFT和无窗FFT的振幅；

图6是两个单频信号频偏估计的RMSE曲线；

图7多频复合信号频偏估计的RMSE曲线；

具体实施方式

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

如图1所示，本发明的本发明的一种改进的相位差频率估计方法，包括以下步骤：

(1)输入M＝2N-1个样本x(-N+1)～x(N-1),其中N为2的正整数次幂，对样本x(-N+1)～x(N-1)作无窗全相位快速傅立叶变换得到频谱Y(k)，同时从样本x(-N+1)～x(N-1)中选取N个点x(0)～x(N-1)作无窗快速傅立叶变换得到频谱X(k)；

(2)分别对频谱X(k)和频谱Y(k)进行搜索，得到其对应的谱峰位置k^*，进而获得谱峰对应的相位和求得两者的差值并进行相位调整得到再除以τ＝(N-1)/2得到频偏估计

(3)对步骤(2)中得到的频偏估计进行大小判断：若(ζ为给定小阈值)，令保存第一次相位差频偏估计的结果，同时构造序列-N+1≤n≤N-1，并用序列s(n)对样本x(-N+1)～x(N-1)进行时域调制，将调制后的样本重新执行步骤(1)和步骤(2)，获得新的频偏估计再与相加得到校正后的频偏估计即若则直接将第一次相位差频率估计得到的频偏估计作为最终的频偏估计即目的是将大频偏自动调小，以此获得更高的频率估计精度。当频偏检测值为较大的正值时，则时域调制操作使得频谱朝负方向搬移，恰好搬到k^*Δf附近；当为较大的负值时，则时域调制操作使得频谱朝正方向搬移，也恰好搬到k^*Δf附近——无论哪种情况，频移后的信号频偏都趋于0，使得无窗apFFT和无窗FFT谱峰都获得大幅值，抵御噪声干扰的能力比频移前都会提升，再经频偏补偿后，可获得更高精度。

(4)将提前设定的N、谱峰位置k^*以及步骤(3)得到的最终的频偏估计代入公式求出归一化频率估计

上述步骤(1)至步骤(4)可以实现对某一单频信号的频率估计，需指出，本方法同样适合于多频信号的频率估计，即当所述步骤(2)中对频谱X(k)和频谱Y(k)搜索得到m个谱峰位置k^*时(m＝2,3……)，重复步骤(2)至步骤(4)逐个算出各个谱峰位置对应的频率估计即可。本发明以DSP(Digital Signal Processor，数字信号处理器)为核心器件，如图2所示。

如图3所示为apFFT谱分析流程图，主要分为两个步骤：

步骤一，用长为(2N-1)的卷积窗w_c对输入数据x(n)加权，然后将间隔为N的数据两两叠加(中间元素除外)而形成N个数据y(0),y(-1),...,y(-N+1)；

步骤二，将步骤一得到的数据y(0),y(-1),...,y(-N+1)做FFT得到离散谱Y(k)；其中，步骤一中的卷积窗w_c是由长为N的前窗f和翻转的长为N的后窗b卷积得到的，即

w_c(n)＝f(n)*b(n),n∈[-N+1,N-1] (1)

令R_N(n)为长度为N的矩形窗，如果f＝b＝R_N，则称图2为无窗全相位FFT，如果f＝b≠R_N，则称图2为双窗全相位FFT。

将式(1)中的前窗f的傅立叶变换写为幅度和相位分离的形式

F(jω)＝F_g(jω)·e^-j(N-1)ω/2 (2)

对于单频的正弦信号0≤n≤N-1，其归一化频率可表示为

ω₀＝2πf₀＝2π(k^*+δ)/N (3)

其中，k^*代表谱峰位置故取值为正整数，δ为频偏估计，取值范围-0.5≤δ<0.5.则信号x(n)的传统加窗FFT谱X(k)为

文献^[17]指出，对于正弦信号-N+1≤n≤N-1的双窗apFFT的幅度谱Y(k)与传统FFT的幅度谱X(k)的平方成正比，其相位谱与频率无关，等于输入中间样点的瞬时相位，即

在峰值谱k＝k^*处，则其传统加窗FFT为

其apFFT的峰值谱为

令和分别表示峰值谱线上观测到的FFT和apFFT相位值，则取式(6)与式(7)的相位差，即可得频率偏离值

在实际应用中，要考虑到峰值谱观测相位和的取值范围都是有限的，即满足故

因而式(8)的观测相位差有可能存在整周模糊问题。为消除模糊，考虑频偏值满足-0.5≤δ<0.5，需对-0.5≤δ<0.5做如下调整

显然经过调整的进而频偏估计和频率估计分别表达式为

从上述分析中可以看出，apFFT/FFT相位差法与基于谱线内插的估计器有所不同。通常，内插谱校正频率估计器除了需要用到峰值谱外，还必须利用附近1～2根旁谱线信息，而相位差法仅需单根峰值谱线(k＝k^*处)的apFFT/FFT相位信息，无需旁谱线参与。因此，为了增强经典apFFT/FFT相位差估计器的抗噪性能和频偏估计精度，我们应该尽量使得能量集中在单根主峰值谱线上。而apFFT和FFT的谱峰幅度及泄漏程度，既与频偏估计δ有关，又与加窗类型有关。下面，我们通过实施例一说明频偏大小与加窗类型对估计精度的影响。

实施例一：

设幅值为1的复信号x(n)＝exp[j(ω₀n+θ₀)],-N+1≤n≤N-1,其中N＝32，ω₀＝(k^*+δ)Δω，k^*＝3，θ₀＝π/3。为了研究频偏估计δ对估计精度的影响，我们让δ分别取0、0.1、0.4。如图4所示，4(a)给出归一化后的双窗apFFT振幅谱(由两个汉明窗卷积而得w_c)，4(b)给出归一化后的加汉明窗的FFT振幅谱。如图5所示，5(a)给出归一化后的无窗apFFT振幅谱，5(b)给出归一化后的无窗FFT振幅谱。表1列出了不同加窗情况下的谱峰幅值。

表1

通过图4、图5以及表1，我们可以总结出二个规律：

(1)无论频偏估计δ如何取值，双窗apFFT峰值谱幅度|Y(k^*)|和加窗FFT峰值谱幅度|X(k^*)|都取不到很大的值；而无窗apFFT和无窗FFT的峰值谱幅度在各个频偏的情况均高于加窗情况；

(2)当频偏估计δ较小时(δ＝0，δ＝0.1)，无窗FFT谱泄漏较小，而无窗apFFT几乎不存在谱泄漏，无窗apFFT峰值谱幅度|Y(k^*)|和无窗FFT峰值谱幅度|X(k^*)|都接近1；当频偏估计δ较大时(δ＝0.4)，无论哪种加窗情况，apFFT和FFT的谱峰峰值都很小，都会有很大一部分能量泄漏到旁谱线上。

需要注意的是，经典的apFFT/FFT相位差法是基于双窗apFFT和加窗FFT做频率估计的，虽然加窗操作可以在大频偏时缩小旁谱线的泄漏范围，但这对于仅需峰值谱信息的相位差频率估计器来说，起不到改善作用；相反，这还要以降低谱峰幅度作为代价，限制其抵御噪声能力，降低频率估计精度。

因此，结合规律(2)中小频偏时谱泄漏情况更小的情况，我们将引入如下两个措施以此提高频率估计精度：1)在选择加窗模式上，采用无窗apFFT和无窗FFT；2)尽可能在δ较小时取相位差。

仿真实验

令N＝32，Δω＝2π/N，考虑两个频率ω₁＝3.1Δω，ω₂＝6.3Δω，对两个单频信号{x_i(n)＝exp(jω_in)+w_i(n)，-N+1≤n≤N-1}，i＝1,2，(w_i(n)为零均值高斯复噪声)，以及复合信号x(n)＝x₁(n)+x₂(n)分别用经典双窗apFFT/FFT相位差法、本发明提出的无窗apFFT/FFT频移补偿相位差法(阈值设为0.1)和Tsui内插估计法[10]，在不同信噪比条件下分别做频率估计。对于每个测试信号的每种SNR情况，做1000次Monte-Carlo测频试验，并统计均方根误差(root-mean-square error,RMSE)。如图6所示，6(a)给出了ω₁＝3.1Δω单频信号测试的RMSE曲线，6(b)给出了ω₂＝6.3Δω单频信号测试的RMSE曲线。如图7所示，7(a)给出了用以上方法对复合信号做测试得到的ω₁＝3.1Δω的RMSE曲线，7(b)给出了用以上方法对复合信号做测试得到的ω₂＝6.3Δω的RMSE曲线，测试过程中对两个频率做了独立的谱峰搜索，并与克拉美-罗限CRB做对照。

同时，文献^[9]给出给定N个样本、信噪比为ρ(单位为能量比值，SNR＝10lg^ρ)的频率估计克拉美-罗限为

由于每次试验各估计器都耗费等长的M＝2N-1个样本，因而相应的CRB为

从图6可以看出，对于单频信号，经过频移补偿和无窗处理这两项改进措施，本发明提出的无窗apFFT/FFT相位差法(‘▽’标记)，相比于经典双窗apFFT/FFT相位差法(‘○’标记)，其频率精度均得以提升，其RMSE曲线更逼近CRB的开方曲线(‘——’标记)。从图7可以看出，对于多频率成分信号，在较低信噪比区域，本发明提出的频率估计方法(‘▽’标记)，相比于经典双窗apFFT/FFT相位差法(‘○’标记)具有更高的精度，其中，‘——’表示CRB的开方曲线，因此本发明所提出的改进方法具有更强的抵抗大噪声能力，这也印证了上面所说的无窗模式小频偏情况更可以凸显单根谱峰幅度的结论。

参考文献

[1]曹文杰,张磊,杜兰,等.基于瞬时频率估计的进动锥体目标微多普勒频率提取方法[J].电子与信息学报,2015,37(5):1091-1096.

[2]Wang Wen-Qin and H.C.SO.Transmit subaperturing for range and angleestimation in frequency diverse array radar[J].IEEE Transactions on SignalProcessing,2014,62(8):2000-2011.

[3]Ansari M,Esmailzadeh E,and Jalili N.Exact frequency analysis of arotating cantilever beam with tip mass subjected to torsional-bendingvibrations[J].Journal of Vibration and Acoustics,2011,133(4):041003.

[4]Aboutanios E and Mulgrew B.Iterative frequency estimation byinterpolation on Fourier coefficients[J].IEEE Transactions on SignalProcessing,2005,53(4):1237-1242.

[5]Candan C.A method for fine resolution frequency estimation fromthree DFT samples[J].IEEE Signal Processing Letters,2011,18(6):351-354.

[6]Duda K.DFT interpolation algorithm for Kaiser–Bessel and Dolph–Chebyshev windows[J].IEEE Transactions on Instrumentation and Measurement,2011,60(3):784-790.

[7]Candan C.Analysis and further improvement of fine resolutionfrequency estimation method from three DFT samples[J].IEEE Signal ProcessingLetters,2013,20(9):913-916.

[8]邓振淼,刘渝.正弦波频率估计的牛顿迭代方法初始值研究[J].电子学报,2007,35(1):104-7.

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

[10]Tsui J.Digital techniques for wideband receivers[M].USA:SciTechPublishing,2004.

[11]黄翔东,王兆华.基于全相位频谱分析的相位差频谱校正法[J].电子与信息学报,2008,30(2):293-297.

[12]谭思炜,任志良,孙常存.全相位FFT相位差频谱校正法改进[J].系统工程与电子技术,2013,35(1):34-39.

[13]黄翔东,孟天伟,丁道贤,等.前后向子分段相位差频率估计法[J].物理学报,2014,63(21):214304.

[14]Huang Xiang-dong and Xia Xiang-gen.A fine resolution frequencyestimator based on double sub-segment phase difference[J].IEEE SignalProcessing Letters,2015,22(8):1055-1059.

[15]曹浩,刘得军,冯叶,等.全相位时移相位差法在电力谐波检测中的应用[J].电测与仪表,2012,49(7):24-28.

[16]姜成昊，杨进华，张丽娟，等.基于激光拍频高准确度相位式测距方法[J].光子学报,2014,43(9):912006.

[17]王兆华,黄翔东.数字信号全相位谱分析与滤波技术[M].北京:电子工业出版社,2009.

[18]方汉方,黄勇,蔡艺剧，等.超声波传输时间的高精度测量[J].信号处理,2012,28(4):595-600.

[19]孙向前，李晴，范展.全相位频谱校正技术在水声通信中的应用研究[J].声学技术,2015,34(2):127-133.

尽管上面结合附图对本发明的功能及工作过程进行了描述，但本发明并不局限于上述的具体功能和工作过程，上述的具体实施方式仅仅是示意性的，而不是限制性的，本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨和权利要求所保护的范围情况下，还可以做出很多形式，这些均属于本发明的保护之内。

Claims (3)

1.一种改进的相位差频率估计方法，其特征在于，包括以下步骤：

(1)输入M＝2N-1个样本x(-N+1)～x(N-1),N为2的正整数次幂，对样本x(-N+1)～x(N-1)作无窗全相位快速傅立叶变换得到频谱Y(k)，同时从样本x(-N+1)～x(N-1)中选取N个点x(0)～x(N-1)作无窗快速傅立叶变换得到频谱X(k)；

(2)分别对频谱X(k)和频谱Y(k)进行搜索，得到其对应的谱峰位置k^*，进而获得谱峰对应的相位和求得两者的差值进行相位调整得到再除以τ＝(N-1)/2得到频偏估计

(3)若(ζ为给定小阈值)，令同时构造序列-N+1≤n≤N-1，并用序列s(n)对样本x(-N+1)～x(N-1)进行时域调制，将调制后的样本重新执行步骤(1)和步骤(2)，重新获得频偏估计再与相加得到校正后的频偏估计若则频偏估计

(4)将步骤(3)得到的频偏估计代入公式求出频率估计

2.根据权利要求1所述的一种改进的相位差频率估计方法，其特征在于，所述步骤(2)中对频谱X(k)和频谱Y(k)搜索得到m个谱峰位置k^*时(m＝2,3……)，重复步骤(2)至步骤(4)算出各个谱峰位置对应的频率估计(i＝0,1,2....)。

3.根据权利要求1所述的一种改进的相位差频率估计方法，其特征在于，所述步骤(2)中相位调整后得到为