一种基于并行FFT的快速精频捕获方法与流程

本发明属于卫星信号处理领域，具体涉及一种基于并行fft的快速精频捕获方法。

背景技术：

卫星导航定位技术提供高精度的定位、导航和授时服务，广泛应用于大众消费、智慧城市、精密机械控制、高空探测、电力、金融计时、交通、公安、减灾救灾、农业和渔业等工作和生活的方方面面。卫星接收机属于用户端，主要由射频前端、捕获、跟踪、导航信息解调、伪距和位置解算五部分组成。对于资源有限的多通道接收机，接收机的每个通道包含一个捕获模块和一个跟踪模块。通过降低各通道采集模块的计算量，可以为接收机节省大量的计算空间。

卫星信号捕获的目的是对观测卫星进行检测，估计出多普勒频移和测距码延迟。对于锁相环(pll，phaselockingloop)来说，来自粗捕获的频率精度可能太粗糙，且捕获得到的多普勒频移越接近实际值，pll达到稳定所需的时间越短。为了使pll模块更快地进入跟踪状态，必须保证由精细捕获模块或锁频环(fll，frequencylockedloop)模块估计的多普勒频移误差在几十赫兹以内。

现有精频捕获方法大致分为三类，第一类是通过补零提高快速傅里叶变换(fft，fastfouriertransformalgorithm)的分辨率进而提升载波多普勒频移的估计精度。然而，要获取10hz的估计精度fft需要处理100ms的数据，假设采样频率为10mhz，则意味着fft过程需要处理1,000,000个数据，这将耗费巨大的计算资源。第二类方法是通过频率牵引实现精频捕获，如fll，这种方法首先通过较短数据实现多普勒频移的粗略估计，然后将多普勒频移的估计精度牵引到pll的带宽范围内。然而牵引过程需要不断地迭代，使多普勒频移的估计精度满足pll进入稳定跟踪状态要求。迭代过程既消耗了计算资源又增加了捕获时间。第三类方法是两步法，通常称为粗-精捕获。第一步使用较短数据长度实现粗捕，第二步，通过缩短搜索步长、载波相位或fft输出可以用sinc函数表示特性等实现载波多普勒频移的精确捕获。然而，如果采用较短的搜索步长依然会消耗大量的计算资源，在信号很弱的时候无法采用载波相位或者fft输出可以用sinc函数表示的特性。

技术实现要素：

发明目的：为了解决上述背景技术中提出的技术问题，本发明旨在提供一种基于并行fft的快速精频捕获方法，利用fft输出幅度谱的特性，设计一种并行结构，实现对载波多普勒频移快速精确估计。

技术方案：本发明所述的一种基于并行fft的快速精频捕获方法，包括以下步骤：

(1)根据测距码的自相关特性，通过将中频信号与本地测距码序列相乘，得到剥离测距码序列后的中频信号s^j(n)；

(2)获取四个不同频率的载波信号，并分别与剥离测距码序列后的中频信号s^j(n)进行混频，得到四个低频信号

(3)对四个低频信号分别进行降采样处理，得到降采样后的四个信号xcn(k)；

(4)将降采样后的四个低频信号分别进行fft处理，得到四个幅度谱序列

(5)分别提取四个序列幅度谱的最大值，记为mfft1,mfft2,mfft3和mfft4；引入参数mdi，明确每个fft输出幅度谱的最大值之间的数值关系，并对这四个数进行排序，将排序后的mfft1,mfft2,和mfft4进行二次曲线拟合。

进一步地，所述步骤(1)通过以下公式实现：

其中，i是复数符号；fif是中频频率，为系统已知常数；表示第j个信号发射端相对接收机的载波多普勒频移；ni和nq是一对正交的等功率的高斯白噪声；r^j表示信号测距码与本地测距码间的相关值；kn表示搜索步长编号；τ表示信号测距码与本地测距码间的相位差；表示调制在第j颗卫星b1i信号上的导航电文序列；ab1i是b1i信号的幅值；表示第j颗卫星的初始载波相位。

进一步地，所述步骤(2)的实现过程如下：

采用复变频技术对码剥离后的信号s^j(n)进行下变频：

其中，下标cn对应四个不同频率的载波信号，分别表示为1,2,3,4，fcn表示cn相关器的中心频率。

进一步地，所述在步骤(3)实现过程如下：

通过积分过程降低信号的采样频率：

其中，ts表示初始采样周期；m表示累加点数；xcn(k)为通过降采样获得的四个信号，该类信号的采样周期变为m·ts。

进一步地，所述步骤(4)实现过程如下：

公式(5)经过矩形窗截断和dtft处理后，可以用公式(6)表示：

有限长度的离散频域信号可以用公式(7)表示：

其中，*表示卷积运算符号，δ(·)表示狄拉克函数，sinc(·)表示sinc函数，xcn(f)表示经过dtft过程后获取的四组连续幅度谱信号，rect(·)为矩形窗函数，表示有限长的离散幅度谱序列，tl表示矩形窗长度，kf表示离散频率量，其相邻两值间的间隔为1/tl，表示对连续幅度谱xcn(f)截断和离散后的信号。

