一种基于频率相对偏差预估的分段综合单频信号频率估计方法与流程

本发明涉及信号处理技术领域，尤其是一种基于频率相对偏差预估的分段综合单频信号频率估计方法。

背景技术：

对含噪声单频信号进行频率估计信号是处理领域的一个重要研究课题，在声纳、雷达以及电子战等领域具有重要的理论和应用价值，尤其在雷达和水声信号信号处理中所占的地位更为突出。

目前国内外学者提出了很多单频信号频率估计方法，主要有最大似然估计方法(maximumlikelihood，ml)、基于时域相位估计的方法、基于参数模型的现代谱估计方法和基于离散傅里叶变换(discretefouriertransform，dft)频谱校正的方法。其中ml估计方法频率估计方差最小，在高信噪比条件下可接近频率估计的克拉美罗下界(cramer-raolowbound，crlb)，但该方法是一种非线性最小二乘拟合问题，算法复杂度高，运算量大；基于时域相位的方法如线性回归频率估计方法，在高信噪比下可以达到crlb，计算量小，但该方法对信噪比要求高；基于参数模型的现代谱估计方法最大的优点是频率分辨率高，而且频率分辨率不依赖于数据长度，但对该类方法对参数谱模型误差敏感，且运算量大。

基于dft频谱校正的方法，是一种非参数模型方法，物理意义十分明确，可使用fft算法快速实现、实时性好，而且该方法还具有高处理增益且对参数不敏感等优点，是综合考虑算法复杂度、稳健性以及估计精度三个方面性能最佳的一种方法，因此得到了广泛的应用和研究。基于dft频谱校正的方法目前主要有能量重心法、相位差法和插值法等方法。其中能量重心法在高信噪比下，估计精度较高，但是信噪比低时，估计精度降低；相位差法利用同一信号两段相同或不同长度时间序列的傅里叶变换相位结果之间的数量关系来估计信号频率，相当于降低了dft长度，降低了信号处理增益。

插值法是国内外学者研究最多的一种参数估计方法，其中最具代表性的是rife插值法。rife插值法利用单频信号dft谱主瓣内两根谱线的幅值比进行插值，从而估计出信号频率。该方法实现简单，得到的应用和研究最为广泛，但当信号频率与dft谱峰之间的相对频率偏差较小时很容易出现插值方向错误的问题，从而导致频率估计误差偏大。

技术实现要素：

本发明所要解决的技术问题在于，提供一种基于频率相对偏差预估的分段综合单频信号频率估计方法，该方法依据预估频率相对偏差的分布范围和设定的门限，自动选择综合相对偏差的取值，从而进一步利用插值公式估计出单频信号的频率参数，在运算量相当的条件下，频率估计性能较常规的rife插值法有明显改进。

为解决上述技术问题，本发明提供一种基于频率相对偏差预估的分段综合单频信号频率估计方法，包括如下步骤：

(1)获取待处理的单频信号采样数据序列x(n),n＝0,1,…,n-1；从传感器接收n个采样点的实时采集数据作为待处理的数据序列x(n),n＝0,1,…,n-1，或从存储器中提取从检测到信号时刻起始的n个采样点的数据作为待处理的数据序列x(n),n＝0,1,…,n-1，所述的n为检测到的单频信号脉宽长度所对应的采样点个数，取值为2的整数次幂；

(2)对所述数据序列x(n)做快速傅里叶变换，计算得到数据序列的离散傅里叶变换x(l)和功率谱p(k)；

(3)搜索功率谱p(k)最大值所对应的离散频率索引k0，并计算x(l)位于离散频率索引k0-1，k0和k0+1处的模值al，am和ar；

(4)预估有向频率相对偏差δr和加权频率相对偏差δw；

(5)利用al，am,ar，δr和δw计算综合频率相对偏差δc；

(6)将δc带入插值公式估计出单频信号频率即

其中δf为长度为n的离散傅里叶变换的频率分辨率，δf＝fs/n，fs为采样频率。

优选的，步骤(2)中，计算得到数据序列的离散傅里叶变换x(l)和功率谱p(k)，具体计算过程如下：x(n)的离散傅里叶变换为

