简介面向北斗导航学习与研究者的 MATLAB 北斗信号模拟器详解包内含可直接运行的信号生成与仿真脚本可帮助用户理解 BDS 信号生成、BOC 调制、主码/子码构造及功率谱密度分析等核心环节。压缩包共 17 个文件其中 11 个 m 脚本为主要代码分别对应 B1C、B1I、B2a、B2b、B3I 等多个频点的信号生成、子码构造与功率谱绘图另有 3 个 zbak 代码备份、1 个知识拓展 zip、1 份 README 说明文档和 1 张 B1 功率谱示意图。整个压缩包仅 6.55MB便于快速下载和本地调试。目前已有 63 人学习浏览特别适合刚接触北斗信号处理、希望从代码层面剖析导航信号生成流程的个人学习者。借助这些脚本读者可以逐步追踪从伪码生成、BOC 调制叠加到频谱绘制的完整链路结合示意图直观验证各频点信号特征为后续定位算法验证与导航应用开发打下基础。1. 北斗信号模拟器 MATLAB没有硬件时如何验证捕获算法北斗信号模拟器 MATLAB 详解重点其实落在“模拟器”三个字上。硬件信号源在导航实验室里动辄几十万但当你想验证的只是捕获算法在某个采样率、某个多普勒频偏下能否出峰用 MATLAB 自己生成一段 B1I 中频信号是完全够用的而且参数完全可控。搜 MATLAB 教程时能搜到大量工具箱用法但真正把伪码、载波、电文、NH 码和量化器串成一条完整链路的现成例子很少缺的那段通常是信号模型到采样序列之间的桥。后面按“信号模型 → 可运行实现 → 接收端验证 → 工程化校验”的顺序展开。直接做基带仿真的人可以跳到第三章抄链路还在搭仿真平台的人建议把第二章也看完因为采样率、中频和动态范围这几项在项目初期定错后面所有结果都要返工。写出来的代码在 R2023b 之后版本都可以直接跑不依赖任何额外工具箱。2. 北斗信号模拟器的信号模型与关键参数设计2.1 从 B1I 信号体制到仿真模型码、电文与载波北斗信号模拟器要做的事情是把卫星发射的射频信号折算到接收机天线之后的模拟前端采样序列。射频载波在 MATLAB 里直接仿真不现实1561 MHz 载波在几十 MHz 采样率下根本无法直接表示所以通常先搬移到中频或者基带再生成采样点。这也是绝大多数软件接收机模拟器的做法。B1I 是公开民用信号仿真时用到的核心参数如下表参数B1I 取值说明载波频率1561.098 MHz仿真时折算为中频测距码速率2.046 McpsB1I 主码速率主码长度2046 chips周期 1 ms主码周期1 ms捕获相干积分的基本单元导航电文速率50 bpsD1 电文20 ms 一个比特二次编码NH 码1 kbps 符号率与电文逐比特模二加中频仿真信号的表达式写出来是这个形式s(t) A · C(t) · D(t) · cos(2π·(f_IF f_d)·t φ0)C(t) 是测距码D(t) 是导航电文与 NH 码复合后的符号序列f_IF 是接收机中频f_d 是多普勒频偏φ0 是初始载波相位。多普勒由卫星与接收机之间的视线相对运动产生静止场景下通常几 kHz高动态场景下可能到几十 kHz。这个模型里最容易忽略的是 D(t)。很多初版模拟器只生成测距码和载波捕获阶段看不到任何问题但一旦做到比特同步或数据解调电文缺失导致的结果就是环路锁定后解不出数据。所以即便你的目标是捕获算法验证也建议把电文层从最开始就留好位置。2.2 用移位寄存器生成 B1I 测距码Gold 码结构与实现B1I 测距码属于 Gold 码族由两个 11 级 m 序列生成器模二和产生码长 2046对应 1 ms。公开文献中 G1 和 G2 的特征多项式分别为G1: x^11 x^10 x^9 x^7 x^5 x^4 x^2 1 G2: x^11 x^10 x^9 x^8 x^7 x^6 x^5 x^3 x^2 1标准实现里不同 PRN 的差异来自 G2 寄存器的相位选择。完整的相位选择表在北斗 ICD 里需要和真实卫星一一对应时可以直接查表。下面给出的函数用 PRN 编号做简化相位映射保证不同 PRN 之间有正交性可以用于捕获链路闭环验证。function ca genB1ICode(prn) % genB1ICode 生成北斗 B1I 测距码长度 2046周期 1 ms。 % 输入 prn卫星号这里用简化相位映射区分不同 PRN。 % 输出 ca双极性 1/-1 序列。 g1 ones(1, 11); g2 ones(1, 11); ca zeros(1, 2046); % 简化相位选择表标准实现应替换为 ICD 中的 G2 相位赋值 phaseTab [1, 3, 6, 10, 13, 17, 21, 25, 29, 33, 37, 41]; pShift mod(phaseTab(mod(prn - 1, 12) 1), 11); for k 1:2046 g1Out g1(11); g2Out g2(11); % G1 反馈抽头对应特征多项式中除最高次外的项 g1Feed xor(xor(xor(xor(g1(2), g1(4)), xor(g1(5), g1(7))), ... xor(g1(9), g1(10))), g1(11)); % G2 反馈抽头 g2Feed xor(xor(xor(xor(g2(2), g2(3)), xor(g2(5), g2(6))), ... xor(g2(7), g2(8))), xor(g2(9), g2(10))); s1 g1Out; s2 xor(g2Out, g2(mod(k pShift - 1, 11) 1)); ca(k) mod(s1 s2, 2); g1 [g1Feed, g1(1:10)]; g2 [g2Feed, g2(1:10)]; end ca 2 * ca - 1; % 0/1 转 1/-1 end代码里 g1Feed 和 g2Feed 的抽头位置来自特征多项式异或结果反馈到寄存器第一位同时最后一级作为输出。简化相位选择通过循环移位实现虽然与 ICD 不完全一致但每个 PRN 得到的码序列是稳定的、自洽的。如果你要严格对齐真实北斗星座只需要把 phaseTab 换成 ICD 中对应卫星的 G2 相位值即可寄存器结构和循环逻辑不用改。2.3 中频、采样率与多普勒参数怎么选采样率选择直接决定模拟信号占用的内存和处理速度也影响捕获精度。常见配置有三种采样率 fs与码速率关系适用场景常用中频10.23 MHz5 倍教学验证、小数据量1.023 MHz20.46 MHz10 倍精度均衡4.092 MHz40.92 MHz20 倍高码相位精度16.5 MHz 左右选采样率时有三个约束。第一是奈奎斯特条件中频加最大多普勒后必须低于 fs/2。第二是码相位量化误差fs 越接近码速率每个采样点对应的码片跨度越大捕获得到的码相位分辨率越差。第三是计算量1 秒 10.23 MHz 采样是 1023 万点double 类型占 80 MBfs 翻倍内存和耗时都翻倍。多普勒范围要根据载体动态预估。地面低速场景 ±5 kHz 足够机载场景要放到 ±15 kHz低轨用户可能要 ±50 kHz。这个范围在第三章会直接影响捕获频率搜索网格的设计所以模拟器入口参数最好把最大多普勒范围单独拉出来而不是写死在信号生成函数里。3. 用 MATLAB 把 B1I 中频信号仿真出来最小可运行链路3.1 测距码、载波与 1 ms 信号生成主函数1 ms 是 B1I 信号最基本的时间片也是捕获并行码相位搜索的积分长度。下面的函数以 1 ms 为单位生成中频信号所有参数由外部输入方便后续自动化扫描。function sig genB1I_IF(prn, fs, fc, dop, navBit, amp) % genB1I_IF 生成 1 ms 北斗 B1I 中频信号。 % prn : 卫星号 % fs : 采样率单位 Hz % fc : 中频频率单位 Hz % dop : 多普勒频偏单位 Hz % navBit : 当前 1 ms 内的导航电文/NH 复合符号取 0 或 1 % amp : 信号幅度 codeRate 2.046e6; % B1I 码速率 T 1e-3; % 1 ms 码周期 N round(fs * T); % 1 ms 内采样点数 t (0:N-1) / fs; code genB1ICode(prn); % 2046 个码片 codeIdx mod(round((0:N-1) * codeRate / fs), 2046) 1; cSeq code(codeIdx); % 重采样为采样率序列 phase 2 * pi * (fc dop) * t; sig amp * cSeq(:) .* (2 * navBit - 1) .* cos(phase); end代码的核心是 codeIdx 这一行把 2046 个码片按采样位置映射成采样序列。round 操作会引入不超过半个采样周期的码相位量化误差10.23 MHz 采样率下这一点对应约 0.1 个码片在捕获精度接受的范围内。相位累加直接用双精度浮点做不要把 cos 放在多层循环里反复算MATLAB 的向量化写法在数据量上更合适。3.2 把导航电文和 NH 码串成秒级信号单看 1 ms 信号D(t) 只是一个常数。真实信号里50 bps 导航电文先与 1 kbps NH 码逐比特模二加得到码速率为 1 kbps 的复合符号序列每个复合符号持续 1 ms正好覆盖一个主码周期。所以生成秒级信号时每个 1 ms 片段的 navBit 都要重新计算。fs 10.23e6; fc 1.023e6; prn 5; dop 1000; amp 1; bitStream randi([0 1], 50, 1); % 50 bps1 秒 50 个电文比特 nhCode randi([0 1], 20, 1); % 占位用标准 NH 码序列替换 sigLen fs * 1; % 1 秒信号 sig zeros(sigLen, 1); for ms 0:999 bit bitStream(floor(ms / 20) 1); nh nhCode(mod(ms, 20) 1); d mod(bit nh, 2); % 电文与 NH 码模二加 idx ms * round(fs * 1e-3) 1 : (ms 1) * round(fs * 1e-3); sig(idx) genB1I_IF(prn, fs, fc, dop, d, amp); end外层循环按毫秒跑1 秒信号就是 1000 次函数调用生成 1023 万点在普通笔记本上约十几秒。nhCode 这里用随机序列占位后续接入真实 NH 码时只需要把这一行替换成 20 bit 的标准序列。要注意 bit 和 nh 都是 0/1 序列模二加等价于异或结果仍为 0/1传入 genB1I_IF 后转成 1/-1。生成完成后可以用 plot 直接看时域包络或者用周期图确认载波频点。这个阶段常见的错误是采样率与码速率比例不对导致信号看起来是“花”的先画出 20 个码片左右的局部波形能明显看到 BPSK 相位翻转即可。3.3 加噪声与 ADC 量化模拟真实前端动态范围真实接收机在 ADC 之前的信号是模拟的量化位数和 AGC 增益会直接影响捕获灵敏度。模拟器加噪声时建议在信号域加而不是在相关系数上折算这样才能验证接收链路各环节的真实损耗。function q quantizeFrontEnd(x, nbit, scale) % quantizeFrontEnd 模拟接收机 ADC 量化。 % x : 输入模拟信号含噪声 % nbit : 量化位数 % scale : AGC 归一化参考电平输入平均幅度与该值比较 x x / scale; x max(-1, min(1, x)); % 限幅 q round(x * (2^(nbit - 1) - 1)); q(q -(2^(nbit - 1))) -(2^(nbit - 1)); q(q (2^(nbit - 1) - 1)) (2^(nbit - 1) - 1); end x sig sqrt(0.5) * randn(size(sig)); % 每采样点加复噪声 q1 quantizeFrontEnd(x, 1, 2 * std(x)); % 1 bit 量化 q4 quantizeFrontEnd(x, 4, 2 * std(x)); % 4 bit 量化量化位数对载噪比损耗有明确的理论经验值1 bit 量化大约损失 1.96 dB2 bit 约 0.55 dB4 bit 以上基本可以忽略。scale 参数对应 AGC 的参考电平scale 太小会让信号频繁削顶太大则量化台阶覆盖不到信号幅度两种情况都会让捕获峰值下降。调试时固定信号幅度、扫描 scale看捕获峰值的变化就能找到当前配置下的最优工作点。4. 接收端捕获验证与常见问题排查4.1 用并行码相位搜索验证模拟器输出模拟器输出是否正确最直接的验证方式是把信号送进一个标准的捕获算法看能否在预设的多普勒频率和码相位处找到峰值。并行码相位搜索是软件接收机里最常用的做法对输入信号做 FFT与本地码的 FFT 共轭相乘再反变换得到所有码相位的相关结果。function [freqPeak, codeChip, grid] acqSearch(sig, fs, fc, prn, fSearch) % acqSearch 并行码相位搜索。 % sig : 1 ms 中频信号 % fs, fc : 采样率、中频 % fSearch : 频率搜索向量单位 Hz % freqPeak : 估计的多普勒频率 % codeChip : 估计的码相位单位码片 N length(sig); t (0:N-1) / fs; code genB1ICode(prn); % 重采样到信号采样率 codeRef code(mod(round((0:N-1) * 2.046e6 / fs), 2046) 1); C conj(fft(codeRef)); grid zeros(length(fSearch), N); for k 1:length(fSearch) xb sig .* exp(-1j * 2 * pi * (fc fSearch(k)) * t); R ifft(fft(xb) .* C); grid(k, :) abs(R); end [~, fIdx] max(max(grid, [], 2)); [~, pIdx] max(grid(fIdx, :)); freqPeak fSearch(fIdx); codeChip (pIdx - 1) / fs * 2.046e6; % 采样点换算为码片 endfreqPeak 应该接近生成信号时设置的 dop 值codeChip 在码相位偏移设置为 0 时也应该接近 0。注意捕获是在 1 ms 数据上做的所以传入 sig 前要先用 genB1I_IF 生成一个完整的 1 ms 片段。频率搜索出来的峰值如果和设置值偏差在一个搜索步长以内链路就算打通了。4.2 多普勒搜索范围与相干积分时间怎么匹配多普勒搜索网格的设计有个固定套路先定最大动态范围再按相干积分时间定步长最后计算频率点数量。第一步定范围。模拟器设置的最大多普勒是 ±D_max接收端搜索范围至少覆盖这个区间。第二步定步长理论上频率误差导致的损耗近似为 sin(πΔfT)/(πΔfT)1 ms 相干积分下 250 Hz 步长边界损耗约 0.9 dB500 Hz 步长边界损耗接近 4 dB所以工程上常用 250 Hz。第三步算网格数±5 kHz 范围、250 Hz 步长需要 41 个频率点。这里有一个常见误配模拟器多普勒设 2 kHz接收端搜索范围只设 ±1 kHz捕获结果当然找不到峰。模拟器和接收端共用一个配置头文件比两边手写参数要稳得多。4.3 模拟器常见的四个坑和排查方法现象可能原因排查方法捕获不到峰值采样率与码速率映射错位检查 codeIdx 计算打印第一个码片对应的采样点峰值出现在 0 Hz 附近电文位或 NH 符号翻转改用同一符号的 1 ms 数据峰值幅度明显偏低量化 scale 不合适或限幅过重扫描 scale观察峰值变化频率估计偏差超过步长载波相位不连续确认相位累加用双精度不每次重置提示捕获失败时先不要怀疑算法把模拟器输出的 1 ms 信号存成 .mat 文件用同一份数据去跑不同版本的捕获代码能快速定位问题在生成端还是接收端。5. 把模拟器推向工程化多通道、加速与闭环校验5.1 多通道合路的功率控制真实场景下接收机同时看到多颗卫星各卫星信号以不同功率叠加。多通道模拟的核心是先把各颗卫星的 1 ms 信号分别生成再叠加归一化numSV 4; prnList [1 5 8 12]; dopList [0 1200 -800 500]; ampList [1 0.8 1.1 0.9]; sig zeros(round(fs * 1e-3), 1); for k 1:numSV sig sig genB1I_IF(prnList(k), fs, fc, dopList(k), 1, ampList(k)); end sig sig / max(abs(sig));合路后全部卫星的功率叠加再经过 AGC 归一化单颗卫星的信噪比会比单独生成时下降这是正常现象。仿真多通道时不要把每颗星幅度都设成 1实际星座的功率差异是接收机灵敏度测试的重要输入。5.2 从脚本到 Simulink 和 C 代码的改造仿真做完后如果要进硬件在环一般把 genB1I_IF 和 genB1ICode 抽成纯函数去掉全局变量和随机流然后交给 MATLAB Coder 生成 C 代码。纯函数的要求是输入输出类型明确循环次数固定codeIdx 这类中间变量不依赖动态分配的数组。改造之后可以在 Simulink 里用 S-Function 封装作为基带信号源模块复用。这样做的收益不只是速度更重要的是把仿真参数和算法实现从脚本里拆开。5.3 用多普勒与码率偏移做交叉验证最后一个校验技巧是拿多普勒和码率偏移互相验证。多普勒频偏 f_d 与视线方向伪距变化率的关系是 f_d -v·f_carrier/c码率偏移则等于伪距变化率除以光速再乘码速率。把这两组独立计算结果放在一起比对c 299792458; fCarrier 1561.098e6; lambdaC c / fCarrier; vLine -freqPeak * lambdaC; % 视线速度单位 m/s codeRateDoppler vLine / c * 2.046e6; % 码率偏移单位 Hz expectedCodeRate dop / fCarrier * 2.046e6; % 由设置的多普勒反推上面 expectedCodeRate 和 codeRateDoppler 应该一致误差在捕获步长换算范围内。这个交叉验证能一次性检查载波频率设置、码速率设置和捕获频率估计三个环节是否自洽也是每次改动模拟器参数后建议跑一遍的回归项。本文还有配套的精品资源点击获取