进一步地，所述步骤(5)的实现过程如下：

四个fft输出幅度谱最大值间的数值关系用公式(8)表示：

其中，mfft1,mfft2,mfft3和mfft4分别表示fft1,fft2,fft3和fft4输出幅度谱的最大值，l(d)表示集合d中元素的个数，l(a)表示集合中元素a的个数，a表示集合u中的最小值，b表示集合cua中的最小值，cu表示余集，d表示集合u中的最大值，c表示集合cud中的最小值；

根据a和b的数值关系，四个fft输出幅度谱最大值可以分为2类，共8种情况：第一类为a远大于b，包括：情况1，a＝mfft1；情况3，a＝mfft4；情况5，a＝mfft3；情况7，a＝mfft2；第二类为a与b近似相等，包括：情况2，(a＝mfft1&b＝mfft4)|(a＝mfft4&b＝mfft1)；情况4，(a＝mfft3&b＝mfft4)|(a＝mfft4&b＝mfft3)；情况6，(a＝mfft2&b＝mfft3)|(a＝mfft3&b＝mfft2)；情况8，(a＝mfft1&b＝mfft2)|(a＝mfft2&b＝mfft1)；

设计参数mdi来明确每个fft输出的最大值之间的数值关系，如等式(9)所示：

其中，mdi表示模型定义指标，为本专利设计的一个参数；

二次曲线拟合方法如公式(10)所示：

其中，p是二次多项式的系数矩阵；y是由fft1,fft2,fft3和fft4输出幅度谱序列的最大值组成的矩阵，根据最小二乘原理，矩阵用式(11)求解：

p＝a⁺y(11)

其中，a⁺是a矩阵的广义逆，为常数矩阵；y矩阵如式(12)：

由式(11)推导二次曲线最大点的位置，如式(13)所示：

其中，p(2)＝v，表示二次多项式的一次项系数；p(1)＝u，表示二次多项式的二次项系数；fco表示二次多项式极值点位置；

二次曲线的最大点位置可以作为fft估计的多普勒频移的补偿，用fco表示，获得多普勒频移用式(14)表示：

其中，ffft1是fft1估计的多普勒频移，ffft4是fft4估计的多普勒频移，fd为多普勒频移的估计值。

有益效果：与现有技术相比，本发明的有益效果：1、本发明基于fft输出幅度谱可以用二次函数拟合的特性，设计了四个并行的结构，通过二次拟合后，可以获得更可靠更高精度的多普勒频移估计，同时由于本发明采用基于并行结构，而不是采用迭代结构，使得捕获过程更加迅速；2、本发明适合卫星信号的捕获，同时也适用于其他采用码分多址(cdma，codedivisionmultipleaccess，)技术调制的信号的捕获。

附图说明

图1是本发明的流程图；

图2是四个fft输出幅度谱最大值随多普勒频移变化的曲线图；

图3是四个fft输出幅度谱最大值间的8种数值关系图；

图4是四个fft输出幅度谱最大值间的数值关系判别图。

具体实施方式

下面结合附图对本发明作进一步详细说明。

本发明提出一种基于并行fft的快速精频捕获方法，如图1所示，具体包括以下步骤：

步骤1：根据测距码的自相关特性，通过将中频信号与本地测距码序列相乘，得到剥离测距码序列后的中频信号s^j(n)。

通过射频前端的下变频和离散化之后，到达接收机捕获模块的信号可用公式(1)表示。

上式中，表示接受到的中频卫星信号；j表示卫星号；ab1i是b1i信号的幅值；表示第j颗卫星的测距码序列；表示调制在第j颗卫星b1i信号上的导航电文序列；i是复数符号；fif是中频频率，这是系统已知常数；表示第j个信号发射端相对接收机的载波多普勒频移；表示第j颗卫星的初始载波相位；ni和nq是一对正交的等功率的高斯白噪声。

保持1ms的本地产生的测距码不动，逐步移动接收到的2ms中频信号相关过程在交叠区域完成，该过程可以用公式(2)表示。

上式中，s^j(n)表示载波剥离后的中频信号；表示本地产生的第j颗卫星的测距码序列；r^j(τ)表示第j颗卫星b1i信号测距码的自相关函数，kn表示搜索步长编号，τ表示信号测距码与本地测距码间的相位差。当τ＝0时，在一个测距码周期内，公式(2)可以用公式(3)表示。

步骤2：设计四个不同频率的载波信号，并分别与剥离测距码序列后的中频信号s^j(n)进行混频，得到四个低频信号

采用复变频技术对码剥离后的信号s^j(n)进行下变频，其过程如公式(4)所示：

上式中，下标cn表示相关器的序列号，可以用1,2,3,4表示；fcn表示cn相关器的中心频率，本地产生的4组混频序列是固定的可以预先存储，当执行下变频操作时调用；表示通过下变频获得的4组低频信号。

步骤3：对这四个低频信号分别进行降采样处理，得到降采样后的四个信号xcn(k)，减小该算法的计算量。

