简介：面向北斗三号信号处理学习与科研实践的 MATLAB 开源工程，聚焦 B1C、B2a 频点的软件接收机实现，适合具备一定 GNSS 基础、希望打通从中频数据到定位结果完整链路的高年级本科生、研究生与工程人员使用。压缩包共 87 个文件，以 82 个 .m 脚本为核心，辅以 md 说明、txt 许可与元数据、1 份 pdf 文档，整体约 3.72MB；代码按初始化、参数配置、捕获、跟踪、导航电文解码、伪距测量与 PVT 解算等模块组织，并包含通用功能与 B1C/B2a 分目录实现。已有 1936 人学习，说明其在同类开源接收机资料中具有一定参考价值。读者可据此复现数据与导频双通道跟踪流程，结合 IF 数据集及元数据完成离线后处理，也可参考文档理解各 SDR 接收机之间的共性设计与测试思路，用于课程实验、算法验证和二次开发。

1. B1C/B2a 双频 SDR 接收机：一段中频数据能跑到定位吗

去年帮同事调试一段 B1C 中频采集数据，用几十行 MATLAB 代码就找到了相关峰，他以为剩下的只是按老办法堆代码。结果一上手就卡住了：B1C 是 BOC 调制，主码周期 10 ms，导频和数据分量还分开，捕获那套老参数直接套用会失效；B2a 码率 10.23 Mcps，采样率一低，相关峰就会被削平。这套 BDS-3 B1C/B2a SDR 接收机的开源 MATLAB 实现，把从 IQ 采样到 PVT 解算的整条链路完整写了出来，适合两类人：一类是想弄明白双频接收机内部到底怎么运转的学生和刚入行的工程师，另一类是手里有采集数据、想拿它当基线做二次开发的人。

2. 工程结构与信号体制：B1C/B2a 双频链路全景

2.1 从 IQ 文件到定位结果的六个环节

这套代码的骨架是一条单线程、文件驱动的处理链，没有实时性要求，也不依赖任何硬件驱动，喂一段复基带采样就能从头跑到尾。整条链路拆成六段：读采样文件并把中频搬移到零中频、按 PRN 做二维捕获、对捕获到的卫星开跟踪通道、比特同步与帧同步、电文解析出星历和钟差、观测量组帧后做最小二乘定位。每一段都是一个独立函数，段与段之间只通过结构体传递，这样改其中一段不会牵连其他部分。

前两段最容易被轻视。采样文件通常是 int16 交织的 I/Q 或者 complex float，读进来第一件事是确认字节序和交织方式，第二件事是确认采样率与标称值是否一致。很多采集设备标称 8.184 MHz，实际晶振偏几十 ppm，落到 B1C 的 10 ms 主码上就是几十纳秒的漂移，短时间看不出来，跑几分钟跟踪就慢慢滑掉。代码里把采样率、中频、信号带宽、码率全部收在一个rx_config.m里，就是为了让这种系统性偏差只有一个源头可以改。

捕获和跟踪的边界也要划清楚。捕获只负责给出「这颗星大概在哪个码相位、多普勒大概多少」，精度不需要很高，码相位误差半个码片、频率误差几十赫兹都能被跟踪环拉回来。跟踪才是把精度做上去的地方，它每个积分周期输出一次载波频率、码相位和 I/Q 累加值，后面的比特同步、电文解析、伪距生成全部建立在这串观测量上。把这两段的职责混在一起写，是新手最容易犯的结构性错误。

最后一段是观测量到 PVT。伪距的物理含义是「接收时刻减去发射时刻乘光速」，发射时刻要从电文的参考时间加上码相位折算，接收时刻来自采样点索引，两个时间必须落在同一个时间系统里。代码里用了一个t_rx数组统一管理每个采样点对应的接收时刻，所有通道共用，避免每个通道各自算一套时间导致伪距之间互相不一致。

2.2 B1C 与 B2a 的信号参数对比

双频接收机不是把两个单频接收机简单堆在一起，两个频点的信号体制差别很大，代码里必须分开处理。下面这张表是我按公开的接口文件整理的关键参数，写代码时对着它改常数基本不会错：

