简介本资源是一套面向卫星导航系统学习者与MATLAB信号处理初学者的北斗基带信号仿真工程聚焦BDS B1频段C/A码生成与数字调制建模解决导航信号原理验证、基带波形生成及软件无线电前端设计等实践问题。压缩包共8个文件含6个核心MATLAB函数.m、1份说明文档.md和1个文本说明.txt总大小仅9KB轻量易部署其中主函数main.m完整实现10秒内含导航数据±1符号、20ms/bit、载波调制与C/A码扩频的基带信号合成配套generateCACode、digitalCA等模块支持码序列生成、采样量化与卫星编号映射。已有142人学习下载代码结构清晰、注释完备提供从C/A码产生→数字化查表→基带信号合成的全流程可运行脚本适合作为导航原理课程实验、GNSS接收机算法入门或MATLAB通信系统仿真实践的可靠参考。 做北斗基带信号MATLAB仿真这件事我从一开始的“拿GPS代码改改参数”到后来自己重新写了整套流程中间绕了不少弯路。很多人以为北斗B1I和GPS L1差别不大直接换频率、换码率就行真跑起来才发现PRN码生成逻辑、NH码二次调制、码长截断这些细节全是坑。这篇就把我完整跑通的一套北斗基带信号仿真代码讲清楚从协议参数到PRN码生成、基带调制、捕获验证每一步都给出可复现的MATLAB实现和我在实际调试中的经验教训适合正在做卫星导航基带信号仿真、软件接收机开发或者刚开始接触北斗信号处理的朋友参考。1. 北斗B1I信号到底长什么样先把仿真的输入参数敲定1.1 从协议层看B1I信号结构北斗B1I信号是北斗二号和北斗三号系统在B1频点播发的公开服务信号标称载波频率1561.098 MHz这个频率本身就和GPS的1575.42 MHz不一样所以做仿真时本地载波NCO的频率字、混频链路的参数都要重新设计直接拿GPS那套移植过来肯定会偏。B1I信号的码速率为2.046 Mcps码长为2046个码片也就是说一个完整码周期刚好是1 ms这个整毫秒结构对后续捕获、跟踪的积分时间设计非常友好。B1I信号最容易被忽略的是它采用了二级码调制也就是在PRN码之上还叠加了一个NH码。NH码的全称是Neumann-Hoffman码B1I使用的NH码长度为20 bit码片速率1 kbps每个NH码片持续1 ms正好对应一个PRN码周期。NH码的作用是给导航电文提供比特同步和子帧同步的辅助信息同时对信号的频谱进行一定程度的扩展。仿真如果漏掉NH码这一层生成出来的信号在捕获阶段可能看不出大问题但一旦进入比特同步或者导航电文解调就会完全对不上时间关系。信号的大致组成可以写成s(t) A × C(t) × D(t) × cos(2πf_c t φ₀)其中A是信号幅度C(t)是扩频码序列包含PRN码和NH码的乘积关系D(t)是导航电文比特流f_c是载波频率1561.098 MHzφ₀是初始载波相位。仿真是要把这四个部分全部搭起来然后通过上采样、混频输出中频采样信号最后丢给捕获算法去验证。1.2 仿真链路设计的总体思路我把整个仿真链路拆成几层来设计。最底层是PRN码生成模块负责产生指定卫星号的2046码片扩频码第二层是NH码与导航电文的组合模块把20 ms的导航比特与NH码进行异或调制再调制到PRN码上第三层是信号生成模块完成从基带序列到采样信号的转换包括码片上采样、脉冲成形、载波混频和加噪第四层是接收端验证模块用捕获算法去确认生成信号的码相位和多普勒频率是可检测的。这套分层设计的好处在于每一层都可以独立测试。我在实际调试中经常遇到的问题是信号生成出来以后捕获不到峰值这时如果是整套代码一起写根本不知道是码生成错误还是调制链路错误分层以后就能顺着链路逐级检查。下面我从PRN码开始讲这是整个仿真最基础也最容易出错的地方。2. 北斗PRN码生成从移位寄存器到2046码片的完整实现2.1 G1/G2移位寄存器原理北斗B1I的PRN码是基于两个11级线性反馈移位寄存器G1和G2生成的。两个寄存器不仅反馈多项式不同而且G2输出还需要经过相位抽头选择后与G1输出做异或得到最终码。G1寄存器的反馈多项式在B1I协议里是固定的简单说是11级全部参与反馈形成的序列本质是一个周期为2047的m序列。G2寄存器同样也是11级但它的输出是指定两个抽头位置做异或以后的结果不同卫星号对应不同的抽头选择方式。这里有个非常关键的细节2047是m序列的自然周期但B1I的PRN码长度是2046比2047少了一个码片。这是因为G1寄存器在产生到第2046个码片后会被强制复位然后重新开始下一轮周期。这个“截断一位”的操作在协议里是有明确规定的如果仿真时忽略这一点出来的PRN码会比真实信号多一位导致整个码周期长度变成2047/2046倍的关系与1 ms对齐结构错位。给大家一个直观类比2047就像一圈楼梯有2047级台阶但北斗要求每一秒钟的采样点必须和整毫秒对齐所以拿走一级台阶让一圈变成2046级这样每圈刚好对应1 ms后续信号处理的时间关系才能保证。2.2 MATLAB代码实现我写了一个函数来生成指定卫星号的B1I PRN码function prn generateB1IPRN(prnNum) % 生成北斗B1I PRN码输出长度2046 % prnNum: 卫星PRN号 1~37 % 寄存器初始化全1 g1 ones(1, 11); g2 ones(1, 11); % 相位抽头查表这里仅示意 % 实际需要完整的抽头对照表格式为[抽头1, 抽头2] tapTable [ 2, 7; % PRN 1 3, 4; % PRN 2 % ... 继续补全到37号 ]; taps tapTable(prnNum, :); prn zeros(1, 2046); for idx 1:2046 % 当前输出码片 g1out g1(end); g2out xor(g2(taps(1)), g2(taps(2))); prn(idx) xor(g1out, g2out); % G1反馈所有位参与异或 fb1 mod(sum(g1), 2); % 等价于异或所有位 g1 [fb1, g1(1:10)]; % G2反馈所有位参与异或 fb2 mod(sum(g2), 2); g2 [fb2, g2(1:10)]; end % 最后一个状态强制复位由下一次循环自然完成 % 每周期正好输出2046码片 end关于代码有几点说明。第一G1反馈我直接用了mod(sum(g1), 2)这样写比连续嵌套xor要清晰得多也方便日后改成其他多项式结构。第二G2的抽头选择不同卫星不一样示例代码里的表格示意了两个PRN号对应的抽头组合完整的相位抽头关系表在我的实现中是单独保存在一个函数里实际使用时应该按照完整协议表补齐。第三寄存器的移位方向是左移还是右移本身不影响码序列的性质只要输出抽头选取与移位方向一致即可。我看到一些文献里用右移写法这没问题关键是确保输出码与标准PRN码一致完成后用相关峰值验证即可。2.3 怎么验证生成的码是对的生成码以后不验证就直接往后面链路接是我见过的最多的错误做法。验证方法其实很简单同两个PRN码做自相关和互相关分析。% 生成PRN 1和PRN 2 prn1 generateB1IPRN(1); prn2 generateB1IPRN(2); prn1(prn1 0) -1; prn2(prn2 0) -1; % 自相关 autoCorr xcorr(prn1, prn1, 1000); figure; plot(-1000:1000, autoCorr); title(PRN1 自相关); % 互相关 crossCorr xcorr(prn1, prn2, 1000); figure; plot(-1000:1000, crossCorr); title(PRN1-PRN2 互相关);正常情况下的自相关应该是一个在零延迟处出现明显尖峰、其他地方接近零的波形互相关则整体接近零没有明显尖峰。我第一次跑的时候自相关出现了多个副峰后来排查发现是G1反馈多项式写错了把所有位异或改成了只有最后一位反馈这种错误从波形上一下子就能看出来。所以验证这一步一定不能跳过它虽然花不了多少时间却能避免后续所有问题都堆在一起无从下手。另外还要注意码的极性映射。PRN码生成出来是0/1序列但基带信号处理当中要用双极性±1来表示0映射为1、1映射为-1或者反过来也可以关键是发射端和接收端要一致。我在代码里统一用0转1、1转-1的映射方式捕获相关计算基于这个约定来写避免符号翻转引起的相关性反相问题。3. 基带调制链路NH码、导航电文与上采样3.1 NH码调制B1I的特殊之处B1I的NH码长度是20 bit对应的码序列在不同文档里写法略有差异我按自己项目中的实现作为示例。NH码的作用是把50 bps的导航电文扩成1 kbps的数据流也就是每个50 bps比特被20个NH码片调制每个NH码片持续1个PRN码周期。在实现中我先把NH码转成双极性序列然后与导航电文比特逐个异或等价于相乘再复制到PRN码序列上% NH码B1I的20bit二次码按实际协议填入 nhCode [0 0 0 0 0 1 0 1 1 0 1 0 0 1 1 1 0 0 0 1]; nhCode(nhCode 0) -1; % 导航电文比特这里以随机序列模拟 navBits randi([0 1], 1, 10); navBits(navBits 0) -1; % 生成一个码周期的PRN prn generateB1IPRN(1); prn(prn 0) -1; % 基带码片序列生成 baseband []; for bIdx 1:length(navBits) for nhIdx 1:20 % 当前码片块 PRN码 × NH码片 × 导航比特 chipBlock prn * nhCode(nhIdx) * navBits(bIdx); baseband [baseband, chipBlock]; end end这段代码运行后baseband的长度是10个导航比特 × 20个NH码片 × 2046个PRN码片也就是409200个码片对应时长200 ms。这里有一个性能隐患用[baseband, chipBlock]这种拼接方式在循环里会不断重建数组数据量大了以后非常慢。我实际处理更长的数据时改用了预分配内存的方式numBits 10; totalChips numBits * 20 * 2046; baseband zeros(1, totalChips); idxStart 1; for bIdx 1:numBits for nhIdx 1:20 block prn * nhCode(nhIdx) * navBits(bIdx); baseband(idxStart:idxStart2045) block; idxStart idxStart 2046; end end这一点对于初学者来说可能觉得无关紧要但如果你的仿真数据量到几十秒甚至上百秒循环里动态拼接会让代码跑到怀疑人生预分配后速度能提升好几个数量级。3.2 导航电文不再只有随机数上面用了随机数模拟导航电文实际项目中通常需要更贴近真实场景。北斗B1I的D1导航电文包含帧同步码、子帧计数、星历参数、电离层修正等结构如果是完整做接收机当然要按协议一层层拼。但如果做的是基带信号仿真和捕获验证用随机比特或固定的伪随机序列完全可以满足需求重点是把信号的时间增益和码结构做对。我在自己项目中用一个简单的函数生成固定模式的导航电文方便定位问题function navBits createNavBits(numBits) % 生成固定的导航比特便于仿真调试 pattern [1 0 1 1 0 0 1 0]; navBits repmat(pattern, 1, ceil(numBits/length(pattern))); navBits navBits(1:numBits); navBits(navBits 0) -1; end固定模式的好处是当你在接收端解调出数据以后一眼就能看出比特顺序对不对、有没有发生滑位。用纯随机序列当然也能验证但出了问题很难判断是解调错误还是本来就预期范围内的随机数据。3.3 从基带到中频采样率设计与混频实现基带码片序列要变成可处理的采样信号中间要经过两个核心步骤上采样和载波混频。这里最关键的参数是采样率它决定了整个仿真的计算量和模拟精度。我常用的做法是先确定每个码片的采样点数再反推采样率。B1I码速率2.046 Mcps如果每个码片取4个采样点采样率就是8.184 MHz。之所以不用整数倍的2.046 MHz而是取4倍采样是为了匹配后续捕获算法的FFT处理同时兼顾频谱观察的直观性。上采样直接用repelem实现它会把每个码片复制成指定数量的采样点sps 4; % 每个码片采样点数 fs 2.046e6 * sps; % 采样率 8.184 MHz sigBase repelem(baseband, sps); % 混频到中频 t (0:length(sigBase)-1) / fs; ifFreq 4.092e6; % 中频选为码速率2倍附近便于观察频谱 lo exp(1j * 2 * pi * ifFreq * t); sigIF real(sigBase .* lo);中频选择4.092 MHz也就是采样率的一半这样频谱上信号位于正半轴负半轴是对称镜像。要注意lo用了复数形式这样信号会被搬移到中频正频率处取实部以后得到实中频信号。如果你希望信号留在基带也可以不做混频直接用sigBase加噪声取决于你的验证目标。捕获算法通常在中频或者基带都能工作但混到中频更接近真实接收机的处理流程。加噪声时有个容易犯错的地方很多人直接用awgn函数随手加一个信噪比但awgn默认按信号功率计算而信噪比的定义要明确是码片信噪比还是采样点信噪比。我更习惯手动计算噪声功率% 按载噪比添加噪声C/N0单位dB-Hz CN0 45; % 典型GNSS载噪比 Tcoh 0.001; % 相干积分时间1ms signalPower mean(sigIF.^2); SNR CN0 - 10*log10(1/Tcoh); % 转换为1ms带宽内的信噪比(dB) noisePower signalPower / (10^(SNR/10)); noise sqrt(noisePower) * randn(size(sigIF)); rxSignal sigIF noise;我习惯把载噪比作为加噪的原始输入而不是直接写信噪比因为GNSS仿真中大家聊的都是C/N0这样你的代码参数和其他接收机软件接口对得上不会出现“我说的是45 dB你说的是15 dB”这种鸡同鸭讲的情况。4. 捕获验证如何用自己生成的信号自证清白4.1 基于FFT的并行码相位搜索信号生成了但不能自说自话就认为它是对的。最有效的自验证手段是让仿真信号通过一个捕获模块如果能捕获到正确的码相位和多普勒频率就说明信号链路是通的。我用的是基于FFT的并行码相位搜索算法这是目前软件接收机中最经典也最有效的捕获方案。核心思想是对接收信号做FFT、对本地PRN码做FFT再取共轭两者在频域相乘后做IFFT结果就是时域的循环相关。这样一次IFFT就能同时计算出所有码相位的相关值避免了逐码相位滑动的巨大计算量。对应的MATLAB实现如下function [doppler, codePhase, peakMetric] acquisition(x, prn, fs, sps, dopplerRange) % 基于FFT的并行码相位捕获 % x 接收中频信号 % prn 本地PRN码0/1序列 % fs 采样率 % sps 每码片采样点数 % dopplerRange搜索的多普勒频率范围 prnCode prn; prnCode(prnCode 0) -1; localCode repelem(prnCode, sps); % 1ms信号的采样点数 N length(localCode); peakMetric 0; doppler 0; codePhase 0; for fd dopplerRange % 补偿多普勒后的接收信号 t (0:N-1) / fs; xComp x(1:N) .* exp(-1j * 2 * pi * fd * t); % FFT相关 Xf fft(xComp); Cf conj(fft(localCode)); R ifft(Xf .* Cf); R abs(R); % 寻找当前多普勒下的相关峰 [peak, idx] max(R); % 记录最大相关峰 if peak peakMetric peakMetric peak; codePhase idx; doppler fd; end end % 计算噪声底 noiseFloor mean(R(R peakMetric * 0.5)); peakMetric peakMetric / noiseFloor; % 输出峰值与噪声底比值 end这里的多普勒搜索范围选择很关键。如果是静止场景搜索范围可以缩小到±5 kHz频率步进取500 Hz整个搜索不到几十次FFT速度非常快。如果模拟高动态场景多普勒范围扩展到±20 kHz步进可以相应增粗到1 kHz牺牲一点精度换取速度。要注意代码中计算噪声底时用了mean(R(R peakMetric * 0.5))这个方式这种简单阈值分离是工程上常用的做法如果相关峰太强会把整个均值抬高影响噪声底的估计准确性加一个峰值排除条件会更稳健。4.2 捕获结果判读捕获完成后结果可以通过三维图直观展示surf(dopplerRange, 1:N, RMatrix); view(0, 90); xlabel(多普勒频率 (Hz)); ylabel(码相位 (采样点));正常的捕获结果应该看到明显的“图钉”状峰值峰值所在位置的X坐标对应多普勒频率Y坐标对应码相位。如果出现峰值不明显或者出现多峰并列的情况通常有几个原因。第一是本地码和信号码相位差超过了一个码片这时候相关值本身就很低需要在更大的范围内搜索或者考虑使用多普勒补偿后再跑一次捕获。第二是多普勒频率没有对齐导致信号能量被分散到多个频点相关峰会变缓变矮。第三是码相位步进与采样点比例不对如果codePhase对应到芯片级别的换算出现偏差会看到峰值虽然存在但出现在意料之外的位置。我调试时注意到一个有意思的现象当生成的信号没有加NH码调制时捕获也能成功峰值甚至更高加了NH码后因为NH码片翻转的原因1 ms积分的能量有一半可能被抵消峰值会有所下降。这不是错误而是真实信号本来就要这样处理。如果要提升捕获灵敏度可以改用半码片长度的相干积分或者在捕获前先做NH码剥离这是另一个话题但在仿真验证时不要因为峰值下降了就怀疑代码出错。5. 仿真中容易踩的坑采样率、码相位与频谱泄漏5.1 采样率不是随便定的我知道很多初学者会习惯性地把采样率设成一个整数值比如10 MHz、20 MHz这样看起来整齐但和码速率2.046 Mcps的匹配关系就变得很别扭。如果采样率与码速率不是整数倍关系每个码片的采样点数就会因为四舍五入而抖动导致信号时间基准出现漂移。尤其在长数据仿真中这种漂移会不断累积最终使信号和本地码的时间对齐完全错乱。我强烈建议把采样率设成码速率的整数倍最常用的是4倍采样也就是8.184 MHz。这样每个码片严格对应4个采样点码相位计算时只需要做一个简单映射不需要处理小数采样偏移。5.2 码相位对齐最容易出错的“最后一米”捕获输出的是采样点索引这个索引要换算成码片偏移才能判断是否正确。以4倍采样为例采样点索引除以4取整才能得到码片偏移但这里的“取整”方式需要注意。我见过有人直接用floor或者round在边界情况下会引进一个码片左右的偏差。正确做法是先明确输出索引的计数起点。如果MATLAB的max返回的索引从1开始那么码片偏移应该是codePhaseSample idx - 1; % 转为0基 codePhaseChip codePhaseSample / sps; % 码片级偏移然后和构造信号时已知的码相位真值比较看是否一致。这个在校验闭环时非常关键我在自己的代码中会打印出真值和捕获值方便一眼看出误差是多少。5.3 频谱泄漏与做窗函数生成中频信号后我习惯先观察一下频谱再跑捕获。观察频谱时如果不加窗矩形截断带来的频谱泄漏会让谱峰变得很宽边瓣看起来像孪生峰造成误判。建议在频谱分析时加上汉宁窗或者布莱克曼窗win hann(length(sigIF)); sigWin sigIF .* win; freqAxis linspace(-fs/2, fs/2, length(sigWin)); plot(freqAxis, fftshift(abs(fft(sigWin))));做频谱分析时加窗做捕获相关时不要去加窗这是两码事。捕获算法处理的是真实信号加窗会破坏信号的恒包络特性反而降低相关峰性能。我在最初调试时有一次把加窗逻辑误用到了捕获路径结果相关峰从600掉到200多排查了很久才发现是这个低级错误。6. 从仿真到实用这套代码还能怎么往前延伸6.1 多通道与多普勒上面讲的是单星信号生成实际接收机面对的是多颗卫星叠加的混合信号。仿真多通道信号时每个通道生成独立的中频信号然后直接相加再加噪声。多普勒频率不同的卫星信号在频谱上占据不同位置捕获算法只要能扫完预设的频率范围就能把多个卫星分别捕获出来。实现多通道时要注意信号功率的分配。所有通道加上噪声后总功率会被推高如果每个通道都按单通道噪声功率加噪信噪比会被严重压低。我通常先把所有通道叠加再根据叠加后的总功率计算噪声功率保证最终载噪比符合设定值。6.2 对接软件接收机这套仿真生成的基带采样数据可以直接写入文件喂给标准的软件接收机处理。数据格式通常是int8或int16的IQ交织采样我按照实际接收机的数据格式定义输出文件% 将生成的信号写入文件格式为IQ交替的int16 sigI int16(real(rxSignal) * 32767); sigQ int16(imag(rxSignal) * 32767); sigIQ reshape([sigI; sigQ], [], 1); fid fopen(sim_bf_b1i.bin, wb); fwrite(fid, sigIQ, int16); fclose(fid);这样生成的数据文件就能被人家的接收机代码直接读取在验证接收机算法时非常方便。我在自己的项目里经常用MATLAB生成一组已知参数的“标准数据”然后交给接收机程序跑看它捕获到的多普勒、码相位这两个参数与设定值是否一致这样能快速定位接收机的问题而不需要拿着真实信号在户外反复测试。6.3 硬件在环之前的最后一道检查从纯仿真走向硬件在环或者实际射频测试之前建议在MATLAB里做一次完整的信号级验证。我关注的指标包括捕获峰值与噪声底的比值、跟踪环路的码相位误差曲线、载波多普勒的估计残差。这些仿真结果先要稳定再去接前端射频板卡否则硬件调试时信号有问题根本分不清是射频链路还是算法链路的问题。在具体的工程经验里我建议把PRN码生成函数单独抽出来做一个单元测试每次修改相关代码后重新跑一遍自相关和互相关验证这个习惯帮我省了无数排查时间。NH码的序列值在协议文档中可能由于版本差异略有不同务必以所使用的ICD文件为准不要照抄网络上某一篇博客的序列这也是我曾经踩过的坑。这套代码做完了以后如果再想往前扩展可以加入多径信道模型、电离层延迟模型或者改成北斗B1C或者B2a的现代化新信号体制。核心的仿真架构保持一致换掉扩频码和调制参数就可以复用这也是当初我花时间把每一层解耦的回报。最后再分享一个小技巧。整个仿真过程中我在关键节点都加了变量检查PRN码生成后检查长度和相关峰、NH码调制后检查20个PRN周期内是否有倒相导致的码片翻转、混频后检查频谱峰值是否落在设定的中频上。这三个检查点只要有一个没通过就停下来排查不要继续往后跑。这套检查习惯让我的仿真从“看起来对”变成了“每一步都能证明对”建议你也照着搭一遍。本文还有配套的精品资源点击获取
北斗B1I基带信号MATLAB仿真：PRN码生成与捕获验证完整实现
简介本资源是一套面向卫星导航系统学习者与MATLAB信号处理初学者的北斗基带信号仿真工程聚焦BDS B1频段C/A码生成与数字调制建模解决导航信号原理验证、基带波形生成及软件无线电前端设计等实践问题。压缩包共8个文件含6个核心MATLAB函数.m、1份说明文档.md和1个文本说明.txt总大小仅9KB轻量易部署其中主函数main.m完整实现10秒内含导航数据±1符号、20ms/bit、载波调制与C/A码扩频的基带信号合成配套generateCACode、digitalCA等模块支持码序列生成、采样量化与卫星编号映射。已有142人学习下载代码结构清晰、注释完备提供从C/A码产生→数字化查表→基带信号合成的全流程可运行脚本适合作为导航原理课程实验、GNSS接收机算法入门或MATLAB通信系统仿真实践的可靠参考。 做北斗基带信号MATLAB仿真这件事我从一开始的“拿GPS代码改改参数”到后来自己重新写了整套流程中间绕了不少弯路。很多人以为北斗B1I和GPS L1差别不大直接换频率、换码率就行真跑起来才发现PRN码生成逻辑、NH码二次调制、码长截断这些细节全是坑。这篇就把我完整跑通的一套北斗基带信号仿真代码讲清楚从协议参数到PRN码生成、基带调制、捕获验证每一步都给出可复现的MATLAB实现和我在实际调试中的经验教训适合正在做卫星导航基带信号仿真、软件接收机开发或者刚开始接触北斗信号处理的朋友参考。1. 北斗B1I信号到底长什么样先把仿真的输入参数敲定1.1 从协议层看B1I信号结构北斗B1I信号是北斗二号和北斗三号系统在B1频点播发的公开服务信号标称载波频率1561.098 MHz这个频率本身就和GPS的1575.42 MHz不一样所以做仿真时本地载波NCO的频率字、混频链路的参数都要重新设计直接拿GPS那套移植过来肯定会偏。B1I信号的码速率为2.046 Mcps码长为2046个码片也就是说一个完整码周期刚好是1 ms这个整毫秒结构对后续捕获、跟踪的积分时间设计非常友好。B1I信号最容易被忽略的是它采用了二级码调制也就是在PRN码之上还叠加了一个NH码。NH码的全称是Neumann-Hoffman码B1I使用的NH码长度为20 bit码片速率1 kbps每个NH码片持续1 ms正好对应一个PRN码周期。NH码的作用是给导航电文提供比特同步和子帧同步的辅助信息同时对信号的频谱进行一定程度的扩展。仿真如果漏掉NH码这一层生成出来的信号在捕获阶段可能看不出大问题但一旦进入比特同步或者导航电文解调就会完全对不上时间关系。信号的大致组成可以写成s(t) A × C(t) × D(t) × cos(2πf_c t φ₀)其中A是信号幅度C(t)是扩频码序列包含PRN码和NH码的乘积关系D(t)是导航电文比特流f_c是载波频率1561.098 MHzφ₀是初始载波相位。仿真是要把这四个部分全部搭起来然后通过上采样、混频输出中频采样信号最后丢给捕获算法去验证。1.2 仿真链路设计的总体思路我把整个仿真链路拆成几层来设计。最底层是PRN码生成模块负责产生指定卫星号的2046码片扩频码第二层是NH码与导航电文的组合模块把20 ms的导航比特与NH码进行异或调制再调制到PRN码上第三层是信号生成模块完成从基带序列到采样信号的转换包括码片上采样、脉冲成形、载波混频和加噪第四层是接收端验证模块用捕获算法去确认生成信号的码相位和多普勒频率是可检测的。这套分层设计的好处在于每一层都可以独立测试。我在实际调试中经常遇到的问题是信号生成出来以后捕获不到峰值这时如果是整套代码一起写根本不知道是码生成错误还是调制链路错误分层以后就能顺着链路逐级检查。下面我从PRN码开始讲这是整个仿真最基础也最容易出错的地方。2. 北斗PRN码生成从移位寄存器到2046码片的完整实现2.1 G1/G2移位寄存器原理北斗B1I的PRN码是基于两个11级线性反馈移位寄存器G1和G2生成的。两个寄存器不仅反馈多项式不同而且G2输出还需要经过相位抽头选择后与G1输出做异或得到最终码。G1寄存器的反馈多项式在B1I协议里是固定的简单说是11级全部参与反馈形成的序列本质是一个周期为2047的m序列。G2寄存器同样也是11级但它的输出是指定两个抽头位置做异或以后的结果不同卫星号对应不同的抽头选择方式。这里有个非常关键的细节2047是m序列的自然周期但B1I的PRN码长度是2046比2047少了一个码片。这是因为G1寄存器在产生到第2046个码片后会被强制复位然后重新开始下一轮周期。这个“截断一位”的操作在协议里是有明确规定的如果仿真时忽略这一点出来的PRN码会比真实信号多一位导致整个码周期长度变成2047/2046倍的关系与1 ms对齐结构错位。给大家一个直观类比2047就像一圈楼梯有2047级台阶但北斗要求每一秒钟的采样点必须和整毫秒对齐所以拿走一级台阶让一圈变成2046级这样每圈刚好对应1 ms后续信号处理的时间关系才能保证。2.2 MATLAB代码实现我写了一个函数来生成指定卫星号的B1I PRN码function prn generateB1IPRN(prnNum) % 生成北斗B1I PRN码输出长度2046 % prnNum: 卫星PRN号 1~37 % 寄存器初始化全1 g1 ones(1, 11); g2 ones(1, 11); % 相位抽头查表这里仅示意 % 实际需要完整的抽头对照表格式为[抽头1, 抽头2] tapTable [ 2, 7; % PRN 1 3, 4; % PRN 2 % ... 继续补全到37号 ]; taps tapTable(prnNum, :); prn zeros(1, 2046); for idx 1:2046 % 当前输出码片 g1out g1(end); g2out xor(g2(taps(1)), g2(taps(2))); prn(idx) xor(g1out, g2out); % G1反馈所有位参与异或 fb1 mod(sum(g1), 2); % 等价于异或所有位 g1 [fb1, g1(1:10)]; % G2反馈所有位参与异或 fb2 mod(sum(g2), 2); g2 [fb2, g2(1:10)]; end % 最后一个状态强制复位由下一次循环自然完成 % 每周期正好输出2046码片 end关于代码有几点说明。第一G1反馈我直接用了mod(sum(g1), 2)这样写比连续嵌套xor要清晰得多也方便日后改成其他多项式结构。第二G2的抽头选择不同卫星不一样示例代码里的表格示意了两个PRN号对应的抽头组合完整的相位抽头关系表在我的实现中是单独保存在一个函数里实际使用时应该按照完整协议表补齐。第三寄存器的移位方向是左移还是右移本身不影响码序列的性质只要输出抽头选取与移位方向一致即可。我看到一些文献里用右移写法这没问题关键是确保输出码与标准PRN码一致完成后用相关峰值验证即可。2.3 怎么验证生成的码是对的生成码以后不验证就直接往后面链路接是我见过的最多的错误做法。验证方法其实很简单同两个PRN码做自相关和互相关分析。% 生成PRN 1和PRN 2 prn1 generateB1IPRN(1); prn2 generateB1IPRN(2); prn1(prn1 0) -1; prn2(prn2 0) -1; % 自相关 autoCorr xcorr(prn1, prn1, 1000); figure; plot(-1000:1000, autoCorr); title(PRN1 自相关); % 互相关 crossCorr xcorr(prn1, prn2, 1000); figure; plot(-1000:1000, crossCorr); title(PRN1-PRN2 互相关);正常情况下的自相关应该是一个在零延迟处出现明显尖峰、其他地方接近零的波形互相关则整体接近零没有明显尖峰。我第一次跑的时候自相关出现了多个副峰后来排查发现是G1反馈多项式写错了把所有位异或改成了只有最后一位反馈这种错误从波形上一下子就能看出来。所以验证这一步一定不能跳过它虽然花不了多少时间却能避免后续所有问题都堆在一起无从下手。另外还要注意码的极性映射。PRN码生成出来是0/1序列但基带信号处理当中要用双极性±1来表示0映射为1、1映射为-1或者反过来也可以关键是发射端和接收端要一致。我在代码里统一用0转1、1转-1的映射方式捕获相关计算基于这个约定来写避免符号翻转引起的相关性反相问题。3. 基带调制链路NH码、导航电文与上采样3.1 NH码调制B1I的特殊之处B1I的NH码长度是20 bit对应的码序列在不同文档里写法略有差异我按自己项目中的实现作为示例。NH码的作用是把50 bps的导航电文扩成1 kbps的数据流也就是每个50 bps比特被20个NH码片调制每个NH码片持续1个PRN码周期。在实现中我先把NH码转成双极性序列然后与导航电文比特逐个异或等价于相乘再复制到PRN码序列上% NH码B1I的20bit二次码按实际协议填入 nhCode [0 0 0 0 0 1 0 1 1 0 1 0 0 1 1 1 0 0 0 1]; nhCode(nhCode 0) -1; % 导航电文比特这里以随机序列模拟 navBits randi([0 1], 1, 10); navBits(navBits 0) -1; % 生成一个码周期的PRN prn generateB1IPRN(1); prn(prn 0) -1; % 基带码片序列生成 baseband []; for bIdx 1:length(navBits) for nhIdx 1:20 % 当前码片块 PRN码 × NH码片 × 导航比特 chipBlock prn * nhCode(nhIdx) * navBits(bIdx); baseband [baseband, chipBlock]; end end这段代码运行后baseband的长度是10个导航比特 × 20个NH码片 × 2046个PRN码片也就是409200个码片对应时长200 ms。这里有一个性能隐患用[baseband, chipBlock]这种拼接方式在循环里会不断重建数组数据量大了以后非常慢。我实际处理更长的数据时改用了预分配内存的方式numBits 10; totalChips numBits * 20 * 2046; baseband zeros(1, totalChips); idxStart 1; for bIdx 1:numBits for nhIdx 1:20 block prn * nhCode(nhIdx) * navBits(bIdx); baseband(idxStart:idxStart2045) block; idxStart idxStart 2046; end end这一点对于初学者来说可能觉得无关紧要但如果你的仿真数据量到几十秒甚至上百秒循环里动态拼接会让代码跑到怀疑人生预分配后速度能提升好几个数量级。3.2 导航电文不再只有随机数上面用了随机数模拟导航电文实际项目中通常需要更贴近真实场景。北斗B1I的D1导航电文包含帧同步码、子帧计数、星历参数、电离层修正等结构如果是完整做接收机当然要按协议一层层拼。但如果做的是基带信号仿真和捕获验证用随机比特或固定的伪随机序列完全可以满足需求重点是把信号的时间增益和码结构做对。我在自己项目中用一个简单的函数生成固定模式的导航电文方便定位问题function navBits createNavBits(numBits) % 生成固定的导航比特便于仿真调试 pattern [1 0 1 1 0 0 1 0]; navBits repmat(pattern, 1, ceil(numBits/length(pattern))); navBits navBits(1:numBits); navBits(navBits 0) -1; end固定模式的好处是当你在接收端解调出数据以后一眼就能看出比特顺序对不对、有没有发生滑位。用纯随机序列当然也能验证但出了问题很难判断是解调错误还是本来就预期范围内的随机数据。3.3 从基带到中频采样率设计与混频实现基带码片序列要变成可处理的采样信号中间要经过两个核心步骤上采样和载波混频。这里最关键的参数是采样率它决定了整个仿真的计算量和模拟精度。我常用的做法是先确定每个码片的采样点数再反推采样率。B1I码速率2.046 Mcps如果每个码片取4个采样点采样率就是8.184 MHz。之所以不用整数倍的2.046 MHz而是取4倍采样是为了匹配后续捕获算法的FFT处理同时兼顾频谱观察的直观性。上采样直接用repelem实现它会把每个码片复制成指定数量的采样点sps 4; % 每个码片采样点数 fs 2.046e6 * sps; % 采样率 8.184 MHz sigBase repelem(baseband, sps); % 混频到中频 t (0:length(sigBase)-1) / fs; ifFreq 4.092e6; % 中频选为码速率2倍附近便于观察频谱 lo exp(1j * 2 * pi * ifFreq * t); sigIF real(sigBase .* lo);中频选择4.092 MHz也就是采样率的一半这样频谱上信号位于正半轴负半轴是对称镜像。要注意lo用了复数形式这样信号会被搬移到中频正频率处取实部以后得到实中频信号。如果你希望信号留在基带也可以不做混频直接用sigBase加噪声取决于你的验证目标。捕获算法通常在中频或者基带都能工作但混到中频更接近真实接收机的处理流程。加噪声时有个容易犯错的地方很多人直接用awgn函数随手加一个信噪比但awgn默认按信号功率计算而信噪比的定义要明确是码片信噪比还是采样点信噪比。我更习惯手动计算噪声功率% 按载噪比添加噪声C/N0单位dB-Hz CN0 45; % 典型GNSS载噪比 Tcoh 0.001; % 相干积分时间1ms signalPower mean(sigIF.^2); SNR CN0 - 10*log10(1/Tcoh); % 转换为1ms带宽内的信噪比(dB) noisePower signalPower / (10^(SNR/10)); noise sqrt(noisePower) * randn(size(sigIF)); rxSignal sigIF noise;我习惯把载噪比作为加噪的原始输入而不是直接写信噪比因为GNSS仿真中大家聊的都是C/N0这样你的代码参数和其他接收机软件接口对得上不会出现“我说的是45 dB你说的是15 dB”这种鸡同鸭讲的情况。4. 捕获验证如何用自己生成的信号自证清白4.1 基于FFT的并行码相位搜索信号生成了但不能自说自话就认为它是对的。最有效的自验证手段是让仿真信号通过一个捕获模块如果能捕获到正确的码相位和多普勒频率就说明信号链路是通的。我用的是基于FFT的并行码相位搜索算法这是目前软件接收机中最经典也最有效的捕获方案。核心思想是对接收信号做FFT、对本地PRN码做FFT再取共轭两者在频域相乘后做IFFT结果就是时域的循环相关。这样一次IFFT就能同时计算出所有码相位的相关值避免了逐码相位滑动的巨大计算量。对应的MATLAB实现如下function [doppler, codePhase, peakMetric] acquisition(x, prn, fs, sps, dopplerRange) % 基于FFT的并行码相位捕获 % x 接收中频信号 % prn 本地PRN码0/1序列 % fs 采样率 % sps 每码片采样点数 % dopplerRange搜索的多普勒频率范围 prnCode prn; prnCode(prnCode 0) -1; localCode repelem(prnCode, sps); % 1ms信号的采样点数 N length(localCode); peakMetric 0; doppler 0; codePhase 0; for fd dopplerRange % 补偿多普勒后的接收信号 t (0:N-1) / fs; xComp x(1:N) .* exp(-1j * 2 * pi * fd * t); % FFT相关 Xf fft(xComp); Cf conj(fft(localCode)); R ifft(Xf .* Cf); R abs(R); % 寻找当前多普勒下的相关峰 [peak, idx] max(R); % 记录最大相关峰 if peak peakMetric peakMetric peak; codePhase idx; doppler fd; end end % 计算噪声底 noiseFloor mean(R(R peakMetric * 0.5)); peakMetric peakMetric / noiseFloor; % 输出峰值与噪声底比值 end这里的多普勒搜索范围选择很关键。如果是静止场景搜索范围可以缩小到±5 kHz频率步进取500 Hz整个搜索不到几十次FFT速度非常快。如果模拟高动态场景多普勒范围扩展到±20 kHz步进可以相应增粗到1 kHz牺牲一点精度换取速度。要注意代码中计算噪声底时用了mean(R(R peakMetric * 0.5))这个方式这种简单阈值分离是工程上常用的做法如果相关峰太强会把整个均值抬高影响噪声底的估计准确性加一个峰值排除条件会更稳健。4.2 捕获结果判读捕获完成后结果可以通过三维图直观展示surf(dopplerRange, 1:N, RMatrix); view(0, 90); xlabel(多普勒频率 (Hz)); ylabel(码相位 (采样点));正常的捕获结果应该看到明显的“图钉”状峰值峰值所在位置的X坐标对应多普勒频率Y坐标对应码相位。如果出现峰值不明显或者出现多峰并列的情况通常有几个原因。第一是本地码和信号码相位差超过了一个码片这时候相关值本身就很低需要在更大的范围内搜索或者考虑使用多普勒补偿后再跑一次捕获。第二是多普勒频率没有对齐导致信号能量被分散到多个频点相关峰会变缓变矮。第三是码相位步进与采样点比例不对如果codePhase对应到芯片级别的换算出现偏差会看到峰值虽然存在但出现在意料之外的位置。我调试时注意到一个有意思的现象当生成的信号没有加NH码调制时捕获也能成功峰值甚至更高加了NH码后因为NH码片翻转的原因1 ms积分的能量有一半可能被抵消峰值会有所下降。这不是错误而是真实信号本来就要这样处理。如果要提升捕获灵敏度可以改用半码片长度的相干积分或者在捕获前先做NH码剥离这是另一个话题但在仿真验证时不要因为峰值下降了就怀疑代码出错。5. 仿真中容易踩的坑采样率、码相位与频谱泄漏5.1 采样率不是随便定的我知道很多初学者会习惯性地把采样率设成一个整数值比如10 MHz、20 MHz这样看起来整齐但和码速率2.046 Mcps的匹配关系就变得很别扭。如果采样率与码速率不是整数倍关系每个码片的采样点数就会因为四舍五入而抖动导致信号时间基准出现漂移。尤其在长数据仿真中这种漂移会不断累积最终使信号和本地码的时间对齐完全错乱。我强烈建议把采样率设成码速率的整数倍最常用的是4倍采样也就是8.184 MHz。这样每个码片严格对应4个采样点码相位计算时只需要做一个简单映射不需要处理小数采样偏移。5.2 码相位对齐最容易出错的“最后一米”捕获输出的是采样点索引这个索引要换算成码片偏移才能判断是否正确。以4倍采样为例采样点索引除以4取整才能得到码片偏移但这里的“取整”方式需要注意。我见过有人直接用floor或者round在边界情况下会引进一个码片左右的偏差。正确做法是先明确输出索引的计数起点。如果MATLAB的max返回的索引从1开始那么码片偏移应该是codePhaseSample idx - 1; % 转为0基 codePhaseChip codePhaseSample / sps; % 码片级偏移然后和构造信号时已知的码相位真值比较看是否一致。这个在校验闭环时非常关键我在自己的代码中会打印出真值和捕获值方便一眼看出误差是多少。5.3 频谱泄漏与做窗函数生成中频信号后我习惯先观察一下频谱再跑捕获。观察频谱时如果不加窗矩形截断带来的频谱泄漏会让谱峰变得很宽边瓣看起来像孪生峰造成误判。建议在频谱分析时加上汉宁窗或者布莱克曼窗win hann(length(sigIF)); sigWin sigIF .* win; freqAxis linspace(-fs/2, fs/2, length(sigWin)); plot(freqAxis, fftshift(abs(fft(sigWin))));做频谱分析时加窗做捕获相关时不要去加窗这是两码事。捕获算法处理的是真实信号加窗会破坏信号的恒包络特性反而降低相关峰性能。我在最初调试时有一次把加窗逻辑误用到了捕获路径结果相关峰从600掉到200多排查了很久才发现是这个低级错误。6. 从仿真到实用这套代码还能怎么往前延伸6.1 多通道与多普勒上面讲的是单星信号生成实际接收机面对的是多颗卫星叠加的混合信号。仿真多通道信号时每个通道生成独立的中频信号然后直接相加再加噪声。多普勒频率不同的卫星信号在频谱上占据不同位置捕获算法只要能扫完预设的频率范围就能把多个卫星分别捕获出来。实现多通道时要注意信号功率的分配。所有通道加上噪声后总功率会被推高如果每个通道都按单通道噪声功率加噪信噪比会被严重压低。我通常先把所有通道叠加再根据叠加后的总功率计算噪声功率保证最终载噪比符合设定值。6.2 对接软件接收机这套仿真生成的基带采样数据可以直接写入文件喂给标准的软件接收机处理。数据格式通常是int8或int16的IQ交织采样我按照实际接收机的数据格式定义输出文件% 将生成的信号写入文件格式为IQ交替的int16 sigI int16(real(rxSignal) * 32767); sigQ int16(imag(rxSignal) * 32767); sigIQ reshape([sigI; sigQ], [], 1); fid fopen(sim_bf_b1i.bin, wb); fwrite(fid, sigIQ, int16); fclose(fid);这样生成的数据文件就能被人家的接收机代码直接读取在验证接收机算法时非常方便。我在自己的项目里经常用MATLAB生成一组已知参数的“标准数据”然后交给接收机程序跑看它捕获到的多普勒、码相位这两个参数与设定值是否一致这样能快速定位接收机的问题而不需要拿着真实信号在户外反复测试。6.3 硬件在环之前的最后一道检查从纯仿真走向硬件在环或者实际射频测试之前建议在MATLAB里做一次完整的信号级验证。我关注的指标包括捕获峰值与噪声底的比值、跟踪环路的码相位误差曲线、载波多普勒的估计残差。这些仿真结果先要稳定再去接前端射频板卡否则硬件调试时信号有问题根本分不清是射频链路还是算法链路的问题。在具体的工程经验里我建议把PRN码生成函数单独抽出来做一个单元测试每次修改相关代码后重新跑一遍自相关和互相关验证这个习惯帮我省了无数排查时间。NH码的序列值在协议文档中可能由于版本差异略有不同务必以所使用的ICD文件为准不要照抄网络上某一篇博客的序列这也是我曾经踩过的坑。这套代码做完了以后如果再想往前扩展可以加入多径信道模型、电离层延迟模型或者改成北斗B1C或者B2a的现代化新信号体制。核心的仿真架构保持一致换掉扩频码和调制参数就可以复用这也是当初我花时间把每一层解耦的回报。最后再分享一个小技巧。整个仿真过程中我在关键节点都加了变量检查PRN码生成后检查长度和相关峰、NH码调制后检查20个PRN周期内是否有倒相导致的码片翻转、混频后检查频谱峰值是否落在设定的中频上。这三个检查点只要有一个没通过就停下来排查不要继续往后跑。这套检查习惯让我的仿真从“看起来对”变成了“每一步都能证明对”建议你也照着搭一遍。本文还有配套的精品资源点击获取