通过积分过程降低信号的采样频率，由于信号的采样频率fs远大于载波多普勒频移，降采样过程可以用公式(5)表示：

上式中，ts表示初始采样周期；m表示累加点数；xcn(k)为通过降采样获得的四个信号，该类信号的采样周期变为m·ts。

步骤4：将降采样后的四个低频信号分别进行fft处理，得到四个幅度谱序列

公式(5)经过矩形窗截断和离散时间傅里叶变换(dtft，discrete-timefouriertransform)处理后，可以用公式(6)表示：

有限长度的离散频域信号可以用公式(7)表示。

上式中，*表示卷积运算符号，δ(·)表示狄拉克函数，sinc(·)表示sinc函数，xcn(f)表示经过dtft过程后获取的四组连续幅度谱信号，rect(·)为矩形窗函数，表示有限长的离散幅度谱序列，tl表示矩形窗长度，kf表示离散频率量，其相邻两值间的间隔为1/tl，表示对连续幅度谱xcn(f)截断和离散后的信号。

步骤5：分别提取四个序列幅度谱的最大值，为：mfft1,mfft2,mfft3和mfft4。根据本专利设计的模型定义指标(mdi,modeldefinitionindicators)对这四个数进行排序，将排序后的mfft1,mfft2,mfft3和mfft4进行二次曲线拟合。

在该步骤中，中心频率fcn被分配为：f1＝-500hz,f2＝-250hz,f3＝0hz和f4＝250hz。理论中，fft3估计得多普勒频移表示真实的载波多普勒频移，但是fft3的频率分辨率仍然是1khz。为了使得多普勒频移的估计精度更加精确，对四个fft输出幅度谱的最大值进行二次曲线拟合。为了使得这四个fft输出幅度谱的最大值坐落在同一条二次曲线上，应该得知这四个fft输出幅度谱最大值间的数值关系。所以，首先从理论上研究这四个fft输出幅度谱最大值的数值关系；然后，研究噪声对理论关系的影响；最后，获取有噪声情况下这四个fft输出幅度谱最大值间的数值关系。为了描述简单，这四个fft输出幅度谱最大值间的数值关系用公式(8)表示：

上式中，mfft1,mfft2,mfft3和mfft4分别表示fft1,fft2,fft3和fft4输出幅度谱的最大值。l(d)表示集合d中元素的个数，l(a)表示集合中元素a的个数，a表示集合u中的最小值，b表示集合cua中的最小值，cu表示余集，d表示集合u中的最大值，c表示集合cud中的最小值。

根据公式(6)和公式(7)可以得出无噪声情况下的四个fft输出幅度谱与真实多普勒频移之间的关系图，如图3所示。

如图2所示，这四个fft输出幅度谱最大值间的数值关系随着多普勒频移变化而改变，改变周期为1khz。根据a和b的数值关系，这四个fft输出幅度谱的最大值可以分为2类，共8种情况，如图3所示。第一类为a远大于b，包括：情况1，a＝mfft1；情况3，a＝mfft4；情况5，a＝mfft3；情况7，a＝mfft2。第二类为a与b近似相等，包括：情况2，(a＝mfft1&b＝mfft4)|(a＝mfft4&b＝mfft1)；情况4，(a＝mfft3&b＝mfft4)|(a＝mfft4&b＝mfft3)；情况6，(a＝mfft2&b＝mfft3)|(a＝mfft3&b＝mfft2)；情况8，(a＝mfft1&b＝mfft2)|(a＝mfft2&b＝mfft1)。

由于噪声的参与，很难获得与图3类似的清晰分类方法。因此，设计一个参数mdi来明确每个fft输出幅度谱的最大值之间的数值关系，如等式(9)所示。

通过实验，mdi的门限值可设置为0.1227。

通过mdi将四个fft输出幅度谱最大值间的数值关系分为两大类，如图4所示。

二次曲线拟合方法如公式(10)所示。

上式中，p是二次多项式的系数矩阵，u表示二次多项式的二次项系数，v表示二次多项式的一次项系数，r表示二次多项式的常数项系数；y是由fft1,fft2,fft3和fft4输出幅度谱序列的最大值组成的矩阵，a是一个与fcn有关的常数矩阵。根据最小二乘原理，矩阵p可用式(11)求解。

p＝a⁺y(11)

上式中，a⁺是a矩阵的广义逆，为常数矩阵。接下来，根据图4中的8种情况，可以得到y矩阵，并在式(12)中示出。

因此，二次曲线最大点的位置可由式(11)推导，如式(13)所示。

上式中，p(2)＝v，表示二次多项式的一次项系数；p(1)＝u，表示二次多项式的二次项系数；fco表示二次多项式极值点位置。

二次曲线的最大点位置可以作为fft估计的多普勒频移的补偿，用fco表示。

由该算法获得的多普勒频移可用式(14)表示。

上式中，ffft1是fft1估计的多普勒频移，ffft4是fft4估计的多普勒频移，ffft1和ffft4的频率分辨率均为1khz，fd为多普勒频移的估计值。

实施例仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明保护范围之内。