B1C 的难点在调制。BOC(1,1) 的功率谱主瓣被搬到了 ±1.023 MHz，BOC(6,1) 更远，合成之后的 QMBOC 相关峰比 BPSK 窄得多，主峰宽度只有 BPSK 的一半左右。窄主峰意味着测距精度更高，代价是对码相位误差更敏感，跟踪环的 DLL 间隔如果还按 BPSK 的 1 个码片去设，鉴别器曲线会直接跑出线性区。代码里超前减滞后间隔默认取 0.5 个码片，B1C 通道可以再压到 0.3 左右。

B2a 的难点在采样率。码率 10.23 Mcps，主瓣宽度 ±10.23 MHz，采样率低于 40.92 MHz 时高次谐波被砍掉，相关峰展宽、峰值下降，捕获的归一化门限要相应调低，否则明明有信号却判决不出来。另外 B2a 主码周期只有 1 ms，相干积分时长最多堆到 1 ms 的整数倍，想靠长时间积分换信噪比，就得走导频分量加次级码那条路。

两个频点的电文结构也不一样，B-CNAV1 和 B-CNAV2 的编码方式、帧长、纠错方式都不同，解析函数必须分开写。这里不要图省事写一个通用解析器，两套电文的字段位置和编码规则差异太大，强行合并只会让调试变成玄学。

2.3 目录结构与各模块职责

代码按功能分目录，一个目录只干一件事。常见做法是下面这种组织方式，我按它拆过几个类似的接收机，改起来最顺手：

BDS3-SDR/ ├── main_acq.m % 捕获入口，输出 acqResult.mat ├── main_track.m % 跟踪入口，读捕获结果，输出 trkResult.mat ├── main_pvt.m % 定位入口，读跟踪结果，输出 PVT ├── config/ │ └── rx_config.m % 采样率、中频、带宽、环路带宽、门限 ├── acq/ │ ├── acquire_b1c.m │ ├── acquire_b2a.m │ └── gen_local_code.m % 主码 + 次级码生成 ├── trk/ │ ├── track_channel.m % 单通道跟踪主循环 │ ├── dll_disc.m │ └── pll_disc.m ├── nav/ │ ├── bcnav1_decode.m │ ├── bcnav2_decode.m │ └── eph2pos.m % 星历转卫星位置 ├── pvt/ │ └── ls_pvt.m ├── data/ % 采样文件放这里 └── utils/ % FFT、坐标转换、绘图

三个main_*.m各自独立可跑，中间结果存成.mat，好处是调跟踪的时候不用每次重跑捕获。捕获一段 100 ms 的数据在普通笔记本上要几十秒，B2a 因为采样率高还要更久，每次改跟踪参数都重跑捕获，一天下来大部分时间都耗在等结果上。

gen_local_code.m是共用模块，两个频点的主码生成都从这里出。B1C 和 B2a 的主码都是 10230 长，但生成多项式和初相不同，函数入参里带一个频点标识，内部走不同分支。次级码单独放在这个文件里，因为导频分量的次级码周期长，跟踪阶段要用它做电文剥离，和主码不是一回事。

config/rx_config.m建议只放常数，不要放逻辑。我见过有人把环路带宽的自适应切换写在配置里，结果调试时改一个参数要翻三个文件才能确认最终生效的是哪个值，血泪经验。

3. 捕获模块拆解：并行码相位搜索的 MATLAB 实现

3.1 捕获主函数逐段讲

捕获做的是二维搜索：码相位一维，多普勒一维。最直接的做法是每个多普勒频点上都做一次时域相关，计算量随多普勒网格数线性增长，B2a 这种 40 MHz 采样率的数据根本跑不动。代码里用的是频域并行码相位搜索，把本地码的 FFT 只算一次，每个频点只做一次 FFT 和一次 IFFT，把码相位那一维的搜索全部压进频域乘法里。

