简介面向卫星导航与MATLAB信号处理学习者这份压缩包完整实现了北斗B1C信号从捕获、跟踪到定位解算的闭环流程覆盖信号仿真、FFT捕获、PLL/DLL跟踪、BPSK解调与最小二乘/卡尔曼定位等关键环节可直接对照学习或作为工程改造基础。包内共51个文件以42个m脚本为主辅以9个mat数据文件整体约25.56MBm脚本按捕获、跟踪、帧同步、解码、PVT解算等功能模块划分mat文件用于保存中间信号与结果便于分步验证。已有530人学习下载。通过逐行运行代码读者不仅能理解B1C扩频码生成、BOC调制及多星定位原理还能掌握从原始中频数据到最终位置输出的完整处理链是入门北斗接收机算法的高价值参考资料。1. 北斗 B1C 信号的从捕获到定位这套 MATLAB 流程能干嘛从射频采样数据直接走到经纬度中间要跨越 BOC 副载波解调、扩频码同步、电文帧解析和最小二乘解算四道坎。这套 MATLAB 工程把四道坎全部跑通文件按my_acquisition、B1C_tracking、B1C_frame_syn、B1C_PVT的链路组织既能看到每颗卫星的捕获峰也能在plotVisiblestars.m里画出可见星和最终定位点。适合两类人一类是刚接触北斗三号 B1C 信号处理、想用真实测距码结构验证捕获跟踪算法的学生另一类是在 GPS L1 C/A 流程上已经很熟、想快速迁移到北斗 B1C 的工程师。拿到代码后建议按signal_generate.m生成仿真数据再对照test0111.m里保存的中间结果逐段执行避免一上来就被多普勒搜索范围吓住。2. B1C 信号模型与测距码构造BOC 副载波与 Weil 码生成在跑B1C_acquisition.m之前先要把信号模型说清楚。北斗三号 B1C 的中心频率是 1575.42 MHz码率 1.023 Mcps主码长度 10230这些参数比 GPS L1 C/A 的 1023 码片整整长了 10 倍直接决定了捕获搜索窗的大小。B1C 不是单一 BPSK 调制而是把数据分量和导频分量分开调制工程包里Boc11.m、Boc61.m是副载波发生器ScWeil.m、PWeil.m负责扩频码legendre.m提供 Legendre 序列基础。2.1 B1C 的数据分量和导频分量B1C 的一个关键设计是数据通道和导频通道共享同一个载波但功率分配不一样。导频分量没有电文调制可以连续跟踪能量更高所以B1C_tracking.m里主要锁导频数据分量只在帧同步之后用来解电文。对于捕获来说用导频分量做相关峰检测更稳因为不会出现电文比特翻转导致的相关损失。代码包里的my_acquisition_data.m和my_acquisition_data_inac.m分别对应数据路和导频路的采集处理后者通常更容易出峰。2.2 Boc11.m、Boc61.m 与副载波生成BOC(1,1) 的含义是副载波频率 1.023 MHz、码率 1.023 Mcps每个码片内副载波翻转一次BOC(6,1) 则是副载波频率 6 × 1.023 MHz每个码片内翻转 12 次用来把导频频谱推向两侧从而与数据分量实现频域复用。一个典型的副载波生成函数长这样function sc boc_subcarrier(fs, fc, m, n) % fs: 采样率(Hz), fc: 码率(Hz), m/n 表示 BOC(m,n) % 返回一个码片内的副载波采样序列 fsc m * 1.023e6; Ts 1/fs; Tchip 1/fc; Nchip round(Tchip / Ts); t (0 : Nchip-1) * Ts; sc sign(sin(2 * pi * fsc * t)); end这段代码用正弦符号函数近似方波副载波Boc11.m调用时传m1, n1Boc61.m传m6, n1。注意fsc和fc的关系决定了每个码片内副载波翻转次数如果采样率不是两个频率的公倍数Nchip会四舍五入导致频谱上出现轻微非对称工程上更稳妥的办法是先生成一个整数码片周期的副载波再在相关器里按小数时延插值。2.3 Weil 码legendre.m、PWeil.m 与 ScWeil.mB1C 的主码不是 m 序列而是 Weil 码。Weil 码由长度为素数的 Legendre 序列与其循环移位序列模 2 相加得到B1C 用长度 10223 的 Legendre 序列作为基础再补足到 10230 码片。legendre.m负责生成 Legendre 符号序列PWeil.m根据 Weil 索引挑选两个移位相位并做异或ScWeil.m处理 B1C 子码拼接。下面的简化代码展示了核心逻辑% 理解 Weil 码生成的最小实现 p 10223; % Legendre 序列长度 leg zeros(1, p); for k 1 : p-1 if mod(k, 2) 1 % 实际用二次剩余判断这里仅为示意 leg(k1) 1; end end WeilIndex 120; % 不同 PRN 对应不同 Weil 索引 seq xor(leg, circshift(leg, WeilIndex));这里circshift的延迟量WeilIndex不能随便选。北斗公开 ICD 会给每个 PRN 指定两组 Weil 索引一组用于数据分量、另一组用于导频分量直接影响自相关旁瓣。工程包里test_legendre.m应该就是用来验证这批索引是否正确的跑完看主峰和最大旁瓣比正常情况下旁瓣抑制要明显优于 m 序列。2.4 参数表与 B1C_code_sample.m 自相关验证关键参数可以整理成表方便跟踪环节直接引用参数数值说明载波频率1575.42 MHzB1C 中心频点码率1.023 Mcps主码码率主码长度102301 ms 码长数据分量调制BOC(1,1)含电文导频分量调制QMBOC(6,1,4/33)导频连续分量子码结构主码 子码ScWeil.m负责拼接使用B1C_code_sample.m可以直接生成指定 PRN 的测距码然后用自相关函数看波形code B1C_code_sample(1, P); % PRN 1 导频主码 [acf, lags] xcorr(code, code, 200); [~, mid] max(acf); acf_norm acf / acf(mid); plot(lags, acf_norm);这段代码里xcorr只截取 ±200 码片主要看主峰周围旁瓣是否平坦。B1C_code_sample.m第二个参数传P表示导频分量传D就是数据分量。若旁瓣超过 0.3先怀疑ScWeil.m里面的截断位置错了不要急着进捕获流程。3. 捕获阶段的 FFT 相关实现my_acquisition 脚本与门限确定捕获的目的是给出每颗可见卫星的粗略码相位和多普勒频率。B1C 主码长度 10230直接在时域滑动相关在 MATLAB 下会非常慢所以工程包用 FFT 循环相关把复杂度从 O(N²) 降到 O(N log N)。my_acquisition.m、my_acquisition2.m是两个版本的实现my_acquisition_data_inac.m专门处理导频数据。运行顺序推荐test_data.m先生成中频数据再跑任意一个 acquisition 函数跑完之后用Final_result.mat里的伪距粗值核对一下能快速判断捕获峰是否选偏。3.1 捕获前先确认信号功率和采样率不要一上来就搜多普勒先把中频采集数据的频谱画出来确认目标信号落在预设频带内。B1C 的双边带宽约 30 MHz采样率如果低于 30 Msps副载波高频分量会混叠。以下代码用于观察频谱fs 62e6; % 采样率工程包常见设置 fIF 12e6; % 中频频率 spec abs(fftshift(fft(signal))); f linspace(-fs/2, fs/2, length(spec)); plot(f/1e6, 20*log10(spec eps));参数说明fs决定能表示的频率范围fIF决定频谱中心位置。如果signal是复数中频FFT 之后只会在fIF附近出现一个明显的峰如果看到镜像对称的两个峰说明数据被当成实数处理了会损失一半信噪比。对 B1C 来说频谱形状应呈现 BOC 调制特有的双峰而不是像 GPS C/A 那样的单三角峰。3.2 FFT 相关搜索码相位和多普勒捕获的核心是把本地产生的 B1C 导频码与输入信号做循环相关。多普勒频率从 -5 kHz 扫到 5 kHz步进一般取 500 Hz。一个最小实现如下function [code_phase, doppler_freq] acquire_B1C(signal, code, fs, fIF) N length(signal); codeFreq fft(code, N); % 本地码 FFT f_step 500; for df -5000 : f_step : 5000 n 0 : N-1; carrier exp(1j * 2 * pi * (fIF df) * n / fs); x signal .* carrier; corr ifft(fft(x) .* conj(codeFreq)); [peak, idx] max(abs(corr)); % 记录本次多普勒下的峰值和位置 doppler_freq df; code_phase idx; end这段代码的循环变量df是相对中频的多普勒偏移。codeFreq用N点 FFT要求code能循环扩展到N个采样点当输入信号长于一个码周期时最好先截取 1 ms 数据避免跨码周期相关把峰值抹平。idx是采样点索引换算成码片时要用idx / fs * 1.023e6这一步漏掉的话跟踪环节的初始码相位会差几十个码片。3.3 门限设计与峰值确认相关峰不能只看最大值还要看它比噪声底高多少。常用做法是取相关结果绝对值的均值乘一个系数作为门限峰值超过门限才判定为捕获成功。下面是一段门限判断逻辑noise_floor mean(abs(corr)); threshold 2.5 * noise_floor; if peak threshold % 记录候选卫星和对应多普勒、码相位 end参数表如下参数推荐范围说明多普勒步进250~500 Hz1 ms 相干积分下 500 Hz 损失约 0.9 dB相干积分长度1~10 ms导频分量无电文可做 10 ms非相干累加次数5~20 次弱信号场景取大值门限倍数2~3 倍噪声均值过低会大量虚警ACF.m在捕获阶段用来画相关函数曲线。正常情况主峰旁边至少有两个副峰间距等于 BOC 副载波相关间距如果副峰高度接近主峰多半是副载波相位没对齐而不是信号本身有问题。捕获完成后把码相位存下来作为跟踪环路的初始值。3.4 捕获失败时先查这三处第一查本地码是否从码片边界开始。B1C_code_sample.m输出的是完整码序列但被fft扩展后必须保证第一个采样点对应码片起始沿否则相关峰会偏一个采样点。第二查多普勒搜索范围。静止接收机 ±5 kHz 通常够用但signal_generate.m如果设置了高动态多普勒可能超过 ±8 kHz这时步进取 250 Hz 才能保证捕获到主峰。第三查信号直流偏置。signal里如果混入直流分量FFT 相关会在零频位置产生一个假峰且所有多普勒搜索都指向同一个df0位置。先做去直流处理再重复捕获。4. 跟踪阶段的 PLL/FLL/DLL 环路B1C_tracking.m 与 lock_detector捕获给出的码相位和多普勒只够开环使用要想稳定解出电文位和伪距必须进入闭环跟踪。工程包里的B1C_tracking.m和B1C_tracking2.m是两版跟踪主循环B1C_PLL.m、B1C_FLL.m、B1C_DLL.m分别实现锁相环、锁频环和码延迟环lock_detector.m输出锁定指示。运行顺序建议先看test0111.m它会加载ch0111.mat的通道数据逐历元更新环路状态跟踪过程中如果lock_detector一直为 0不要调滤波器带宽先检查初始码相位是否来自捕获结果。4.1 环路组合和各自职责载波环用 FLL 先拉大频偏再用 PLL 精确锁相位码环用 DLL 跟踪码相位延迟。B1C 跟踪里还有个特殊点导频通道没有电文PLL 可以直接用二象限反正切鉴相不会出现 180° 相位模糊。B1C_tracking.m通常把 FLL 输出作为 PLL 的预滤波DLL 则独立使用早迟相关器。三个环路的带宽和更新率不一样分开调优于一起调。环路输入输出推荐带宽FLL同相信号 I、正交信号 Q多普勒频率误差10~20 HzPLLI/Q 积分结果载波相位误差5~15 HzDLLE/P/L 三路相关值码相位误差0.5~2 Hz4.2 载波环路实现B1C_PLL.m 和 B1C_FLL.mPLL 常见做法是取即时支路的Ip和Qp用atan(Qp/Ip)得到相位差。代码如下function phi B1C_PLL(Ip, Qp) % Ip/Qp: 即时支路相干积分各一个标量 phi atan(Qp / Ip); end用二象限反正切而不是atan2(Qp, Ip)是为了避免电文比特翻转导致相位跳变B1C 导频通道没有数据但数据通道跟踪时仍需保留这一约束。phi超过 ±π/4 就说明环路失锁需要重新回到捕获阶段。FLL 的鉴频器一般用交叉积。给出代码function w B1C_FLL(cross, dot) % cross I_k * Q_{k-1} - Q_k * I_{k-1} % dot I_k * I_{k-1} Q_k * Q_{k-1} w cross / (dot eps); end这里cross和dot来自前后两个历元的积分值w是归一化频率误差。FLL 环路滤波器输出直接纠正载波 NCOPLL 输出只负责残余相位。工程包里B1C_tracking2.m会把 FLL 和 PLL 串联成二阶环路B1C_PLL.m返回的相位误差经过滤波器后加到B1C_FLL.m输出的频率上。4.3 码环实现B1C_DLL.m 与码相位误差DLL 使用超前、即时、滞后三路相关器。超前码和滞后码相对即时码偏移半个码片BOC 调制的 B1C 信号建议偏移更小例如 0.1 个码片因为 BOC 自相关主峰更窄。误差公式为function error B1C_DLL(E, L) % E: 超前支路相关值, L: 滞后支路相关值 error (E - L) / (E L eps); end这个归一化误差不受信号幅度影响适合直接输入环路滤波器。如果不做归一化信号变强时误差也会变大环路增益会随载噪比抖动。B1C_DLL.m里还可以加入早迟间隔参数跟踪开始用 0.3 码片稳定后缩到 0.1 码片。4.4 环路滤波器参数与初始化跟踪启动时载波 NCO 初始频率设为捕获到的中频 多普勒码 NCO 初始相位设为捕获到的码相位。带宽参数直接影响动态响应和噪声FLL 带宽取 12 Hz 时能容忍较大动态但噪声也大PLL 带宽取 8 Hz 时热噪声引起的相位抖动约 1° 量级。调节时先固定码环带宽再单独调载波环。lock_detector.m通常用相干积分后的信号功率估算载噪比大于 32 dB-Hz 才认为锁定。4.5 跟踪结果验证E_data.mat 和 E.mat跟踪完成后E_data.mat和E.mat保存了每个历元的即时支路积分结果。检查里面的Ip符号是否稳定如果符号跳变说明 B1C 数据通道的电文位在翻转这是正常现象如果Qp始终较大而Ip接近零说明 PLL 没有锁相应查看B1C_PLL.m的鉴相器输出是否一直在边界跳。test0111.m跑完可以打印lock_detector结果用mean(locked)统计锁定率低于 90% 时优先检查环路滤波器系数而不是信号功率。5. 帧同步、电文解码与伪距提取B1C_frame_syn.m 到 B1C_decode.m跟踪环路稳定之后数据通道输出的比特流还只是一串符号必须找到帧头、解出导航电文才能得到卫星星历和时钟参数。工程包里的B1C_frame_syn.m负责帧同步B1C_decode.m负责按 ICD 的电文格式解析crc_correct.m做校验bin_comp2dec.m把二进制补码转成十进制。N_data.mat和N.mat保存了解码后的导航数据B1C_time.m负责维护 B1C 系统时间。5.1 B1C 电文帧结构与同步头B1C 电文采用 B-CNAV1 格式基本帧长 1800 符号符号速率 100 sps等效 18 秒一帧。同步头固定为 ICD 规定的码型B1C_frame_syn.m做的事情就是在比特流里滑动搜索这个固定序列。一个简化版实现如下function offset B1C_frame_syn(bits) % bits: 解调后的 0/1 序列 syn [1 0 0 0 1 0 1 1]; % 示例同步头实际以 ICD 为准 len length(syn); for i 1:length(bits)-len if isequal(bits(i:ilen-1), syn) offset i; return end end offset -1;逻辑说明syn序列必须与发射端严格一致出现一个比特翻转就无法同步。offset是同步头在比特流中的起始位置后续所有字段都从offset开始按位偏移读取。如果返回 -1先检查B1C_tracking2.m输出的符号是否发生了极性反转B1C 数据通道的符号极性错误会导致同步头永远找不到。5.2 电文解码与 CRC 校验B-CNAV1 电文里大量字段是有符号整数用二进制补码表示。bin_comp2dec.m的典型实现如下function val bin_comp2dec(bits) % bits: 0/1 行向量最高位为符号位 b double(bits); if b(1) 1 raw bi2de(b(2:end), left-msb); val raw - 2^(length(b)-1); else val bi2de(b(2:end), left-msb); end end参数说明bi2de是通信工具箱函数没有工具箱时可用bin2dec(num2str(bits(2:end)))替代。B1C_decode.m读字段时要注意位宽例如星历参数crs是 16 位toc是 11 位位宽读错会造成数量级错误。crc_correct.m对整帧做 CRC 校验校验不过的电文不能进定位解算否则B1C_satpos.m算出的卫星位置会出现几十公里的跳跃。5.3 从码相位到伪距的提取伪距由两部分组成跟踪环路的码相位给出小数码片部分帧同步和子帧计数给出整毫秒部分。简化写法如下c 299792458; code_freq 1.023e6; % code_phase_unit: 码片integer_code: 整码周期数 pseudo_range (integer_code * 10230 code_phase_unit) * c / code_freq;这里integer_code是接收时刻与信号发射时刻之间的整码周期数code_phase_unit由B1C_DLL.m输出的码相位误差累加得到。两组结果相加前要统一单位integer_code没有小数部分时伪距误差可能达 300 km这也是为什么B1C_frame_syn.m的同步结果必须紧跟在跟踪后面中间不能隔历史数据。B1C_time.m把信号接收时刻换算出 B1C 系统时间避免各通道之间出现毫秒级偏差。5.4 中间数据文件与调试顺序Final_result.mat、base0111.mat、ch0111.mat、result0111.mat是不同阶段的存档。建议调试顺序先跑test0110.m生成基带数据再跑test0111.m完成捕获跟踪最后用pvt0111.m解算位置。如果某颗卫星的电文解不出来直接打开E_data.mat看对应通道的比特流定位问题在前端还是解调环节。GPS_CA_receiver.m可以当作对照组同样结构跑一遍 GPS L1 C/A能把 B1C 特有的 Weil 码和 BOC 副载波问题隔离出来。下面这张表对应关系很关键文件内容用途E_data.mat跟踪历元数据检查环路的即时支路输出N_data.mat解码后的导航数据给B1C_satpos.m使用Final_result.mat捕获和跟踪最终结果定位解算输入result0111.mat解算后的位置结果与最终结果对比6. 定位解算与结果验证leastSquarePos、坐标转换和可见星图伪距和星历都齐了以后本体位置由B1C_PVT.m和pvt0111.m完成。工程包里自带leastSquarePos.m按加权最小二乘迭代解四元方程组。B1C_satpos.m和satpos.m负责卫星位置计算iono_correct.m、tropo.m修正电离层和对流层误差e_r_corr.m处理地球自转改正。最后用cart2geo.m、togeod.m把 ECEF 坐标转成经度、纬度和高度。6.1 加权最小二乘 PVT 解算观测方程组是y G * x ε其中x是位置增量和接收机钟差。核心迭代如下G(:,1:3) -(sat_pos - pos_est) ./ range; G(:,4) 1; x (G * W * G) \ (G * W * y); pos_est pos_est x(1:3); clock_bias clock_bias x(4);G是几何矩阵range是估计几何距离W是伪距噪声方差的倒数。卫星仰角越低伪距噪声越大W的对应权重就越小如果不加权低仰角卫星会把定位结果拉偏好几米。6.2 坐标转换与误差修正cart2geo.m和togeod.m都是 ECEF 转大地坐标区别在于后者用迭代算法支持椭球参数不同。iono_correct.m通常用 Klobuchar 模型tropo.m用 Hopfield 或 Saastamoinen 模型参数取自B1C_decode.m解出的星历和电文辅助参数。e_r_corr.m补偿信号传播期间地球旋转造成的卫星位置偏移幅度在几百米量级忽略它定位结果会在东向上出现系统性偏置。6.3 用可见星图和中间结果验证一致性跑完pvt0111.m后用plotVisiblestars.m画出参与解算卫星的仰角和方位角。正常情况下可见星会分布在多个方向且都在仰角 10° 以上。再与result0111.mat比较位置差小于 1 m 说明链路连贯如果差了几个公里重点查B1C_satpos.m的星历参考时间是否已经过了有效期或者B1C_time.m给出的系统时间和伪距之间不是同一个观测历元。对比时注意伪距单位要先统一成米B1C_time.m返回的钟差单位是秒直接拿秒去减伪距会造成 1e-7 量级的偏差。本文还有配套的精品资源点击获取
