简介本资源是一套面向卫星导航信号处理学习者与北斗系统开发者的MATLAB仿真代码包聚焦北斗B1I路信号扩频码生成与基带处理全流程解决扩频序列建模、BPSK调制、AWGN信道仿真及接收端同步解扩等核心实践问题。压缩包共3个文件均为.m脚本涵盖BD2B1_Acquisition捕获模块、BD2B1_SoftReceiver软接收机实现及GPS_Tracking可复用的跟踪框架参考总大小仅7KB轻量易读适合嵌入式导航算法验证与教学实验快速上手。已有575人学习下载代码严格依据北斗B1I信号规范设计完整实现Gold长码与短码生成、扩频调制、载波加载、噪声叠加及基于相关峰检测的码相位同步附带清晰注释与模块化结构便于理解扩频原理、调试接收链路或拓展多系统联合仿真。 去年我拿到一份北斗B1I中频采集数据本来以为之前写过GPS L1 C/A码的捕获程序改成北斗也就是换几个参数的事。结果本地码一生成相关峰怎么都搜不出来折腾了一整晚。后来把B1I的扩频码生成机制从头捋了一遍才发现这东西和GPS C/A码的差别远比想象中大码长不是1023而是10230周期不是1ms而是5msG2移位寄存器的初相也不是全1卫星号是靠G2序列的不同相位抽头来区分的。这篇就把我从“看ICD文档”到“在MATLAB里生成B1I扩频码、验证、再对接SDR捕获”的整个过程写出来。里面包括完整可跑的MATLAB代码、自相关互相关验证方法以及我踩过的几个坑。适合刚接触北斗信号处理、想在MATLAB里跑通B1I测距码生成和捕获仿真的朋友。1. 北斗B1I扩频码的三个关键参数先把最容易错的地方说清楚1.1 码速率2.046Mcps码长10230周期5ms先看一组对比数字这是B1I和GPS L1 C/A码最容易混淆的地方。参数GPS L1 C/A北斗B1I载波频率1575.42 MHz1561.098 MHz码速率1.023 Mcps2.046 Mcps码长1023 chip10230 chip码周期1 ms5 ms调制方式BPSK-R(1)BPSK-R(2)注意看B1I的码速率是2.046Mcps码长是10230所以周期是10230 / 2.046e6 5ms不是1ms。这个坑特别隐蔽。很多人写代码的时候看到“B1I”就默认码速率是1.023Mcps的倍数直接把GPS代码里的码长改成了10230但相关处理窗口还是按1ms来切结果捕获灵敏度掉了大约7dB——因为每5ms码周期里你只用了20%的码能量。B1I的D1导航电文速率是50bps一个电文比特20ms正好包含4个完整的码周期这个关系在设计捕获累加策略时会用到。1.2 Gold码不是普通的m序列G1/G2两个多项式必须配套B1I测距码属于Gold码族由两个m序列模二和异或生成标准里管这两个序列叫G1和G2。G1生成多项式1 x x³ x⁴ x⁸ G2生成多项式1 x² x³ x⁴ x⁸两个都是8级移位寄存器所以每个m序列的周期是2⁸ − 1 255。G1和G2都是周期255的序列模二和之后得到的Gold码序列理论上周期是255的倍数。这里有一个关键数字10230 255 × 40。也就是说一个完整的B1I码周期内G1和G2各自完整循环了40次。这不是巧合B1I码不是对Gold码做截短而是让两个m序列连续跑40个完整周期凑出10230个码片。GPS C/A码的码长1023恰好等于m序列周期102310级寄存器所以GPS是一整个m序列周期。B1I是40个m序列周期的拼接在移位寄存器实现上其实一样简单——不需要中途做任何复位操作让寄存器一直跑就行。1.3 卫星号是靠G2相位选择区分的不是换个初相就行这是B1I和GPS在码生成上最大的区别。GPS L1 C/A码区分卫星是通过改变G2移位寄存器的初相状态实现的。比如PRN1的G2初相是1001001PRN2是0101001每个卫星一套初相值一共32套。北斗B1I反过来所有卫星的G1、G2初始相位完全一样不同卫星号靠的是G2输出端的相位选择。也就是说每颗卫星并不是重新设定寄存器状态而是在G2移位寄存器不同级之间做异或组合取不同的抽头作为G2序列的输出。实际实现中每个PRN对应G2寄存器的一对抽头索引。以PRN1为例G2输出 第2级 异或 第7级。这个抽头选择表在ICD文档里有完整列表我后面代码里会给出部分。正是这个机制决定了如果你拿GPS那套“每颗卫星一个初相”的思路来写B1I或者把抽头表和初相搞混生成的码序列看起来像模像样但实际上没有一个卫星是正确的。2. 扩频码生成器原理G1/G2移位寄存器是怎么拧出码序列的2.1 移位寄存器的工作方式理解B1I扩频码生成核心就是搞清楚两级LFSR线性反馈移位寄存器怎么跑。以G1为例8级寄存器分别为r1到r8每个时钟周期做两步操作按反馈抽头计算反馈值feedback r1 ⊕ r3 ⊕ r4 ⊕ r8寄存器集体右移一位r1变为feedback原来的r1变成新的r2以此类推。抽头从哪来直接从生成多项式里读G1多项式是1 x x³ x⁴ x⁸去掉常数项1剩下的x¹、x³、x⁴、x⁸对应寄存器第1、3、4、8级。所以反馈抽头是[1 3 4 8]。G2多项式是1 x² x³ x⁴ x⁸所以G2的反馈抽头是[2 3 4 8]。这里有个很容易搞混的点不同教程里寄存器的排列方向可能不一样有人从左到右是r8到r1有人是r1到r8。抽头索引在不同排列下会变但生成的序列本质上是一样的。关键是初相和抽头必须配套使用不能前面G2初相按A方案的写法后面抽头又按B方案的索引来那必然错。2.2 为什么一个码周期要循环40个m序列周期每次移位寄存器跑255个时钟G1和G2就各自完成一个m序列周期。B1I码长10230刚好是255的40倍所以从第1个码片到第10230个码片硬件上不需要任何中间处理寄存器从头到尾不停跑就行。但有一点必须注意每个5ms码周期开始前G1和G2的初相要重新设置。因为下一个码周期仍然从同一个起点开始这是信号发射端的行为接收端生成本地码时也要遵守同样的规则否则相关处理无法对齐。B1I标准定义的初相G1初相11111111全1G2初相0101010注意G2初相在ICD文档里写的是7位“0101010”对应到8级移位寄存器的第2到第8级第1级是0。所以MATLAB数组初始化写出来是G1 ones(1, 8); G2 [0 0 1 0 1 0 1 0];这个[0 0 1 0 1 0 1 0]的每一位从左到右对应寄存器第1到第8级第2到第8级就是0、1、0、1、0、1、0正好是ICD里的“0101010”。如果你在网上找到的代码G2写的是[0 1 0 1 0 1 0 0]不用急着怀疑大概率是寄存器方向定义不同但你的抽头表要做相应调整。2.3 初相和相位选择表的正确使用方式既然所有卫星的初相一样那卫星号体现在哪里体现在G2的输出端。具体做法是每个码片周期G1输出固定取第8级G2不直接取第8级而是根据当前PRN对应的抽头把G2寄存器里两个级的状态异或后作为输出。打个比方G2序列是一条首尾相接的环形伪随机码带每个卫星相当于在这条带上不同位置插一根探针取信号。抽头不一样取到的G2序列相位就不一样与固定G1模二和之后得到的Gold码就不同。部分PRN的抽头表完整表需要查北斗ICDPRNG2抽头PRNG2抽头1(2,7)5(2,3)2(3,4)6(4,5)3(5,6)7(6,7)4(1,8)8(1,4)这里还要提醒一点有些资料里不是给出抽头对而是给出“G2延迟量”比如PRN1对应G2序列相位移为1。这两种表示本质上是一回事因为从G2寄存器的第a级和第b级异或取输出等价于取G2序列的某个循环移位后的相位。但如果你把延迟量直接当成抽头索引来用生成的码就完全不对了。所以拿到资料先确认写的是抽头对还是延迟量。3. MATLAB实现从函数到码表的完整代码与逐段解释3.1 生成函数单星码生成器下面这个函数是我在MATLAB里用的B1I单星码生成器逻辑直接对应上一节的移位寄存器原理。function code genB1ICode(prn) % 生成北斗B1I单颗卫星测距码 % 输入prn卫星号取1~37 % 输出code1x10230的0/1序列 % G1、G2初相 G1 ones(1, 8); G2 [0 0 1 0 1 0 1 0]; % 第2~8级为0101010 % 相位选择抽头表按ICD节选前16个PRN tap_table [ 2 7; 3 4; 5 6; 1 8; 2 3; 4 5; 6 7; 1 4; 2 5; 3 6; 4 7; 5 8; 1 2; 3 5; 4 6; 7 8 ]; if prn size(tap_table, 1) error(当前抽头表只覆盖PRN1~%d请补充完整表, size(tap_table, 1)); end a tap_table(prn, 1); b tap_table(prn, 2); code zeros(1, 10230); for i 1:10230 % G2按抽头输出 g2out xor(G2(a), G2(b)); % 当前码片G1第8级 xor G2抽头输出 code(i) xor(G1(8), g2out); % 计算反馈值 fb1 mod(sum(G1([1 3 4 8])), 2); % G1多项式 1xx^3x^4x^8 fb2 mod(sum(G2([2 3 4 8])), 2); % G2多项式 1x^2x^3x^4x^8 % 移位 G1 [fb1, G1(1:7)]; G2 [fb2, G2(1:7)]; end end代码逻辑不复杂但几个细节值得展开说明。反馈计算用的是mod(sum(...), 2)这是异或的向量化写法。sum里几个0/1相加后取模2等价于连续做xor。在MATLAB里循环10230次每次做两次sum取模速度完全能接受不是性能瓶颈。G2的抽头索引a和b是在每个码片周期内、移位之前直接作用于当前寄存器的状态。这和“G2延迟量”的关系是G2每移一位抽头取到的就相当于G2序列的下一个相位。所以抽头表一旦定下来整颗卫星的G2序列相位就固定了。3.2 相位选择表与输出逻辑代码里的tap_table我给了前16个PRN的抽头对。实际项目里如果只做几颗卫星的捕获实验这个表够用如果要做全星座仿真建议把ICD里的完整表整理成一个CSV文件用readmatrix读进来不要硬编码在函数里。整理表的时候有个建议把PRN、抽头a、抽头b三列存好生成函数直接按prn索引取行不要用手工在代码里复制粘贴避免抄错。输出序列是0/1二值序列对应BPSK调制的逻辑电平。做相关处理的时候要转成±1形式code * 2 - 1。0变成-11变成1。如果不转直接拿0/1序列和接收信号相乘等于把一个BPSK信号当成了OOK信号在处理相关结果会很难看。3.3 码表缓存与批量生成单次生成一颗卫星的码在MATLAB里大约几十毫秒看起来不慢。但如果你在捕获程序里要对多颗卫星、多个多普勒频点做搜索每次都重新调用genB1ICode生成码整体耗时就会变得不可接受。我的做法是用persistent变量做缓存function code getB1ICode(prn) persistent cache if isempty(cache) cache containers.Map(KeyType, double, ValueType, any); end if cache.isKey(prn) code cache(prn); else code genB1ICode(prn); code code * 2 - 1; % 转成±1 cache(prn) code; end end这样同一颗卫星的码只生成一次。如果想把所有PRN一次性预生成还可以拼成一个37 × 10230的矩阵后续做批量相关的时候直接用矩阵运算比循环快很多。4. 验证生成的码对不对自相关、互相关与FFT检查法4.1 自相关峰与旁瓣检查生成完码之后第一件事不是急着去捕获信号而是先验证码本身对不对。最直接的手段是自相关。c1 getB1ICode(1); % 循环移位自相关检查0~15个码片延迟 corr_vals zeros(1, 16); for shift 0:15 c_shift circshift(c1, [0, shift]); corr_vals(shift 1) sum(c1 .* c_shift); end stem(0:15, corr_vals);这里必须用circshift循环移位不能用普通索引平移。B1I码是周期序列相关特性是基于循环相关的定义来的。用线性移位去做旁瓣会很高容易误判码质量。期望的结果是shift0处峰值等于10230码长其他位置的值要远小于10230。正常B1I码的自相关旁瓣在几百到一千多这个量级不会出现第二个接近峰值的尖峰。如果自相关峰值明显不对比如峰值变成了255的整数倍相关值那就说明初相或者抽头表错了G1和G2的序列关系没有达到Gold码的优选对特性。4.2 不同卫星互相关性能Gold码的另一个优势是互相关有理论上限。8级m序列构成的Gold码族非归一化互相关最大值为33在一个255码片m序列周期内。B1I码实际是40个m序列周期拼接所以不同卫星之间的互相关峰值会比33大一些但相对自相关峰10230来说通常能低10到20dB。实际验证互相关的代码c2 getB1ICode(2); cross_corr zeros(1, 21); for shift -10:10 c_shift circshift(c2, [0, shift]); cross_corr(shift 11) sum(c1 .* c_shift); end max(abs(cross_corr))如果你算出来某两颗卫星的互相关峰值接近甚至超过了自相关峰值的一半那就是抽头表或者初相配置有问题。这个检查在批量生成卫星码的时候特别好用能快速定位哪几颗卫星的码写错了。4.3 与ICD参考码对拍还有一种更严格的验证方式拿ICD文档里给出的某颗卫星前几十个码片去对。B1I ICD里会给出观察时刻的测距码相位但文档翻起来比较费劲。更实用的做法是对照国外一些开源GNSS库里的B1I码表比如某些GPS-SDR仿真项目里已经有北斗B1I码生成模块把它的输出和你的MATLAB函数输出对比一致就说明伪距尺度和码相位没有系统性偏差。如果做真实信号捕获最直接的验证就是拿你生成的码和采集到的中频数据做捕获找到相关峰然后核对多普勒频率是否符合该卫星的预期范围。这个能出峰基本就说明码生成正确了。5. 从码片序列到SDR采样信号让扩频码真正参与捕获5.1 码片序列重采样成采样率对应的本地码扩频码生成出来是10230个码片但SDR采集到的信号采样率通常不是2.046Mcps。常见的采集配置比如16.369MHz采样率5ms一帧内有81845个采样点而码片只有10230个平均每个码片对应约8.00098个采样点。要生成本地码的采样序列不能简单用kron(code, ones(1,8))凑数因为16.369MHz和2.046MHz不是严格8倍关系。最稳妥的做法是索引映射fs 16.369e6; fcode 2.046e6; N round(fs * 5e-3); % 一帧采样点数 code getB1ICode(1); % ±1形式 local zeros(1, N); % 每个采样点对应一个码片索引 for k 1:N chipIdx floor((k - 1) * fcode / fs) 1; p a hrefhttps://download.csdn.net/download/wouderw/87444338 stylecolor:#ec7500;font-size:14px; 本文还有配套的精品资源点击获取 /a img altmenu-r.4af5f7ec.gif srchttps://csdnimg.cn/release/wenkucmsfe/public/img/menu-r.4af5f7ec.gif stylewidth:16px;margin-left:4px;vertical-align:text-bottom;cursor:text; /p 北斗B1I扩频码生成原理与MATLAB实现及捕获验证
简介本资源是一套面向卫星导航信号处理学习者与北斗系统开发者的MATLAB仿真代码包聚焦北斗B1I路信号扩频码生成与基带处理全流程解决扩频序列建模、BPSK调制、AWGN信道仿真及接收端同步解扩等核心实践问题。压缩包共3个文件均为.m脚本涵盖BD2B1_Acquisition捕获模块、BD2B1_SoftReceiver软接收机实现及GPS_Tracking可复用的跟踪框架参考总大小仅7KB轻量易读适合嵌入式导航算法验证与教学实验快速上手。已有575人学习下载代码严格依据北斗B1I信号规范设计完整实现Gold长码与短码生成、扩频调制、载波加载、噪声叠加及基于相关峰检测的码相位同步附带清晰注释与模块化结构便于理解扩频原理、调试接收链路或拓展多系统联合仿真。 去年我拿到一份北斗B1I中频采集数据本来以为之前写过GPS L1 C/A码的捕获程序改成北斗也就是换几个参数的事。结果本地码一生成相关峰怎么都搜不出来折腾了一整晚。后来把B1I的扩频码生成机制从头捋了一遍才发现这东西和GPS C/A码的差别远比想象中大码长不是1023而是10230周期不是1ms而是5msG2移位寄存器的初相也不是全1卫星号是靠G2序列的不同相位抽头来区分的。这篇就把我从“看ICD文档”到“在MATLAB里生成B1I扩频码、验证、再对接SDR捕获”的整个过程写出来。里面包括完整可跑的MATLAB代码、自相关互相关验证方法以及我踩过的几个坑。适合刚接触北斗信号处理、想在MATLAB里跑通B1I测距码生成和捕获仿真的朋友。1. 北斗B1I扩频码的三个关键参数先把最容易错的地方说清楚1.1 码速率2.046Mcps码长10230周期5ms先看一组对比数字这是B1I和GPS L1 C/A码最容易混淆的地方。参数GPS L1 C/A北斗B1I载波频率1575.42 MHz1561.098 MHz码速率1.023 Mcps2.046 Mcps码长1023 chip10230 chip码周期1 ms5 ms调制方式BPSK-R(1)BPSK-R(2)注意看B1I的码速率是2.046Mcps码长是10230所以周期是10230 / 2.046e6 5ms不是1ms。这个坑特别隐蔽。很多人写代码的时候看到“B1I”就默认码速率是1.023Mcps的倍数直接把GPS代码里的码长改成了10230但相关处理窗口还是按1ms来切结果捕获灵敏度掉了大约7dB——因为每5ms码周期里你只用了20%的码能量。B1I的D1导航电文速率是50bps一个电文比特20ms正好包含4个完整的码周期这个关系在设计捕获累加策略时会用到。1.2 Gold码不是普通的m序列G1/G2两个多项式必须配套B1I测距码属于Gold码族由两个m序列模二和异或生成标准里管这两个序列叫G1和G2。G1生成多项式1 x x³ x⁴ x⁸ G2生成多项式1 x² x³ x⁴ x⁸两个都是8级移位寄存器所以每个m序列的周期是2⁸ − 1 255。G1和G2都是周期255的序列模二和之后得到的Gold码序列理论上周期是255的倍数。这里有一个关键数字10230 255 × 40。也就是说一个完整的B1I码周期内G1和G2各自完整循环了40次。这不是巧合B1I码不是对Gold码做截短而是让两个m序列连续跑40个完整周期凑出10230个码片。GPS C/A码的码长1023恰好等于m序列周期102310级寄存器所以GPS是一整个m序列周期。B1I是40个m序列周期的拼接在移位寄存器实现上其实一样简单——不需要中途做任何复位操作让寄存器一直跑就行。1.3 卫星号是靠G2相位选择区分的不是换个初相就行这是B1I和GPS在码生成上最大的区别。GPS L1 C/A码区分卫星是通过改变G2移位寄存器的初相状态实现的。比如PRN1的G2初相是1001001PRN2是0101001每个卫星一套初相值一共32套。北斗B1I反过来所有卫星的G1、G2初始相位完全一样不同卫星号靠的是G2输出端的相位选择。也就是说每颗卫星并不是重新设定寄存器状态而是在G2移位寄存器不同级之间做异或组合取不同的抽头作为G2序列的输出。实际实现中每个PRN对应G2寄存器的一对抽头索引。以PRN1为例G2输出 第2级 异或 第7级。这个抽头选择表在ICD文档里有完整列表我后面代码里会给出部分。正是这个机制决定了如果你拿GPS那套“每颗卫星一个初相”的思路来写B1I或者把抽头表和初相搞混生成的码序列看起来像模像样但实际上没有一个卫星是正确的。2. 扩频码生成器原理G1/G2移位寄存器是怎么拧出码序列的2.1 移位寄存器的工作方式理解B1I扩频码生成核心就是搞清楚两级LFSR线性反馈移位寄存器怎么跑。以G1为例8级寄存器分别为r1到r8每个时钟周期做两步操作按反馈抽头计算反馈值feedback r1 ⊕ r3 ⊕ r4 ⊕ r8寄存器集体右移一位r1变为feedback原来的r1变成新的r2以此类推。抽头从哪来直接从生成多项式里读G1多项式是1 x x³ x⁴ x⁸去掉常数项1剩下的x¹、x³、x⁴、x⁸对应寄存器第1、3、4、8级。所以反馈抽头是[1 3 4 8]。G2多项式是1 x² x³ x⁴ x⁸所以G2的反馈抽头是[2 3 4 8]。这里有个很容易搞混的点不同教程里寄存器的排列方向可能不一样有人从左到右是r8到r1有人是r1到r8。抽头索引在不同排列下会变但生成的序列本质上是一样的。关键是初相和抽头必须配套使用不能前面G2初相按A方案的写法后面抽头又按B方案的索引来那必然错。2.2 为什么一个码周期要循环40个m序列周期每次移位寄存器跑255个时钟G1和G2就各自完成一个m序列周期。B1I码长10230刚好是255的40倍所以从第1个码片到第10230个码片硬件上不需要任何中间处理寄存器从头到尾不停跑就行。但有一点必须注意每个5ms码周期开始前G1和G2的初相要重新设置。因为下一个码周期仍然从同一个起点开始这是信号发射端的行为接收端生成本地码时也要遵守同样的规则否则相关处理无法对齐。B1I标准定义的初相G1初相11111111全1G2初相0101010注意G2初相在ICD文档里写的是7位“0101010”对应到8级移位寄存器的第2到第8级第1级是0。所以MATLAB数组初始化写出来是G1 ones(1, 8); G2 [0 0 1 0 1 0 1 0];这个[0 0 1 0 1 0 1 0]的每一位从左到右对应寄存器第1到第8级第2到第8级就是0、1、0、1、0、1、0正好是ICD里的“0101010”。如果你在网上找到的代码G2写的是[0 1 0 1 0 1 0 0]不用急着怀疑大概率是寄存器方向定义不同但你的抽头表要做相应调整。2.3 初相和相位选择表的正确使用方式既然所有卫星的初相一样那卫星号体现在哪里体现在G2的输出端。具体做法是每个码片周期G1输出固定取第8级G2不直接取第8级而是根据当前PRN对应的抽头把G2寄存器里两个级的状态异或后作为输出。打个比方G2序列是一条首尾相接的环形伪随机码带每个卫星相当于在这条带上不同位置插一根探针取信号。抽头不一样取到的G2序列相位就不一样与固定G1模二和之后得到的Gold码就不同。部分PRN的抽头表完整表需要查北斗ICDPRNG2抽头PRNG2抽头1(2,7)5(2,3)2(3,4)6(4,5)3(5,6)7(6,7)4(1,8)8(1,4)这里还要提醒一点有些资料里不是给出抽头对而是给出“G2延迟量”比如PRN1对应G2序列相位移为1。这两种表示本质上是一回事因为从G2寄存器的第a级和第b级异或取输出等价于取G2序列的某个循环移位后的相位。但如果你把延迟量直接当成抽头索引来用生成的码就完全不对了。所以拿到资料先确认写的是抽头对还是延迟量。3. MATLAB实现从函数到码表的完整代码与逐段解释3.1 生成函数单星码生成器下面这个函数是我在MATLAB里用的B1I单星码生成器逻辑直接对应上一节的移位寄存器原理。function code genB1ICode(prn) % 生成北斗B1I单颗卫星测距码 % 输入prn卫星号取1~37 % 输出code1x10230的0/1序列 % G1、G2初相 G1 ones(1, 8); G2 [0 0 1 0 1 0 1 0]; % 第2~8级为0101010 % 相位选择抽头表按ICD节选前16个PRN tap_table [ 2 7; 3 4; 5 6; 1 8; 2 3; 4 5; 6 7; 1 4; 2 5; 3 6; 4 7; 5 8; 1 2; 3 5; 4 6; 7 8 ]; if prn size(tap_table, 1) error(当前抽头表只覆盖PRN1~%d请补充完整表, size(tap_table, 1)); end a tap_table(prn, 1); b tap_table(prn, 2); code zeros(1, 10230); for i 1:10230 % G2按抽头输出 g2out xor(G2(a), G2(b)); % 当前码片G1第8级 xor G2抽头输出 code(i) xor(G1(8), g2out); % 计算反馈值 fb1 mod(sum(G1([1 3 4 8])), 2); % G1多项式 1xx^3x^4x^8 fb2 mod(sum(G2([2 3 4 8])), 2); % G2多项式 1x^2x^3x^4x^8 % 移位 G1 [fb1, G1(1:7)]; G2 [fb2, G2(1:7)]; end end代码逻辑不复杂但几个细节值得展开说明。反馈计算用的是mod(sum(...), 2)这是异或的向量化写法。sum里几个0/1相加后取模2等价于连续做xor。在MATLAB里循环10230次每次做两次sum取模速度完全能接受不是性能瓶颈。G2的抽头索引a和b是在每个码片周期内、移位之前直接作用于当前寄存器的状态。这和“G2延迟量”的关系是G2每移一位抽头取到的就相当于G2序列的下一个相位。所以抽头表一旦定下来整颗卫星的G2序列相位就固定了。3.2 相位选择表与输出逻辑代码里的tap_table我给了前16个PRN的抽头对。实际项目里如果只做几颗卫星的捕获实验这个表够用如果要做全星座仿真建议把ICD里的完整表整理成一个CSV文件用readmatrix读进来不要硬编码在函数里。整理表的时候有个建议把PRN、抽头a、抽头b三列存好生成函数直接按prn索引取行不要用手工在代码里复制粘贴避免抄错。输出序列是0/1二值序列对应BPSK调制的逻辑电平。做相关处理的时候要转成±1形式code * 2 - 1。0变成-11变成1。如果不转直接拿0/1序列和接收信号相乘等于把一个BPSK信号当成了OOK信号在处理相关结果会很难看。3.3 码表缓存与批量生成单次生成一颗卫星的码在MATLAB里大约几十毫秒看起来不慢。但如果你在捕获程序里要对多颗卫星、多个多普勒频点做搜索每次都重新调用genB1ICode生成码整体耗时就会变得不可接受。我的做法是用persistent变量做缓存function code getB1ICode(prn) persistent cache if isempty(cache) cache containers.Map(KeyType, double, ValueType, any); end if cache.isKey(prn) code cache(prn); else code genB1ICode(prn); code code * 2 - 1; % 转成±1 cache(prn) code; end end这样同一颗卫星的码只生成一次。如果想把所有PRN一次性预生成还可以拼成一个37 × 10230的矩阵后续做批量相关的时候直接用矩阵运算比循环快很多。4. 验证生成的码对不对自相关、互相关与FFT检查法4.1 自相关峰与旁瓣检查生成完码之后第一件事不是急着去捕获信号而是先验证码本身对不对。最直接的手段是自相关。c1 getB1ICode(1); % 循环移位自相关检查0~15个码片延迟 corr_vals zeros(1, 16); for shift 0:15 c_shift circshift(c1, [0, shift]); corr_vals(shift 1) sum(c1 .* c_shift); end stem(0:15, corr_vals);这里必须用circshift循环移位不能用普通索引平移。B1I码是周期序列相关特性是基于循环相关的定义来的。用线性移位去做旁瓣会很高容易误判码质量。期望的结果是shift0处峰值等于10230码长其他位置的值要远小于10230。正常B1I码的自相关旁瓣在几百到一千多这个量级不会出现第二个接近峰值的尖峰。如果自相关峰值明显不对比如峰值变成了255的整数倍相关值那就说明初相或者抽头表错了G1和G2的序列关系没有达到Gold码的优选对特性。4.2 不同卫星互相关性能Gold码的另一个优势是互相关有理论上限。8级m序列构成的Gold码族非归一化互相关最大值为33在一个255码片m序列周期内。B1I码实际是40个m序列周期拼接所以不同卫星之间的互相关峰值会比33大一些但相对自相关峰10230来说通常能低10到20dB。实际验证互相关的代码c2 getB1ICode(2); cross_corr zeros(1, 21); for shift -10:10 c_shift circshift(c2, [0, shift]); cross_corr(shift 11) sum(c1 .* c_shift); end max(abs(cross_corr))如果你算出来某两颗卫星的互相关峰值接近甚至超过了自相关峰值的一半那就是抽头表或者初相配置有问题。这个检查在批量生成卫星码的时候特别好用能快速定位哪几颗卫星的码写错了。4.3 与ICD参考码对拍还有一种更严格的验证方式拿ICD文档里给出的某颗卫星前几十个码片去对。B1I ICD里会给出观察时刻的测距码相位但文档翻起来比较费劲。更实用的做法是对照国外一些开源GNSS库里的B1I码表比如某些GPS-SDR仿真项目里已经有北斗B1I码生成模块把它的输出和你的MATLAB函数输出对比一致就说明伪距尺度和码相位没有系统性偏差。如果做真实信号捕获最直接的验证就是拿你生成的码和采集到的中频数据做捕获找到相关峰然后核对多普勒频率是否符合该卫星的预期范围。这个能出峰基本就说明码生成正确了。5. 从码片序列到SDR采样信号让扩频码真正参与捕获5.1 码片序列重采样成采样率对应的本地码扩频码生成出来是10230个码片但SDR采集到的信号采样率通常不是2.046Mcps。常见的采集配置比如16.369MHz采样率5ms一帧内有81845个采样点而码片只有10230个平均每个码片对应约8.00098个采样点。要生成本地码的采样序列不能简单用kron(code, ones(1,8))凑数因为16.369MHz和2.046MHz不是严格8倍关系。最稳妥的做法是索引映射fs 16.369e6; fcode 2.046e6; N round(fs * 5e-3); % 一帧采样点数 code getB1ICode(1); % ±1形式 local zeros(1, N); % 每个采样点对应一个码片索引 for k 1:N chipIdx floor((k - 1) * fcode / fs) 1; p a hrefhttps://download.csdn.net/download/wouderw/87444338 stylecolor:#ec7500;font-size:14px; 本文还有配套的精品资源点击获取 /a img altmenu-r.4af5f7ec.gif srchttps://csdnimg.cn/release/wenkucmsfe/public/img/menu-r.4af5f7ec.gif stylewidth:16px;margin-left:4px;vertical-align:text-bottom;cursor:text; /p