function acq = acquire_b1c(x, fs, prn, cfg) % x : 复基带采样，列向量 % fs : 采样率 Hz % prn : 卫星 PRN 号 % cfg : 捕获配置（多普勒范围、步进、相干积分时长、判决门限） codeRate = 1.023e6; % B1C 主码速率 Nc = 10230; % 主码长度 % 1) 生成本地主码，并按采样率重采样到与信号同一时间栅格 code = gen_local_code(prn, 'B1C'); % 1 x 10230，取值 ±1 splPerCode = round(fs * Nc / codeRate); % 每个主码周期的采样点数 tLocal = (0 : splPerCode-1) / fs; localCode = interp1((0:Nc-1)/codeRate, code, mod(tLocal, Nc/codeRate), 'linear'); localCode = localCode(:); % 2) 相干积分长度：B1C 导频分量可做 10 ms 整数倍，数据分量压到 1 ms 量级 Ncoh = round(fs * cfg.cohMs * 1e-3); nBlk = floor(length(x) / Ncoh); % 3) 多普勒搜索网格 dopGrid = cfg.fMin : cfg.fStep : cfg.fMax; % 典型 -6000:250:6000 % 4) 本地码 FFT 只算一次，长度取线性卷积不混叠的下界 L = Ncoh + splPerCode - 1; codeF = fft(localCode, L); best = struct('peak', 0, 'codePhase', 0, 'doppler', 0, 'metric', []); for k = 1:nBlk seg = x((k-1)*Ncoh + (1:Ncoh)); for d = 1:numel(dopGrid) % 载波剥离：把多普勒假设搬掉 n = (0:Ncoh-1).'; baseband = seg .* exp(-1j*2*pi*dopGrid(d)*n/fs); % 频域并行码相位搜索 corr = ifft(fft(baseband, L) .* conj(codeF)); corr = corr(1:splPerCode); % 只取一个主码周期内的相关值 [v, idx] = max(abs(corr)); if v > best.peak best.peak = v; best.codePhase = idx - 1; % 单位：采样点 best.doppler = dopGrid(d); best.metric = abs(corr); end end end acq = best; acq.prn = prn; acq.splPerCode = splPerCode; acq.detRatio = best.peak / mean(abs(best.metric)); acq.detected = acq.detRatio > cfg.thrRatio; end

代码里有几个点值得展开说。第一，本地码重采样用的是interp1线性插值，对 ±1 的方波码来说插值会引入一点幅度损失，但比起码相位量化的误差，这点损失可以接受。如果想做得更干净，可以按采样率直接生成码表，省掉插值。第二，FFT 长度L取Ncoh + splPerCode - 1，保证循环卷积不产生混叠，取完corr的前splPerCode点才是有效的线性相关结果，这段长度算错是新手最常见的翻车点，相关峰会莫名其妙地分裂成两个。第三，多普勒网格是外层循环、FFT 是内层，这样本地码 FFT 可以复用，但如果多普勒范围很大，可以把 IFFT 那一维也并行化，用pagefun或者直接循环展开。

detRatio用的是峰值比均值，不是峰值比次峰。峰值比均值对单颗强星效果好，但在强信号压制弱信号的时候，均值会被强星拉高，弱星就漏了。更稳的做法是峰值比次峰，把主峰附近 ±2 个码片、±1 个频点以外的区域取最大值当次峰，代码里可以加一个掩膜。这一步多十几行，但能明显降低虚警。

3.2 多普勒网格、相干积分与门限怎么定

多普勒范围取决于应用场景。地面静态接收，B1C 的最大多普勒在 ±5 kHz 量级，B2a 因为频率低一些，量级相近但符号可能不同；车载场景要放到 ±10 kHz，机载还要更大。宁可在捕获阶段把范围给宽一点，也不要在跟踪阶段因为初值偏太远拉不回来。代价是网格点变多，计算时间线性增长。

频率步进和相干积分时长是一对。二阶环的频率牵引范围大约是1/(2*Tcoh)的量级，10 ms 相干积分对应 50 Hz 左右的半带宽，步进取 250 Hz 时相邻网格之间会有重叠，不会漏掉真正的峰。如果为了省时间把步进取到 1 kHz，牵引范围覆盖不住，跟踪环大概率失锁。我一般按fStep <= 1/(2*Tcoh)先算一遍，再往上放一档取整，250 Hz 和 500 Hz 是最常用的两个值。

