简介这是一套面向 GNSS 信号处理与 Matlab 开发者的代码生成工具能够生成 GPSL1CA/L2CM/L2CL/L5I/L5Q、GalileoE1B/E1C/E5aI/E5aQ/E5bI/E5bQ以及北斗 B1I 的 PRN 码、二级码与无数据信号波形并输出对应频谱和自相关函数可帮助简化接收机算法验证中的信号构造环节。资源含 PDF 格式的 GNSS 信号原理摘要、按系统整理的 PRN 码 .mat 数据容器以及一段 L2 频段 GPS L2C 实际采集数据方便学习者对照真实捕获结果测试自己的处理流程。压缩包共 56 个文件以 Matlab 脚本.m和 .mat 数据文件为主另有 TEX/PDF 文档、PNG 示意图及开源许可说明整体仅 6.36MB目录划分清楚、易于检索。已有 921 人学习下载。通过阅读源码可以快速掌握 BOC 调制、频谱生成与二级码构造的具体实现适合正在做 GNSS 接收机仿真、课程设计或相关毕业设计的科研人员和学生。1. 生成一套能用的GNSS PRN和频谱比想象中多三步路做GNSS接收机或者信号模拟器的人迟早都要面对同一个问题手里没有真实的中频采样数据软件接收机的捕获、跟踪、定位全都没法验证。这个标题给的出路是用Matlab代码生成GPS的PRN码、二级代码和无数据信号再拉出频谱确认信号没毛病。它能解决的是从零造出一段“看起来像卫星发出来”的基带信号用来喂给相关器、测捕获灵敏度、验证多普勒搜索策略。适合谁正在写捕获算法、做信号级仿真、给学生搭GNSS实验平台的人。反直觉的是生成PRN码本身不是查表而是用移位寄存器抽头组合出来的二级码也不是另一个PRN是叠加在导频通道上的二次编码。等把这三层拆开你会发现做出正确频谱的关键根本不在FFT而在载波频率和采样率的关系上。2. 用Matlab复现GPS C/A码移位寄存器、抽头表与最小生成函数2.1 C/A码的生成原理为什么是G1和G2两个寄存器GPS L1 C/A码的本质是Gold码序列码长1023码率1.023 Mcps周期1 ms。生成结构是固定的一个G1寄存器生成多项式 (1x^3x^{10})一个G2寄存器生成多项式 (1x^2x^3x^6x^8x^9x^{10})。两颗卫星的区别不在于G1或G2的反馈多项式而在于从G2寄存器哪些抽头位置引出信号与G1输出做模二加。抽头组合定了PRN号就定了。所以生成流程分三步初始化两个10级寄存器为全1按各自多项式做反馈移位把G2的指定抽头异或后与G1异或输出。这里有个新手容易搞反的点反馈多项式里最高次项对应的是移位寄存器最右侧那一级不是最左侧。我见过有人把抽头序号按从左往右数生成的码序列和真实C/A码长得完全不一样。下面是最小可运行的Matlab生成函数。2.2 最小生成函数一段可以直接保存的代码function ca generateCAcode(prn) % 生成GPS L1 C/A码返回长度1023的±1序列 % 输入prn卫星号(1~32)输出ca1023x1的double元素为1或-1 % G2抽头表第1列为主抽头第2列为辅抽头定义见IS-GPS-200 tapTable [2,6; 3,7; 4,8; 5,9; 1,9; 2,10; 1,8; 2,9; 3,10; 2,3; 3,4; 5,6; 6,7; 7,8; 8,9; 9,10; 1,4; 2,5; 3,6; 4,7; 5,8; 6,9; 1,3; 4,6; 5,7; 6,8; 7,9; 8,10; 1,6; 2,7; 3,8; 4,9]; if prn 1 || prn 32 error(PRN必须在1到32之间); end % 寄存器初始化全1这是Gold码的标准起始状态 g1 ones(1, 10); g2 ones(1, 10); ca zeros(1, 1023); tap1 tapTable(prn, 1); tap2 tapTable(prn, 2); for i 1:1023 % 输出是g1和g2抽头异或后再异或 ca(i) xor(g1(10), xor(g2(tap1), g2(tap2))); % G1反馈第3级和第10级异或 g1Feedback xor(g1(3), g1(10)); % G2反馈第2、3、6、8、9、10级异或 g2Feedback xor(g2(2), xor(g2(3), xor(g2(6), xor(g2(8), xor(g2(9), g2(10)))))); % 右移一位把反馈值放到最左侧 g1 [g1Feedback, g1(1:9)]; g2 [g2Feedback, g2(1:9)]; end % 转成±1逻辑1映射为1逻辑0映射为-1 ca 2 * ca - 1; end代码逻辑不复杂但有几处必须强调。抽头表里用的是G2寄存器的级号级号从左侧输入端开始数第1级是刚移入的反馈位。循环里每轮先取当前状态输出一个码片再移位顺序不能反反了会让整个序列错位一个码片。反馈多项式的实现方式是把对应级数异或起来G2的反馈项比较多我拆成了嵌套xor如果你觉得可读性差可以先用中间变量累加异或再赋值。调用方式也简单ca generateCAcode(3)就得到PRN 3的本地码。要验证生成对不对一个土办法是数正负号的游程长度C/A码里最长连“1”是17个码片最长连“0”是18个如果这个特征不对寄存器级数或反馈多项式基本有问题。2.3 验证生成结果用自相关和互相关确认码族正交性生成完不能直接拿去用要先验证码的性质。Gold码族最核心的特性是互相关旁瓣被抑制在-65 dB以下自相关峰是尖锐的单峰。验证方法是用FFT做循环互相关避免直接写双重循环1023点用FFT几百微秒就出来了。% 验证脚本检查PRN 1到5的自相关峰和互相关旁瓣 fsCode 1.023e6; % 码率这里只用于时间轴标注 prnList 1:5; caCodes zeros(1023, length(prnList)); for idx 1:length(prnList) caCodes(:, idx) generateCAcode(prnList(idx)); end % 循环互相关用FFT实现等效于圆周移位相关 fftLen 2048; % 补零到2倍长度避免圆周混叠干扰 corrMat zeros(1023, length(prnList)); for idx 1:length(prnList) ref caCodes(:, idx); corrSeq ifft(fft(ref, fftLen) .* conj(fft(caCodes(:, 1), fftLen))); corrMat(:, idx) abs(corrSeq(1:1023)); end % 输出互相关旁瓣峰值单位dBc for idx 2:length(prnList) sidelobe max(corrMat(2:end, idx)) / max(corrMat(:, idx)); fprintf(PRN %d vs PRN 1 互相关旁瓣峰值: %.2f dB\n, prnList(idx), 20*log10(sidelobe)); endFFT相关运算里fftLen取2048而不是1023是为了让圆周卷积退化成线性卷积避免码首尾互相污染。取相关序列前1023个点看主瓣和旁瓣主瓣位置应该在第一个点。这里有个判断标准旁瓣峰值低于主瓣24 dB以上说明码族性质基本正确更好的理论值在-65 dB左右但短序列加非整周期观测会抬噪底不必死磕。如果你发现某些PRN之间的旁瓣明显高于其他组合先查抽头表复制粘贴时有没有错行。2.4 关于通信工具箱的现成对象什么时候不用自己写Matlab通信工具箱里有comm.GoldSequence对象可以直接配置两路移位寄存器的反馈多项式和抽头索引。它适合快速原型验证不需要自己维护抽头表。但我个人的习惯是保留自己的生成函数原因有两个一是接收机仿真里经常要控制码相位偏移、码率误差、多普勒频移自己写的函数改动成本低二是工具箱对象输出的逻辑电平格式和后续基带调制代码有时存在类型不匹配来回转换反而容易出bug。如果你打算用工具箱对象配置思路是GoldSequence的FirstPolynomial设为x^10x^31SecondPolynomial设为x^10x^9x^8x^6x^3x^21再用PrivateSequence设置两颗卫星的抽头。注意这里多项式字符串的写法是高位在左和前面移位寄存器代码里从右往左数的级号是反的别对着IS-GPS-200的抽头表直接填容易差一位。3. 二级码是什么、无数据信号怎么合成把导航电文和导频信号分开看3.1 二级码的两层含义数据比特流和overlay code标题里的“二级代码”在GNSS语境下有两种常见指向。第一种是导航电文的数据比特流GPS L1 C/A的导航电文速率是50 bps一个比特持续20 ms对应20个C/A码周期。这是信号体制里的“第二层编码”扩频码码片是第一层数据比特是叠加其上的第二层。第二种指向是新一代信号比如GPS L1C、北斗B1C的导频通道overlay code它是在扩频码之上再调制的一段低速码序列用于辅助跟踪和数据通道分离。如果做的是L1 C/A信号模拟二级代码通常指导航电文比特。无数据信号对应的就是“导频信号”的概念——不调制任何导航电文数据比特恒为1。这种信号常用于接收机捕获阶段的验证因为没有比特翻转相关积分不会跨比特边界损失能量。3.2 生成导航电文比特流帧结构、奇偶校验与速率匹配生成“无数据信号”前先把带数据的信号做出来数据模式可以自定义。最简单的做法是生成50 bps的随机比特映射成±1然后每个比特重复20次去匹配1 ms的C/A码周期。% 生成1秒长度的GPS基带信号含导航数据 C/A码 numMs 1000; % 仿真时长1000ms 1s codeLen 1023; % 每周期码片数 fsCode 1.023e6; % 码片速率1.023Mcps samplesPerChip 2; % 每个码片采2个点采样率2.046MHz fs fsCode * samplesPerChip; % 取PRN 1的C/A码 caCode generateCAcode(1); caCodeUp repelem(caCode(:), samplesPerChip); % 上采样到采样率 % 导航电文这里用随机比特模拟实际可用真帧头替换 navBits randi([0 1], numMs/20, 1); navSymbols 2 * navBits - 1; % 0→-1, 1→1 navUp repelem(navSymbols, 20 * codeLen * samplesPerChip); % 基带扩频数据比特乘以C/A码 baseband navUp .* repmat(caCodeUp, numMs, 1);导航比特重复20次的逻辑是1个数据比特 20 ms 20个C/A码周期 20×102320460个码片。上采样后每个码片2个点所以每个数据比特对应40920个采样点。如果你换采样率比如改成8.184 MHz每码片8个点那navUp的重复次数就要相应改成20 * codeLen * 8。这里最常见的错误是只改采样率忘了同步改重复因子结果数据比特速率对不上。需要说明的是真实GPS导航电文有帧同步头比如遥测字的前8位是10001011和奇偶校验码汉明码32位中30位是数据、2位是校验如果做的是接收机算法验证建议把电文按真实格式生成否则跟踪环路里比特同步那一步会找不到帧边界。但单纯做频谱分析和捕获验证随机比特就够了。3.3 无数据信号合成把数据位焊死在1专注扩频码和载波无数据信号在Matlab里生成反而更简单因为少了一层调制。直接把扩频码重复N个周期不做数据比特映射。它的价值在于频谱测量没有数据比特翻转时频谱主瓣形状更规整旁瓣结构清晰可辨适合验证扩频码本身的频谱特征。% 无数据基带信号纯C/A码不含导航电文 numSpreads 100; % 重复100个码周期 pureCode repmat(caCodeUp, numSpreads, 1); % 加上载波变成中频信号 fc 4.092e6; % 中频4.092MHz常见接收机中频 t (0:length(pureCode)-1) / fs; ifSignal pureCode .* cos(2 * pi * fc * t);载波频率的选择会影响后续频谱图的观感。4.092 MHz是1.023 MHz的4倍与1码片8采样点的配置刚好整数倍关系FFT观测时不会出现频谱泄漏带来的“假峰”。如果你把中频选成3.5 MHz这种与码率非整数倍的值频谱图里的谱线会呈现梳状分布那不是信号有问题是采样率与载波频率不匹配导致的栅栏效应。这个点在第4章避坑里再展开。无数据信号还有个用途是测相位一致性。用纯导频信号做载波相位测量相位误差的标准差比带数据信号低一个量级因为消除了比特翻转对相位估计的扰动。如果你在写高精度测距算法这一步值得单独验证。4. 频谱分析验证信号质量从时域波形看到1.023 MHz主瓣宽度4.1 理想频谱长什么样sinc包络、主瓣宽度与零点位置扩频码的功率谱密度是sinc函数的平方包络主瓣宽度由码率决定第一零点在码率的整数倍处。对C/A码码率1.023 MHz所以主瓣零点在±1.023 MHz。副瓣峰值比主瓣低约13.26 dB这是sinc函数的基本性质。判断一个生成的GPS信号对不对第一条就是看主瓣零点位置是否精确落在±1.023 MHz上。如果主瓣宽度不对几乎可以断定码率设置错误。比如你意外把码率设成了1.25 Mbps主瓣零点就跑到了±1.25 MHz。这种错误在时域上几乎看不出差别但在频谱上一眼就能抓出来。4.2 用pwelch和FFT两种方式拉频谱参数选择与归一化% 方法一周期图法适合快速看一眼 N length(pureCode); spectrum fftshift(fft(pureCode .* hann(N))); freqAxis (-N/2 : N/2-1) / N * fs; plot(freqAxis / 1e6, 20*log10(abs(spectrum) / max(abs(spectrum)))); xlim([-3 3]); ylim([-60 0]); xlabel(频率 (MHz)); ylabel(归一化幅度 (dB)); title(无数据C/A码基带信号频谱周期图法); grid on; % 方法二平均周期图适合细看旁瓣结构 windowLen 2048; overlap 512; nfft 2048; [p, f] pwelch(pureCode, windowLen, overlap, nfft, fs); plot(f/1e6, 10*log10(p / max(p))); xlim([-3 3]); ylim([-60 0]); xlabel(频率 (MHz)); ylabel(相对功率谱密度 (dB)); grid on;两种方法的取舍周期图法频率分辨率高但方差大谱线毛刺明显pwelch平均次数多曲线平滑但频率分辨率下降。看主瓣零点位置用周期图法更精确看旁瓣整体趋势用pwelch更清楚。hann窗会轻微展宽主瓣零点位置移动不超过一个频率分辨率单元可接受。参数里有几个细节值得说。hann窗的长度N如果正好等于信号长度那频域主瓣宽度会略微偏大因为窗函数本身的频谱和信号频谱卷积了。想得到锐利主瓣可以让信号长度是窗长的整数倍或者干脆不加窗直接做FFT。pwelch的windowLen取2048对应的时间是1 ms的整数倍2048点在2.046 MHz采样率下恰好是1 ms这样每个窗段都包含至少一个完整C/A码周期频谱估计不会因为截断不完整周期而失真。4.3 频谱异常怎么定位多根谱线、鼓包和偏移实际生成中常见的频谱异常有三种。第一种是主瓣零点位置偏离±1.023 MHz原因多半是采样率设错或者码率与采样率的关系换算错误。第二种是主瓣内有周期性起伏呈现“波纹状”这通常是数据比特翻转造成的——数据比特率的1/T_b50 Hz调制在信号上频谱上会产生离散谱线。第三种是频谱整体平移比如基带信号显示中心频率不在0 Hz那是因为混入了残余载波分量多半是cos调制时初始相位不为0导致载波泄漏。% 载波泄漏检查计算零频附近的功率占比 carrierLeak sum(abs(spectrum(floor(N/2)-10 : floor(N/2)10)).^2); totalPower sum(abs(spectrum).^2); fprintf(载波泄漏功率占比: %.2f dB\n, 10*log10(carrierLeak/totalPower));这个检查脚本会告诉你残留载波分量有多强。正常无数据信号应该完全抑制载波占比低于-40 dB。如果高于-20 dB回去检查cos的初始相位和采样时间向量是否对齐常见错误是t从1开始而不是从0开始导致相位偏了2πfc/fs。5. 生成PRN码的避坑记录五个最容易“翻车”的细节5.1 抽头表乱序导致PRN错号现象生成的PRN 1的码序列在相关器里和真实PRN 1对不上但却能和PRN 3的捕获峰对上号。原因IS-GPS-200里的抽头表是按卫星编号排列的第1行是PRN 1的抽头第2行是PRN 2的抽头。有人从网上下载的代码里抽头表顺序错位或者编程语言数组下标从0开始时没做偏移导致取出来的是相邻卫星的抽头组合。解决生成后立刻做自检——打印G2抽头索引和标准表人工核对一遍。同时用一个已知的PRN码序列做相关验证比如用PRN 1的码和信号模拟器输出的PRN 1码做滑动相关相关峰超过25 dB才算通过。5.2 寄存器初值不是全1现象生成的自相关峰不是单一尖峰而是有两个等高的旁瓣且最大旁瓣只比主瓣低10 dB左右。原因Gold码的初值约定是G1和G2都置全1。如果初值是全0移位寄存器就永远输不出序列这一步会直接报错但如果初值是其他随机值生成的是同一个Gold码族的循环移位版本自相关旁瓣结构被破坏。解决初始化一律用ones(1,10)别用zeros也别用randi。如果要生成不同码相位的C/A码正确做法是初始状态不变把输出序列做循环移位而不是改寄存器初值。5.3 采样率不是码率整数倍导致的频谱“假峰”现象频谱主瓣零点位置看着对但±1.023 MHz之外多出几根高谱线且位置随FFT点数变化漂移。原因采样率如果对码率不是整数倍关系码片边沿的采样位置在相邻码片里不一致等效于给信号叠了一层周期性的采样抖动频谱上表现为额外谱线。解决把采样率设置为码率的整数倍最常用的是2、4、8倍。如果实际接收机的中频采样率确实不是整数倍在仿真里宁可先加一个重采样步骤也不要直接用非整数倍采样率去模拟。5.4 代码长度截断导致相关峰展宽现象捕获阶段相关峰不是1个采样点的尖峰而是出现三角形平台峰值宽度横跨好几个采样点。原因C/A码周期是1023个码片但你在生成基带信号时如果按repmat(caCodeUp, numSpreads, 1)而caCodeUp长度被repelem扩成了2046个采样点信号总长度对应的码周期数不是整数时最后的截断破坏了一个完整码周期的自相关特性。解决repmat的次数要和码周期精确对齐。比如每码片2个采样点一个码周期是2046个采样点总采样长度必须是2046的整数倍。如果仿真时长需要1000 ms而1 ms一个周期那长度就是2046×1000不多不少。5.5 pwelch窗口长度与码周期错位现象pwelch画出的频谱主瓣边缘模糊旁瓣底部“鼓包”明显但周期图法画出来是正常的。原因pwelch的每个窗段如果截断在码周期的中间段内信号不是整周期的循环平稳过程频谱估计里混入了截断效应。尤其当windowLen不是码周期的整数倍时这种影响持续存在。解决把windowLen设成码周期采样长度的整数倍。在2.046 MHz采样率下码周期长度2046所以windowLen取4092、8184或2046都行。我习惯取4092这样每段包含2个完整码周期频率分辨率足够看主瓣细节。6. 用相关捕获验证生成的GPS信号让整条链路自证清白6.1 一个捕获验证脚本把PRN、基带和频谱串起来生成PRN码和基带信号之后最后一步验证是模拟一次粗捕获把生成的信号当作“接收信号”用本地生成的PRN码做二维搜索码相位×多普勒频移确认在正确的码相位上出现尖锐相关峰。这一步能同时验证PRN码生成的正确性、扩频调制的正确性以及信号幅度的合理性。% 捕获验证对无数据基带信号做码相位搜索 rxSignal pureCode(1:fs*0.01); % 取10ms信号 localCode generateCAcode(1); localUp repelem(localCode, samplesPerChip); % 用FFT实现循环相关在码相位域搜索 corrLen length(rxSignal); corrSeq ifft(fft(rxSignal) .* conj(fft(localUp, corrLen))); corrPower abs(corrSeq).^2; % 归一化并输出峰值位置 [peakVal, peakIdx] max(corrPower); fprintf(捕获峰值位置: %d 采样点, 归一化峰值: %.2f\n, ... peakIdx, 10*log10(peakVal / mean(corrPower)));如果一切正确峰值位置应该在第一个采样点附近因为信号和本地码从同一相位开始归一化峰均比在25 dB以上。如果峰值位置偏移了几个采样点说明repelem和repmat的维度对错了。还有个值得一试的测试把接收信号换成PRN 2生成的信号本地码仍用PRN 1捕获峰会消失这个负向验证能确认码族隔离度是真的。6.2 参数选择建议采样率、信号长度和频偏范围做接收机捕获算法仿真时我常用这样的参数组合采样率4.092 MHz码率4倍信号长度10 ms多普勒搜索范围±5 kHz步进500 Hz。这个组合在普通笔记本上跑完一次全二维搜索大约需要几十秒适合反复调参。如果只验证PRN生成信号1 ms就够如果要测捕获灵敏度至少10 ms如果要测比特同步需要1秒以上包含完整导航电文的数据。我在实际做信号模拟器时踩过的最深一个坑是把所有代码都调通之后才发现中频频率和采样率的关系没留裕量导致频谱图里镜像分量和主瓣重叠。后来养成了一个习惯每写一段信号生成代码先跑频谱看主瓣零点再跑一次捕获看相关峰两个检查过了才继续往下走。这个方法帮我在后续的多径仿真里省了大量排查时间。希望帮到你。本文还有配套的精品资源点击获取