Matlab生成GNSS PRN码：Gold码、二级代码与频谱分析全解析

简介：面向Matlab开发者的GNSS信号仿真工具包，聚焦GPS、Galileo、北斗二号三大卫星导航系统，提供L1CA、L2CM、L2CL、L5I、L5Q、E1B、E1C、E5aI、E5aQ、E5bI、E5bQ、B1I等多种频点信号的PRN码、二级码、无数据信号与频谱生成能力，满足不同频点与码型的仿真需求。资源共56个文件，以.m脚本和.mat数据容器为主，附带PDF/TeX格式的技术文档、许可证及文本说明，压缩后仅6.36MB，目录划分清晰，便于按文档、码表、真实数据分块检索。文档部分收录GNSS代码、信号、频谱和自相关函数摘要，并注明官方ICD文件获取途径；prn_codes子目录保存各系统尚未重采样的码序列Matlab数据容器，可直接加载用于仿真验证；real_data子目录提供一段L2频段实际采集数据，采样率为2.5MSps、时长为1.5秒，其中已知含有GPS L2C 26号PRN，适合检验捕获算法与信号处理流程，也可用于教学演示。已有921人学习下载，对从事GNSS接收机算法验证、通信仿真、教学实验的研究人员具有实用价值。

1. 用Matlab手搓GNSS PRN：从码生成、二级代码到频谱分析，一次讲透

做GNSS接收机或者信号级仿真的人早晚会遇到同一个坎：手里没有采集板、没有射频前端，但需要用Matlab生成一段“看起来像GPS卫星发出来的基带信号”。不管是验证捕获算法、测距精度，还是给课程设计交差，第一关都是PRN码。真正动手之后会发现，PRN码不是随机数，而是Gold码，靠G1/G2两个移位寄存器就能生成；GPS L1 C/A码只是第一步，L5这类新信号还要叠加二级代码（secondary code），构造无数据信号（pilot），再配上频谱分析确认带宽和码率对不对。这篇文章按“码生成 → 叠加二级代码 → 频谱验证 → 避坑 → 接入捕获”这条线，把一套能直接跑的Matlab方案拆给你。

2. 生成GPS PRN码：G1/G2移位寄存器与1023个码片的完整实现

2.1 为什么PRN码不能直接查表，要先搞懂Gold码

GPS L1 C/A码的长度是1023个码片，码速率1.023 Mchip/s，周期1ms。网上能查到现成序列，甚至有人直接把1023个码片存成数组。但我不建议这么做，原因有两个：一是查表得到的序列一旦贴错、截断或者极性反了，你连错在哪都发现不了；二是后续做L5、B1I、E1的时候，同样是Gold码结构，只是多项式不同，查表法换一个信号就要重新找数据，改动成本很高。

正确做法是用生成多项式去推。C/A码由两个10级移位寄存器G1和G2生成，本质上是两个最大长度序列（m序列）的异或。G1的特征多项式是G1(x)=1+x³+x¹⁰，G2的特征多项式是G2(x)=1+x²+x³+x⁶+x⁸+x⁹+x¹⁰。两个寄存器初始状态全为1，每个时钟周期做一次反馈移位，G1和G2的输出异或之后就是C/A码。

这里有一个容易被忽略的关键点：G2不能直接拿来和G1异或，必须对G2的输出做“相位选择”。不同PRN编号对应G2寄存器不同抽头位置的异或结果。标准里叫code delay，实现上通常取两个抽头，比如PRN 1取G2的第2级和第6级异或，PRN 2取第3级和第7级。IS-GPS-200里给了32颗卫星的完整相位选择表，网上流传的版本可能有细微差别，我一般以gnss-sdr源码里那张表为准。

2.2 生成函数代码与相位选择表

下面这段是我自己工程里在用的版本，返回单颗卫星一个完整周期的C/A码，取值是0/1逻辑电平，没有映射到±1，方便后面做BPSK调制时自己决定极性。