相干积分时长受电文符号率限制。B1C 导频分量没有电文翻转，可以把相干积分堆到 10 ms 甚至 20 ms；B1C 数据分量和 B2a 数据分量有电文符号，符号率分别是 1000 sps 和 500 sps 量级，相干积分超过符号长度就会因为符号跳变而能量对消。所以捕获阶段如果拿不准用的是哪个分量，保守取 1 ms 相干积分，用非相干累加补信噪比，累加次数 10 到 20 次。

判决门限thrRatio没有普适值，和积分时长、采样率、噪声带宽都有关。经验值是 2.5 到 3.5 之间，先用一段确认有信号的数据把detRatio打出来，看真峰和噪声底的实际差距，再定门限。B2a 因为采样率高、噪声带宽大，同样信噪比下detRatio会比 B1C 低一截，门限要相应下调。

3.3 捕获结果怎么验证

捕获结果不能只看detected是 true 就放心。最直接的办法是把acq.metric画成二维图，横轴码相位、纵轴多普勒，看主峰是不是一个干净的尖峰。如果主峰旁边有明显旁瓣，或者主峰被削平，多半是采样率不够或者本地码重采样出了问题。

第二个验证方法是对同一时刻的多颗卫星做码相位对比。同一段数据里捕获到的几颗星，码相位折算成时间之后，两两之差应该和它们到接收机的几何距离差对得上，量级在毫秒以内。差得离谱，说明码相位到时间的换算系数写错了，常见错误是把splPerCode当成码周期用了。

第三个方法是反算伪距率。捕获给出的多普勒是载波多普勒，乘上光速除以载波频率就是伪距率，和卫星运动引起的距离变化率应该在同一个量级。如果某颗星的伪距率符号反了，说明载波剥离时用了+j而不是-j，这个错误在跟踪阶段会表现为载波频率一直往一个方向漂。

4. 跟踪与电文解调：从环路鉴别器到比特流

4.1 跟踪主循环的完整代码

跟踪是整条链路里最容易出问题的一段。捕获给的码相位和多普勒只是粗值，跟踪要在每个积分周期里把这个粗值收紧，同时输出载波相位和码相位的小数部分，供后面的伪距计算使用。代码按标准的二阶环写，DLL 和 PLL 各一路，共用同一个积分周期。

function trk = track_channel(x, fs, acq, cfg) % 单通道跟踪：输入已捕获的码相位与多普勒，输出每积分周期的观测量 codeRate = 1.023e6; Nc = 10230; splPerCode = acq.splPerCode; code = gen_local_code(acq.prn, 'B1C'); T = splPerCode / fs; % 积分周期，B1C 为 10 ms [w_dll, a_dll] = loop_gain(cfg.Bn_dll, T, 0.707); [w_pll, a_pll] = loop_gain(cfg.Bn_pll, T, 0.707); codePhase = acq.codePhase; % 单位：采样点，可带小数 carrFreq = acq.doppler; % Hz carrPhase = 0; nIter = floor((length(x) - 3*splPerCode) / splPerCode); trk = struct('prn', acq.prn, 'I_P', zeros(nIter,1), 'Q_P', zeros(nIter,1), ... 'carrFreq', zeros(nIter,1), 'codePhase', zeros(nIter,1), ... 'CN0', zeros(nIter,1)); for k = 1:nIter % ---- 1) 本地三路码，超前/即时/滞后间隔 0.5 码片 ---- offE = -0.5 * fs / codeRate; offL = 0.5 * fs / codeRate; idxE = mod(round(codePhase + offE) + (0:splPerCode-1), Nc) + 1; idxP = mod(round(codePhase) + (0:splPerCode-1), Nc) + 1; idxL = mod(round(codePhase + offL) + (0:splPerCode-1), Nc) + 1; cE = code(idxE); cP = code(idxP); cL = code(idxL); % ---- 2) 载波剥离与相关累加 ---- n = (0:splPerCode-1).'; localCarr = exp(1j*(2*pi*carrFreq*n/fs + carrPhase)); seg = x((k-1)*splPerCode + (1:splPerCode)); yE = sum(seg .* localCarr .* cE.'); yP = sum(seg .* localCarr .* cP.'); yL = sum(seg .* localCarr .* cL.'); I_E = real(yE); Q_E = imag(yE); I_P = real(yP); Q_P = imag(yP); I_L = real(yL); Q_L = imag(yL); % ---- 3) 鉴别器 ---- pE = I_E^2 + Q_E^2; pL = I_L^2 + Q_L^2; dCode = (pE - pL) / max(pE + pL, eps); % 归一化超前减滞后功率 dCarr = atan2(Q_P, I_P); % 四象限相位鉴别器 % ---- 4) 环路滤波与 NCO 更新 ---- codeFreq = codeRate + w_dll * dCode; carrFreq = carrFreq + w_pll * dCarr; carrPhase = mod(carrPhase + 2*pi*carrFreq*T, 2*pi); codePhase = mod(codePhase + codeFreq * T, Nc); % ---- 5) 记录观测量与载噪比 ---- trk.I_P(k) = I_P; trk.Q_P(k) = Q_P; trk.carrFreq(k) = carrFreq; trk.codePhase(k) = codePhase; trk.CN0(k) = 10*log10(max(I_P^2 + Q_P^2, eps) / ...); end end