其中l为x(l)的离散频率索引，j表示虚数单位，即则该数据序列的功率谱为

其中k为p(k)的离散频率索引，||代表取模值运算，|x(l)|为数据序列的幅度谱。

优选的，步骤(2)中，式(1)所示的x(n)的离散傅里叶变换x(l)通过快速傅里叶变换实现。

优选的，步骤(3)中，计算x(l)位于离散频率索引k0-1，k0和k0+1处的模值al，am和ar，计算过程如下：搜索功率谱p(k)最大值所对应的索引k0

其中表示在1≤k≤n/2-1范围内搜索p(k)的最大值所对应的离散频率索引；

取x(l)在离散频率索引k0-1，k0和k0+1处的模值，分别记为al，am和ar，即

al＝|x(k0-1)|,am＝|x(k0)|,ar＝|x(k0+1)|(4)。

优选的，步骤(4)中，预估有向频率相对偏差δr和加权频率相对偏差δw，计算过程如下：利用所述模值al，am和ar计算有向频率相对偏差δr和加权频率相对偏差δw，即

优选的，步骤(5)中，利用al，am,ar，δr和δw计算综合频率相对偏差δc，计算过程如下：首先，计算初始综合相对偏差δm

其中tl为相对频率偏差阈值下限，th为相对频率偏差阈值上限，要求tl和th同时满足：

然后，对δm进行修正得到综合相对偏差δc

优选的，步骤(5)中，当相对频率偏差阈值下限tl＝0.07，相对频率偏差阈值上限th＝0.14，估计效果较好。

本发明的有益效果为：(1)本发明的估计方法通过预估计有向频率相对偏差δr和加权频率相对偏差δw的大小，依据预估频率相对偏差分布范围的不同采用不同的估计方法：当|δr|较大时，即|δr|>th，信号dft主瓣内两根谱线的幅值均较大，抗噪声能力强，且主瓣内次大谱线的幅值远大于旁瓣谱线的幅值，利用rife插值法出现插值方向错误的概率很小，因此比对所有相对偏差δ均利用rife插值法的估计精度高；(2)本发明的估计方法通过预估计有向频率相对偏差δr和加权频率相对偏差δw的大小，依据预估频率相对偏差分布范围的不同采用不同的估计方法：当δ较小时，即δw<th且|δr|<th，受噪声影响，信号dft主瓣内次大谱线的幅值与旁瓣谱线幅值接近，直接利用两者的大小来确定插值方向，容易出现插值方向的错误，而同时利用离散频率索引为k0-1，k0和k0+1处的三根离散线谱模值al，am和ar进行插值，可以较好的避免插值方向的错误，降低频率估计误差；(3)本发明的估计方法首先计算初始综合相对偏差δm，再对δm进行修正得到综合相对偏差δc，可以有效避免如式(7)所示的综合相对偏差估计数学模型与信号dft主瓣模型失配，从而导致频率估计误差急剧增大的问题；(4)本发明的估计方法与现有的rife插值法相比，可以在不增加计算量的前提下，提高频率估计的精度，工程实用性强，适合对信号进行实时处理。

附图说明

图1为本发明的方法流程示意图。

图2为本发明的实施例1中仿真接收单频信号的功率谱示意图。

图3为本发明的实施例1中仿真接收单频信号的幅度谱示意图。

图4为本发明的实施例2中仿真接收单频信号的功率谱示意图。

图5为本发明的实施例2中仿真接收单频信号的幅度谱示意图。

图6为本发明的实施例3中仿真接收单频信号的功率谱示意图。

图7为本发明的实施例3中仿真接收单频信号的幅度谱示意图。

具体实施方式

如图1所示，一种基于频率相对偏差预估的分段综合单频信号频率估计方法，包括如下步骤：

第一步，获取待处理的单频信号采样数据序列x(n),n＝0,1,…,n-1：从传感器接收n个采样点的实时采集数据作为待处理的数据序列x(n),n＝0,1,…,n-1，或从存储器中提取从检测到信号时刻起始的n个采样点的数据作为待处理的数据序列x(n),n＝0,1,…,n-1，所述的n为检测到的单频信号脉宽长度所对应的采样点个数，取值为2的整数次幂；