function ca = generateCA(prn) % generateCA 生成GPS L1 C/A码 % 输入prn为卫星编号(1-32)，输出ca为1x1023的0/1序列 % G1: 1 + x^3 + x^10 % G2: 1 + x^2 + x^3 + x^6 + x^8 + x^9 + x^10 g1 = ones(1, 10); % G1移位寄存器，初始全1 g2 = ones(1, 10); % G2移位寄存器，初始全1 taps = phaseTable(prn); % 查G2相位选择表 ca = zeros(1, 1023); % 预分配，1周期1023个码片 for k = 1:1023 % 当前码片 = G1末级 xor (G2两个抽头异或) ca(k) = xor(g1(10), xor(g2(taps(1)), g2(taps(2)))); % G1反馈：G1的第10级与第3级异或 fb1 = xor(g1(10), g1(3)); % G2反馈：G2的第10、9、8、6、3、2级异或 fb2 = xor(xor(xor(xor(g2(10), g2(9)), g2(8)), g2(6)), ... xor(g2(3), g2(2))); % 移位：新反馈值放入第1级，原有内容向右移 g1 = [fb1, g1(1:9)]; g2 = [fb2, g2(1:9)]; end end function taps = phaseTable(prn) % phaseTable GPS L1 C/A码G2相位选择表（部分） % 完整32项见IS-GPS-200或gnss-sdr源码，这里列前16颗 table = { 1, [2 6]; 2, [3 7]; 3, [4 8]; 4, [5 9]; 5, [1 9]; 6, [2 10]; 7, [1 8]; 8, [2 9]; 9, [3 10]; 10, [2 3]; 11, [3 4]; 12, [5 6]; 13, [6 7]; 14, [7 8]; 15, [8 9]; 16, [9 10] }; taps = table{prn, 2}; end

这段代码的核心就是把标准里“G2延迟N个码片”的写法，换成了“两个抽头异或”的等价实现。需要注意G2的反馈项顺序：特征多项式从高次到低次写是x¹⁰+x⁹+x⁸+x⁶+x³+x²+1，所以反馈异或的寄存器下标是10、9、8、6、3、2，不是抽头表里的那些下标。很多人第一次写会把反馈项和抽头表搞混，反馈算错会导致整个序列变成噪声。

代码里的phaseTable目前只写了16颗。如果你的工程用到PRN 17以后，把标准表剩余项补进去就行，结构不变。我不建议用查表方式在运行时去读这32项以外的内容，因为民用GPS卫星PRN范围就是1到32，再多也不需要。

2.3 先验证自相关/互相关，再往下走

生成完码不能直接拿去叠加信号，第一步永远是把码的统计特性验一遍。C/A码的数学性质就两条：自相关在零偏移处等于1023，其余位置接近0；不同PRN码之间互相关接近0。下面这段代码做两个检查：

% 验证PRN 1的自相关性 ca1 = generateCA(1); % 做周期相关，偏移范围取±200个码片即可 [acor, lag] = xcorr(ca1, ca1, 200, 'normalized'); figure; plot(lag, acor); title('PRN 1 自相关'); xlabel('码片偏移'); ylabel('相关值'); grid on; % 验证PRN 1和PRN 2的互相关性 ca2 = generateCA(2); [ccor, clag] = xcorr(ca1, ca2, 200, 'normalized'); figure; plot(clag, ccor); title('PRN 1 vs PRN 2 互相关'); xlabel('码片偏移'); ylabel('相关值'); grid on;

xcorr的normalized选项会把峰值归一化到1，C/A码周期相关在零偏移处应当是1，非零偏移处绝对值大约在0.03到0.08这个量级。如果你的结果出现了副峰超过0.5的情况，基本可以断定抽头表写错了，不要继续往下做。

另外一个容易出问题的细节是极性约定。generateCA返回的是0/1序列，你在做BPSK之前要把0映射成+1、1映射成-1（或者反过来）。两种约定在捕获结果里只差一个符号，但如果你同时用了外部数据源（比如采集的IQ文件），极性不一致会导致相关峰变成负峰。我习惯统一用“0→+1、1→-1”，并且在代码注释里写清楚，避免隔一个月自己都忘了。

