优化的PFA算法在GPS软件接收机中的应用

姚相振，崔绍龙，方金云

（1.中国科学院 计算技术研究所，北京 100190；2.中国科学院 研究生院，北京 100049）

1 引言

卫星信号捕获是GPS接收机里的第一个必要操作，它用来检测来自天线接收的信号中包含有哪些导航卫星的信息，捕获的最终目的是为了得到载波的多普勒频移以及每颗卫星C/A码的初始相位。在此过程中执行C/A码相位和多普勒频移的二维搜索，具体的操作是用本地生成的码跟接收到的信号做交叉相关，然后把相关的结果累加起来，设定一个门限，超过此门限则判定相应的卫星被捕获。

GPS信号采用扩频调制技术，导航数据首先被C/A码或者P码进行扩频调制，然后对1575.42MHz的L1载波进行BPSK调制。同相分量上含有C/A码，正交分量上含有P码。因此，利用Gold码的相关特性进行相关调制是 GPS接收机中最基本的运算。相关器的性能和C/A码的捕获成功率，往往决定了接收机的性能。在噪声环境，尤其是弱信号的时候，利用相干累计和非相干累计技术可以提高接收机的增益。但是无论是相干累计还是非相干累计都包含有大量的相关运算，这对于软件接收机来说无疑是一个沉重的负担[1]，Van Nee在文献[2]中提出了用FFT运算来代替相关运算，此方法大大提高了运算的效率，也奠定了软件无线电设备的基础，但是在软件导航接收机中要实现实时信号处理，传统的FFT算法还是无法胜任，另外，FFT中存在的大量浮点运算也使得在嵌入式导航接收机上更是难以实现实时性，因为市场上现有的嵌入式设备对于浮点运算的支持还不是很完善。

为了减少捕获阶段的运算复杂度以及处理时间，大量的尝试和实验已经被实施，在文献[3]中提出一种平均相关算法，此算法对采样值求平均，以降低参与 FFT/IFFT的运算点数。在文献[4,5]中，作者提出一种自我调整捕获方案，它能够根据信噪比的大小自动的调节捕获门限的大小，文献[6]中提到了一种利用特殊软件相关器的方法，使用并行位操作来实现 32个点同时工作来模拟硬件相关器，结果显示其效率有很大程度的提高。

本文中提出一种优化 PFA算法，首先，提出一种基于位操作的基带混合算法，从RF前端得到的数字信号被存储为一种特殊的2bit形式，本地载波体也被量化为2bit形式，为2组数据创建一个数据表，通过索引来实现基带混合的过程，避免了乘法运算，在FFT运算过程中，将旋转因子预存为16bit的补码形式，并将其映射为一个较大的整数，中间数据以及最终结果都用一个32bit的整数来存储，在每一次蝶形运算后进行溢出检测，以确保数据的正确性，另外 PFA算法过程中不必要的乘法也将会被省略，因此在整个运算过程中只存在整数运算，此方法在浮点运算支持不利嵌入式设备上的运行取得良好的运行效果。

2 信号捕获

2.1 循环相关到FFT的演化

所有的GPS卫星都发射1.57542GHz(L1波段)和1.2276GHz（L2波段）2种频率的信号，GPS的信号发射采用的是CDMA模式[7]，低频的导航数据被高频的伪随机码进行调制，每颗卫星的信号通过调制不同的伪随机码来实现身份识别，GPS信号采用了 2种不同的伪随机码:一种是 C/A码（Gold码的一种），码片率为 1.023Mchip/s, 另一种是 P码，码片率为 10.23Mchip/s。P码被调制在 L1和L2载波上，而C/A码仅调制在L1载波上且与P码相位相差90°。C/A码被广泛用于民用接收机应用，而P码被加密为一个被称为P(Y)码的编码，只用于军事装备。在此本文重点讨论C/A码的应用模式。

从RF前端接收到的信号可以表示为

其中，Ns代表可见卫星的个数，Ai代表信号的幅度，Di代表导航数据，Ci代表第i颗卫星的C/A码，tτi代表码时延，fdi代表多普勒频移，fL1代表 L1载波的频率，其值为1575.42MHz，fIF是RF下变频后的中频频率，φi代表卫星i的初始相位，n代表一个均值为零方差为的高斯白噪声。用于计算多普勒频移造成的码元长度变化。