第二步：对所述数据序列x(n)做快速傅里叶变换，计算得到数据序列的离散傅里叶变换x(l)和功率谱p(k)，计算过程如下：

x(n)的离散傅里叶变换为

其中l为x(l)的离散频率索引，j表示虚数单位，即则该数据序列的功率谱为

其中k为p(k)的离散频率索引，||代表取模值运算，|x(l)|为数据序列的幅度谱；

在第二步中，x(n)的离散傅里叶变换即式(1)，是通过快速傅里叶变换实现的，利用快速傅里叶变换可以降低算法的运算量，提高算法的计算效率；式(2)中k＝0,1,2…,n/2是因为实数据序列的离散傅里叶变换关于中心对称，因此k可以只取前n/2+1个点，可以降低一半的计算量。

第三步：搜索功率谱p(k)最大值所对应的离散频率索引k0，并计算x(l)位于离散频率索引k0-1，k0和k0+1处的模值al，am和ar，计算过程如下：

搜索功率谱p(k)最大值所对应的索引k0

其中表示在1≤k≤n/2-1范围内搜索p(k)的最大值所对应的离散频率索引；

取x(l)在离散频率索引k0-1，k0和k0+1处的模值，分别记为al，am和ar，即

al＝|x(k0-1)|,am＝|x(k0)|,ar＝|x(k0+1)|式(4)

在第三步中，分两步实现的：第一步搜索功率谱p(k)最大值所对应的离散频率索引k0；第二步依据搜索的最大值所对应的离散频率索引k0，计算x(l)位于离散频率索引k0-1，k0和k0+1处的模值al，am和ar。

第四步：预估有向频率相对偏差δr和加权频率相对偏差δw，计算过程如下：

利用所述模值al，am和ar计算有向频率相对偏差δr和加权频率相对偏差δw，即

在第三步中，利用式(5)计算有向频率相对偏差δr，例如《电子学报》2004年，32(4)的第625页至第628页中公开的rife插值方法，rife插值方法的频率相对偏差为式(5)中前两行所示的表达式。本文发明的方法对rife插值方法的频率相对偏差做了修正，即当al＝ar时将δr置0，这是因为当al＝ar时说明，信号频率相对偏差为0，无需再依据次大线谱去计算相对频率偏差。

第五步：利用al，am，ar，δr和δw计算综合频率相对偏差δc，计算过程如下：

首先，计算初始综合相对偏差δm

其中tl为相对频率偏差阈值下限，th为相对频率偏差阈值上限，要求tl和th同时满足：

然后，对δm进行修正得到综合相对偏差δc

在第五步中，分两步实现的：第一步利用式(8)计算初始综合相对偏差δm；第二步利用式(9)对δm进行修正得到综合相对偏差δc。

在第五步中，相对频率偏差阈值下限tl和上限th只要同时满足式(8)所示的三个条件即可，但是通过蒙特卡洛仿真实验得到的优选值为tl＝0.07，th＝0.14。

第六步：将δc带入插值公式估计出单频信号频率即

其中δf为长度为n的离散傅里叶变换的频率分辨率，δf＝fs/n，fs为采样频率；

本发明的一种基于频率相对偏差预估的分段综合单频信号频率估计方法：首先依据预估频率相对偏差的分布范围和设定的门限，自动选择综合相对偏差的取值，从而进一步利用插值公式估计出单频信号的频率参数。

本发明的实施例中，仿真接收单频信号模型为：

其中a为信号幅度，为初始相位，n为信号采样点数，f0为信号频率，fs为采样频率，w(n)为均值为0，方差为σ²的高斯白噪声，方差σ²的大小由信噪比snr决定：snr＝10log10[a²/(2σ²)]。

实施例1：

仿真信号参数分别设置为：信号幅度a＝1，初始相位采样点数n＝1024，采样频率fs＝4096hz，则离散傅里叶变换的频率分辨率δf＝fs/n＝4hz，信号频率f0＝301hz，信噪比snr＝0db。设定相对频率偏差阈值下限tl＝0.07，相对频率偏差阈值上限th＝0.14。

