簡介本資源是一套面向衛(wèi)星導(dǎo)航信號處理初學(xué)者與通信專業(yè)學(xué)生的MATLAB仿真實(shí)踐材料聚焦北斗B1I頻段擴(kuò)頻碼生成、信號調(diào)制及接收解調(diào)全流程建模解決理論學(xué)習(xí)中缺乏可運(yùn)行實(shí)例、難以理解擴(kuò)頻同步與抗干擾機(jī)制的痛點(diǎn)。壓縮包為7KB的ZIP文件內(nèi)含MATLAB腳本與Simulink模型如BD2B1_SoftReceiver等核心文件涵蓋擴(kuò)頻碼序列生成、正弦載波調(diào)制、加性高斯白噪聲信道模擬、滑動相關(guān)碼同步、匹配濾波解擴(kuò)等關(guān)鍵模塊代碼結(jié)構(gòu)清晰、注釋完整便于逐層調(diào)試與原理驗(yàn)證。已有264人下載學(xué)習(xí)適合課程設(shè)計(jì)、畢業(yè)設(shè)計(jì)或北斗系統(tǒng)入門實(shí)驗(yàn)使用讀者可直接運(yùn)行復(fù)現(xiàn)B1I信號從生成到解調(diào)的全過程掌握偽隨機(jī)碼特性、擴(kuò)頻增益計(jì)算、時域同步實(shí)現(xiàn)等核心知識點(diǎn)并為后續(xù)多徑信道建模與定位算法開發(fā)奠定基礎(chǔ)。 寫這個項(xiàng)目的初衷其實(shí)很簡單北斗B1I信號的捕獲和跟蹤是所有做衛(wèi)星導(dǎo)航接收機(jī)的人繞不開的第一道門檻。B1I是北斗二號和北斗三號都在播發(fā)的公開服務(wù)信號頻率1561.098 MHz碼速率2.046 Mcps碼長2046。相比GPS的C/A碼B1I測距碼更長、碼率更高捕獲難度也略大一些。我最初在MATLAB里做這個仿真時最大的感受就是“資料不少但能直接跑通的完整鏈路不多”——很多文章只講了碼生成或者只講了捕獲很少把發(fā)射端擴(kuò)頻調(diào)制、信道加噪、接收端解擴(kuò)這一整條鏈路串起來說清楚。所以我決定把自己調(diào)試通過的一套流程整理出來從擴(kuò)頻碼生成到接收端相關(guān)解擴(kuò)逐步拆解順便把踩過的坑也一并交代清楚。這篇內(nèi)容適合三類人剛接觸衛(wèi)星導(dǎo)航信號處理、想搞明白B1I測距碼到底怎么生成的初學(xué)者已經(jīng)在做導(dǎo)航信號仿真、但在捕獲或解擴(kuò)環(huán)節(jié)卡殼的研究生以及需要用MATLAB快速搭建信號級仿真鏈路的工程師。讀完你至少能實(shí)現(xiàn)三件事生成任意PRN號的北斗B1I測距碼、完成B1I信號的擴(kuò)頻調(diào)制并驗(yàn)證頻譜特性、實(shí)現(xiàn)一個完整的相關(guān)解擴(kuò)接收流程并正確還原導(dǎo)航數(shù)據(jù)。1. 項(xiàng)目背景與整體設(shè)計(jì)思路1.1 北斗B1I信號結(jié)構(gòu)解析先把這個信號的基本盤說清楚。北斗B1I信號的載波頻率是1561.098 MHz這個頻率不是隨便定的它落在L頻段適合電離層延遲修正和與GPS L1頻段的兼容性設(shè)計(jì)。信號結(jié)構(gòu)上B1I由I支路和Q支路組成但民用上我們最關(guān)心的是I支路因?yàn)閷?dǎo)航電文和測距碼都在I支路上調(diào)制而Q支路沒有公開的測距碼。本項(xiàng)目聚焦的“B1I路”指的就是這個I支路。I支路的信號表達(dá)式可以寫成[ S_{B1I}(t) A \cdot D(t) \cdot C(t) \cdot \cos(2\pi f_0 t \varphi) ]其中 (D(t)) 是導(dǎo)航電文數(shù)據(jù)碼碼速率50 bps(C(t)) 是測距碼序列也就是我們要生成的擴(kuò)頻碼(f_01561.098) MHz是載波頻率(\varphi) 是初始載波相位。整個發(fā)射端的任務(wù)本質(zhì)上就是把 (D(t)) 和 (C(t)) 做模二加再調(diào)制到載波上。這里有個容易忽略的細(xì)節(jié)B1I的測距碼碼長是2046個碼片chip碼速率2.046 Mcps所以一個完整碼周期持續(xù)1毫秒。而導(dǎo)航電文比特率50 bps一個電文比特持續(xù)20毫秒也就是一個數(shù)據(jù)比特內(nèi)包含20個完整碼周期。這個“20倍關(guān)系”在后面設(shè)計(jì)解擴(kuò)累加器時非常關(guān)鍵——每20個碼周期做一次相干累加正好對應(yīng)一個數(shù)據(jù)比特的判決。1.2 為什么選B1I做擴(kuò)頻通信仿真很多教材用GPS C/A碼做例子因?yàn)榇a長1023、結(jié)構(gòu)簡單、資料多。但B1I碼其實(shí)更適合作為擴(kuò)頻通信的教學(xué)案例原因有三第一B1I碼序列更長2046互相關(guān)特性更接近隨機(jī)序列能更直觀地體現(xiàn)擴(kuò)頻碼“用戶區(qū)分”和“抗干擾”的本質(zhì)。用B1I做互相關(guān)分析時不同PRN號之間的隔離度明顯比C/A碼更容易觀察到。第二B1I碼的生成結(jié)構(gòu)更有代表性。它采用Gold碼衍生結(jié)構(gòu)由兩個m序列異或生成主碼寄存器為11級但截?cái)嗔藘蓚€序列的末位得到了2046長的碼。這種“截短”的設(shè)計(jì)思路在工程上很常見——通過調(diào)整相位選擇網(wǎng)絡(luò)可以從同一組寄存器得到不同PRN號的碼序列非常巧妙。第三實(shí)際工程中B1I的捕獲接收機(jī)設(shè)計(jì)和GPS有差異把仿真做扎實(shí)了你在處理真實(shí)中頻數(shù)據(jù)或硬件平臺時能少踩很多坑。GPS C/A碼的捕獲算法直接套到B1I上往往要改相關(guān)器的積分長度和碼頻率參數(shù)。1.3 系統(tǒng)方案選型與MATLAB實(shí)現(xiàn)路線整套仿真在我的設(shè)計(jì)中劃分為五個模塊參數(shù)初始化、擴(kuò)頻碼生成、發(fā)射端擴(kuò)頻調(diào)制、信道加噪、接收端解擴(kuò)。這五個模塊不復(fù)雜但每個模塊都有可展開的細(xì)節(jié)。參數(shù)初始化定義采樣率 (f_s)、碼率 (f_c)、載波頻率 (f_0)、數(shù)據(jù)率、PRN號、信噪比等。擴(kuò)頻碼生成根據(jù)PRN號查表得到G2相位選擇用MATLAB生成2046長度的B1I測距碼。發(fā)射端擴(kuò)頻調(diào)制將導(dǎo)航電文比特復(fù)制擴(kuò)展與擴(kuò)頻碼異或再與載波相乘得到B1I中頻信號。信道加噪在信號上疊加高斯白噪聲模擬接收機(jī)前端的熱噪聲環(huán)境。接收端解擴(kuò)完成載波剝離、碼相關(guān)、累加判決恢復(fù)出導(dǎo)航數(shù)據(jù)比特。關(guān)于采樣率的選擇我的實(shí)際做法是碼率是2.046 MHz仿真中采樣率取4.092 MHz也就是每個碼片2個采樣點(diǎn)。這個折中方案兼顧了仿真速度和捕獲精度的需求。如果你要追求更精細(xì)的碼相位分辨率可以取12.276 MHz每個碼片6個采樣點(diǎn)但仿真時間會明顯增加。在學(xué)習(xí)和驗(yàn)證階段4.092 MHz完全夠用。2. 擴(kuò)頻碼生成原理與MATLAB實(shí)現(xiàn)2.1 B1I測距碼的生成機(jī)制B1I測距碼的本質(zhì)是Gold碼由兩個11級m序列 (G1) 和 (G2) 模二加生成。生成多項(xiàng)式分別是(G1)(1 x x^3 x^4 x^{11})(G2)(1 x x^2 x^3 x^4 x^5 x^8 x^9 x^{11})兩個序列的長度都是 (2^{11} - 1 2047)模二加后的Gold碼也是2047長。但B1I測距碼的實(shí)際碼長是2046怎么來的關(guān)鍵在“截短”——G2序列的相位選擇網(wǎng)絡(luò)輸出時從特定的兩個寄存器抽頭出來并且特意少了一個狀態(tài)使得組合后碼長為2046。嚴(yán)格來說B1I碼的生成方式相當(dāng)講究G1序列直接從第一個寄存器輸出G2序列則通過一個相位選擇器選擇兩個抽頭位置異或后輸出。選擇哪兩個抽頭由衛(wèi)星的PRN號決定。每個PRN號對應(yīng)一組唯一的抽頭組合這樣用同一組移位寄存器就能生成互相關(guān)性良好的多個不同碼序列。在MATLAB里實(shí)現(xiàn)時不需要真的構(gòu)建完整的Gold碼生成器——可以用相位選擇表直接構(gòu)造G2序列的初始狀態(tài)然后一次性生成整段碼。下面是具體的實(shí)現(xiàn)代碼。function b1i_code generateB1ICode(prn) % generateB1ICode 生成北斗B1I測距碼 % 輸入prn - 衛(wèi)星PRN號1~63 % 輸出b1i_code - 1x2046 雙極性碼1/-1 % 相位選擇表每個PRN對應(yīng)G2序列的兩個抽頭位置 % 這里列出部分PRN的抽頭完整表可參考北斗ICD文件 phase_sel [ 1, (1, 2); % PRN 1 2, (1, 3); % PRN 2 % ... 省略實(shí)際要查表 ]; % 根據(jù)PRN獲取抽頭 tap1 table_value(prn, 1); tap2 table_value(prn, 2); % G1初始狀態(tài)全1 g1 ones(1, 11); % G2初始狀態(tài)全1 g2 ones(1, 11); code zeros(1, 2046); for n 1:2046 % G1輸出為第1個寄存器 g1_out g1(1); % G2輸出為抽頭位置異或 g2_out xor(g2(tap1), g2(tap2)); % 碼輸出為G1和G2異或 code(n) xor(g1_out, g2_out); % G1反饋 fb1 xor(g1(11), xor(g1(4), xor(g1(3), g1(1)))); g1 [fb1, g1(1:10)]; % G2反饋 fb2 xor(xor(xor(xor(xor(xor(g2(11), g2(10)), g2(9)), g2(8)), g2(5)), g2(4)), g2(1)); g2 [fb2, g2(1:10)]; end % 轉(zhuǎn)雙極性 b1i_code 1 - 2 * code; end這段代碼中G1反饋的計(jì)算來自生成多項(xiàng)式 (1 x x^3 x^4 x^{11})所以抽的是第1、第4、第11個寄存器做異或。注意MATLAB里寄存器向量順序是左到右第一個元素是 (x^1)最后一個是 (x^{11})。G2反饋多項(xiàng)式是 (1 x x^2 x^3 x^4 x^5 x^8 x^9 x^{11})對應(yīng)的抽頭是1、2、3、4、5、8、9、11位。這個細(xì)節(jié)很多人會搞反生成出來的碼相關(guān)性不對排查問題時會很痛苦。2.2 G2序列相位選擇的關(guān)鍵細(xì)節(jié)相位選擇表是整個B1I碼生成的核心。每個PRN號對應(yīng)G2序列的兩個抽頭位置這個關(guān)系不是任意指定的而是經(jīng)過設(shè)計(jì)保證了不同PRN碼之間的互相關(guān)性能。完整表格在北斗ICD文件里我這里列出前幾個常用的PRN方便你驗(yàn)證代碼是否正確PRN抽頭位置1抽頭位置2112213314415516627注意這個抽頭位置的索引是按G2寄存器從左到右排列來數(shù)的第1位置即 (x^1) 的寄存器。實(shí)際應(yīng)用中你可以從ICD或開源代碼庫如GNSS-SDR里把完整表拷貝出來放到MATLAB常量數(shù)組里避免手輸錯誤。我當(dāng)年第一版代碼就是手敲表格結(jié)果第23號PRN的抽頭輸錯了一個數(shù)字導(dǎo)致那個衛(wèi)星的碼始終捕獲不上排查了大半天。所以強(qiáng)烈建議直接從可信源拷貝完整相位表不要手工錄入。驗(yàn)證碼生成是否正確的黃金標(biāo)準(zhǔn)是自相關(guān)/互相關(guān)函數(shù)。B1I碼理想的自相關(guān)函數(shù)在主峰處等于2046其余位置接近0互相關(guān)函數(shù)處處接近0。用MATLAB的xcorr函數(shù)可以快速驗(yàn)證。prn1_code generateB1ICode(1); prn2_code generateB1ICode(2); % 自相關(guān) [autocorr_val, lags] xcorr(prn1_code, prn1_code); center (length(autocorr_val) - 1) / 2 1; main_peak autocorr_val(center); % 應(yīng)為 2046 side_peak max(abs(autocorr_val([1:center-10, center10:end]))); % 旁瓣 % 互相關(guān) [mutual_corr, ~] xcorr(prn1_code, prn2_code); cross_peak max(abs(mutual_corr));正常情況下主峰2046旁瓣和互相關(guān)峰值都遠(yuǎn)小于主峰一般不超過100。如果算出來旁瓣接近主峰那基本可以斷定碼生成邏輯有誤優(yōu)先檢查反饋抽頭和相位選擇表。2.3 周期擴(kuò)展與多碼周期拼接實(shí)現(xiàn)生成2046碼片長度的碼后發(fā)射端還需要把碼周期擴(kuò)展因?yàn)橐粋€導(dǎo)航電文比特要持續(xù)20個碼周期。在MATLAB里這用repmat函數(shù)一行就能搞定。% 生成PRN1的碼 b1i_code generateB1ICode(1); % 擴(kuò)展20個碼周期 code_20ms repmat(b1i_code, 1, 20);然后和導(dǎo)航電文比特做異或前的處理這里有個數(shù)據(jù)類型轉(zhuǎn)換的細(xì)節(jié)測距碼是雙極性1/-1導(dǎo)航電文也轉(zhuǎn)成雙極性1代表數(shù)據(jù)位1-1代表位0那么擴(kuò)頻調(diào)制的過程就是兩個雙極性序列逐元素相乘。乘法和異或本質(zhì)等效1對0、-1對1但雙極性乘法在處理加噪、解擴(kuò)時更直觀因?yàn)楹罄m(xù)混頻和積分的公式都是基于實(shí)數(shù)乘法的。所以整個仿真鏈路我建議統(tǒng)一采用雙極性表示避免0/1和1/-1來回切換導(dǎo)致混亂。3. 發(fā)射端擴(kuò)頻調(diào)制仿真3.1 B1I信號生成的整體流程發(fā)射端的流程可以總結(jié)為四步生成導(dǎo)航電文序列。仿真中我經(jīng)常用隨機(jī)比特或者一段有規(guī)律的二進(jìn)制序列比如交替的0/1作為電文方便后續(xù)核對解調(diào)結(jié)果。將每一個電文比特重復(fù)20次對應(yīng)20個碼周期擴(kuò)展成長序列。電文擴(kuò)展序列與擴(kuò)展后的測距碼逐元素相乘得到擴(kuò)頻序列。擴(kuò)頻序列調(diào)制到載波上生成B1I中頻信號。用代碼實(shí)現(xiàn)就是% 參數(shù)配置 fs 4.092e6; % 采樣率 4.092 MHz fc 2.046e6; % 碼率 2.046 Mcps sps round(fs / fc); % 每碼片采樣點(diǎn)數(shù)這里為2 num_bit 10; % 仿真10個導(dǎo)航電文比特 total_ms num_bit * 20; % 總時間 200ms % 導(dǎo)航電文隨機(jī)數(shù)據(jù)比特 data_bits randi([0, 1], 1, num_bit); data_bipolar 1 - 2 * data_bits; % 測距碼 b1i_code generateB1ICode(1); % 擴(kuò)展到總碼片長度 code_chips repmat(b1i_code, 1, total_ms); % 電文擴(kuò)展每比特重復(fù)20次碼周期 data_expanded repelem(data_bipolar, 20 * length(b1i_code)); % 擴(kuò)頻序列 spread_seq data_expanded .* code_chips; % 上采樣到采樣域 spread_signal repelem(spread_seq, sps);這段代碼里有個細(xì)節(jié)要重點(diǎn)解釋repelem和repmat的區(qū)別。repmat是把整個序列整體復(fù)制比如repmat([1,2], 1, 3)得到[1,2,1,2,1,2]而repelem([1,2], 3)得到[1,1,1,2,2,2]。電文擴(kuò)展要用repelem因?yàn)槊總€比特要連續(xù)持續(xù)20個碼周期而不是交替跳躍。碼序列擴(kuò)展用repmat因?yàn)榇a周期本來就是連續(xù)的。搞反這兩個函數(shù)信號結(jié)構(gòu)就全錯了。3.2 載波調(diào)制與中頻信號生成有了擴(kuò)頻序列的采樣域表示接下來疊加載波。中頻頻率的選擇有個講究仿真里不需要太高中頻但要保證載波頻率和采樣率滿足奈奎斯特采樣定理。我經(jīng)常選用中頻 (f_{IF} 1.023) MHz這樣在4.092 MHz的采樣率下信號頻譜不會混疊且方便后續(xù)數(shù)字下變頻處理。注意這里數(shù)字中頻和真實(shí)射頻前端的模擬中頻不同仿真中純粹是為了模擬接收機(jī)的處理流程不必照搬真實(shí)硬件的中頻頻率。要實(shí)現(xiàn)載波調(diào)制先在時間軸上生成采樣點(diǎn)索引再構(gòu)造載波信號% 時間軸 t (0 : length(spread_signal) - 1) / fs; % 中頻載波 fIF 1.023e6; % 中頻 1.023 MHz carrier cos(2 * pi * fIF * t); % B1I中頻信號 b1i_if spread_signal .* carrier;生成后建議做一個關(guān)鍵檢查用pwelch或fft觀察信號的頻譜。B1I擴(kuò)頻信號的特點(diǎn)很明顯BPSK(2)調(diào)制的擴(kuò)頻信號頻譜是 (\text{sinc}^2) 形狀主瓣寬度為2倍的碼率也就是約4.092 MHz兩個第一零點(diǎn)位于 (f_{IF} \pm f_c)?？匆谎垲l譜就能確認(rèn)擴(kuò)頻調(diào)制是否正確——如果頻譜主瓣寬度不對大概率是采樣率和碼率的關(guān)系設(shè)錯了。3.3 信噪比定義與高斯噪聲疊加很多初學(xué)者在這里會疑惑擴(kuò)頻通信里信噪比到底怎么定義是每個碼片的信噪比還是每個數(shù)據(jù)比特的信噪比這個必須明確。在擴(kuò)頻系統(tǒng)中接收端解擴(kuò)的過程會把信號能量累積起來噪聲能量通過相關(guān)處理被抑制所以處理增益體現(xiàn)在碼片到數(shù)據(jù)比特的累積上。仿真中我習(xí)慣按“每個數(shù)據(jù)比特的信噪比 (E_b/N_0)”來設(shè)置噪聲然后用公式反算噪聲功率。具體做法是先計(jì)算信號功率 (P_s)通常是1因?yàn)槲覀冇玫氖请p極性±1序列然后根據(jù)給定的 (E_b/N_0) 計(jì)算噪聲功率譜密度[ N_0 \frac{E_b}{10^{(E_b/N_0)/10}} \frac{P_s \cdot T_b}{10^{(E_b/N_0)/10}} ]其中 (T_b) 是一個數(shù)據(jù)比特的持續(xù)時間20 ms。噪聲功率為 (P_n N_0 \cdot f_s)噪聲標(biāo)準(zhǔn)差就是 (\sqrt{P_n})。疊加噪聲的MATLAB代碼EbN0_dB 5; % 每比特信噪比 5 dB Tb 20e-3; % 比特持續(xù)時間 Ps mean(b1i_if.^2); % 信號功率 Eb Ps * Tb; N0 Eb / (10^(EbN0_dB / 10)); Pn N0 * fs; noise_std sqrt(Pn); % 加噪 rx_signal b1i_if noise_std * randn(size(b1i_if));這里randn生成的是標(biāo)準(zhǔn)正態(tài)分布噪聲乘上標(biāo)準(zhǔn)差就是所需功率的高斯白噪聲。有人會問為什么不直接用awgn函數(shù)awgn函數(shù)很方便但它需要指定信號功率和信噪比的相對關(guān)系在擴(kuò)頻場景下容易算錯“信號功率”到底是碼片級的還是比特級的。手動設(shè)置噪聲功率雖然多寫兩行代碼但每一步物理意義都清楚調(diào)參也直觀。4. 接收端解擴(kuò)仿真4.1 接收機(jī)處理流程框架接收機(jī)的任務(wù)是拿到含噪的B1I中頻信號后恢復(fù)出導(dǎo)航電文比特。整體流程為載波剝離、解除擴(kuò)頻碼相關(guān)、積分累加、判決。其中載波剝離需要本地復(fù)現(xiàn)一個與信號同頻同相的載波解除擴(kuò)頻需要本地生成一個碼相位對齊的測距碼。在仿真中我假設(shè)接收機(jī)已經(jīng)完成了信號捕獲和跟蹤——也就是說載波頻率、載波相位、碼相位都已知。這樣的簡化能讓我們聚焦于解擴(kuò)本身不用陷入捕獲和跟蹤環(huán)路的復(fù)雜細(xì)節(jié)。當(dāng)然實(shí)際接收機(jī)不可能事先知道這些參數(shù)需要靠捕獲和跟蹤來估計(jì)。在文章后半部分我會專門討論在沒有先驗(yàn)信息時如何用二維搜索實(shí)現(xiàn)捕獲以及跟蹤環(huán)的基本考慮。4.2 載波剝離與低通濾波載波剝離的核心是將接收信號乘以本地產(chǎn)生的cos分量和sin分量得到I路和Q路兩路基帶信號。這種結(jié)構(gòu)對應(yīng)的是正交解調(diào)它可以同時處理載波相位不確定性問題。具體做法% 本地載波假設(shè)頻率已知初始相位為0 local_carrier_cos cos(2 * pi * fIF * t); local_carrier_sin -sin(2 * pi * fIF * t); % 正交下變頻 i_baseband rx_signal .* local_carrier_cos; % I支路 q_baseband rx_signal .* local_carrier_sin; % Q支路下變頻之后信號包含基帶分量和 (2f_{IF}) 的高頻分量。一定要用低通濾波器把高頻分量濾掉否則后續(xù)相關(guān)累加時高頻分量會在碼周期內(nèi)積分出非零的殘留干擾判決。濾波器的設(shè)計(jì)我用的是最省事的方式——直接用一個平均濾波器或者簡單的等紋波低通截止頻率設(shè)在碼率以下。% 設(shè)計(jì)一個簡單的FIR低通濾波器 fcut 2.046e6 / 2; % 截止頻率 1.023 MHz fir_coeff designfilt(lowpassfir, PassbandFrequency, 0.8e6, ... StopbandFrequency, 1.2e6, PassbandRipple, 1, ... StopbandAttenuation, 40, SampleRate, fs); i_filtered filter(fir_coeff, 1, i_baseband); q_filtered filter(fir_coeff, 1, q_baseband);提示濾波器會引入群延遲如果后續(xù)要做碼相位精確對齊需要對濾波后的信號做相應(yīng)的延遲補(bǔ)償。仿真里如果本地碼相位按理想情況對齊這個延遲會造成幾碼片偏差需要在相關(guān)前手動對齊。4.3 相關(guān)解擴(kuò)與積分累加解擴(kuò)的核心操作就是相關(guān)運(yùn)算。將濾波后的基帶信號與本地測距碼相乘如果碼相位對齊則乘積中數(shù)據(jù)比特的極性會被完整保留如果碼相位沒對齊則乘積正負(fù)交替積分后趨近于0。這個特性正是擴(kuò)頻系統(tǒng)抗干擾和測距的基礎(chǔ)。先對基帶信號做碼片的匹配一個碼片采樣2點(diǎn)需要分成2點(diǎn)求和再映射到碼片% 將采樣點(diǎn)轉(zhuǎn)換為碼片 i_chips reshape(i_filtered, sps, []); i_chip sum(i_chips, 1); q_chips reshape(q_filtered, sps, []); q_chip sum(q_chips, 1);然后按碼周期整形逐周期與本地碼做相關(guān)code_len length(b1i_code); % 2046 num_cycle floor(length(i_chip) / code_len); i_chip_crop i_chip(1 : num_cycle * code_len); q_chip_crop q_chip(1 : num_cycle * code_len); % 按碼周期排列 i_mat reshape(i_chip_crop, code_len, []); q_mat reshape(q_chip_crop, code_len, []); % 逐周期與本地碼相關(guān) i_corr i_mat. * b1i_code.; % 每行是一個碼周期的相關(guān)結(jié)果 q_corr q_mat. * b1i_code.;相關(guān)結(jié)果 (i_corr) 和 (q_corr) 的長度等于碼周期數(shù)每個值對應(yīng)一個1毫秒的相關(guān)累積。如果碼相位對齊這些值應(yīng)該集中在復(fù)平面的某個方向——具體方向取決于載波相位和導(dǎo)航電文比特的乘積。4.4 數(shù)據(jù)比特判決與性能驗(yàn)證現(xiàn)在要把20個碼周期的相關(guān)結(jié)果累加得到一個數(shù)據(jù)比特的積分值。這里要非常注意如果載波相位沒有完全對齊到零I路和Q路的能量如何合成最穩(wěn)妥的做法是求相關(guān)結(jié)果的模值平方和% 把20個碼周期累加成1個比特 bit_len 20; i_bit reshape(i_corr, bit_len, []); q_bit reshape(q_corr, bit_len, []); % 相干累加 I_sum sum(i_bit, 1); Q_sum sum(q_bit, 1); % 判決量 decision I_sum .* sign(I_sum) Q_sum .* sign(Q_sum); % 或簡單地用幅度判決 amplitude sqrt(I_sum.^2 Q_sum.^2);數(shù)據(jù)符號的判決有兩種經(jīng)典方式若載波相位已經(jīng)跟蹤得很穩(wěn)殘余相位接近0判決量直接用 (I_sum) 的符號即可。(I_sum 0) 判為數(shù)據(jù)比特1否則判為0。若殘余載波相位不可忽略比如還沒有進(jìn)入相位鎖定就利用 IQ 兩路能量合成用 (\text{sign}(I_sum)) 判決同時計(jì)算 (\sqrt{I_{sum}^2 Q_{sum}^2}) 輔助判斷信號質(zhì)量。為了驗(yàn)證解擴(kuò)效果我通常畫三張圖第一張是接收信號的頻譜確認(rèn)加噪后信號仍可辨識第二張是相關(guān)結(jié)果的幅值序列可以看到每個碼周期相關(guān)峰的高度第三張是恢復(fù)的數(shù)據(jù)比特與原始數(shù)據(jù)的對比。最后直接統(tǒng)計(jì)誤碼率recovered_bits double(I_sum 0); ber sum(recovered_bits ~ data_bits) / num_bit;在 (E_b/N_0) 大于0 dB時B1I這類BPSK擴(kuò)頻系統(tǒng)應(yīng)該能達(dá)到很低的誤碼率。仿真中如果把 (E_b/N_0) 設(shè)為5 dB理論上誤碼率在 (10^{-3}) 量級以下但由于我們只有10個比特的仿真量可能一個錯誤都沒有如果要看統(tǒng)計(jì)特性需要增加比特?cái)?shù)到1000個以上。5. 常見問題與排查技巧實(shí)錄5.1 碼相位對齊錯誤導(dǎo)致的“鬼影相關(guān)峰”這是整個仿真里最容易踩的坑。我見過很多人修改了碼速率或采樣率后發(fā)現(xiàn)相關(guān)峰消失就開始懷疑濾波器、懷疑載波頻率、懷疑代碼邏輯。其實(shí)大多數(shù)情況下就是碼相位和采樣點(diǎn)之間對不齊了。原因在于(f_s / f_c 2) 時每個碼片正好2個采樣點(diǎn)這個整數(shù)倍關(guān)系保證了碼片邊界恰好落在采樣點(diǎn)上。但有些仿真里會把采樣率設(shè)成 (f_s 5) MHz為了滿足其他約束此時每個碼片不是整數(shù)個采樣點(diǎn)碼片邊界和采樣點(diǎn)之間就會有漂移必須通過插值或調(diào)整本地碼采樣策略來對齊。我的建議是如果暫時不需要精細(xì)的碼相位分辨率就把采樣率設(shè)為碼率的整數(shù)倍把問題先繞過去。排查方法在本地碼和接收信號做相關(guān)之前先把接收信號按 I/Q 解調(diào)后的基帶信號打印出來手工檢查碼片邊界是否和本地碼的起始位置對齊。如果發(fā)現(xiàn)偏差調(diào)整本地碼的起始索引即可。5.2 載波頻率偏差導(dǎo)致的相關(guān)能量衰減仿真中假設(shè)本地載波和信號載波同頻同相但實(shí)際接收機(jī)中多普勒頻移會破壞這個假設(shè)。當(dāng)載波頻率有偏差 (\Delta f) 時相關(guān)積分的輸出會附加一個 (\text{sinc}(\Delta f \cdot T_{\text{int}})) 的衰減因子。假設(shè)相干積分時間為1個碼周期1 ms頻率偏差1 kHz時(\text{sinc}(1)) 已經(jīng)明顯衰減。所以如果在仿真中加入多普勒頻移卻不修改接收端處理解擴(kuò)性能會急劇下降。最直觀的驗(yàn)證方式是以1 kHz為步進(jìn)畫一條相關(guān)峰高度隨頻率偏差變化的曲線你會看到典型的sinc形狀——這就是捕獲階段需要做頻率搜索的物理原因。學(xué)習(xí)階段可以先不加多普勒等基本流程跑通后再嘗試加入頻率偏移并做頻率搜索。5.3 誤碼率曲線與理論值的偏差調(diào)優(yōu)如果你把仿真比特?cái)?shù)增加畫出誤碼率隨 (E_b/N_0) 變化的曲線可能會發(fā)現(xiàn)仿真結(jié)果和理論BPSK誤碼率曲線有偏差。偏差來源無非三個低通濾波器引入了群延遲導(dǎo)致碼相位偏移相關(guān)增益下降。載波相位殘余不為零IQ合成策略不當(dāng)導(dǎo)致能量損失。碼片匹配之前沒有做正確的采樣點(diǎn)分配導(dǎo)致碼片能量被兩個采樣點(diǎn)“切碎”后沒有完整累加。我的排查順序是先用無噪聲信號跑一遍看能否完美恢復(fù)如果能再逐步加噪聲如果無噪聲都有誤碼那一定是處理鏈路存在邏輯錯誤絕口不提噪聲的影響。這是調(diào)試仿真系統(tǒng)的通用原則——先保證系統(tǒng)確定性正確再引入隨機(jī)性。5.4 常見問題速查表問題現(xiàn)象可能原因解決辦法自相關(guān)旁瓣接近主峰G2抽頭或寄存器方向錯誤對照ICD核對相位選擇表、檢查寄存器左到右順序相關(guān)峰忽大忽小碼相位未對齊或采樣非整數(shù)倍使用整數(shù)倍采樣率、調(diào)整本地碼起始索引加噪后完全無法解調(diào)信噪比設(shè)置錯誤或載波頻率偏差手動計(jì)算噪聲功率確認(rèn) (E_b/N_0) 定義低通濾波后波形畸變?yōu)V波器階數(shù)過高導(dǎo)致延遲過大檢查群延遲并補(bǔ)償或改用簡單平均濾波恢復(fù)數(shù)據(jù)整體翻轉(zhuǎn)載波相位約180度用IQ合成判決或加入差分編碼多PRN同時存在時串?dāng)_碼間互相關(guān)未被平均器平均掉增加相干積分時間或引入碼跟蹤環(huán)5.5 MATLAB數(shù)值精度與性能優(yōu)化建議最后給幾條工程層面的建議。第一盡量用向量化操作代替循環(huán)。上面代碼里碼生成用了循環(huán)這是正確的因?yàn)橐莆患拇嫫鞯姆答伿侵鹞灰蕾嚨臎]法向量化。但發(fā)射端的擴(kuò)頻、解擴(kuò)、積分累加都盡量用矩陣運(yùn)算MATLAB在矩陣化運(yùn)算上的效率遠(yuǎn)高于循環(huán)。第二大數(shù)據(jù)量時注意內(nèi)存管理。如果仿真1000個數(shù)據(jù)比特總信號長度是 (1000 \times 20 \times 2046 \times 2) 個采樣點(diǎn)約8100萬點(diǎn)單個雙精度數(shù)組就是640 MB。這種規(guī)模的數(shù)據(jù)一次性生成很容易把內(nèi)存打爆。對策是把數(shù)據(jù)分段處理——生成一段、處理一段、釋放一段。我常用timer或手動分段每段100個比特循環(huán)處理。第三濾波器設(shè)計(jì)要控制階數(shù)。designfilt默認(rèn)可能生成很高階的濾波器處理8100萬點(diǎn)數(shù)據(jù)時極慢?？梢杂胒ir1(64, Wn)這種固定階數(shù)的設(shè)計(jì)方式精度足夠速度也快得多。6. 進(jìn)一步擴(kuò)展信號捕獲與多衛(wèi)星場景6.1 從“已知碼相位”到“二維搜索捕獲”前面整個鏈路假設(shè)接收機(jī)知道碼相位和載波頻率這在驗(yàn)證解擴(kuò)原理時沒問題但它回避了接收機(jī)最核心的問題接收機(jī)一開始什么都不知道怎么確定碼相位和多普勒頻率這個問題的答案就是捕獲。捕獲的本質(zhì)是在二維平面上搜索碼相位維度上搜索2046個碼片頻率維度上以一定的步進(jìn)搜索可能的載波頻率偏差。搜索指標(biāo)是相關(guān)峰的大小。當(dāng)某個碼相位-頻率組合下相關(guān)峰顯著高于其他組合時就認(rèn)為捕獲到了信號。在MATLAB里實(shí)現(xiàn)常用的并行碼相位搜索算法function [peak_freq_idx, peak_code_idx, corr_plane] acquisition_parallel(rx_signal, prn, fIF, fs, fc, freq_search_range) % 并行碼相位捕獲 code generateB1ICode(prn); sps round(fs / fc); code_sampled repelem(code, sps); freq_steps -freq_search_range : 250 : freq_search_range; % 250Hz步進(jìn) corr_plane zeros(length(freq_steps), length(code_sampled)); for k 1:length(freq_steps) % 本地載波含頻率偏差 t (0:length(rx_signal)-1) / fs; local_osc exp(1j * 2 * pi * (fIF freq_steps(k)) * t); % 下變頻到復(fù)數(shù)基帶 baseband rx_signal .* local_osc; % 與本地碼做循環(huán)相關(guān)FFT法 spec fft(baseband) .* conj(fft(code_sampled)); corr ifft(spec); corr_plane(k, :) abs(corr); end end這里的核心是利用FFT實(shí)現(xiàn)循環(huán)相關(guān)復(fù)雜度從 (O(N^2)) 降到 (O(N \log N))這在2046個碼片長度上效果非常明顯。但要注意FFT相關(guān)是循環(huán)相關(guān)當(dāng)碼相位超出信號長度時會產(chǎn)生環(huán)繞效應(yīng)所以實(shí)際使用需要將信號長度擴(kuò)展為碼周期的整數(shù)倍并處理好邊界。6.2 多衛(wèi)星信號疊加的仿真驗(yàn)證北斗接收機(jī)從來不會只收一顆衛(wèi)星的信號。在真實(shí)場景中可見的北斗衛(wèi)星通常在8顆以上它們的信號在接收機(jī)前端被混疊在一起。擴(kuò)頻系統(tǒng)的多址能力使得接收機(jī)能夠用不同的本地碼分別捕獲不同的衛(wèi)星。在仿真中模擬多衛(wèi)星場景很簡單生成多路B1I信號不同PRN號、不同碼相位、不同功率疊加后送入接收機(jī)然后用不同PRN的本地碼分別做捕獲和解擴(kuò)。你會發(fā)現(xiàn)只要各衛(wèi)星信號功率沒有懸殊的差異解擴(kuò)后的誤碼率依然能保持。如果有一顆衛(wèi)星功率特別強(qiáng)近衛(wèi)星而另一顆特別弱遠(yuǎn)衛(wèi)星弱信號會被強(qiáng)信號的互相關(guān)噪聲淹沒這就是遠(yuǎn)近效應(yīng)。做這個擴(kuò)展實(shí)驗(yàn)?zāi)茏尅皵U(kuò)頻碼區(qū)分用戶”這個抽象概念變得特別具象。我當(dāng)時跑通兩顆衛(wèi)星疊加時才真正理解為什么GPS/北斗要精心設(shè)計(jì)擴(kuò)頻碼的互相關(guān)特性。6.3 從仿真走向硬件平臺的落地點(diǎn)軟件仿真驗(yàn)證通過后往硬件平臺遷移的過程中以下幾個點(diǎn)需要特別留意一是采樣率的調(diào)整。硬件平臺的中頻采樣率往往是16.368 MHz或更高和仿真里的4.092 MHz差很多。碼環(huán)和載波環(huán)的參數(shù)帶寬、增益都是和采樣率強(qiáng)相關(guān)的不能直接照搬。二是整型化處理。MATLAB里默認(rèn)用雙精度浮點(diǎn)硬件平臺為了速度和功耗通常用定點(diǎn)數(shù)實(shí)現(xiàn)。定點(diǎn)化的過程中要重點(diǎn)關(guān)注相關(guān)累加器的位寬——2046個碼片的累加如果再疊加20個碼周期的相干積分累加值動態(tài)范圍很大位寬不足會直接導(dǎo)致飽和失真。三是實(shí)時性約束。軟件仿真可以逐塊處理數(shù)據(jù)硬件接收機(jī)必須在一個碼周期內(nèi)完成相關(guān)運(yùn)算這對算法復(fù)雜度和實(shí)現(xiàn)效率提出了更高要求。好在MATLAB仿真時的參數(shù)和結(jié)構(gòu)驗(yàn)證能為硬件實(shí)現(xiàn)提供可靠的參考基準(zhǔn)。7. 個人的一點(diǎn)實(shí)操心得整套仿真做下來我最大的體會是北斗B1I的擴(kuò)頻仿真不是一個“跑通代碼就行”的事情每一個參數(shù)背后都有信號體制的約束。比如2046這個數(shù)它來自Gold碼截短的巧妙設(shè)計(jì)20倍的關(guān)系它對應(yīng)導(dǎo)航電文速率與碼周期之間的除法關(guān)系2倍采樣率的選擇則是接收機(jī)資源與性能的折中。把這些問題想清楚代碼只是它們自然的表達(dá)。調(diào)試過程中我建議你準(zhǔn)備一個“信號正確性檢查板”每完成一個模塊就驗(yàn)證一次碼生成后檢查自相關(guān)、擴(kuò)頻后檢查頻譜、加噪后檢查信噪比、解擴(kuò)后檢查誤碼率。模塊級驗(yàn)證遠(yuǎn)比整體調(diào)試高效——整體調(diào)試的問題往往一環(huán)套一環(huán)定位成本極高。寧可每個模塊多花十分鐘確認(rèn)也不要最后一口氣面對幾十個報錯或錯誤結(jié)果。最后再分享一個小技巧在仿真中故意制造一些極端情況比如把 (E_b/N_0) 拉到 -5 dB觀察解擴(kuò)前后的信號變化。你會直觀理解擴(kuò)頻增益到底是怎么把信號從噪聲里“撈”出來的——這個感覺比任何教科書描述都深刻。接下來如果你要繼續(xù)深入可以從捕獲算法優(yōu)化、碼環(huán)跟蹤實(shí)現(xiàn)、或者載波環(huán)輔助著手每一步都有新的坑等著你去踩但也有對應(yīng)的快感。本文還有配套的精品資源點(diǎn)擊獲取