C/A码一个周期长度为1023个码元，码率为1.023MHz，一个周期的C/A码可以用式(3)表示为

p(t)表示以原点为中心，具有单位高度和单位宽度的脉冲波形。xn为C/A码的第n个码元。可以用+1 和−1 来表示，Tc=1/(1.023MHz)= (1/1.023)μs，是C/A码的码元宽度，C/A码的周期Tcode= 1023Tc=1ms。

在硬件GPS接收机中，由于硬件相关器的应用使得处理循环相关的速度非常快，而在软件接收机中则是用FFT来代替相干运算以达到提速的目的。对于N点数据的离散傅里叶变换（DFT）表示如下:

x(n) 和y(n) 2个N点有限时域信号的循环相关表达式为

基于式(4)和式(5)推出x和y序列循环相关的DFT公式:

其中，* 代表复数的复共轭形式，X和Y是x和y信号的频域表达形式。

从式(6)中可以看到x和y序列循环相关的DFT可以用X和Y的乘积来表示，见文献[8]，其中一个是复共轭的形式。只需要对其结果做一次逆傅里叶变换(IFFT)就可以转化为时域表示形式:

z(i)max代表 z( n)集合中最大的一个,i暗示了GPS信号C/A码的初始相位。

本地载波信号分为同向和正交2种，载波信号跟本地C/A码信号相乘作为整体的本地数据:

Cj代表j号卫星的C/A码，是卫星j的估计初始相位。

2.2 循环相关过程

接收机接收到的数据跟本地数据相乘得:

3 优化PFA算法

本文将RF前端得到的中频信号以及本地载波信号进行量化，用特定的位格式存储，通过移位元操作来代替原有的乘法操作，因此所有的浮点运算都被转化为了位运算，在大多数的微处理器上，移位操作比数学运算要快得多。在PFA算法之前首先要进行基带混合运算，将接收到的导航数据跟本地载波数据相乘；PFA算法是FFT算法的变种[9]，为了解决非基二点的数据输入，但是对于输入点个数也有相应的要求，必须为 m×2n的形式，其中，m为质数，n为正整数。本文根据实际情况设计了一种适用于PFA算法的插值方案。

3.1 带通采样定理和插值参数的选择

对于中心频率在0f，信号带宽为B的中频模拟信号，如果用采样率为sf进行采样量化，采样后的信号频谱，必然是以sf为周期，对原始模拟信号的频谱进行周期扩展的结果。如果要求采样后的频谱不发生混迭，必然要求

式(11)就是带通采样定理的主要内容[10]。

因此原始采样频率sf=5.714MHz是满足带通采样定理的。

在捕获时，对原始数据进行线性插值，必然会改变采样频率，从上面的讨论可以看出，在5.332MHz ≤ fs≤ 6.572MHz 这个范围内的采样频率，不会使得一个C/A码周期（1ms）的样点数正好等于2n。由于3×211=6144，对于这个组合数，有快速的 PFA算法，而与此相对应的采样频率fs= 6.144MHz ，也能满足带通采样定理的要求。

3.2 基带混合中的数据类型转换

本文采用的RF前端为ZARLINK的GP2015，信号的十进制值为+1，−1，+3，−3，每个值存放在一个字节中，有用的信息存放在一个字节的低2bit中，本文将其先提取出来每4个存储到一个字节中，然后以字节读取的方式一次性读取 4个导航数据进行处理，数字信号的 4个十进制值表示方式见表1。本地载波的同相跟正交两路信号也被量化为2bit形式，见表2，在这种定义方式中假设本地载波信号的幅度为2.5，导航数据跟载波信号的量化值被预先存起来，用(Sign, Mag)数组表示，根据系统特定的要求，如果要进行高灵敏度处理，则可以通过增加本地载波的量化级数来实现。对于前面所述的理论推导过程，量化只是改变了数据的表示精度从而简化运算过程以提高效率，满足原有的理论，只是需要在PFA算法上进行运算方式的调整。

表1 导航数据的2bit量化

表2 本地载波的2bit量化