3. 二级代码与无数据信号：L5导频支路的搭建方法

3.1 二级代码是什么，为什么L1 C/A没有而L5有

二级代码（secondary code）也叫叠加码或NH码（Neuman-Hofman码），是叠在PRN码之上的又一层扩频序列。GPS L1 C/A信号没有二级代码，一个PRN周期就是1ms，直接调制50bps的导航数据。到了GPS L5和北斗B1I，信号结构变了：码速率更高，数据率也做了分层，于是引入二级代码来改善谱线分布和相关性能。

L5信号的完整结构分两个支路：I支路（In-phase）在PRN码之上叠加NH10码，再调制50bps的CNAV导航数据；Q支路（Quadrature）同样叠加NH20码，但不加导航数据，这是一个真正的无数据信号，也叫pilot分量。接收机通常优先用pilot支路做载波跟踪，因为没有数据比特翻转，环路不会因为位同步错误而失锁。

不同系统的二级代码参数差异很大，做仿真前先对一张表：

一个常见的误区是以为二级代码会改变PRN码本身的码片序列。实际上二级代码是以“符号级”叠加在PRN序列外层的：每个二级码片持续若干个PRN周期，PRN码内部结构完全不变。L5 Q支路里，NH20的一个码片恰好等于10.23 MHz码速率下的1ms，也就是10230个PRN码片，一个完整NH20周期20ms，包含204600个码片。

3.2 L5 Q支路无数据信号的生成代码

下面这段代码生成L5 Q支路一个完整NH20周期内的无数据信号，采样率取10.23 MHz，也就是每个PRN码片采一个点：

% 生成L5 Q支路无数据信号（pilot） % 参数：采样率 = 码速率 = 10.23 MHz，1个采样点 = 1个码片 fs = 10.23e6; codeRate = 10.23e6; samplesPerChip = round(fs / codeRate); % 此处 = 1 % NH20二级代码，以IS-GPS-705为准 nh20 = [0 0 0 0 1 0 1 0 0 1 1 0 1 1 0 1 0 1 1 1]; % 取PRN 1的L5码（L5的PRN码和L1不同，这里用L1码代替占位） % 实际工程中L5使用自己的PRN生成多项式，下文说明 ca = generateCA(1); % 1023码片，实际L5 PRN周期是10230码片 pri = [ca, ca]; % L5的PRN码长正好是L1的10倍，此处拼接占位 % 构造完整的20ms导频序列 signal = zeros(1, codeRate * 0.02); % 20ms * 10.23MHz = 204600个采样点 for nhIdx = 1:20 % 当前NH码片极性：0映射为+1，1映射为-1 polarity = 1 - 2 * nh20(nhIdx); startIdx = (nhIdx - 1) * 10230 + 1; endIdx = nhIdx * 10230; signal(startIdx:endIdx) = polarity * pri; end % 检查信号能量和基础统计量 fprintf('信号长度: %d 采样点\n', length(signal)); fprintf('信号能量: %.2f\n', sum(abs(signal).^2));

这段代码里最需要解释的是polarity = 1 - 2 * nh20(nhIdx)这一行，它把0/1序列映射成±1。映射关系直接决定了后续频谱相关分析里峰值是正是负，建议固定住这个约定。另外这里有一个偷懒的处理：L5真实PRN码是10230个码片，不是L1的1023个码片。为了把流程跑通，我把L1码拼接了10倍，严格来说L5有自己的生成多项式，L5 PRN序列不能由L1序列拼接得到。如果你只做频谱形态验证，这种方法够用；但如果要做接收机捕获级别的仿真，必须换成L5真正多项式生成的10230码片序列。

3.3 数据支路与导频支路的功率分配和叠加顺序

GPS L5把I、Q两支路放在正交载波上，I支路是数据信号，Q支路是导频信号。实际发送时，I支路功率和Q支路功率相等或者按约1:1分配，因此完整的L5基带信号可以写成：

