简介面向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)1x³x¹⁰G2的特征多项式是G2(x)1x²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选项会把峰值归一化到1C/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码速率二级代码二级码速率周期数据支路GPS L1 C/A1.023 Mchip/s无——50 bpsGPS L5 I10.23 Mchip/sNH101 kbps10 msCNAV 50 bpsGPS L5 Q10.23 Mchip/sNH201 kbps20 ms无北斗B1I2.046 Mchip/sNH码20bit1 kbps20 ms50 bps一个常见的误区是以为二级代码会改变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 MHz1个采样点 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映射为11映射为-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周期边界没有对齐仿真出来的信号相位会在整体上产生随机跳变这种跳变不是错误但会让频谱分析和捕获结果看起来不稳定。我一般在仿真里把第一个导航数据位的起始点定义为t0所有序列都从这个时间零点开始排。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重叠2048FFT点数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的窗函数选hammingFFT点数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周期至少取10msL5码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信号的前端带宽有限、有多普勒频偏你生成的理想基带信号直接做相关当然对不上。更隐蔽的问题是信号起始时间没对齐仿真信号从t0开始而采集数据里有任意码相位偏移。解决先做码相位搜索而不是指望零偏移对齐。用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码之后第一件事永远是打一张自相关图然后打一张功率谱图确认码长和码速率两个数字都正确再谈捕获和测距。这两个图只要花两分钟但能省掉后面一整天的排查时间。希望帮到你。本文还有配套的精品资源点击获取