loop_gain按二阶环的标准形式写，自然角频率wn = 2*Bn/(zeta + 1/(4*zeta))，两个滤波系数分别是a = 4*zeta*wn*T/(1+2*zeta*wn*T+(wn*T)^2)和与之配套的另一个系数。这里不需要手推，只要保证传入的带宽Bn和阻尼系数zeta是物理量纲，系数由公式算出来，改带宽时只改一个数。

三路码的间隔用固定的 0.5 码片，B1C 因为相关峰窄，可以改到 0.3。间隔越小，鉴别器增益越低，但线性区越窄，跟踪精度反而可能更好。实际调的时候先把间隔设 0.5 跑通，再往下压，每压一次看CN0有没有掉，掉了就说明压过头了。

4.2 鉴别器选型和环路带宽的经验值

DLL 用的是归一化超前减滞后功率，这个鉴别器在 ±0.5 码片内近似线性，超过之后曲线会翻转，所以跟踪一旦失锁，鉴别器输出会乱跳，环路自己拉不回来。判断失锁的简单办法是看CN0，连续几个周期掉到 25 dB-Hz 以下就可以认为这个通道废了，标记一下不再输出观测量，别让它污染后面的定位解算。

PLL 用四象限atan2，线性区 ±90 度，比二象限的atan稳。B1C 导频分量没有电文翻转，可以直接用纯 PLL；数据分量因为符号跳变，PLL 鉴别器会在符号边界上产生 ±180 度的跳变，这时候要么先用 FLL 把频率拉稳再切 PLL，要么用判决反馈把符号剥掉。代码里最简单的做法是前 20 个积分周期只用 FLL，频率误差收敛到几十赫兹以内再切 PLL，切的时候相位从 0 开始，不要沿用 FLL 阶段的相位。

环路带宽的经验值：静态场景 PLL 取 10 到 15 Hz，DLL 取 0.5 到 1 Hz；车载场景 PLL 放到 25 到 30 Hz，DLL 放到 2 Hz；高动态场景还要再宽，但宽到一定程度噪声就压不住了。带宽改窄会让跟踪更平滑，但牵引范围变小，捕获给出的多普勒初值如果偏了 200 Hz 以上，窄带宽拉不回来，表现为频率一直朝一个方向缓慢漂。

积分周期也不是越长越好。B1C 主码周期 10 ms，数据分量的电文符号率是 1000 sps，一个主码周期里有 10 个符号，跨符号积分会互相抵消。导频分量没有这个问题，可以老老实实积 10 ms。代码里如果两个分量都跟踪，数据通道的积分周期要压到 1 ms 量级，导频通道保持 10 ms，最后在观测量层面把两个通道的结果合起来。