% 构造L5 I支路数据信号（示意，忽略导航数据内容） % 假设导航数据全部为+1，叠加NH10 nh10 = [0 0 0 0 1 0 1 0 0 1]; dataI = zeros(size(signal)); for nhIdx = 1:20 % NH10周期10ms，一个NH10周期内dataI每1ms切换一次极性 nhPos = mod(nhIdx - 1, 10) + 1; polarity = 1 - 2 * nh10(nhPos); startIdx = (nhIdx - 1) * 10230 + 1; endIdx = nhIdx * 10230; dataI(startIdx:endIdx) = polarity; end % 完整L5基带：I支路 + j*Q支路 l5Signal = dataI + 1j * signal;

叠加顺序要特别注意：数据支路的导航数据速率是50bps，一个导航数据位持续20ms，恰好等于NH10的2个周期；导频支路没有数据，NH20直接以20ms为周期循环。如果导航数据位边界和NH10周期边界没有对齐，仿真出来的信号相位会在整体上产生随机跳变，这种跳变不是错误，但会让频谱分析和捕获结果看起来不稳定。我一般在仿真里把第一个导航数据位的起始点定义为t=0，所有序列都从这个时间零点开始排。

4. 无数据信号的频谱分析：用pwelch验证主瓣与码速率

4.1 频谱要回答的三个问题

把无数据信号生成出来之后，频谱分析是验证“这个信号像不像真东西”的最直接手段。对着功率谱图，我一般只看三件事：主瓣宽度对不对、零点位置对不对、带外能量有没有异常。

BPSK-R调制信号的功率谱包络是sinc²函数形状，主瓣零点之间的宽度是2倍码速率。以L5 Q支路为例，码速率10.23 MHz，主瓣两个零点之间的宽度应当是20.46 MHz。如果看到的主瓣零点宽度只有一半或者翻倍，要么码速率设置错了，要么采样率和码速率的对应关系不对。

L1 C/A信号同理：码速率1.023 MHz，零点宽度2.046 MHz。接收机带宽通常取2 MHz，就是因为信号主瓣正好占这么宽。这个数字可以直接从频谱图上读出来，不需要看任何文档。

还有一个容易被忽略的点：二级代码叠加后，频谱包络不会变。NH码的速率是1 kbps，远低于PRN码速率，它对频谱的影响只体现在极窄频率范围内的谱线结构，不会改变sinc²的包络形状。如果你发现加上二级代码之后主瓣变宽了，说明你的叠加方式根本不是“信号级叠加”，而是把PRN序列本身给改了。

4.2 功率谱估计代码与参数设置

用pwelch做功率谱估计是最稳的选择，直接fft也行但方差太大。下面这段代码对前面生成的L5导频信号做频谱分析：

% 对L5 Q无数据信号做功率谱估计 % 注意：上一节的signal信号是204600点，时长为20ms fs = 10.23e6; % pwelch参数：窗长4096，重叠2048，FFT点数4096 [pxx, f] = pwelch(signal, hamming(4096), 2048, 4096, fs, 'centered'); % 转为dB pxxDb = 10 * log10(pxx + eps); % 画频谱 figure; plot(f / 1e6, pxxDb, 'LineWidth', 1); xlabel('频率 (MHz)'); ylabel('功率谱密度 (dB)'); title('L5 Q无数据信号功率谱'); grid on; % 标注主瓣零点位置 hold on; xline(-10.23, 'r--', '左零点'); xline(10.23, 'r--', '右零点'); hold off; % 输出主瓣宽度 fprintf('主瓣零点宽度: %.2f MHz (理论值20.46 MHz)\n', 2 * 10.23);

pwelch的窗函数选hamming，FFT点数4096，在10.23 MHz采样率下频率分辨率约2.5 kHz，足够分辨NH20码产生的1 kHz级谱线影响。如果你的电脑性能足够，把FFT点数提升到16384，谱线会更平滑，但是计算时间会明显增加。204600个点对个人电脑来说没有任何压力，全程跑下来应该在一秒以内。