图2所示是频率为301hz的仿真接收单频信号的功率谱p(k)；图3所示是频率为301hz的仿真接收单频信号的幅度谱|x(l)|。

下面对x(n)进行频率估计。

首先谱搜索功率谱p(k)最大值所对应的索引k0，由图2可知k0＝75；

然后取x(l)在离散频率索引k0-1，k0和k0+1处的模值al，am和ar，由图3可知al，am和ar分别为

al＝|x(k0-1)|＝0.2008，am＝|x(k0)|＝0.8844，ar＝|x(k0+1)|＝0.2836

由此可知al<ar，因此依据式(5)可以计算出有向频率相对偏差δr

依据式(6)可以计算出加权频率相对偏差δw

由此可知|δr|>th，因此依据第五中的式(7)计算得到初始综合相对偏差δm

δm＝δr＝0.2428

由此可知|δm|<0.5，进一步依据式(9)对δm进行修正得到综合相对偏差δc

δc＝δm＝0.2428

最后依据第六步的式(10)插值估计出正弦波脉冲信号频率

估计频率相对误差为

实施例2：

仿真信号参数分别设置为：信号幅度a＝2，初始相位采样点数n＝1024，采样频率fs＝4096hz，则离散傅里叶变换的频率分辨率δf＝fs/n＝4hz，信号频率f0＝500hz，信噪比snr＝3db。设定相对频率偏差阈值下限tl＝0.07，相对频率偏差阈值上限th＝0.14。

图4所示是频率为500hz的仿真接收单频信号的功率谱p(k)；图5所示是频率为500hz的仿真接收单频信号的幅度谱|x(l)|。

下面对x(n)进行频率估计。

首先谱搜索功率谱p(k)最大值所对应的索引k0，由图4可知k0＝125；

然后取x(l)在离散频率索引k0-1，k0和k0+1处的模值al，am和ar，由图5可知al，am和ar分别为

al＝|x(k0-1)|＝0.01663，am＝|x(k0)|＝2.001，ar＝|x(k0+1)|＝0.04218

由此可知al<ar，因此依据式(5)可以计算出有向频率相对偏差δr

依据式(6)可以计算出加权频率相对偏差δw

由此可知δw<tl且|δr|<th，因此依据第五中的式(7)计算得到初始综合相对偏差δm

由此可知|δm|<0.5，进一步依据式(9)对δm进行修正得到综合相对偏差δc

δc＝δm＝0.000094

最后依据第六步的式(10)插值估计出正弦波脉冲信号频率

估计频率相对误差为

实施例3：

仿真信号参数分别设置为：信号幅度a＝2，初始相位采样点数n＝1024，采样频率fs＝4096hz，则离散傅里叶变换的频率分辨率δf＝fs/n＝4hz，信号频率f0＝799.6hz，信噪比snr＝6db。设定相对频率偏差阈值下限tl＝0.07，相对频率偏差阈值上限th＝0.14。

图6所示是频率为799.6hz的仿真接收单频信号的功率谱p(k)；图7所示是频率为799.6hz的仿真接收单频信号的幅度谱|x(l)|。

下面对x(n)进行频率估计。

首先谱搜索功率谱p(k)最大值所对应的索引k0，由图6可知k0＝200；

然后取x(l)在离散频率索引k0-1，k0和k0+1处的模值al，am和ar，由图7可知al，am和ar分别为

al＝|x(k0-1)|＝0.2355，am＝|x(k0)|＝1.9465，ar＝|x(k0+1)|＝0.1540

由此可知al>ar，因此依据式(5)可以计算出有向频率相对偏差δr

依据式(6)可以计算出加权频率相对偏差δw

由此可知tl<δw<th且|δr|<th，因此依据第五中的式(7)计算得到初始综合相对偏差δm

由此可知|δm|<0.5，进一步依据式(9)对δm进行修正得到综合相对偏差δc

δc＝δm＝-0.0238

最后依据第六步的式(10)插值估计出正弦波脉冲信号频率

估计频率相对误差为

尽管本发明就优选实施方式进行了示意和描述，但本领域的技术人员应当理解，只要不超出本发明的权利要求所限定的范围，可以对本发明进行各种变化和修改。