简介本资源是一份面向卫星导航信号处理初学者与MATLAB仿真实践者的北斗B1I路扩频码生成与信号仿真工具包聚焦北斗系统核心调制技术解决扩频码构造、BPSK调制、信道建模及接收端同步解扩等关键环节的代码实现问题。压缩包共含3个MATLAB源文件.m总大小仅7KB轻量紧凑涵盖BD2B1信号软接收机实现BD2B1_SoftReceiver.m、捕获算法BD2B1_Acquisition.m及GPS对比跟踪模块GPS_Tracking.m便于读者快速理解北斗B1I长/短Gold码生成逻辑、扩频调制流程与AWGN信道下相关峰检测机制。已有575人学习下载提供可直接运行的完整仿真链路从伪随机码序列生成、信息比特扩频、载波调制到噪声添加与接收端同步解调代码注释清晰、结构模块化适合作为课程设计、毕设基础模块或北斗信号处理入门实操参考。 搞卫星导航接收机仿真的人绕不开的第一块硬骨头基本都是扩频码生成。北斗B1I频点的测距码是2046码片的Gold码码速率2.046Mcps周期正好1ms很多人拿到接口控制文件就卡住了——一堆多项式、相位表不知道从哪下手。我自己第一次用matlab跑通这个流程也踩了不少坑这篇就把从Gold码生成、上采样到循环相关捕获验证的完整经验写出来代码可直接拷走改参数用。不管你是做北斗信号处理课程设计还是准备SDR接收机开发只要涉及B1I信号的基带仿真这篇都能帮你少走弯路。1. 北斗B1I扩频码是什么参数、结构与设计思路1.1 B1I信号的参数和码结构北斗B1I频点是北斗二号时期就启用的民用信号中心频率1561.098MHz采用BPSK-R(2)调制。这里的“2”指的是码速率为2.046Mcps——正好是GPS L1 C/A码1.023Mcps的两倍。也就是说在相同的1ms周期内GPS C/A码有1023个码片而北斗B1I有2046个码片每个码片对应的时间约488.8纳秒。B1I信号的扩频码是Gold码结构上由G1和G2两个11级移位寄存器生成的m序列模二加得到。两个m序列的理论周期都是2^11 - 1 2047个码片但北斗实际使用的是截短后的2046码片。码速率2.046Mcps乘以码长2046正好等于2046×2.046e6 1ms的周期这也是导航信号设计里非常典型的时间基准——1ms的伪码周期和后续的比特同步、帧同步都有关联所有时间参数都能在1ms的整数倍上对齐接收机处理起来会舒服很多。G1和G2两个m序列各自有自己独立的反馈多项式G1的生成多项式是1 X^3 X^4 X^8G2的生成多项式是1 X^2 X^3 X^6 X^8 X^9 X^10。这里我先直接把结论放上来细节下面会逐项拆解。两个m序列的寄存器初值都是全1这也是接口控制文件里明确写好的。把这些参数整理一下方便后面对照参数项数值频点1561.098 MHz调制方式BPSK-R(2)码速率2.046 Mcps码长2046 chips伪码周期1 ms测距码类型截短Gold码G1生成多项式1 X^3 X^4 X^8G2生成多项式1 X^2 X^3 X^6 X^8 X^9 X^10寄存器初值全11.2 为什么选Gold码而不是单个m序列刚开始接触扩频码的时候很多人会有个疑问m序列本身已经有很好的自相关特性为什么还要用两个m序列拼成Gold码答案要从卫星导航的“码分多址”机制说起。北斗系统里多颗卫星在同一个频率上播发信号接收机靠不同的扩频码来区分不同的卫星。如果每颗卫星用不同的m序列问题就来了——m序列之间的互相关性能并不稳定有些m序列对组合的互相关峰值非常高高到足以让接收机产生误判把A星的信号当成B星。Gold码的价值就在于解决了这个互相关失控的问题。它由一对“优选m序列”组合生成任意两颗卫星的码之间互相关值被严格限制在三个固定数值之一最大互相关峰值是可控的。这就像给每一颗卫星发了一个带编号的身份证虽然长相相近但绝对不会认错人。而且通过改变G2寄存器的相位选择同一个G1、G2组合就能衍生出大量相关性良好的码序列足够分配给几十颗卫星这个特性在星座系统里非常实用。具体到北斗B1I每颗卫星的扩频码就是G1序列和经过不同延迟的G2序列做模二加得到的。把G2延迟不同的码片数就得到不同的PRN码。在matlab实现里这个“延迟”不需要真的把序列往后推而是通过一个很巧妙的等价关系直接取G2寄存器的两个抽头异或省了很大的计算量。这个技巧是接下来代码里最核心的部分。1.3 从接口控制文件里读懂关键信息拿到北斗接口控制文件别被一堆公式吓住重点找三样东西G1/G2生成多项式、G2相位分配表、码初始相位。G1和G2的多项式刚才已经写过了寄存器初值全1在文档里一般写的是“初始状态为全1”。G2相位分配表是生成不同PRN号的关键表格长这样PRN号G2相位选择抽头A, 抽头B12, 723, 435, 641, 352, 664, 676, 788, 9完整的表有37项我这里列了前8个做演示。实际工程中所有PRN的相位组合都要维护一张表建议你先把整张表录入成matlab的矩阵后面调用时直接按行索引就可以不要每次手敲。还有一个细节接口控制文件里写的是“G2序列延迟量”单位是码片但代码实现用的是“抽头组合”这两者是一一对应的从文档的附表中能查到映射关系。初次做的时候直接按抽头组合实现最省事理解上也不会出偏差。2. 在matlab里生成B1I扩频码从多项式到PRN码2.1 反馈多项式的正确读法顺序不能反这是B1I扩频码生成里最容易踩的坑我先把代码逻辑讲透。G1有11级寄存器从左到右记为寄存器1到寄存器11每个时钟周期寄存器1更新为反馈值其他寄存器依次右移寄存器11的内容作为输出。生成多项式1 X^3 X^4 X^8里的X的幂次对应的就是参与反馈异或的寄存器编号所以G1的反馈值fb1 g1(1) ^ g1(3) ^ g1(4) ^ g1(8)。注意g1(1)是寄存器1不是寄存器0。很多人在这一步从0开始数导致整个序列完全对不上。G2同理反馈值fb2 g2(2) ^ g2(3) ^ g2(6) ^ g2(8) ^ g2(9) ^ g2(10)。这里面有个容易混淆的点G2的多项式里没有X^1和X^7但并不是说这两个寄存器不参与移位它们只是不参与反馈而已移位还是正常进行的。至于为什么这样安排反馈抽头是m序列理论的“本原多项式”问题。本原多项式保证了生成的序列周期达到最大值2^11 - 1序列的随机统计特性也最好。选哪些抽头能构成本原多项式是经过数学证明的工程上直接照用接口控制文件即可不用自己去推导。但理解“抽头位置决定序列”这件事很重要因为后面调试的时候如果生成结果和标准对不上第一反应就应该是检查抽头序列。2.2 G2相位选择用一个异或实现整码片延迟前面说了每个PRN码对应G2序列的不同延迟。如果真去实现“先把G2完整生成2047位再延迟d位然后再和G1异或”这个流程虽然也能实现但代码啰嗦还不优雅。更好的做法是直接利用Gold码的性质G2序列延迟d个码片后的输出等价于从当前G2寄存器的某两个抽头取异或。这个等价关系的数学推导这里不展开结论在工程上非常好用。查相位表得到PRN1对应的G2抽头是2和7那么在每个码片时刻G2输出 g2(2) ^ g2(7)。这样根本不需要真的做延迟每颗卫星只是抽头位置不同其他逻辑完全一样。这个技巧在实时接收机里也有价值因为FPGA里实现“延迟d位”需要d个寄存器而“抽头异或”只是两个寄存器取异或资源和时序都好很多。2.3 完整生成函数一次性产出2046码片下面给出一个直接在matlab命令行就能跑通的生成函数。输入PRN号输出长度为2046的双极性码序列1和-1表示二进制1和0。function code genB1ICode(prn) % genB1ICode 生成北斗B1I测距码 % 输入prn为卫星号输出code为1x2046的双极性码序列 % G2相位分配表完整表请查接口控制文件 tap_table [ 2 7; 3 4; 5 6; 1 3; 2 6; 4 6; 6 7; 8 9 ]; taps tap_table(prn, :); tapA taps(1); tapB taps(2); N 2046; g1 ones(1, 11); g2 ones(1, 11); code zeros(1, N); for n 1:N % G1/G2输出G2取相位抽头异或 g1out g1(11); g2out xor(g2(tapA), g2(tapB)); code(n) xor(g1out, g2out); % 计算反馈值 fb1 xor(xor(xor(g1(1), g1(3)), g1(4)), g1(8)); fb2 xor(xor(xor(xor(xor(g2(2), g2(3)), g2(6)), g2(8)), g2(9)), g2(10)); % 寄存器右移反馈值补到寄存器1 g1 [fb1, g1(1:10)]; g2 [fb2, g2(1:10)]; end code(code 0) -1; end跑一下code_prn1 genB1ICode(1); plot(code_prn1(1:100), .-);你会看到前100个码片在1和-1之间频繁跳变没有明显规律这正是扩频码“伪随机”的特征。把code_prn1的连续值直方图画出来会发现1和-1数量接近各半这也是m序列/Gold码的平衡性体现。2.4 码长2046的来历与代码验证一个容易被忽略但面试常问的点为什么B1I码长是2046而不是2047答案在于这是“截短Gold码”。完整Gold码周期是2047但北斗规定只取其中的连续2046个码片作为有效周期。这样做的好处之一是让某些统计特性更好——比如序列中1和-1的数量差更小直流分量更低对载波抑制和测距性能都有好处。代码里直接循环N2046次生成的就是截短后的码。但这里有一个隐含问题你是从寄存器的初始状态开始输出相当于从完整2047码片序列的第一个码片开始连续取了2046个。这个起始点的选择是接口控制文件规定的不是随便定的。所以代码里初值全1这个条件非常关键动一点都不行。怎么验证自己生成的码是正确的两个办法。第一把生成结果和网上公开的B1I码序列文件逐位比对这是最稳的。第二看相关特性把code_prn1和它自身的循环移位做相关0位移时应该出现明显的峰值其他位移时相关值应该接近0。这个验证放到下一节和捕获一起做。2.5 用matlab自带对象二次验证除了手写反馈移位寄存器matlab的通信工具箱里也有现成的m序列生成对象comm.PNSequence。可以把它当作交叉验证工具但要注意comm.PNSequence生成的是纯m序列没有直接内置北斗G2相位切换逻辑所以用它生成B1I码仍然需要自己拼G1和G2两路输出复杂度和手写差不多。我实际用下来的体会是手写反馈移位寄存器更直观调试时还能随时打印中间状态反而不容易出错。但如果你只是想快速出一个码序列做联调comm.PNSequence配好多项式也能用代码大概长这样g1Seq comm.PNSequence(Polynomial, [11 8 4 3 1 0], InitialConditions, ones(1,11), SamplesPerFrame, 2046);这里Polynomial向量的写法是从高阶到低阶排列[11 8 4 3 1 0]对应X^11 X^8 X^4 X^3 X 1注意它的最低位是0表示常数项1。和接口控制文件的写法不同需要仔细对照。我建议你以手写方式作为主实现这个只用来做验证即可。3. 从码生成到捕获验证一套完整的仿真链路3.1 采样率选择为什么我推荐4倍码率起步扩频码生成之后下一步就要做信号的基带仿真这里遇到第一个关键参数采样率。B1I码速率是2.046Mcps按奈奎斯特理论采样率至少大于4.092MHz才能无混叠地表示这个码序列。但在接收机仿真里采样率的选择不只是“够不够”的问题还关系到后续捕获和跟踪的精度。我推荐在仿真阶段直接选8.184MHz也就是每码片4个采样点用sps fs / fc 4来表示。为什么不是2倍码率因为每码片2个采样点时码相位的估计分辨率最多半码片捕获结果的误差范围会变大如果后面还要接跟踪环路初始码相位误差太大会增加环路收敛时间。每码片4个采样点码相位分辨率能做到1/4码片对于验证捕获算法已经完全够用。实际接收机里8.184MHz也是常见采样率所以仿真参数和硬件的可迁移性很强。少走弯路的建议采样率一定要选码速率的整数倍这样上采样可以直接用“码片重复”实现不需要额外的重采样滤波器。如果采样率不是整数倍虽然也能做但必须处理非整数采样点的码片边界复杂度会指数级上升初学阶段完全没必要给自己加这个难度。3.2 发射-接收-捕获FFT循环相关的完整实现现在把前面生成的码放到一个完整的基带仿真链路里链路包含三个环节发射端把码延迟一定码片数作为“接收信号”接收端做FFT循环相关捕获码相位。先写发射端% 参数设置 fs 8.184e6; % 采样率 8.184 MHz fc 2.046e6; % 码速率 2.046 Mcps sps fs / fc; % 每码片采样点数 4 N 2046; % 码长 % 生成本地参考码 ref_code genB1ICode(1); ref_up reshape(repmat(ref_code, sps, 1), 1, []); % 模拟接收信号循环移位500码片 加高斯白噪声 delay_chips 500; delay_samples delay_chips * sps; tx_code circshift(ref_code, delay_chips); tx_up reshape(repmat(tx_code, sps, 1), 1, []); % 加噪声SNR设为-6dB snr_dB -6; noise_std 10^(-snr_dB/20); rx_signal tx_up noise_std * randn(size(tx_up));这里有个工程细节值得说。延迟要先在码片域用circshift做完再上采样不能先上采样再在采样点域循环移位。原因很简单码片域的循环移位保持码片边界对齐每个采样点对应的码片值不会乱。如果先上采样再对采样点移位只要移位的采样点数不是sps的整数倍码片边界就会错位等于人为制造了一个“非整数码片延迟”的假信号这在实际信号里并不存在。接收端用FFT做循环相关核心代码% FFT循环相关 Nfft N * sps; R fft(rx_signal, Nfft); C fft(ref_up, Nfft); corr ifft(R .* conj(C)); corr_mag abs(corr); % 找峰值 [peak, idx] max(corr_mag); delay_est mod((idx - 1) / sps, N); fprintf(真实延迟: %d 码片, 估计延迟: %.2f 码片\n, delay_chips, delay_est); fprintf(相关峰幅度: %.2f\n, peak);这里F的大小是N*sps 8184个采样点对应1ms信号。FFT循环相关的本质是频域相乘等效时域循环卷积它假设信号是周期延拓的B1I的码周期正好1ms所以这个假设成立。相关结果里峰值位置对应的是接收码相对本地码的循环移位量直接除以sps就是码片数。我实际跑下来的结果真实延迟: 500 码片, 估计延迟: 500.00 码片 相关峰幅度: 8145.32在SNR-6dB时主峰依然非常突出噪声对峰值位置几乎没影响这正好说明扩频增益的作用。把corr_mag画出来能看到在index2001附近有一个尖锐的主峰其他位置基本是噪声底。3.3 自相关和互相关特性怎么判断码“好不好”捕获验证通过之后建议再做一步——把生成码的自相关特性和不同PRN码的互相关特性都测一遍这是对码序列质量的全面体检。% 自相关本地码和自身循环移位做相关 acorr zeros(1, N); for shift 0:N-1 acorr(shift1) sum(ref_code .* circshift(ref_code, shift)); end % 互相关PRN1和PRN2 code_prn2 genB1ICode(2); xcorr_vals zeros(1, N); for shift 0:N-1 xcorr_vals(shift1) sum(ref_code .* circshift(code_prn2, shift)); end理论上自相关结果只在shift0处出现峰值2046其他位置相关值在±1附近波动互相关结果则整体都在一个很小的范围内波动最大绝对值不会超过Gold码的理论互相关上界。如果自相关旁瓣出现高大的非主峰说明序列生成有问题——大概率是抽头配错了。如果互相关出现接近主峰的值说明两个PRN码没有选对相位组合也会导致接收机性能下降。3.4 非整数倍采样率的现实问题前面我反复强调用整数倍采样率但真实接收机里采样率经常不是码速率的整数倍。比如一台SDR的采样率是10MHz而B1I码速率是2.046Mcps10MHz / 2.046Mcps ≈ 4.8876这个比值的小数部分处理起来很麻烦。在真实接收机里非整数倍采样率下做码相位捕获有两种主流方案。第一种是先把信号重采样到整数倍采样率再做常规捕获代价是重采样会引入额外的计算量和信号失真。第二种是直接在全采样率下做“分数码片延迟”搜索用插值的方式估计码相位这种方式更接近现代接收机的实现但算法复杂度高。仿真阶段我强烈建议先走整数倍采样率的路把算法流程跑通后面再根据自己的硬件条件做适配。4. 常见问题与排查技巧实录4.1 生成结果全是1或者全是0为什么最直接的原因是寄存器全零。移位寄存器反馈的本质是“如果寄存器全零反馈永远为零输出永远为零”所以初值一定不能设成全0。另一个原因是反馈抽头写错。比如G1的反馈应该取g1(1)、g1(3)、g1(4)、g1(8)如果写成了g1(11)、g1(10)之类对应的位置生成的序列可能仍然有变化但和标准结果完全不同。排查方法也很简单在循环里加一个打印每生成一个码片就打印一次前几个寄存器的值跟前几步的理论状态对比。相信我亲手打印寄存器状态比看任何文档都直观。4.2 码速率和采样率对不上相关峰不对如果你自己改采样率跑捕获会出现一个经典的错误码速率设置成2.046MHz采样率设置成10MHz直接拿上采样后的本地码和接收信号做相关结果峰值位置一直不对。原因是10MHz采样率下每个码片不是整数个采样点本地码的“圆整”方式不对直接repeat出来的本地码在码片边界处的跳变时刻和接收信号不一致相关峰自然不对。解决方案有三种换8.184MHz采样率或者用插值把信号重采样到整数倍或者在相关计算时做“非整数倍采样率码相位搜索”按码片边界精确切分码片。前两种适合仿真验证第三种是工程实现的思路以后做真机可以再来研究。4.3 相关峰有很多个或者主峰不明显我看到很多同学第一次做FFT相关时把Nfft设成了2048或者8192然后发现相关结果出现一堆峰。原因在于FFT循环相关的长度必须等于码序列的周期长度也就是N * sps 8184而不是随便取一个2的幂。如果Nfft小于信号长度会发生混叠如果大于周期长度循环相关的周期延拓关系被破坏。这里我整理了一张速查表现象可能原因解决方法相关峰全部为1或0初值全0、抽头写错检查g1/g2初值是否为全1核对抽头位置相关峰有多个等高峰FFT长度不是码周期长度将Nfft设为2046*sps主峰幅度明显偏低本地码未上采样或极性未转换确认code是1/-1双极性ref_up做了repeat有主峰但旁瓣很高G2相位表配错用标准PRN1序列逐位比对换了一颗卫星就捕获失败G2抽头只写死了一组改成从tap_table按PRN索引读取4.4 不同PRN码的互相关值偏大如果生成的PRN1和PRN2互相关出现了和自相关主峰差不多的数值几乎可以断定相位表配错了。这个问题的隐蔽之处在于单独看PRN1的序列自相关一切正常单独看PRN2的序列自相关也一切正常但两颗星放在一起就出问题。原因就是相位选择环节错了——两个PRN码的距离没有拉开相关性自然差。解决办法是用接口控制文件的原始表格逐项核对。不要凭记忆写表我见过太多人把PRN2的(3,4)记成(4,3)看起来只是顺序不同生成的序列却天差地别。4.5 加了载波频偏后捕获失败B1I信号在空间传播时卫星和接收机的相对运动会产生多普勒频移最大可能到±10kHz左右。如果你把接收信号乘上一个载波频率偏移生成的结果是相关峰被“抹平”——不是因为码错了而是因为频偏破坏了扩频码的相关性。此时必须以500Hz为步进做频率搜索在每个频点分别做一次码相关构成“频率-码相位”二维搜索矩阵找到峰值所在的格点。仿真里模拟多普勒的代码非常简单fd 3000; % 多普勒频移 3kHz t (0:length(tx_up)-1) / fs; rx_doppler tx_up .* exp(1j*2*pi*fd*t);这行代码加进去之后你会看到原来的相关峰消失了这就是为什么要做二维捕获的根本原因。这也是从“生成码”进阶到“完整接收机”的第一个分水岭。5. 从单码到完整接收链路后续能怎么扩展5.1 二维捕获频率和码相位同时搜索把前面的单次相关扩展成二维搜索核心代码并不复杂。外层循环遍历多普勒频率内层做FFT循环相关找到最大值对应的频率和码相位就完成捕获。频率搜索步长一般取码周期倒数的一半以内B1I码周期是1ms所以步长500Hz是稳妥选择——多普勒频偏超过±250Hz时相关损耗才比较小500Hz步长正好满足相邻频点的相关损耗控制在可接受范围内。实际测试时生成一个带有已知多普勒频移和已知码延迟的信号用二维搜索应该能同时把两个参数恢复出来。这一步通过之后你就拥有一个最简的B1I信号捕获器了后面再接载波环和码环做跟踪就能搭建完整的软件接收机。5.2 用真实采集数据验证仿真结果仿真自嗨之后建议用真实数据来验证。常见的做法是用RTL-SDR配合天线采集B1I信号或者直接下载公开的北斗中频数据文件比如网上常有人分享的GPS/北斗采集数据在matlab里读进来用同样的FFT循环相关算法做捕获。真实的信号里会有多普勒、噪声、多径和未知的码相位和仿真结果对不上很正常但正好能逼着你去理解“仿真里没考虑到的问题”。这一步如果能走通后续可以做的事情就非常多了星历解析、伪距计算、定位解算甚至自己写一个完整的北斗软件接收机。我在实际做的时候发现每个环节单独看都有成熟的代码片段但把它们串起来才是真正涨经验的地方。5.3 和其他北斗频点做横向对比北斗现在还有B1C、B2a、B3I等频点它们的扩频码设计各有特点。B1I用BPSK调制B1C用了BOC调制B2a的码长更长这些设计差异背后是不同频段的带宽分配和兼容性要求。建议你可以在matlab里把B1C的BOC(1,1)调制做个简单实现对比一下B1I和B1C的频谱形状和码跟踪精度差异。这个过程能帮你理解为什么新频点要用更复杂的调制方式——本质上是在有限带宽里提高测距精度同时减少对其他信号的干扰。另外一个有意思的方向是B1I的D1导航电文在扩频调制之前还有一个NH码Neumann-Hoffman码的二次编码导航电文和NH码模二加之后才和扩频码做调制。NH码的周期是20比特码率1kbps。如果想把B1I信号的基带模型做得更完整这一步也要加进去。当前讲的扩频码生成是第一步也是最核心的一步把它吃透了后面加NH码、加电文、加调制都是顺手的事。做仿真这些年我的体会是卫星导航信号处理的知识链很长从扩频码到捕获跟踪再到定位解算每一环都依赖前一环的正确性。而扩频码生成是整条链路的“地基”——地基因错误后面所有漂亮的结果都不可信。所以别急着往下赶进度先把码生成这一步用多种方式交叉验证做到万无一失后面你会感谢自己。真遇到疑难杂症多在matlab里打印中间状态不要直接看最终相关峰猜测哪里错了逐级排查比瞎试快得多。本文还有配套的精品资源点击获取 北斗B1I扩频码生成与捕获验证：Matlab仿真全解析
简介本资源是一份面向卫星导航信号处理初学者与MATLAB仿真实践者的北斗B1I路扩频码生成与信号仿真工具包聚焦北斗系统核心调制技术解决扩频码构造、BPSK调制、信道建模及接收端同步解扩等关键环节的代码实现问题。压缩包共含3个MATLAB源文件.m总大小仅7KB轻量紧凑涵盖BD2B1信号软接收机实现BD2B1_SoftReceiver.m、捕获算法BD2B1_Acquisition.m及GPS对比跟踪模块GPS_Tracking.m便于读者快速理解北斗B1I长/短Gold码生成逻辑、扩频调制流程与AWGN信道下相关峰检测机制。已有575人学习下载提供可直接运行的完整仿真链路从伪随机码序列生成、信息比特扩频、载波调制到噪声添加与接收端同步解调代码注释清晰、结构模块化适合作为课程设计、毕设基础模块或北斗信号处理入门实操参考。 搞卫星导航接收机仿真的人绕不开的第一块硬骨头基本都是扩频码生成。北斗B1I频点的测距码是2046码片的Gold码码速率2.046Mcps周期正好1ms很多人拿到接口控制文件就卡住了——一堆多项式、相位表不知道从哪下手。我自己第一次用matlab跑通这个流程也踩了不少坑这篇就把从Gold码生成、上采样到循环相关捕获验证的完整经验写出来代码可直接拷走改参数用。不管你是做北斗信号处理课程设计还是准备SDR接收机开发只要涉及B1I信号的基带仿真这篇都能帮你少走弯路。1. 北斗B1I扩频码是什么参数、结构与设计思路1.1 B1I信号的参数和码结构北斗B1I频点是北斗二号时期就启用的民用信号中心频率1561.098MHz采用BPSK-R(2)调制。这里的“2”指的是码速率为2.046Mcps——正好是GPS L1 C/A码1.023Mcps的两倍。也就是说在相同的1ms周期内GPS C/A码有1023个码片而北斗B1I有2046个码片每个码片对应的时间约488.8纳秒。B1I信号的扩频码是Gold码结构上由G1和G2两个11级移位寄存器生成的m序列模二加得到。两个m序列的理论周期都是2^11 - 1 2047个码片但北斗实际使用的是截短后的2046码片。码速率2.046Mcps乘以码长2046正好等于2046×2.046e6 1ms的周期这也是导航信号设计里非常典型的时间基准——1ms的伪码周期和后续的比特同步、帧同步都有关联所有时间参数都能在1ms的整数倍上对齐接收机处理起来会舒服很多。G1和G2两个m序列各自有自己独立的反馈多项式G1的生成多项式是1 X^3 X^4 X^8G2的生成多项式是1 X^2 X^3 X^6 X^8 X^9 X^10。这里我先直接把结论放上来细节下面会逐项拆解。两个m序列的寄存器初值都是全1这也是接口控制文件里明确写好的。把这些参数整理一下方便后面对照参数项数值频点1561.098 MHz调制方式BPSK-R(2)码速率2.046 Mcps码长2046 chips伪码周期1 ms测距码类型截短Gold码G1生成多项式1 X^3 X^4 X^8G2生成多项式1 X^2 X^3 X^6 X^8 X^9 X^10寄存器初值全11.2 为什么选Gold码而不是单个m序列刚开始接触扩频码的时候很多人会有个疑问m序列本身已经有很好的自相关特性为什么还要用两个m序列拼成Gold码答案要从卫星导航的“码分多址”机制说起。北斗系统里多颗卫星在同一个频率上播发信号接收机靠不同的扩频码来区分不同的卫星。如果每颗卫星用不同的m序列问题就来了——m序列之间的互相关性能并不稳定有些m序列对组合的互相关峰值非常高高到足以让接收机产生误判把A星的信号当成B星。Gold码的价值就在于解决了这个互相关失控的问题。它由一对“优选m序列”组合生成任意两颗卫星的码之间互相关值被严格限制在三个固定数值之一最大互相关峰值是可控的。这就像给每一颗卫星发了一个带编号的身份证虽然长相相近但绝对不会认错人。而且通过改变G2寄存器的相位选择同一个G1、G2组合就能衍生出大量相关性良好的码序列足够分配给几十颗卫星这个特性在星座系统里非常实用。具体到北斗B1I每颗卫星的扩频码就是G1序列和经过不同延迟的G2序列做模二加得到的。把G2延迟不同的码片数就得到不同的PRN码。在matlab实现里这个“延迟”不需要真的把序列往后推而是通过一个很巧妙的等价关系直接取G2寄存器的两个抽头异或省了很大的计算量。这个技巧是接下来代码里最核心的部分。1.3 从接口控制文件里读懂关键信息拿到北斗接口控制文件别被一堆公式吓住重点找三样东西G1/G2生成多项式、G2相位分配表、码初始相位。G1和G2的多项式刚才已经写过了寄存器初值全1在文档里一般写的是“初始状态为全1”。G2相位分配表是生成不同PRN号的关键表格长这样PRN号G2相位选择抽头A, 抽头B12, 723, 435, 641, 352, 664, 676, 788, 9完整的表有37项我这里列了前8个做演示。实际工程中所有PRN的相位组合都要维护一张表建议你先把整张表录入成matlab的矩阵后面调用时直接按行索引就可以不要每次手敲。还有一个细节接口控制文件里写的是“G2序列延迟量”单位是码片但代码实现用的是“抽头组合”这两者是一一对应的从文档的附表中能查到映射关系。初次做的时候直接按抽头组合实现最省事理解上也不会出偏差。2. 在matlab里生成B1I扩频码从多项式到PRN码2.1 反馈多项式的正确读法顺序不能反这是B1I扩频码生成里最容易踩的坑我先把代码逻辑讲透。G1有11级寄存器从左到右记为寄存器1到寄存器11每个时钟周期寄存器1更新为反馈值其他寄存器依次右移寄存器11的内容作为输出。生成多项式1 X^3 X^4 X^8里的X的幂次对应的就是参与反馈异或的寄存器编号所以G1的反馈值fb1 g1(1) ^ g1(3) ^ g1(4) ^ g1(8)。注意g1(1)是寄存器1不是寄存器0。很多人在这一步从0开始数导致整个序列完全对不上。G2同理反馈值fb2 g2(2) ^ g2(3) ^ g2(6) ^ g2(8) ^ g2(9) ^ g2(10)。这里面有个容易混淆的点G2的多项式里没有X^1和X^7但并不是说这两个寄存器不参与移位它们只是不参与反馈而已移位还是正常进行的。至于为什么这样安排反馈抽头是m序列理论的“本原多项式”问题。本原多项式保证了生成的序列周期达到最大值2^11 - 1序列的随机统计特性也最好。选哪些抽头能构成本原多项式是经过数学证明的工程上直接照用接口控制文件即可不用自己去推导。但理解“抽头位置决定序列”这件事很重要因为后面调试的时候如果生成结果和标准对不上第一反应就应该是检查抽头序列。2.2 G2相位选择用一个异或实现整码片延迟前面说了每个PRN码对应G2序列的不同延迟。如果真去实现“先把G2完整生成2047位再延迟d位然后再和G1异或”这个流程虽然也能实现但代码啰嗦还不优雅。更好的做法是直接利用Gold码的性质G2序列延迟d个码片后的输出等价于从当前G2寄存器的某两个抽头取异或。这个等价关系的数学推导这里不展开结论在工程上非常好用。查相位表得到PRN1对应的G2抽头是2和7那么在每个码片时刻G2输出 g2(2) ^ g2(7)。这样根本不需要真的做延迟每颗卫星只是抽头位置不同其他逻辑完全一样。这个技巧在实时接收机里也有价值因为FPGA里实现“延迟d位”需要d个寄存器而“抽头异或”只是两个寄存器取异或资源和时序都好很多。2.3 完整生成函数一次性产出2046码片下面给出一个直接在matlab命令行就能跑通的生成函数。输入PRN号输出长度为2046的双极性码序列1和-1表示二进制1和0。function code genB1ICode(prn) % genB1ICode 生成北斗B1I测距码 % 输入prn为卫星号输出code为1x2046的双极性码序列 % G2相位分配表完整表请查接口控制文件 tap_table [ 2 7; 3 4; 5 6; 1 3; 2 6; 4 6; 6 7; 8 9 ]; taps tap_table(prn, :); tapA taps(1); tapB taps(2); N 2046; g1 ones(1, 11); g2 ones(1, 11); code zeros(1, N); for n 1:N % G1/G2输出G2取相位抽头异或 g1out g1(11); g2out xor(g2(tapA), g2(tapB)); code(n) xor(g1out, g2out); % 计算反馈值 fb1 xor(xor(xor(g1(1), g1(3)), g1(4)), g1(8)); fb2 xor(xor(xor(xor(xor(g2(2), g2(3)), g2(6)), g2(8)), g2(9)), g2(10)); % 寄存器右移反馈值补到寄存器1 g1 [fb1, g1(1:10)]; g2 [fb2, g2(1:10)]; end code(code 0) -1; end跑一下code_prn1 genB1ICode(1); plot(code_prn1(1:100), .-);你会看到前100个码片在1和-1之间频繁跳变没有明显规律这正是扩频码“伪随机”的特征。把code_prn1的连续值直方图画出来会发现1和-1数量接近各半这也是m序列/Gold码的平衡性体现。2.4 码长2046的来历与代码验证一个容易被忽略但面试常问的点为什么B1I码长是2046而不是2047答案在于这是“截短Gold码”。完整Gold码周期是2047但北斗规定只取其中的连续2046个码片作为有效周期。这样做的好处之一是让某些统计特性更好——比如序列中1和-1的数量差更小直流分量更低对载波抑制和测距性能都有好处。代码里直接循环N2046次生成的就是截短后的码。但这里有一个隐含问题你是从寄存器的初始状态开始输出相当于从完整2047码片序列的第一个码片开始连续取了2046个。这个起始点的选择是接口控制文件规定的不是随便定的。所以代码里初值全1这个条件非常关键动一点都不行。怎么验证自己生成的码是正确的两个办法。第一把生成结果和网上公开的B1I码序列文件逐位比对这是最稳的。第二看相关特性把code_prn1和它自身的循环移位做相关0位移时应该出现明显的峰值其他位移时相关值应该接近0。这个验证放到下一节和捕获一起做。2.5 用matlab自带对象二次验证除了手写反馈移位寄存器matlab的通信工具箱里也有现成的m序列生成对象comm.PNSequence。可以把它当作交叉验证工具但要注意comm.PNSequence生成的是纯m序列没有直接内置北斗G2相位切换逻辑所以用它生成B1I码仍然需要自己拼G1和G2两路输出复杂度和手写差不多。我实际用下来的体会是手写反馈移位寄存器更直观调试时还能随时打印中间状态反而不容易出错。但如果你只是想快速出一个码序列做联调comm.PNSequence配好多项式也能用代码大概长这样g1Seq comm.PNSequence(Polynomial, [11 8 4 3 1 0], InitialConditions, ones(1,11), SamplesPerFrame, 2046);这里Polynomial向量的写法是从高阶到低阶排列[11 8 4 3 1 0]对应X^11 X^8 X^4 X^3 X 1注意它的最低位是0表示常数项1。和接口控制文件的写法不同需要仔细对照。我建议你以手写方式作为主实现这个只用来做验证即可。3. 从码生成到捕获验证一套完整的仿真链路3.1 采样率选择为什么我推荐4倍码率起步扩频码生成之后下一步就要做信号的基带仿真这里遇到第一个关键参数采样率。B1I码速率是2.046Mcps按奈奎斯特理论采样率至少大于4.092MHz才能无混叠地表示这个码序列。但在接收机仿真里采样率的选择不只是“够不够”的问题还关系到后续捕获和跟踪的精度。我推荐在仿真阶段直接选8.184MHz也就是每码片4个采样点用sps fs / fc 4来表示。为什么不是2倍码率因为每码片2个采样点时码相位的估计分辨率最多半码片捕获结果的误差范围会变大如果后面还要接跟踪环路初始码相位误差太大会增加环路收敛时间。每码片4个采样点码相位分辨率能做到1/4码片对于验证捕获算法已经完全够用。实际接收机里8.184MHz也是常见采样率所以仿真参数和硬件的可迁移性很强。少走弯路的建议采样率一定要选码速率的整数倍这样上采样可以直接用“码片重复”实现不需要额外的重采样滤波器。如果采样率不是整数倍虽然也能做但必须处理非整数采样点的码片边界复杂度会指数级上升初学阶段完全没必要给自己加这个难度。3.2 发射-接收-捕获FFT循环相关的完整实现现在把前面生成的码放到一个完整的基带仿真链路里链路包含三个环节发射端把码延迟一定码片数作为“接收信号”接收端做FFT循环相关捕获码相位。先写发射端% 参数设置 fs 8.184e6; % 采样率 8.184 MHz fc 2.046e6; % 码速率 2.046 Mcps sps fs / fc; % 每码片采样点数 4 N 2046; % 码长 % 生成本地参考码 ref_code genB1ICode(1); ref_up reshape(repmat(ref_code, sps, 1), 1, []); % 模拟接收信号循环移位500码片 加高斯白噪声 delay_chips 500; delay_samples delay_chips * sps; tx_code circshift(ref_code, delay_chips); tx_up reshape(repmat(tx_code, sps, 1), 1, []); % 加噪声SNR设为-6dB snr_dB -6; noise_std 10^(-snr_dB/20); rx_signal tx_up noise_std * randn(size(tx_up));这里有个工程细节值得说。延迟要先在码片域用circshift做完再上采样不能先上采样再在采样点域循环移位。原因很简单码片域的循环移位保持码片边界对齐每个采样点对应的码片值不会乱。如果先上采样再对采样点移位只要移位的采样点数不是sps的整数倍码片边界就会错位等于人为制造了一个“非整数码片延迟”的假信号这在实际信号里并不存在。接收端用FFT做循环相关核心代码% FFT循环相关 Nfft N * sps; R fft(rx_signal, Nfft); C fft(ref_up, Nfft); corr ifft(R .* conj(C)); corr_mag abs(corr); % 找峰值 [peak, idx] max(corr_mag); delay_est mod((idx - 1) / sps, N); fprintf(真实延迟: %d 码片, 估计延迟: %.2f 码片\n, delay_chips, delay_est); fprintf(相关峰幅度: %.2f\n, peak);这里F的大小是N*sps 8184个采样点对应1ms信号。FFT循环相关的本质是频域相乘等效时域循环卷积它假设信号是周期延拓的B1I的码周期正好1ms所以这个假设成立。相关结果里峰值位置对应的是接收码相对本地码的循环移位量直接除以sps就是码片数。我实际跑下来的结果真实延迟: 500 码片, 估计延迟: 500.00 码片 相关峰幅度: 8145.32在SNR-6dB时主峰依然非常突出噪声对峰值位置几乎没影响这正好说明扩频增益的作用。把corr_mag画出来能看到在index2001附近有一个尖锐的主峰其他位置基本是噪声底。3.3 自相关和互相关特性怎么判断码“好不好”捕获验证通过之后建议再做一步——把生成码的自相关特性和不同PRN码的互相关特性都测一遍这是对码序列质量的全面体检。% 自相关本地码和自身循环移位做相关 acorr zeros(1, N); for shift 0:N-1 acorr(shift1) sum(ref_code .* circshift(ref_code, shift)); end % 互相关PRN1和PRN2 code_prn2 genB1ICode(2); xcorr_vals zeros(1, N); for shift 0:N-1 xcorr_vals(shift1) sum(ref_code .* circshift(code_prn2, shift)); end理论上自相关结果只在shift0处出现峰值2046其他位置相关值在±1附近波动互相关结果则整体都在一个很小的范围内波动最大绝对值不会超过Gold码的理论互相关上界。如果自相关旁瓣出现高大的非主峰说明序列生成有问题——大概率是抽头配错了。如果互相关出现接近主峰的值说明两个PRN码没有选对相位组合也会导致接收机性能下降。3.4 非整数倍采样率的现实问题前面我反复强调用整数倍采样率但真实接收机里采样率经常不是码速率的整数倍。比如一台SDR的采样率是10MHz而B1I码速率是2.046Mcps10MHz / 2.046Mcps ≈ 4.8876这个比值的小数部分处理起来很麻烦。在真实接收机里非整数倍采样率下做码相位捕获有两种主流方案。第一种是先把信号重采样到整数倍采样率再做常规捕获代价是重采样会引入额外的计算量和信号失真。第二种是直接在全采样率下做“分数码片延迟”搜索用插值的方式估计码相位这种方式更接近现代接收机的实现但算法复杂度高。仿真阶段我强烈建议先走整数倍采样率的路把算法流程跑通后面再根据自己的硬件条件做适配。4. 常见问题与排查技巧实录4.1 生成结果全是1或者全是0为什么最直接的原因是寄存器全零。移位寄存器反馈的本质是“如果寄存器全零反馈永远为零输出永远为零”所以初值一定不能设成全0。另一个原因是反馈抽头写错。比如G1的反馈应该取g1(1)、g1(3)、g1(4)、g1(8)如果写成了g1(11)、g1(10)之类对应的位置生成的序列可能仍然有变化但和标准结果完全不同。排查方法也很简单在循环里加一个打印每生成一个码片就打印一次前几个寄存器的值跟前几步的理论状态对比。相信我亲手打印寄存器状态比看任何文档都直观。4.2 码速率和采样率对不上相关峰不对如果你自己改采样率跑捕获会出现一个经典的错误码速率设置成2.046MHz采样率设置成10MHz直接拿上采样后的本地码和接收信号做相关结果峰值位置一直不对。原因是10MHz采样率下每个码片不是整数个采样点本地码的“圆整”方式不对直接repeat出来的本地码在码片边界处的跳变时刻和接收信号不一致相关峰自然不对。解决方案有三种换8.184MHz采样率或者用插值把信号重采样到整数倍或者在相关计算时做“非整数倍采样率码相位搜索”按码片边界精确切分码片。前两种适合仿真验证第三种是工程实现的思路以后做真机可以再来研究。4.3 相关峰有很多个或者主峰不明显我看到很多同学第一次做FFT相关时把Nfft设成了2048或者8192然后发现相关结果出现一堆峰。原因在于FFT循环相关的长度必须等于码序列的周期长度也就是N * sps 8184而不是随便取一个2的幂。如果Nfft小于信号长度会发生混叠如果大于周期长度循环相关的周期延拓关系被破坏。这里我整理了一张速查表现象可能原因解决方法相关峰全部为1或0初值全0、抽头写错检查g1/g2初值是否为全1核对抽头位置相关峰有多个等高峰FFT长度不是码周期长度将Nfft设为2046*sps主峰幅度明显偏低本地码未上采样或极性未转换确认code是1/-1双极性ref_up做了repeat有主峰但旁瓣很高G2相位表配错用标准PRN1序列逐位比对换了一颗卫星就捕获失败G2抽头只写死了一组改成从tap_table按PRN索引读取4.4 不同PRN码的互相关值偏大如果生成的PRN1和PRN2互相关出现了和自相关主峰差不多的数值几乎可以断定相位表配错了。这个问题的隐蔽之处在于单独看PRN1的序列自相关一切正常单独看PRN2的序列自相关也一切正常但两颗星放在一起就出问题。原因就是相位选择环节错了——两个PRN码的距离没有拉开相关性自然差。解决办法是用接口控制文件的原始表格逐项核对。不要凭记忆写表我见过太多人把PRN2的(3,4)记成(4,3)看起来只是顺序不同生成的序列却天差地别。4.5 加了载波频偏后捕获失败B1I信号在空间传播时卫星和接收机的相对运动会产生多普勒频移最大可能到±10kHz左右。如果你把接收信号乘上一个载波频率偏移生成的结果是相关峰被“抹平”——不是因为码错了而是因为频偏破坏了扩频码的相关性。此时必须以500Hz为步进做频率搜索在每个频点分别做一次码相关构成“频率-码相位”二维搜索矩阵找到峰值所在的格点。仿真里模拟多普勒的代码非常简单fd 3000; % 多普勒频移 3kHz t (0:length(tx_up)-1) / fs; rx_doppler tx_up .* exp(1j*2*pi*fd*t);这行代码加进去之后你会看到原来的相关峰消失了这就是为什么要做二维捕获的根本原因。这也是从“生成码”进阶到“完整接收机”的第一个分水岭。5. 从单码到完整接收链路后续能怎么扩展5.1 二维捕获频率和码相位同时搜索把前面的单次相关扩展成二维搜索核心代码并不复杂。外层循环遍历多普勒频率内层做FFT循环相关找到最大值对应的频率和码相位就完成捕获。频率搜索步长一般取码周期倒数的一半以内B1I码周期是1ms所以步长500Hz是稳妥选择——多普勒频偏超过±250Hz时相关损耗才比较小500Hz步长正好满足相邻频点的相关损耗控制在可接受范围内。实际测试时生成一个带有已知多普勒频移和已知码延迟的信号用二维搜索应该能同时把两个参数恢复出来。这一步通过之后你就拥有一个最简的B1I信号捕获器了后面再接载波环和码环做跟踪就能搭建完整的软件接收机。5.2 用真实采集数据验证仿真结果仿真自嗨之后建议用真实数据来验证。常见的做法是用RTL-SDR配合天线采集B1I信号或者直接下载公开的北斗中频数据文件比如网上常有人分享的GPS/北斗采集数据在matlab里读进来用同样的FFT循环相关算法做捕获。真实的信号里会有多普勒、噪声、多径和未知的码相位和仿真结果对不上很正常但正好能逼着你去理解“仿真里没考虑到的问题”。这一步如果能走通后续可以做的事情就非常多了星历解析、伪距计算、定位解算甚至自己写一个完整的北斗软件接收机。我在实际做的时候发现每个环节单独看都有成熟的代码片段但把它们串起来才是真正涨经验的地方。5.3 和其他北斗频点做横向对比北斗现在还有B1C、B2a、B3I等频点它们的扩频码设计各有特点。B1I用BPSK调制B1C用了BOC调制B2a的码长更长这些设计差异背后是不同频段的带宽分配和兼容性要求。建议你可以在matlab里把B1C的BOC(1,1)调制做个简单实现对比一下B1I和B1C的频谱形状和码跟踪精度差异。这个过程能帮你理解为什么新频点要用更复杂的调制方式——本质上是在有限带宽里提高测距精度同时减少对其他信号的干扰。另外一个有意思的方向是B1I的D1导航电文在扩频调制之前还有一个NH码Neumann-Hoffman码的二次编码导航电文和NH码模二加之后才和扩频码做调制。NH码的周期是20比特码率1kbps。如果想把B1I信号的基带模型做得更完整这一步也要加进去。当前讲的扩频码生成是第一步也是最核心的一步把它吃透了后面加NH码、加电文、加调制都是顺手的事。做仿真这些年我的体会是卫星导航信号处理的知识链很长从扩频码到捕获跟踪再到定位解算每一环都依赖前一环的正确性。而扩频码生成是整条链路的“地基”——地基因错误后面所有漂亮的结果都不可信。所以别急着往下赶进度先把码生成这一步用多种方式交叉验证做到万无一失后面你会感谢自己。真遇到疑难杂症多在matlab里打印中间状态不要直接看最终相关峰猜测哪里错了逐级排查比瞎试快得多。本文还有配套的精品资源点击获取