【GNSS】GPS P码仿真与实现——从C/A码到精密测距码【含matlab代码】

博客七：GPS P码仿真与实现——从C/A码到精密测距码

本文说明：本系列博客基于作者发表于IEEE的学术论文——Z. Shang, “GPS C/A code simulation analysis and GPS signal capture analysis based on MATLAB,” in 2023 IEEE International Conference on Integrated Circuits and Communication Systems (ICICACS), Xi’an, China, 2023, pp. 1-6. doi: 10.1109/ICICACS57324.2023.10248505。

当年论文中只重点讨论了c/a码

一、为什么需要P码？——从“粗捕获”到“精密测距”

在博客一中，我们详细讨论了C/A码（粗捕获码）的生成原理。C/A码作为民用GPS信号的核心，具有周期短（1ms）、结构简单、易于捕获等优点。但它也有明显的局限性：

P码的“精密”体现在两个层面：

P码被调制在L1载波的同相（I）支路上，与C/A码（正交Q支路）形成正交复用。这种设计使得民用接收机只需处理C/A码即可获得定位服务，而军用/授权接收机则可同时利用P码实现更高精度。

二、P码的生成原理：四个移位寄存器的协同工作

2.1 与C/A码的本质区别

C/A码由两个10级移位寄存器（G1和G2）组合生成。P码的生成结构则复杂得多：

P码由两个12级的M序列发生器组合而成：

$$P(t)=PN1(t)\oplus PN2(t)$$

其中PN1和PN2各自又由两个子序列组合而成，形成四寄存器结构。

2.2 四寄存器结构详解

P码生成器由四个12位移位寄存器组成：

┌─────────────────────────────────────────────────────────────────┐
│ P码生成器架构 │
├─────────────────────────────────────────────────────────────────┤
│ │
│ ┌──────────────┐ ┌──────────────┐ │
│ │ X1A寄存器 │ │ X1B寄存器 │ │
│ │ (12级) │ │ (12级) │ │
│ │ 反馈: 14501 │ │ 反馈: 14501 │ │
│ └──────┬───────┘ └──────┬───────┘ │
│ │ │ │
│ └───────┬───────────┘ │
│ ↓ │
│ PN1(t) = X1A ⊕ X1B │
│ ↓ │
│ ┌───────┴───────┐ │
│ ↓ ↓ │
│ ┌──────────────┐ ┌──────────────┐ │
│ │ X2A寄存器 │ │ X2B寄存器 │ │
│ │ (12级) │ │ (12级) │ │
│ │ 反馈: 17147 │ │ 反馈: 17147 │ │
│ └──────┬───────┘ └──────┬───────┘ │
│ │ │ │
│ └───────┬───────────┘ │
│ ↓ │
│ PN2(t) = X2A ⊕ X2B │
│ ↓ │
│ ┌───────┴───────┐ │
│ ↓ ↓ │
│ P(t) = PN1(t) ⊕ PN2(t + 延时) │
│ │
└─────────────────────────────────────────────────────────────────┘

关键参数：

为什么是4092 chips？

12级移位寄存器的完整周期本应为 $2^{12}-1=4095$。但X1A和X2A被设计为短周期（shortcycled）到4092 chips。这一设计使得PN1和PN2的周期为：

$$T_{PN1}=4092\times 3750=15,345,000\text{ chips}\approx 1.5\text{ 秒}$$

X1A和X2A各运行3750次短周期，覆盖1.5秒。PN1和PN2再组合，最终形成周期约267天的超长P码序列。

2.3 P码的相位分配

与C/A码类似，不同卫星的P码通过不同的时间延迟来区分。GPS ICD文档中定义了C/A码与P码的配对相位分配表。

关键事实：C/A码与P码的相位分配是不可分割的配对——每颗卫星的C/A码相位和P码相位由ICD文档共同定义，不能独立选择。

2.4 为什么P码如此难以生成？

与C/A码相比，P码的生成面临几个核心挑战：

三、MATLAB仿真实现

3.1 设计思路

由于P码的完整周期极长（267天），仿真时通常采用按需生成策略：

3.2 核心数据结构

% P码相位分配表（部分，完整表见IS-GPS-200）
% 每行: [PRN, X2_Delay, P_Relative_Advance]
pCodePhaseTable = [
1, 5, 0; % PRN 1: X2延时5 chips
2, 6, 0; % PRN 2: X2延时6 chips
% … 完整表共32+ entries
];

3.3 12级移位寄存器核心函数

function [seq, state] = lfsr12(init_state, feedback, num_chips, reset_after)
% 12级线性反馈移位寄存器
% 输入: init_state – 12位初始状态 (1×12, 0/1)
% feedback – 反馈多项式（八进制表示）
% num_chips – 需要生成的码片数
% reset_after – 每多少chip重置一次（短周期）
% 输出: seq – 生成的码序列

state = init_state;
seq = zeros(1, num_chips);

% 将八进制反馈多项式转换为抽头位置
taps = octalToTaps(feedback); % 例如 14501 → [1, 3, 5, 10, 12]

for i = 1:num_chips
% 输出最后一级
seq(i) = state(end);