本地载波信号的频率跟接收到的导航数据的载波频率近似相同，然后，将接收到的数据跟本地载波的同相正交两路信号相乘，为基带混合过程，以实现后续的低通滤波，基于以上的2个定义，只通过(Sign, Mag)数组的异或运算就可得到结果，在此定义另一种存放结果的数据格式，见表 3，2个表中的数据相乘时，表3中的High Mag比特为表1中的Mag比特；Low Mag比特是表2中的Mag比特；Sign比特则是前面两表中Sign比特异或运算的结果。此时将原本的浮点数乘法运算转化为了简单的位操作，最终结果只需一次查表便可得到。

表3 基带混合结果查询表

C/A码本身就是二进制的码流，因此不用做转换而直接被预先存储起来。接下来的步骤就是对接收信号与本地信号的乘积以及C/A码分别做快速傅里叶变换。

3.3 整数化与制表法结合的PFA算法

FFT的输入数据在前面的步骤中已经被完全整数化，本系统中采用的FFT算法都是基二时间抽取方式的，见文献[11]，基二FFT的本质是将原始的DFT运算分割成为一系列2点DFT运算，首先将一个原始的N点DFT分裂为2个N/2点DFT运算，然后再将N/2点DFT分裂为2个N/4点DFT，如此操作直到每个DFT的操作点数减少到2，此算法要求FFT的点个数是2的整数幂。一个基二时间抽取方式FFT的蝶形运算器结构见图1，其中，输入数据（A, B）和输出数据（A', B'）都是复数形式，Wn代表在FFT每个运算阶段所需要的旋转因子。本系统采用的是时间抽取的方式，在进行FFT运算之前先进行数据的逆序操作，将排好序的数据序列进行后续的2点DFT运算。

本系统的输入数据只有如表3所示的8个值，并且都是整数形式，不需要转换，而对于一般形式的浮点FFT输入数据，则要先进行整形化，具体的做法可以将其扩大一定倍数转换为整数，因为FFT算法会产生固有的处理增益，对于实数输入，其处理增益是 N/2，对于复数输入，其处理增益是 N，其中，N为FFT的点数，然后在每一级蝶形运算以后对中间数据进行适当缩小以防止数据溢出。这种做法对于点数比较少的FFT应用比较合适，但是当点数很大时，由于每次蝶形运算以后对中间数据的缩小会损失一定的精度，当多次缩小操作后会使得精度大幅降低，因此，本文采取了前述的量化方式，使得输入到FFT的数据值很小，最大的模值仅为6，运算过程中出现整数越界的情况会很少，至多进行一次数据缩小处理，保证了处理精度。

图1 蝶形运算器结构

为了减少运算时间，旋转因子（正弦函数以及余弦函数值）将不会被实时计算，而是预计算出来并且存储在表格中，此表格只存储了1/4周期的正弦函数值，用16bit二进制补码表示，其对应的十进制的整数值代表了−1.0到1.0之间的小数。例如，在表中一个值为16384的数实际代表的值为16384/32768=0.5，其中，32768为16bit二进制补码代表的最大正数，FFT算法需要用到正余弦函数1/2周期的值，本文中为了节省内存开销，只存储了1/4周期的正弦函数值，记为SinTable[]，并且每个数值都有相应的索引，在后续运算用到的正余弦函数值都可以通过此索引从表中找到，计算规则见表4。其中，A代表FFT点个数，B代表正弦函数的索引值。

表4 旋转因子映射公式

SinTable[]的预计算取决于FFT的点数，在本系统中由于采样频率（6.144MHz）是固定的，因此对于C/A码的一个周期（1ms）进行操作，即6144个点，6144不是2的整数幂，但可以拆分为3×211的形式，可以采用一种特殊的FFT算法即PFA算法来进行计算，对于6144点的输入，PFA算法先进行数据分组，分为3组，然后对3组2048点的数据分别进行FFT运算，最后在进行重新组合，因此旋转因子表就只针对固定的2048点进行预计算。以下列举了SinTable[]部分数据。

当在蝶形运算器中需要用到旋转因子的时候则通过查表的方式从 SinTable[]中得到，正弦函数1/4周期外的数据以及余弦函数的数据，通过表 4中的计算规则获得，这样在整个FFT过程中，所有的浮点运算都被转换为了整数运算，并且在可能的情况下避免乘法运算，除法运算用右移操作来替代，使得运算速度大幅提高。本系统采用 32bit有符号数来存储蝶形运算中的中间数据，这样在极端的情况下(32768×6×2048)也不会出现数据越界。