读图的时候有一个技巧：把光标放在主瓣零点的位置，读数应当正好是±10.23 MHz。如果零点位置偏移，先检查fs参数是不是写错了；如果零点位置对但整体有直流尖峰，说明信号有直流偏置，回查极性映射那一步，确认0和1是否真的映射成了对称的±1。

4.3 有/无二级代码的频谱差异在哪里

很多人做完带NH码的信号和纯PRN信号对比，发现两张频谱图几乎一样，然后开始怀疑自己代码写错了。其实这是正常现象：PRN码速率和NH码速率相差一万倍，NH码的频谱贡献被压缩在基带极窄区域，在宽带频谱图上根本看不出来。

想要直观看到二级代码的贡献，需要把频谱分析范围缩小。对L5信号的204600点序列做FFT，观察±5 kHz范围内的频谱细节，会发现带NH码的信号在1 kHz间隔处出现离散谱线，而纯PRN信号只有连续谱。这段代码可以做细看：

% 仔细观察低频段的谱线结构 N = length(signal); faxis = (0:N-1) / N * fs; % 单边频率轴 faxis = faxis(1:N/2); spec = abs(fft(signal)).^2; spec = spec(1:N/2); figure; plot(faxis / 1e3, 10 * log10(spec + eps)); xlim([0 5]); % 只看0-5kHz范围 xlabel('频率 (kHz)'); ylabel('功率 (dB)'); title('低频段谱线：NH20码的1kHz间隔分量'); grid on;

从工程角度说，这个分辨能力意味着你的仿真信号可以用于接收机的精确码相位对齐测试。接收机捕获pilot信号之后，依靠的就是这些在1 kHz间隔上的谱线来辅助比特同步。如果你的信号生成代码不包含二级代码，接收机在实际跟踪时就会丢掉这一层同步信息。

5. 避坑：生成GNSS PRN时最容易翻车的5个细节

5.1 现象：自相关峰矮了一半，副峰抬到0.5

原因：G2抽头表写错了，某几个PRN的抽头下标和标准不一致。我遇到过用网上某篇博客的相位表，PRN 1到10都正常，PRN 11往后自相关就出现了明显副峰。

解决：不要凭记忆写相位表。打开IS-GPS-200或者gnss-sdr的ca_codes_generator.cpp，把32项逐条核对。核对的方法是用已知正确的PRN 1序列做对照，把自相关峰旁边的副峰电平打印出来，如果某颗星副峰超过0.1，立刻查那张表。

5.2 现象：NH码极性搞反，相关峰从1变成负峰

原因：NH码序列有两种极性约定，标准文档给出的是0/1位流，但接收机内部常用的本地序列可能是直接存的±1。如果你从某个资料里抄了NH码数组，没确认映射关系就直接用，很容易整体反相。

解决：先做一个自检。生成10000个码片的导频信号，和一个已知正确的NH码序列做相关，最大化搜索正负两种极性，看哪个方向峰值接近1。以后每次换信号体制都重复这个检查，不需要跑完整链路。

5.3 现象：频谱图全是毛刺，看不出sinc²主瓣

原因：采样点数太少，或者观测时长不是PRN周期整数倍。C/A码是周期信号，如果截取的时长只有1ms但采样率不是1.023 MHz的整数倍，FFT频谱会发生频谱泄漏，主瓣被撕碎。

解决：把采样率固定为码速率整数倍，观测时长取整毫秒。做频谱分析之前先确认信号长度对不上码周期的情况，L1码1ms周期，至少取10ms；L5码1ms周期（10230个码片），同样至少取10个周期以上。如果观测时长受限，就在pwelch里加大窗函数重叠率，别直接对短序列做矩形窗FFT。

5.4 现象：仿真20ms L5信号内存就爆炸了

原因：采样率设置出现数量级错误。有人会把L5的采样率写成102.3 MHz甚至更高，一个20ms的信号就是200多万个点，再做全序列FFT，内存和计算时间都扛不住。

