简介:这是一套面向 GNSS 信号处理与 Matlab 开发者的代码生成工具,能够生成 GPS(L1CA/L2CM/L2CL/L5I/L5Q)、Galileo(E1B/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寄存器,生成多项式 (1+x^3+x^{10});一个G2寄存器,生成多项式 (1+x^2+x^3+x^6+x^8+x^9+x^{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),输出ca:1023x1的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^10+x^3+1',SecondPolynomial设为'x^10+x^9+x^8+x^6+x^3+x^2+1',再用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×1023=20460个码片。上采样后每个码片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_b=50 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秒以上包含完整导航电文的数据。
我在实际做信号模拟器时踩过的最深一个坑,是把所有代码都调通之后才发现中频频率和采样率的关系没留裕量,导致频谱图里镜像分量和主瓣重叠。后来养成了一个习惯:每写一段信号生成代码,先跑频谱看主瓣零点,再跑一次捕获看相关峰,两个检查过了才继续往下走。这个方法帮我在后续的多径仿真里省了大量排查时间。希望帮到你。
本文还有配套的精品资源,点击获取