CN108414833B - 一种信号分量频率的精确估计方法 - Google Patents
一种信号分量频率的精确估计方法 Download PDFInfo
- Publication number
- CN108414833B CN108414833B CN201810206809.8A CN201810206809A CN108414833B CN 108414833 B CN108414833 B CN 108414833B CN 201810206809 A CN201810206809 A CN 201810206809A CN 108414833 B CN108414833 B CN 108414833B
- Authority
- CN
- China
- Prior art keywords
- frequency
- sequence
- signal
- amplitude spectrum
- search
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
- G01R23/00—Arrangements for measuring frequencies; Arrangements for analysing frequency spectra
- G01R23/16—Spectrum analysis; Fourier analysis
Landscapes
- Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- General Physics & Mathematics (AREA)
- Measuring Frequencies, Analyzing Spectra (AREA)
- Radar Systems Or Details Thereof (AREA)
Abstract
一种基于计算机程序对信号分量的频率进行精确估计的方法,包括三个步骤,(1)读取待分析信号的采样序列x(n)n=0,1,…,N‑1,它有频率为f0的感兴趣信号分量,(2)计算序列x(n)幅频谱|X(f)|在全频域的N个等间隔的离散值,|X(k)|k=0,1,…,N‑1,(3)求解一优化问题以确定感兴趣信号分量频率f0的精确估值,其特点是,不求幅频谱本身的峰值点,而求幅频谱与构造信号幅频谱间相关系数的峰值点,不是信号幅值谱线间两两比较,而是整个幅频谱之间两两比较,将加窗和噪声对频率估计的影响降到了最小。
Description
技术领域
本申请涉及一种基于计算机程序对信号分量的频率进行精确估计的方法。实际中,很多场合都涉及基于计算机程序对信号分量的频率进行精确估计的问题,如在电力***、雷达、生物医学信号处理等领域。
在电力***中,***信号分量的频率是电能质量的重要指标之一,是反映电力***运行状态的重要参数。一般地,在电力***稳定运行的状况下,发电机输出功率与***负荷及损耗维持平衡,电力***信号分量的频率表现为标称值,或者随着负荷之波动在小范围内缓慢地变化。如果出现大容量负荷变化或发电机的投切以及控制设备的不完善,则可能导致较大频率偏移,由此给电力***的稳定运行和用户设备的正常工作带来影响。要实现一个完善的智能化的安全稳定运行的电力***,以保障它的电能质量,检测和精确估计***信号分量的频率是一项必须的核心技术。
在雷达领域中,调频连续波(Frequency Modulated Continuous Wave,FMCW)技术是一种在高精度雷达测距中使用的技术。其基本原理如附图6所示:发射波为高频连续波,其频率随时间按照三角波规律变化;雷达接收的回波的频率与发射的频率变化规律相同,都是三角波规律,只是相同频率有一个时间差Δt,利用这个微小的时间差可计算出目标距离d=c·Δt/2(c为电磁波信号在空气中传播的速度)。其中,时间差Δt的计算,首先要精确估计每个时间点发射信号与回波信号的频差,然后根据峰值频差Δf间接计算Δt;按附图6有,Δt=(Δf·T)/(2ΔF),其中T为频率调制信号周期,ΔF为频率调制信号调频范围。
在生物医学领域中,信号分量的频率估计方法也有很多应用。例如,在采集人体生理信号的时候,不可避免的会采集到50Hz或60Hz市电干扰信号,这对生理信号的后续处理有很大的干扰,应想办法在进一步处理前去除,这也涉及到50Hz或60Hz市电干扰信号频率的精确估计问题,它是高质量陷波处理的前提。
背景技术
最常见的信号分量频率估计方法自然是在频域进行的。其一般原理是:设有一个待估计信号分量之频率的实时间信号x(t),它的采样序列为x(n)n=0,1,...,N-1,采样率为fsHz,则其傅里叶变换如果信号x(t)中有一个以f0Hz为频率的信号分量,则其傅里叶变换的模函数,即幅频响应或者说幅频谱|X(f)|,在f=f0处会出现峰值,据此精确估计信号频率f0。
基于计算机程序估计以上信号分量频率f0,表现为求对应峰值的搜索。根据f0所在范围的先验知识的不同,可分别进行全频谱域的搜索或某一频段内的搜索,全频谱域的搜索区间为[0,fs/2]。关于具体搜索算法,可采用等间隔的枚举搜索算法或不等间隔的启发式搜索算法,也可对两种方式进行某种形式的结合来搜索。全频谱域离散N/2个等间隔的搜索算法,即离散傅里叶变换(DFT)的方法,有快速算法,称为快速傅里叶变换(FFT);在某一局部频段内进行若干等间隔的搜索,也有快速算法,叫ZoomFFT;也可将FFT和ZoomFFT相结合,兼顾搜索精度和搜索效率。不等间隔的基于梯度的或启发式的搜索算法,可采用重心法、二分法、梯度法、进化算法等优化算法。
以上精确估计信号分量频率的一大类方法,它们的关键或者说核心,不是采用何种峰值搜索算法以减小栅栏效应、提高搜索精度,而是在原理上幅频谱|X(f)|在f=f0处是否会出现峰值,或者说峰值频率点与真实信号分量频率f0是否存在偏差。
事实上,按照数字信号处理的理论,由于不可避免的加窗,频谱泄漏和频谱混叠是客观存在的,这会造成幅值谱|X(f)|的峰值频率点与真实频率f0多多少少出现一点偏差;另一方面,信号x(t)或者x(n)中的随机噪声和其他干扰分量,也会造成幅值谱|X(f)|的峰值频率点与真实频率f0出现一点偏差。因此,以上方法在理论上就是一种不完美的情况。
衡量信号分量频率估计算法精确度的指标,可采用估计的频率与实际频率的均方误差(mean-square error,MSE)。在一定信噪比水平下,一个信号分量频率估计方法能达到的均方误差有一个理论上的极限,称为克拉美罗界(Cramer-Rao Lower Bound,CRB)。前人尚没有找到一种不同噪声水平下能够逼近克拉美罗界的频率估计方法,因此,信号分量频率的准确估计还有很大的可改进的空间。
参考文献:
[1]B.Boashash,“Estimating and interpreting the instantaneousfrequency of a signal.I.Fundamentals,”Proc.IEEE,vol.80,pp.520-538,Apr.1992.
[2]B.Boashash,“Estimating and interpreting the instantaneousfrequency of a signal.II.Algorithms and applications,”Proc.IEEE,vol.80,pp.540-568,Apr.1992.
[3]D.C.Rife,R.R.Boorstyn,Single-tone parameter estimation fromdiscrete-time observations,IEEE Trans.Inform.Theory 20(1974)591-598.
[4]E.Aboutanios,A modified dichotomous search frequency estimator,IEEE Signal Process.Lett.11(2004)186-188.
[5]Y.V.Zakharov,T.C.Tozer,Frequency estimator with dichotomous searchof periodogram peak,Electron.Lett.35(1999)1608-1609.
[6]B.G.Quinn,Estimating frequency by interpolation using Fouriercoefficients,IEEE Trans.Signal Process.42(1994)1264-1268.
[7]B.G.Quinn,Estimation of frequency,amplitude,and phase from the DFTof a time series,IEEE Trans.Signal Process.45(1997)814-817.
[8]B.G.Quinn,E.J.Hannan,The Estimation and Tracking of Frequency,Cambridge University Press,New York,2001.
[9]M.D.Macleod,Fast nearly ML estimation of the parameters of real orcomplex single tones or resolved multiple tones,IEEE Trans.Signal Process.46(1998)141-148.
[10]E.Aboutanios,B.Mulgrew,Iterative frequency estimation byinterpolation on Fourier coefficients,IEEE Trans.Signal Process.53(2005)1237-1241.
[11]I.Perisa,J.Lindner,Employing simple FFT-interpolation forimproved complex tone detection and fine estimation,in:3rd ISWCS,2006,pp.744-748.
发明内容
发明目的
提出一种信号分量频率的精确估计方法,使频率估计均方误差不论信噪比高低都能够趋近于克拉美罗界,接近信号分量频率估计准确性的极限。
技术方案
提出一种基于计算机程序对信号分量的频率进行精确估计的方法,包括三个步骤,(1)读取待分析信号的采样序列x(n)n=0,1,...,N-1,采样率为fsHz,它有频率为f0的感兴趣信号分量,(2)计算序列x(n)幅频谱|X(f)|在全频域的N个等间隔的离散值,|X(k)|k=0,1,..,N-1,(3)求解一优化问题从而确定感兴趣信号分量频率f0的精确估值,其特征在于,第(3)步求解的优化问题定义为如下最大化问题
式中,f为自变量,R为幅频域两序列|Xf(k)|k=0,1,...,N-1和|X(k)|k=0,1,..,N-1的相关系数函数,|Xf(k)|k=0,1,...,N-1为构造序列的幅频谱序列,为任一常数,f的搜索范围根据先验知识确定,此最大值问题的最大函数值Rmax对应的频率点fmax值即确定为信号分量频率f0的精确估值。
以上方案的原理是:(a)其构造序列的幅频谱序列|Xf(k)|k=0,1,...,N-1,只与构造序列的假定频率f有关,与相位无关;(b)构造序列的假定频率f作为自变量变化时,对应的构造序列幅频谱序列|Xf(k)|与原序列x(n)的幅频谱序列|X(k)|的相关系数作为函数也发生变化;(c)当构造序列的假定频率f等于信号真实频率f0时,其幅频谱序列|Xf(k)|k=0,1,...,N-1与原序列的幅频谱序列|X(k)|k=0,1,..,N-1的相关系数会出现峰值,据此求一个最大化问题确定信号分量频率f0的精确估值。
根据以上提出的一种基于计算机程序对信号分量频率进行精确估计的方法,第(3)步求解一个最大化问题,见式(1),f的搜索范围根据先验知识确定,其特征在于,对|X(f)|的N个等间隔离散值|X(k)|k=0,1,..,N-1,在此全频谱或感兴趣频段内,找出两个最高的相邻谱线,其频率点从左到右分别记为f-1和f1,则将频率搜索区间确定为finterval=[f-1,f1]。finterval的确定方法如附图2流程框图所示。
根据以上提出的一种基于计算机程序对信号分量的频率进行精确估计的方法,第(3)步求解一个最大化问题,其特征在于,在频率搜索区间finterval=[f-1,f1]内,按照二分法进行搜索,(I)在f-1和f1两个频率点上,分别构造序列 计算相应的幅频谱序列并计算它们与原序列x(n)的幅频谱序列|X(k)|的相关系数(II)令 构造序列计算其幅频谱序列并计算它与原序列x(n)的幅频谱序列|X(k)|的相关系数(III)比较相关系数与若则令f-1=f0,f1=f1,否则令f-1=f-1,f1=f0,以缩小频率范围finterval=[f-1,f1];(IV)若deltf=f1-f-1已经缩小到需要的频率精度,则结束搜索、输出频率估值否则回到步骤(II)继续搜索。具体操作步骤如附图3所示。
上述二分法求解最大化问题,基于的原理是:构造序列的假定频率f作为自变量变化时,对应的构造序列幅频谱序列|Xf(k)|与原序列x(n)的幅频谱序列|X(k)|的相关系数作为函数是单峰对称的,如附图4所示,因此,通过二分法寻找相关系数的最大值,相对于通过枚举法寻找,速度提高很多。
有益效果
进一步对附图4说明如下:发明人在频率范围[49.8Hz,50.3Hz]内按0.1Hz的频率间隔构造系列正弦信号,将它们的幅频谱序列逐个与一频率为50Hz的带加性白噪声的正弦信号的幅频谱序列进行相关而得到了此相关系数图。从图4中可以看出,当假定频率取为50Hz时,所构造正弦波的幅频谱序列和带加性白噪声的50Hz正弦信号的幅频谱序列,具有最大的相关系数。此例验证了本发明方案所依据原理的正确性。
尽管本发明方案也是要求解一个最大化问题,但与传统方案的思路和方法大相径庭。传统方案是求幅频谱本身的峰值点,而本方案是求幅频谱间相关系数的峰值点。传统方案求峰值是幅值谱线之间两两比较,而本方案是整个幅频谱之间两两比较。基于信号处理理论,不难理解,本方案将加窗和噪声的影响降到了最小。
为进一步验证本发明的效果,产生不同信噪比水平的仿真信号进行了对比性实验。用MATLAB产生形如的序列模拟待分析信号。其中,f为仿真的真实信号频率,在范围[49.8Hz,50.3Hz]内随机取值;fs=1kHz为采样频率;为初始相位,A为信号幅值,不失一般性,取A=1,ω(n)为均值为零、方差为σ2的白噪声;采样点数取N=1024。
实验中,变化ω(n)的大小以产生不同噪声水平的信号,分别用本发明的方法和另外两种传统估计方法估计这些信号的频率,并比较它们的效果,如附图5所示。其中,传统方法一为文献[5]的二分搜索FFT峰值法,传统方法二为文献[10]的傅里叶插值算法;横坐标为信号的性噪比实验分别在[-5dB,45dB]范围内每隔5dB在该信噪比水平下随机产生10000个信号并估计;纵坐标为各方法的在这10000次实验中表现出来的频率估计误差均方值MSE与克拉美罗界CRB的比MSE/CRB。从仿真结果中可以看出,本技术方案无论在什么信号信噪比水平下,均方误差都能接近克拉美罗界,即接近精确度的极限,而传统方案不能做到。也就是说,本方案在精确度方面相对现有方案有很大的提高。
附图说明
图1本发明信号分量频率估计方法流程图
图2本发明频率搜索区间的确定方法流程图
图3本发明在频率搜索区间二分法求解最大化问题流程图
图4本发明构造信号幅频谱和原信号幅频谱相关系数函数示意图
图5本发明与两种传统方法信号分量频率估计效果对比图
图6 FMCW雷达测距原理图
实施例
结合雷达测距给出本发明的一个实施例。如附图6所示,已知发射信号波动频率变化量ΔF=200MHz,波动频率调制周期T=5ms,实际测量距离为30.000m,则在理想情况下,对应时间为Δt=0.1μs,其回波与发射波的差拍频率f=16kHz。用如下仿真信号来模拟此场景下的差频信号:
其中,信号幅值A=1,差频信号的频率f=16kHz,初始相位采样点数N=2048,采样频率fs=1MHz,ω(n)是均值为0、方差为σ2的白噪声,此处取方差σ2=0.01。利用本发明所述的信号频率估计方法,按照技术方案的步骤,估计此信号的频率。信号频率f估计精度误差界设定为1Hz。
首先,读取信号x(n),对x(n)做1024点的FFT,得到幅频谱序列记为|X(k)|k=0,1,..,1023;然后,找到|X(k)|的最高的两条谱线为f-1=15625Hz、f1=16113Hz,确定接下来二分法求解最大化问题的频率搜索范围。按(除f外各参数取值与上述仿真信号相同)构造频率为f-1和f1的正弦信号并求出相应的幅频谱序列,然后分别与|X(k)|做相关得到相关系数与比较得则确定新的频率搜索范围 f1=f1=16113Hz;因为f1-f-1>1Hz,继续比较与以不断缩小频率范围直到f1-f-1<1Hz。经过39次迭代后,达到估计精度,则结束搜索、输出频率估值
由此方法估计得的频率16001Hz,进而推出距离为30.002m。若直接取频谱的最高点,得到的估计频率为16016Hz,由此推出距离为30.030m。本方法将误差从0.1%降到了0.0067%。
Claims (3)
1.一种基于计算机程序对信号分量的频率进行精确估计的方法,包括三个步骤,(1)读取待分析信号的采样序列x(n)n=0,1,…,N-1,采样率为fsHz,它有频率为f0的感兴趣信号分量,(2)计算序列x(n)幅频谱|X(f)|在全频域的N个等间隔的离散值,|X(k)|,k=0,1,…,N-1,(3)求解一优化问题从而确定感兴趣信号分量频率f0的精确估值,其特征在于,第(3)步求解的优化问题定义为如下最大化问题
2.根据权利 要求1所述的一种基于计算机程序对信号分量的频率进行精确估计的方法,在第(3)步求解式(1)所示最大化问题时,f的搜索范围根据先验知识确定,其特征在于,对|X(f)|的N个等间隔离散值|X(k)|,k=0,1,…,N-1,在此全频谱或感兴趣频段内,找出两个最高的相邻谱线,记相应的两频率点从左到右依次为f-1和f1,将频率搜索的初始区间确定为finterval=[f-1,f1]。
3.根据权利 要求1或权利 要求2所述的一种基于计算机程序对信号分量的频率进行精确估计的方法,在第(3)步求解一个最大化问题,其特征在于,在频率搜索区间finteravl=[f-1,f1]内,按照二分法进行搜索,(I)在f-1和f1两个频率点上,分别构造序列计算相应的幅频谱序列|Xf-1(k)|,|Xf1(k)|,k=0,1,…,N-1,并计算它们与原序列x(n)的幅频谱序列|X(k)|的相关系数R(|Xf-1(k)|,|X(k)|),R(|Xf1(k)|,|X(k)|);(II)令f0=(f-1+f1)/2,构造序列计算其幅频谱序列|Xf0(k)|,并计算它与原序列x(n)的幅频谱序列|X(k)|的相关系数R(|Xf0(k)|,|X(k)|);(III)比较相关系数R(|Xf-1(k)|,|X(k)|)与R(|Xf1(k)|,|X(k)|),若R(|Xf-1(k)|,|X(k)|)<R(|Xf1(k)|,|X(k)|),则令f-1=f0,f1=f1,否则令f-1=f-1,f1=f0,即有新的频率范围finteravl=[f-1,f1];(IV)若Δf=f1-f-1已经缩小到需要的频率精度,则结束搜索、输出频率估值f0=(f-1+f1)/2,否则回到步骤(II)继续搜索。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201810206809.8A CN108414833B (zh) | 2018-03-07 | 2018-03-07 | 一种信号分量频率的精确估计方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201810206809.8A CN108414833B (zh) | 2018-03-07 | 2018-03-07 | 一种信号分量频率的精确估计方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN108414833A CN108414833A (zh) | 2018-08-17 |
CN108414833B true CN108414833B (zh) | 2020-11-10 |
Family
ID=63131379
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201810206809.8A Active CN108414833B (zh) | 2018-03-07 | 2018-03-07 | 一种信号分量频率的精确估计方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN108414833B (zh) |
Families Citing this family (8)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN109374966A (zh) * | 2018-10-23 | 2019-02-22 | 国网重庆市电力公司电力科学研究院 | 一种电网频率估计方法 |
CN109342813B (zh) * | 2018-12-24 | 2020-09-01 | 常州工学院 | 一种基于dft和二分法的正弦信号频率估计方法 |
CN112394223B (zh) * | 2020-11-10 | 2022-04-01 | 南京大学 | 一种信号分量频率和初相位的联合估计方法 |
CN113219243B (zh) * | 2021-04-27 | 2021-12-31 | 电子科技大学 | 一种校正频率搜索间隔量化误差的高精度测频方法 |
CN113204005B (zh) * | 2021-05-07 | 2023-05-02 | 北京理工大学 | 一种提高调频连续波激光雷达距离解算精度的方法及装置 |
CN113219248B (zh) * | 2021-05-07 | 2022-09-20 | 南京大学 | 基于时域波形比对的信号分量估计方法 |
CN114280366B (zh) * | 2021-12-21 | 2023-10-31 | 中国航天科工集团八五一一研究所 | 基于改进型频率插值算法的正弦信号频率估计方法 |
CN118068287A (zh) * | 2024-04-18 | 2024-05-24 | 天津健君科技有限公司 | 一种基于雷达滤波检测的生理数据处理方法 |
Family Cites Families (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
AU2012233133B2 (en) * | 2011-03-30 | 2014-11-20 | Exxonmobil Upstream Research Company | Convergence rate of full wavefield inversion using spectral shaping |
CN103501201B (zh) * | 2013-09-03 | 2016-01-06 | 华南理工大学 | 一种基于线性调频信号的跳频脉位编码水声通信方法 |
-
2018
- 2018-03-07 CN CN201810206809.8A patent/CN108414833B/zh active Active
Non-Patent Citations (2)
Title |
---|
插值FFT估计正弦信号频率的精度分析;齐国清等;《电子学报》;20040425(第04期);全文 * |
改进的DFT插值频率估计算法及其DSP实现;郑威等;《数据采集与处理》;20170515(第03期);全文 * |
Also Published As
Publication number | Publication date |
---|---|
CN108414833A (zh) | 2018-08-17 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN108414833B (zh) | 一种信号分量频率的精确估计方法 | |
Su et al. | Power harmonic and interharmonic detection method in renewable power based on Nuttall double‐window all‐phase FFT algorithm | |
Zygarlicki et al. | A reduced Prony's method in power-quality analysis—parameters selection | |
Belega et al. | Frequency estimation by two-or three-point interpolated Fourier algorithms based on cosine windows | |
CN107102255B (zh) | 单一adc采集通道动态特性测试方法 | |
Zhang et al. | Analysis of white noise on power frequency estimation by DFT-based frequency shifting and filtering algorithm | |
CN105866543B (zh) | 一种消除基波、谐波对间谐波检测干扰的间谐波检测方法 | |
CN105137373B (zh) | 一种指数信号的去噪方法 | |
Ando et al. | An exact direct method of sinusoidal parameter estimation derived from finite Fourier integral of differential equation | |
CN112051583A (zh) | Fmcw距离测量***中拍频信号非线性校正方法 | |
Giarnetti et al. | Non recursive multi-harmonic least squares fitting for grid frequency estimation | |
Ma et al. | Harmonic and interharmonic analysis of mixed dense frequency signals | |
Li et al. | Frequency estimation based on modulation FFT and MUSIC algorithm | |
CN111579867A (zh) | 一种电力***中谐波和间谐波的测量方法及装置 | |
CN108710029B (zh) | 一种信号谐波分量初相位的精确估计方法 | |
CN112394223B (zh) | 一种信号分量频率和初相位的联合估计方法 | |
Zhang et al. | Micro-motion frequency estimation of radar targets with complicated translations | |
CN105372492B (zh) | 基于三条dft复数谱线的信号频率测量方法 | |
CN109239680A (zh) | 一种低截获概率雷达lfm信号的参数估计方法 | |
CN106053936B (zh) | 一种获取电学信号瞬时频率的方法及*** | |
Huang et al. | Estimation of single-tone signal frequency with special reference to a frequency-modulated continuous wave system | |
Dutra et al. | High-precision frequency estimation of real sinusoids with reduced computational complexity using a model-based matched-spectrum approach | |
Li et al. | A compact super resolution phase retrieval method for Fourier measurement | |
Xiong et al. | High-precision frequency estimation for FMCW radar applications based on parameterized de-alternating and modified ICCD | |
Juhlin et al. | Fast gridless estimation of damped modes |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant | ||
TR01 | Transfer of patent right | ||
TR01 | Transfer of patent right |
Effective date of registration: 20220310 Address after: 811-3, 8th floor, building 4, 209 Zhuyuan Road, Suzhou high tech Zone, Jiangsu Province, 215011 Patentee after: Suzhou greede medical sensor technology Co.,Ltd. Address before: College of electronics, Nanjing University, 163 Xianlin Avenue, Qixia District, Nanjing, Jiangsu, 210023 Patentee before: NANJING University |