4.3 比特同步、帧同步与电文解析

跟踪出来的I_P是带符号的累加值，符号里承载着电文比特，但比特边界在哪里还不知道。比特同步就是在若干个可能的边界位置里找出正确的那个。直方图法最省事：把一个比特周期内的I_P值按位置分箱累加，能量最大的那个箱就是比特边界。

% 比特同步：直方图法 % 假设一个电文符号跨 Nsym 个主码周期，每个周期输出一个 I_P Nsym = cfg.symPerCode; % B1C 数据分量约为 1 ms 符号 / 10 ms 主码 acc = zeros(1, Nsym); cnt = zeros(1, Nsym); for k = 1:numel(trk.I_P) bin = mod(k-1, Nsym) + 1; acc(bin) = acc(bin) + abs(trk.I_P(k)); cnt(bin) = cnt(bin) + 1; end acc = acc ./ max(cnt, 1); [~, symBoundary] = max(acc); % 比特边界所在位置

找到边界之后，把每个符号周期内的I_P按边界对齐相加，得到符号判决值，取符号位就是比特。这里有一个容易忽略的点：极性可能整体反了。I_P的正负取决于本地载波和信号载波的相位关系，如果本地载波取了共轭，整串比特会全部翻转。解决办法是在帧同步阶段用已知的帧头或前导码做相关，相关峰为正说明极性对，为负就把整串比特取反。

帧同步之后是电文解析。B-CNAV1 和 B-CNAV2 的帧结构、纠错编码都不一样，解析函数分开写。B2a 的 B-CNAV2 用了 LDPC 编码，解析前要先做译码，译码失败就丢弃这一帧，不要硬着头皮往下解，否则星历参数里会混进垃圾值，最后定位结果飞出几百公里都算正常。星历解析出来之后至少要做一次合理性检查：轨道半长轴应该在 26000 km 量级，倾角在 55 度左右，偏心率是个小数，任何一项对不上就说明帧对齐或者译码出了问题。

电文解析的产出是每颗星的星历、钟差参数、健康标志。健康标志一定要用，不健康的卫星即使观测量质量很好也要剔掉，这是最省事的一道保险。

5. 定位解算与常见问题排查：5 个在实际调试中遇到的坑

5.1 从伪距到 PVT 的最小二乘

伪距的生成依赖两个时刻。接收时刻从采样点索引和采样率算出来，发射时刻从电文里的参考时间加上码相位折算，两者相减乘光速就是伪距。代码里用一个统一的t_rx数组管理接收时刻，所有通道共用同一套时间基准，这是保证通道间一致性的前提。

function [pos, dtr, nIter] = ls_pvt(satPos, pr, satClk, cfg) % satPos : N x 3，ECEF 卫星位置，单位米 % pr : N x 1，伪距，单位米，已扣除卫星钟差 % satClk : N x 1，卫星钟差等效距离，单位米 % cfg : 初始位置与迭代上限 c = 299792458; x = [cfg.initPos(:); 0]; % [x y z c*dt] nIter = 0; for it = 1:cfg.maxIter nIter = it; rho = sqrt(sum((satPos - x(1:3).').^2, 2)); pred = rho + x(4) - satClk; H = [-(satPos - x(1:3).')./rho, ones(size(rho))]; dx = (H.' * H) \ (H.' * (pr - pred)); x = x + dx; if norm(dx(1:3)) < 1e-3 break; end end pos = x(1:3); dtr = x(4) / c; end

这段代码本身不复杂，难的是前置条件。卫星位置必须和伪距在同一时刻，卫星在信号发射时刻的位置要用迭代的方式逼近，直接用接收时刻的位置会引入几十米的误差。另外至少要 4 颗卫星才能解出三维位置加钟差，只有 3 颗的时候要么固定高度，要么放弃。

迭代收敛的判据用位置增量小于 1 mm 就够了，但初始位置不能给得太离谱。初始位置偏差几百公里

本文还有配套的精品资源，点击获取