解决：先想清楚你要验证什么。验证PRN码相关性，用码片级仿真（每码片1个点）；验证频谱包络，用10.23 MHz就够；验证带外辐射或者滤波效果，才需要过采样到30 MHz以上。不要一上来就把采样率拉满，等小规模方案跑通再放大。

5.5 现象：仿真信号拿去和真实硬件数据做捕获，相关峰总对不上

原因：真实GNSS信号的前端带宽有限、有多普勒频偏，你生成的理想基带信号直接做相关当然对不上。更隐蔽的问题是信号起始时间没对齐，仿真信号从t=0开始，而采集数据里有任意码相位偏移。

解决：先做码相位搜索而不是指望零偏移对齐。用FFT循环相关把全部1023个码相位都扫一遍，找到峰值位置后，再去比对码相位和多普勒估计值。和真实数据比对时，先把信号经过一个和硬件前端一致的通带滤波器，再进捕获模块，这样结果才有可比性。

6. 把生成的码接到捕获验证与真实GNSS原始数据上

6.1 用FFT循环相关做捕获峰值验证

前面生成的PRN和导频信号，最终价值要落在捕获验证上。最实用的验证方法是用FFT循环相关，把全部码相位一次性扫完。核心思路是：本地码做FFT，接收信号做FFT，频域相乘后反变换，时域结果就是全相位相关峰。下面这段代码用L1 C/A码演示：

% 构造一个带已知延迟的接收信号 ca = generateCA(1); codeBpsk = 1 - 2 * ca; % 0/1转换为±1 txSignal = [zeros(1, 500), codeBpsk, zeros(1, 523)]; % 人为延迟500个码片 % FFT循环相关 fftSize = length(txSignal); rxFft = fft(txSignal, fftSize); codeFft = conj(fft(codeBpsk, fftSize)); corrResult = ifft(rxFft .* codeFft, fftSize); % 找峰值位置 [peakVal, peakIdx] = max(abs(corrResult)); fprintf('峰值位置: %d (期望: 501)\n', peakIdx);

如果峰值位置不是500附近，先检查周期相关的边界效应。fftSize取的是整个接收信号长度，本地码直接补零到同样长度，这种方法在信号不是整周期时会出现相关平台。工程里更稳的做法是把本地码重复拼接成两倍周期长，再和单周期接收信号做循环相关，峰值位置对应真正延迟。这个差异要自己在仿真里体会一遍，不然会误判捕获模块的代码有问题。

6.2 与手机GNSS原始数据、STK联合仿真的对照边界

用Matlab生成的信号和真实世界之间，最终要对齐一次才算闭环。常见做法是拿手机在空旷地采集一份GNSS原始数据（Android设备的GNSS Logger可以直接输出包含码相位、多普勒、载波相的测量值），把采集数据里的真实信号和你的仿真PRN做互相关比对。这里的边界是：手机采集数据经过了射频前端滤波、自动增益控制和量化，幅度和相位和你理想的±1序列差别很大，互相关峰不会特别尖锐，但只要码相位和多普勒估计正确，峰的位置应当和手机输出的码相位测量值一致。

另外有人会用STK 12.2和Matlab联合仿真，把卫星轨道、钟差、多普勒变化率加进信号生成模型。我建议先在本地把静态信号的PRN和相关峰验证通过，再考虑加多普勒。多普勒会让信号频谱整体搬移，对你的PRN码生成没有影响，但会显著影响捕获搜索维度。先把二维搜索（码相位、多普勒）做通，再去做轨道级的联合仿真，否则排查问题时你会分不清是PRN码的问题还是多普勒模型的问题。

我的习惯是：生成任何一种新信号体制的PRN码之后，第一件事永远是打一张自相关图，然后打一张功率谱图，确认码长和码速率两个数字都正确，再谈捕获和测距。这两个图只要花两分钟，但能省掉后面一整天的排查时间。希望帮到你。

本文还有配套的精品资源，点击获取