% 计算反馈（各抽头的异或）
fb = 0;
for t = taps
fb = xor(fb, state(t));
end

% 移位
state = [fb, state(1:end–1)];

% 短周期重置（如需要）
if reset_after > 0 && mod(i, reset_after) == 0
state = init_state;
end
end
end

3.4 主生成函数结构

function pCode = generatePCode(prn, gps_week, tow_ms, num_chips)
% 生成指定卫星在指定时刻的P码
% 输入: prn – 卫星号 (1-32)
% gps_week – GPS周数
% tow_ms – 周内时（毫秒）
% num_chips – 需要生成的码片数
% 输出: pCode – P码序列 (+1/-1)

% 1. 查表获取该PRN的P码参数
[x2_delay, p_advance] = getPCodeParams(prn);

% 2. 计算四个寄存器的初始状态
% (基于GPS时间和PRN参数的复杂计算)
[init_X1A, init_X1B, init_X2A, init_X2B] = …
computeInitialStates(gps_week, tow_ms, prn, x2_delay);

% 3. 生成PN1序列
pn1 = lfsr12(init_X1A, 14501, num_chips, 4092) …
xor lfsr12(init_X1B, 14501, num_chips, 4092);

% 4. 生成PN2序列（含延时）
pn2_delayed = lfsr12(init_X2A, 17147, num_chips + x2_delay, 4092) …
xor lfsr12(init_X2B, 17147, num_chips + x2_delay, 4092);
pn2 = pn2_delayed(x2_delay + 1 : end);

% 5. 组合为P码
pCode_binary = xor(pn1, pn2);

% 6. 映射为 +1/-1
pCode = 2 * pCode_binary – 1; % 0→-1, 1→+1
end

3.5 一个重要的工程提醒

MATLAB从R2021b开始提供了官方的 gpsPCode 函数，可直接生成GPS卫星的P码：

% MATLAB官方函数（R2021b+）
pCode = gpsPCode(prn, 'NumChips', 10230, 'SampleRate', 10.23e6);

在学术研究和教学场景中，建议先独立实现生成算法以理解原理，再考虑使用官方工具函数提高效率。

四、仿真验证与分析

4.1 P码的时域特征

与C/A码一样，P码也是±1的伪随机序列。但由于码速率是C/A码的10倍（10.23MHz），P码的码片宽度仅为C/A码的1/10：

C/A码: |‾‾‾‾‾‾‾‾‾‾‾|___|‾‾‾‾‾‾‾‾‾‾‾|___| (1 chip ≈ 293m)
P码: |‾‾|__|‾‾|__|‾‾|__|‾‾|__|‾‾|__| (1 chip ≈ 29.3m)

更窄的码片宽度意味着更高的测距分辨率，这是P码“精密”二字的物理基础。

4.2 自相关与互相关特性

P码的自相关/互相关特性与C/A码类似——尖锐的自相关峰和低互相关。但由于P码周期极长，其互相关性能更为优越：

4.3 P码与C/A码的频谱对比

由于码速率相差10倍，两者的频谱主瓣宽度差异显著：

功率谱密度
↑
│ C/A码 (1.023 MHz) P码 (10.23 MHz)
│ ████ ████████████
│ ██████ ██████████████
│ ████████ ████████████████
│ ██████████ ██████████████████
│ ████████████ ████████████████████
│ ██████████████ ██████████████████████
└──────────────────────────────────────────────→ 频率
0 1 2 3 4 5 6 7 8 9 10 (MHz)

P码的主瓣宽度约为C/A码的10倍。更宽的频谱意味着：

• 更强的抗窄带干扰能力
• 更高的处理增益（$10{\log }_{10}(10.23\times 10^6)\approx 70$ dB，远超C/A码的约43dB）

五、P码仿真中的常见问题

5.1 关于“P(Y)码”的说明

严格来说，GPS军用的加密码称为P(Y)码——即P码经过加密后的版本。普通P码（未加密）在理论上是可以生成的，但实际在轨卫星发射的是P(Y)码，没有解密密钥就无法使用。

5.2 状态初始化

P码生成中最复杂的部分是四个寄存器的初始状态计算。这需要：

• 精确的GPS时间（周数 + 周内秒）
• 正确的PRN对应的延时参数
• 短周期（4092 chips）的精确控制

5.3 存储与计算

P码的完整周期约267天，在1.023MHz码速率下相当于约 $2.35\times 10^{14}$ 个码片。切勿尝试一次性生成完整的P码序列——这不仅不现实，也完全没有必要。实际应用中只需生成所需时长的局部序列即可。

六、总结与系列关联

本篇核心收获

与系列博客的关联

系列博客完整目录

引用本文：本系列博客介绍的GPS C/A码与P码生成、信号模拟、捕获、跟踪与PVT解算算法，基于作者发表于IEEE的学术论文——Z. Shang, “GPS C/A code simulation analysis and GPS signal capture analysis based on MATLAB,” in 2023 IEEE International Conference on Integrated Circuits and Communication Systems (ICICACS), Xi’an, China, 2023, pp. 1-6. doi: 10.1109/ICICACS57324.2023.10248505。

（博客七 完）