4 算法分析与测试结果

4.1 运算效率分析及结果

对于一个N点DFT，其中的复数乘法次数为N2，而对于一个N点FFT，其复数运算次数为

当N的值越大的时候，FFT的加速效果越显著，每次复数运算可以被拆分为4个浮点乘和6个浮点加运算，在6144点的常规FFT中，浮点乘运算个数为（3×1024×11×4=135168），浮点加运算个数为（3×1024×11×6=202752），计算旋转因子需要的三角运算个数为（3×1024×2=6144），浮点乘运算个数同样为（3×1024×2=6144）。本算法将 FFT存在的大量浮点运算全部转化为了整数运算，并且通过预存正弦函数表的方法省掉了计算旋转因子的时间，其中，蝶形运算中存在的一些不必要的乘运算可以被省掉，当乘法两边存在零值的时候便直接将结果置为零，当乘法两边存在一的时候则保存原始数据作为结果，具体操作中会对数据两边的值进行检测，以确定时候符合条件，这样复数乘法的次数将减少:

其中，N为FFT点个数。本系统采用2种类型的数据做了测试，一种为采样率4.096MHz的模拟数据，一种为真实接收的导航数据，重采样率为6.144MHz，测试平台为双核1.83GHz PC，操作系统 Ubuntu-Linux，重复 50次实验，对结果进行了平均，将本算法跟定点运算和浮点运算做了比较，算法的性能如图2所示。

文献[6]中采用的方法是利用特殊软件相关器的方法，使用并行位操作来实现 32个点同时工作来模拟硬件相关器，其效率有很大程度的提高，但是其实施方式依然是采用循环相关的方式，与本文中采用的位操作结合整数化PFA算法做了对比，如图3所示。

可以看出由于并行位操作的循环相关运算方式，使得捕获耗时跟处理数据量成正比关系，而优化PFA算法的本质是FFT运算，PFA的运算量不会随着处理数据量的成倍增大而过快的增长，因此在长时间累积的情况下有着明显的优势。

4.2 资源节省效率分析及捕获结果

导航数据采取2bit存储的方式，在捕获阶段，为了提高灵敏度，需要做一个较长时间的累积[12]，本系统采用了 7ms的非相干累积，内存消耗量为（7 × 6144× 2 × 2/8 = 21504）byte，用字节形式存储相比，2bit方式资源占用量减少75%。在基带混合阶段，原始的复数信号是由2个浮点数组成的，内存消耗量为（21× 6144 × 2× 4 = 1032192）byte，以及相应数量的中间结果存储，本系统采用了通过表1和表2这2个表位运算的方式来得到表3的索引，以查表的方式来得到结果，因此省去了所有的浮点数存储空间，内存开销仅为3个表中的数据，可以忽略。在软件接收机的跟踪阶段，需要跟踪并处理相当长的数据，本方法对于资源的节省效率更加突出，尤其在嵌入式设备上，资源的效率显得更为重要。

本文的实验平台采用 Linux-Ubuntu 系统，Intel双核2.8处理器PC机，实验数据采用北京回龙观地区的实地接受信号。对浮点运算和整数运算2种方式得到的捕获结果做了分析，图4为浮点运算结果，图5为本系统采用的整数化方法得到的结果。实验数据采用信噪比大约为−15dB的实地接收信号，由实验结果可以看出，对于载波以及旋转因子的量化并未对捕获的精度造成很大影响，而速度得到了大幅提升。

图4 浮点运算方式的捕获结果

图5 整数化方法的捕获结果

4.3 量化级别与灵敏度分析

本地载波的量化级别会影响到捕获灵敏度，实验中采取2bit、4bit、8bit、16bit的量化方式分别对信噪比不同的导航数据进行了测试，对于捕获结果以及捕获时间进行了对比，实验结果如图6所示。

图6 不同量化级别的捕获结果

