一种低样本数信号中心频率估计方法及系统

1.本发明涉及信号中心频率估计，特别是涉及一种低样本数信号中心频率估计方法及系统。


背景技术：

2.信号中心频率估计在随机信号处理中起到重要作用，尤其在语音信号处理、雷达信号处理、光学图像处理等信号检测和目标识别等领域，信号的中心频率的估计精度对信号特征提取目标属性分析起到至关重要的重要。
3.常规的信号频谱计算和信号中心频率估计都是基于计算机等大型信号处理设备，判断是根据信号频谱或者功率谱曲线峰值对应的频率点作为中心频率，这些设备能处理的信号样本数大，而且能够承受大型fft(快速傅里叶变换)计算，但是在许多物联网前端设备中的信号处理单元使用的是诸如dsp(数字信号处理器)系统等最小系统，计算能力相对较差，为保证信号处理效率及系统稳定性，能够实时fft计算的数据长度有限，频率采样点数影响了信号频谱的频率分辨率，导致信号中心频率偏移，因此，频谱最高峰对应的频率不一定是信息的真实中心频率。


技术实现要素：

4.本发明的目的在于克服现有技术的不足，提供一种低样本数信号中心频率估计方法及系统，能够获得较为准确的中心频率。
5.本发明的目的是通过以下技术方案来实现的：一种低样本数信号中心频率估计方法，包括以下步骤：
6.s1.对输入信号进行采样和自相关计算；
7.s101.设置采样频率为fs，分别对输入信号i/q通道进行采样，并将两通道数据保存在二维数组data[n][2]中，其中n为每个通道的数据采样点数，若i通道数据的保存在二维数组第一列，即data[i][1]，则q通道数据保存在二维数组第二列，即data[i][2]，i＝1,
…
,n；
[0008]
s102.设置序号k，在k＝1,
…
,n时，分别按照如下公式计算当前序号下自相关值；
[0009][0010][0011]
s103.当k＝1,
…
,n中所有取值下的自相关值均计算完毕后，得到自相关二维数组，记为：
[0012]
r[1][1],r[2][1],...,r[n][1]；
[0013]
r[1][2],r[2][2],...,r[n][2]。
[0014]
s2.根据自相关计算结果进行信号功率谱估计；
[0015]
s201.根据自相关二维数组构建复数数组r[k],将其作为自相关序列:
[0016]
r[k]＝r[k][1]+jr[k][2]，k＝1,
…
,n；
[0017]
s202.设置fft变换长度n
fft
，对自相关序列r[k]进行点fft变换计算，得到复数数组 p[k]＝p
real
[k]+jp
image
[k]，其中k＝1,
…
,n
fft
；
[0018]
s203.对于复数数组p[k],k＝1,
…
,n
fft
,计算其幅度：
[0019][0020]
s3.对估计得到的功率谱进行平滑和处理，得到频谱序列；
[0021]
s301.截取功率谱前一半样点，即p
′
[k]＝p[k]，k＝1,
…
,n
fft
/2；
[0022]
s302.设置滑动时窗长度m，
[0023]
在k＝1,
…
,n
fft
/2时，分别按照如下方式计算得到pm[k]：
[0024][0025]
pm[k]＝p
′
[k]，k＝n
fft
/2-m,
…
,n
fft
/2；
[0026]
s303.在k＝1,/,n
fft
/2取值计算完成后，得到频谱序列pm[k]：k＝1,
…
,n
fft
/2。
[0027]
s4.根据频谱序列，实现功率谱主峰采样点范围的计算；
[0028]
s401.获取频谱序列pm[k],k＝1,
…
,n
fft
/2最大值p
max
及其对应的序号k
max
；
[0029]
s402.在序号k＝k
max
,
…
,1时：
[0030]
a1、初始化k＝k
max
；
[0031]
a2、判断是否满足pm[k]≤0.1
×
p
max
，若满足，则k
l
＝k，并执行步骤s403，若不满足，则更新k＝k-1，并继续执行步骤a2；其中，k
l
为功率谱主峰采样点范围的下限值；
[0032]
s403.在序号k＝k
max
,
…
,n
fft
/2时：
[0033]
b1、初始化k＝k
max
；
[0034]
b2、判断是否满足pm[k]≤0.1
×
p
max
，若满足则kr＝k，并执行步骤s404，若不满足，则更新k＝k+1，并继续执行步骤b2；其中，kr为功率谱主峰采样点范围的上限值；
[0035]
s404.得到功率谱主峰采样点范围为[k
l
,kr]。
[0036]
s5.计算信号中心频率：
[0037]
s501.在序号k＝k
l
,
…
,k
max
时；
[0038]
c1、初始化k＝k
l
；
[0039]
c2、判断是否满足pm[k]≥0.707
×
p
max
：
[0040]
若满足，则k
lp
＝k，并执行步骤s502；
[0041]
否则更新k＝k+1，并继续执行步骤c2；
[0042]
s502.在序号k＝k
max
,
…
,kr时；
[0043]
d1、初始化k＝k
max
；
[0044]
d2、判断是否满足pm[k]≥0.707
×
p
max
；
[0045]
若满足，则k
rp
＝k，执行步骤s503；
[0046]
否则更新k＝k+1，继续执行步骤d2；
[0047]
s503.计算信号中心频率f0[0048][0049]
一种低样本数信号中心频率估计系统，包括：
[0050]
自相关计算模块，用于对输入信号进行采样和自相关计算；
[0051]
功率谱估计模块，用于根据自相关计算结果进行信号功率谱估计；
[0052]
频谱序列计算模块，用于对估计得到的功率谱进行平滑和处理，得到频谱序列；
[0053]
主波峰计算模块，用于根据频谱序列，实现功率谱的主波峰计算；
[0054]
中心频率计算模块，用于计算信号中心频率。
[0055]
本发明的有益效果是：本发明以功率谱最高峰为参照，先找到功率谱的主峰范围，排除其他频段信号对有用信号的干扰；然后在主峰范围内确定有用信号能量最集中范围，降低截断效应对信号功率谱的影响；最后在信号能量最集中的频率分为内选择频率中间值作为中心频率，降低频率偏移对信号的影响，能实现在低于采样率数据fft变换下实现信号的中心频率估计；并能满足dsp、单片机等小型系统快速信号处理的需求。
附图说明
[0056]
图1为本发明的方法流程图；
[0057]
图2为本发明的系统原理框图。
具体实施方式
[0058]
下面结合附图进一步详细描述本发明的技术方案，但本发明的保护范围不局限于以下所述。
[0059]
如图1所示，一种低样本数信号中心频率估计方法，包括以下步骤：
[0060]
s1.对输入信号进行采样和自相关计算；在本技术的实施例中，所述步骤s1包括以下子步骤：
[0061]
s101.设置系统数据采样频率为fs，分别对输入信号i/q通道进行采样，并将两通道数据保存在二维数组data[n][2]中，其中n为每个通道的数据采样点数，若i通道数据的保存在二维数组第一列，即data[i][1]，则q通道数据保存在二维数组第二列，即data[i][2]， i＝1,
…
,n；
[0062]
s102.设置序号k，在k＝1,
…
,n时，分别按照如下公式计算当前序号下自相关值；
[0063][0064][0065]
s103.当k＝1,
…
,n中所有取值下的自相关值均计算完毕后，得到自相关二维数组，记为：
[0066]
r[1][1],r[2][1],...,r[n][1]；
[0067]
r[1][2],r[2][2],...,r[n][2]。
[0068]
s2.根据自相关计算结果进行信号功率谱估计；在本技术的实施例中，所述步骤s2包括以下子步骤：
[0069]
s201.根据自相关二维数组构建复数数组r[k],将其作为自相关序列:
[0070]
r[k]＝r[k][1]+jr[k][2]，k＝1,
…
,n；
[0071]
s202.设置fft变换长度n
fft
，对自相关序列r[k]进行点fft变换计算，得到复数数组 p[k]＝p
real
[k]+jp
image
[k]，其中k＝1,
…
,n
fft
；
[0072]
s203.对于复数数组p[k],k＝1,
…
,n
fft
,计算其幅度：
[0073][0074]
s3.对估计得到的功率谱进行平滑和处理，得到频谱序列；在本技术的实施例中，所述步骤s3包括以下子步骤：
[0075]
s301.截取功率谱前一半样点，即p
′
[k]＝p[k]，k＝1,
…
,n
fft
/2；
[0076]
s302.设置滑动时窗长度m，
[0077]
在k＝1,
…
,n
fft
/2时，分别按照如下方式计算得到pm[k]：
[0078][0079]
pm[k]＝p
′
[k]，k＝n
fft
/2-m,
…
,n
fft
/2；
[0080]
s303.在k＝1,
…
,n
fft
/2取值计算完成后，得到频谱序列pm[k]：k＝1,
…
,n
fft
/2。
[0081]
s4.根据频谱序列，实现功率谱主峰采样点范围的计算；在本技术的实施例中，所述步骤s4包括以下子步骤：
[0082]
s401.获取频谱序列pm[k],k＝1,
…
,n
fft
/2最大值p
max
及其对应的序号k
max
；
[0083]
s402.在序号k＝k
max
,
…
,1时：
[0084]
a1、初始化k＝k
max
；
[0085]
a2、判断是否满足pm[k]≤0.1
×
p
max
，若满足，则k
l
＝k，并执行步骤s403，若不满足，则更新k＝k-1，并继续执行步骤a2；其中，k
l
为功率谱主峰采样点范围的下限值；
[0086]
s403.在序号k＝k
max
,
…
,n
fft
/2时：
[0087]
b1、初始化k＝k
max
；
[0088]
b2、判断是否满足pm[k]≤0.1
×
p
max
，若满足则kr＝k，并执行步骤s404，若不满足，则更新k＝k+1，并继续执行步骤b2；其中，kr为功率谱主峰采样点范围的上限值；
[0089]
s404.得到功率谱主峰采样点范围为[k
l
,kr]。
[0090]
步骤s4中，先找到功率谱的幅度最大值p
max
及其对应的序号k
max
，然后以幅度最大值 p
max
为参照，找到离幅度最大值p
max
最近的能量最集中的一段频率范围(用序号表示 k
l
‑‑
kr)，序号k
l
‑‑
kr之间的部分当成主峰采样范围，k
l
和kr分别表示主峰最左边和最右边的序号；因为在低采样信号中，信号的功率谱能量会分散，功率谱最大值对应的频率不一定是信号的真实中心频率。
[0091]
以功率谱最大值p
max
对应的位置k
max
为中心向两边搜索功率谱边界[k
l
,kr]，k
l
为k
max
左侧功率谱第一次降低到0.1p
max
对应的样点位置，kr为k
max
右侧功率谱第一次降低到
0.1p
max
对应的样点位置。
[0092]
s5.计算信号中心频率，在本技术的实施例中，所述步骤s5包括以下子步骤：
[0093]
由于[k
l
,kr]范围太大，因此，需要分别从k
l
和kr向中间搜寻功率谱第一次升到 0.707p
max
对应的样点位置，对应位置分别是k
lp
和k
rp
；这样做的目的是，低采样和信号成分复杂时，信号的功率谱可能是带有一定带宽的，而不是理论的单脉冲，具体地：
[0094]
s501.在序号k＝k
l
,
…
,k
max
时；
[0095]
c1、初始化k＝k
l
；
[0096]
c2、判断是否满足pm[k]≥0.707
×
p
max
：
[0097]
若满足，则k
lp
＝k，并执行步骤s502；
[0098]
否则更新k＝k+1，并继续执行步骤c2；
[0099]
s502.在序号k＝k
max
,
…
,kr时；
[0100]
d1、初始化k＝k
max
；
[0101]
d2、判断是否满足pm[k]≥0.707
×
p
max
；
[0102]
若满足，则k
rp
＝k，执行步骤s503；
[0103]
否则更新k＝k+1，继续执行步骤d2；
[0104]
在本技术的实施例中，使用步骤s401-s404先确定大范围，防止过早陷入虚假带宽边界；然后使用步骤s501-s502确定功率谱能量最集中的部分，这段能量最集中的部分可以看作的功率谱主峰。
[0105]
s503.计算信号中心频率f0[0106][0107]
如图2所示，一种低样本数信号中心频率估计系统，包括：
[0108]
自相关计算模块，用于对输入信号进行采样和自相关计算；
[0109]
功率谱估计模块，用于根据自相关计算结果进行信号功率谱估计；
[0110]
频谱序列计算模块，用于对估计得到的功率谱进行平滑和处理，得到频谱序列；
[0111]
主波峰计算模块，用于根据频谱序列，实现功率谱的主波峰计算；
[0112]
中心频率计算模块，用于计算信号中心频率。
[0113]
以上所述是本发明的优选实施方式，应当理解本发明并非局限于本文所披露的形式，不应该看作是对其他实施例的排除，而可用于其他组合、修改和环境，并能够在本文所述构想范围内，通过上述教导或相关领域的技术或知识进行改动。而本领域人员所进行的改动和变化不脱离本发明的精神和范围，则都应在本发明所附权利要求的保护范围内。