由实验可以看出当量化等级高的时候在同一信噪比下捕获到的卫星数量比较多，量化等级低的时候随着信噪比的降低，会出现丢星的情况，但是在信噪比为−25dB的情况下2bit量化仍然能够捕获到6颗卫星，达到了位置解算的要求。在实际的应用中，常常设定一个参数来衡量卫星是否被捕获到，此参数的值为循环累积的最大值与次大值的比值，比值越大说明相关峰越明显，再根据此参数与阈值的比较来判断捕获结果，在本实验中发现，2bit与4bit量化方式下，捕获参数随着信噪比的下降而降低的比较明显，信噪比每下降 5dB，捕获参数大约降低 22.5%，当信噪比非常低的情况下，参数值衰减比较大从而低于阈值而被判为未被捕获到，而高量化等级的方式下参数值随信噪比降低而衰减的幅度比较小。因此，可以根据实际情况来选取量化的级别，信噪比高的时候选择2bit或者4bit的方式，信噪比低的时候选取8bit或者16bit或者更高级别的量化方式，这样在提高捕获效率的同时保证了系统灵敏度的需求。

5 结束语

本文的重点在于提高软件接收机捕获阶段的效率，考虑到GPS数据的特殊存储结构，提出了一种使用2bit进行操作的PFA算法。在此方法中包括GPS数据，载波频率以及FFT算法中的旋转因子都将被量化为比特类型，FFT中的浮点运算被转化为整数运算。实验表明此算法使得存在大量FFT运算的捕获阶段效率得到大幅提升，另外，采用多种量化方法也大大减少了资源开销，尤其在后续的信号跟踪阶段效果更加明显。嵌入式软件接收机的实时化实现一直是一个难题，本算法消除了系统中的绝大部分浮点运算，并且进行了大量的量化操作，在性能以及资源都处于劣势的嵌入式平台上将会发挥更大的作用，为嵌入式软件接收机的实时化奠定了基础。

[1]AKOS D M.Software Radio Approach to Global Navigation Satellite System Receiver Design[D].Ohio:Ohio University, 1997.

[2]VAN NEE D, COENEN A.A new fast GPS code-acquisition technique using FFT[J].Electronics Letters, 1991, 27(2):158-160.

[3]STARZYK J A, ZHU Z.Averaging correlation for C/A code acquisition and tracking in frequency domain[A].IEEE MWSCAS 2001[C].Dayton, Ohio, 2001.905-908

[4]KWONHUE C, KYUNGWHOON C, TAEJIN J.Adaptive PN code acquisition using instantaneous power-scaled detection threshold under rayleigh fading and pulsed gaussian noise jamming[J].IEEE Trans Commun, 2002, 50(8):1232-1235.

[5]CHUANG M Y, FENG K T.Adaptive GPS acquisition technique in weak signal environment[A].IEEE Vehicular Technology Conference[C].Melbourne, 2006.1612-1616.

[6]BRENT M, LEDVINA, MARK L, et al.Bit-wise parallel algorithms for efficient software correlation applied to a GPS software receiver[J].IEEE Transactions on Wireless Communications, 2004, 3(5):1469-1473.

[7]ANDREW J, VITERB I.CDMA Principles of Spread Spectrum Communication[M].MA:Addison-Wesley, 1995.120-135.

[8]OPPENHEIM A V, SCHAFER R W.Discrete-time Signal Processing,2nd Ed[M].New Jersey:PrenticeHall, Inc, 1999.541-669.

[9]KOLBA D, PARKS T.A prime factor FFT algorithm using high-speed convolution[J].IEEE Transactions on Acoustics, Speech and Signal Processing, 1977,25(4):281-294.

[10]曹志刚, 钱亚生.现代通信原理[M].北京:清华大学出版社, 2006.108-110.CAO Z G, QIAN Y S.The Principle of Modern Communication[M].Beijing:Tsinghua University Press,2006.108-110.

[11]BURRUS C, ESCHENBACHER P.An in-place, in-order prime factor FFT algorithm[J].IEEE Transactions on Acoustics, Speech, and Signal Processing, 1981, 29(4):806-816.

[12]ZIEDAN N I, GARRISON J L.Unaided acquisition of weak GPS signals using circular correlation or double-block zero padding[A].IEEE PLANS,2004 Position Location and Navigation Symposium[C].2004.461-470.