自定義模塊開發(fā)：從原理到QPSK通信鏈路實現(xiàn))
1. 從“搭積木”到“造積木”為什么我們需要自定義S函數(shù)模塊如果你用過Simulink肯定熟悉那些琳瑯滿目的庫模塊增益、積分、傳遞函數(shù)、各種信號源和示波器。用這些現(xiàn)成的模塊搭建系統(tǒng)就像玩高級樂高效率很高。但當你需要實現(xiàn)一個標準庫里沒有的算法或者需要與外部硬件、特定數(shù)據(jù)格式、甚至是C/C代碼深度交互時你就會發(fā)現(xiàn)這些“標準積木”不夠用了。比如你想在Simulink里仿真一個自定義的通信調(diào)制解調(diào)算法或者實現(xiàn)一個復雜的非線性控制器標準模塊要么組合起來非常繁瑣要么根本無法實現(xiàn)。這時你就需要自己“造積木”。在Simulink中自定義模塊的終極武器就是S函數(shù)System-Function。它不是一個具體的模塊而是一個編程接口規(guī)范。你可以用MATLAB語言、C、C甚至Fortran來編寫S函數(shù)描述一個動態(tài)系統(tǒng)的輸入、輸出、狀態(tài)以及它們隨時間演化的關(guān)系。然后Simulink就能像調(diào)用內(nèi)置模塊一樣在仿真的每個時間步調(diào)用你的代碼。這相當于你為Simulink擴展了一套專屬的、功能強大的新“積木”。本次我們要“造”的就是一個實現(xiàn)QPSK正交相移鍵控調(diào)制與解調(diào)的通信鏈路“積木”。QPSK是數(shù)字通信中非常經(jīng)典的一種調(diào)制方式它通過改變載波的相位來傳遞信息一個符號可以攜帶2個比特的信息頻譜效率是BPSK的兩倍。在Simulink中雖然有Communications Toolbox提供專業(yè)的通信模塊但通過S函數(shù)親手實現(xiàn)一遍其價值遠超“仿真出一個波形”。你能徹底理解數(shù)據(jù)如何在仿真中流動、狀態(tài)如何更新、離散系統(tǒng)如何與連續(xù)仿真環(huán)境交互這是深入掌握Simulink仿真引擎和數(shù)字通信原理的絕佳路徑。2. 解剖S函數(shù)理解Simulink仿真引擎的“回調(diào)機制”在深入代碼之前我們必須先理解S函數(shù)是如何與Simulink仿真引擎協(xié)同工作的。Simulink的仿真可以看作一個按時間步推進的循環(huán)。在每個時間步引擎需要知道系統(tǒng)當前的狀態(tài)是什么輸入信號的值是多少基于當前狀態(tài)和輸入輸出應該是什么下一個狀態(tài)如果有應該更新為什么仿真初始化、結(jié)束時要做什么S函數(shù)通過一系列預定義好的“回調(diào)函數(shù)”Callback Methods來回答這些問題。你可以把S函數(shù)想象成一個你編寫的“演員”Simulink引擎是“導演”。在仿真的不同階段“導演”會喊出不同的指令調(diào)用對應的回調(diào)函數(shù)“演員”則執(zhí)行相應的動作。主要的回調(diào)函數(shù)包括mdlInitializeSizes: 這是“演員”的第一次亮相。它告訴“導演”我這個系統(tǒng)的基本信息有幾個輸入端口幾個輸出端口是連續(xù)系統(tǒng)、離散系統(tǒng)還是混合系統(tǒng)有沒有內(nèi)部狀態(tài)離散狀態(tài)或連續(xù)狀態(tài)采樣時間是多少這個函數(shù)在仿真一開始只調(diào)用一次用于定義系統(tǒng)的“骨架”。mdlInitializeSampleTimes: 緊接著它告訴“導演”我的“節(jié)奏”。是繼承上一級的采樣時間是固定的離散采樣時間比如0.001秒還是連續(xù)采樣對于我們的QPSK系統(tǒng)通常處理的是離散的數(shù)字符號所以會定義一個固定的離散采樣時間。mdlOutputs: 這是仿真循環(huán)中最重要的回調(diào)之一?！皩а荨眴枴敖o定當前時間、當前輸入和當前狀態(tài)你的輸出是什么”這個函數(shù)負責計算并輸出當前時間步的端口值。對于無狀態(tài)的靜態(tài)系統(tǒng)輸出只依賴于當前輸入對于有狀態(tài)的動態(tài)系統(tǒng)輸出依賴于當前狀態(tài)和/或輸入。mdlUpdate: 對于離散狀態(tài)系統(tǒng)這是另一個核心回調(diào)?！皩а荨眴枴霸谟嬎阃贻敵鲋蠡诋斍拜斎肽愕南乱粋€狀態(tài)在下一個采樣時刻的狀態(tài)應該更新為什么”這個函數(shù)負責更新離散狀態(tài)。注意mdlOutputs使用的是t時刻的狀態(tài)而mdlUpdate計算的是tTs時刻下一個采樣點的狀態(tài)。這個順序保證了仿真的因果性。mdlDerivatives: 如果是連續(xù)系統(tǒng)這個函數(shù)用來計算連續(xù)狀態(tài)的導數(shù)dx/dt供Simulink的數(shù)值積分器如ode45使用。我們的QPSK是離散系統(tǒng)通常用不到這個函數(shù)。mdlTerminate: 仿真結(jié)束時調(diào)用用于執(zhí)行一些清理工作如關(guān)閉文件、釋放內(nèi)存等。對于用MATLAB語言編寫的S函數(shù)即M語言S函數(shù)這些回調(diào)函數(shù)被寫在一個巨大的switch語句里由一個標志flag來區(qū)分。這是M語言S函數(shù)的經(jīng)典結(jié)構(gòu)模板。理解了這個“導演-演員”的交互模型再看代碼就不會覺得是一團亂麻了。注意S函數(shù)有Level-1和Level-2兩種版本。Level-1較為古老使用flag切換Level-2面向?qū)ο蟾鼜姶箪`活支持多實例、更豐富的端口和數(shù)據(jù)類型。對于復雜模塊推薦使用Level-2。但為了清晰理解基本原理我們從經(jīng)典的Level-1 M語言S函數(shù)入手。3. QPSK鏈路核心原理與S函數(shù)設(shè)計映射在動手寫代碼前我們需要把通信原理“翻譯”成S函數(shù)能理解的“系統(tǒng)描述”。一個最簡單的QPSK調(diào)制解調(diào)鏈路可以分解為以下幾個步驟我們將其映射到S函數(shù)的設(shè)計要素上發(fā)射端調(diào)制二進制數(shù)據(jù)源產(chǎn)生隨機的0/1比特流。這可以作為S函數(shù)的輸入也可以直接在S函數(shù)內(nèi)部生成。我們設(shè)計為外部輸入更通用。串并轉(zhuǎn)換將連續(xù)的比特流兩兩一組分成同相I路和正交Q路兩支路。例如比特對[0 0]映射到(I, Q) (-1, -1)。這屬于數(shù)據(jù)處理邏輯將在mdlOutputs中實現(xiàn)。脈沖成型可選但重要為了限制信號帶寬需要將離散的I/Q符號通過一個成型濾波器如升余弦滾降濾波器。這會引入記憶性即當前輸出不僅取決于當前符號還取決于之前的符號。這需要S函數(shù)內(nèi)部有**狀態(tài)State**來存儲濾波器狀態(tài)或之前的符號。狀態(tài)更新在mdlUpdate中完成。載波調(diào)制將基帶I/Q信號分別與載波cos(2πfct)和-sin(2πfct)相乘后相加得到帶通信號。在仿真中我們通常直接處理復基帶信號即I j*Q載波調(diào)制隱含在后續(xù)處理中。因此S函數(shù)的輸出可以是復基帶信號。信道簡化 5.加性高斯白噪聲AWGN在信號中加入噪聲。這可以在另一個S函數(shù)或使用Simulink內(nèi)置的AWGN Channel模塊實現(xiàn)。為了簡化我們可以將噪聲添加集成到調(diào)制S函數(shù)中或者作為獨立的鏈路環(huán)節(jié)。接收端解調(diào)6.匹配濾波/采樣接收端通過匹配濾波器通常與成型濾波器相同最大化信噪比并在最佳時刻采樣。這同樣需要狀態(tài)濾波器狀態(tài)。 7.判決根據(jù)采樣得到的I/Q值判斷它離哪個星座點(±1, ±1)最近并映射回比特對。這屬于無狀態(tài)的決策邏輯在mdlOutputs中實現(xiàn)。我們的S函數(shù)設(shè)計決策為了聚焦S函數(shù)本身我們將實現(xiàn)一個包含成型濾波的QPSK調(diào)制器和一個包含匹配濾波與判決的QPSK解調(diào)器作為兩個獨立的S函數(shù)模塊。它們通過一個簡單的加法器模擬噪聲連接。這樣我們可以清晰地看到調(diào)制器S函數(shù)輸入為二進制比特流輸出為復基帶信號I/Q兩路。內(nèi)部需要狀態(tài)來存儲成型濾波器的延遲線。解調(diào)器S函數(shù)輸入為含噪的復基帶信號輸出為解調(diào)后的二進制比特流。內(nèi)部同樣需要狀態(tài)來存儲匹配濾波器的延遲線。成型濾波器我們選擇簡單的矩形脈沖不整形這樣就沒有碼間串擾也無需狀態(tài)便于入門。在后續(xù)優(yōu)化中我們再引入升余弦濾波器及其狀態(tài)處理。4. 手把手實現(xiàn)QPSK調(diào)制器M語言S函數(shù)讓我們從最簡單的開始一個無脈沖成型即使用矩形脈沖的QPSK調(diào)制器。它沒有記憶性因此不需要內(nèi)部狀態(tài)。第一步創(chuàng)建S函數(shù)模板文件在MATLAB中新建一個.m文件命名為qpsk_modulator_simple.m。開頭是標準的函數(shù)定義function [sys,x0,str,ts,simStateCompliance] qpsk_modulator_simple(t,x,u,flag, Ts, EbN0_dB) % QPSK調(diào)制器S函數(shù) (簡化版矩形脈沖成型) % 輸入?yún)?shù): % t: 當前仿真時間 % x: 當前狀態(tài)向量本例無狀態(tài)為空 % u: 輸入向量本例為單路二進制比特流0/1 % flag: Simulink引擎調(diào)用的回調(diào)標志 % Ts: 符號周期秒通過S函數(shù)參數(shù)傳入 % EbN0_dB: 信噪比Eb/N0 (dB)用于計算噪聲功率本例暫未使用 % 輸出: % sys: 根據(jù)flag不同返回不同的系統(tǒng)信息 % x0: 初始狀態(tài)向量 % str: 保留參數(shù)始終設(shè)為空矩陣 [] % ts: 采樣時間矩陣 [采樣周期, 時間偏移] % simStateCompliance: 保留參數(shù) switch flag case 0 % flag 0: 初始化回調(diào) (mdlInitializeSizes) [sys,x0,str,ts,simStateCompliance]mdlInitializeSizes(Ts); case 3 % flag 3: 輸出計算回調(diào) (mdlOutputs) sysmdlOutputs(t,x,u); case {1, 2, 4, 9} % flag 1,2,4,9: 本例未使用返回空矩陣 sys []; otherwise error([Unhandled flag ,num2str(flag)]); end第二步實現(xiàn)初始化函數(shù)mdlInitializeSizes這個函數(shù)定義了模塊的“規(guī)格說明書”。function [sys,x0,str,ts,simStateCompliance]mdlInitializeSizes(Ts) % 初始化函數(shù) sizes simsizes; % 獲取一個標準的結(jié)構(gòu)體模板 sizes.NumContStates 0; % 連續(xù)狀態(tài)數(shù)量0 sizes.NumDiscStates 0; % 離散狀態(tài)數(shù)量0 (無記憶) sizes.NumOutputs 2; % 輸出數(shù)量2 (I路和Q路我們輸出兩路實數(shù)) sizes.NumInputs 1; % 輸入數(shù)量1 (單路比特流) sizes.DirFeedthrough 1; % 直接饋通是。因為輸出直接由當前輸入決定不依賴狀態(tài)。 sizes.NumSampleTimes 1; % 采樣時間數(shù)量1 sys simsizes(sizes); % 將設(shè)置應用到sys x0 []; % 初始狀態(tài)為空 str []; % 保留 ts [Ts 0]; % 采樣時間[采樣周期, 時間偏移]。偏移為0表示在整數(shù)倍Ts時刻采樣。 simStateCompliance UnknownSimState;這里的關(guān)鍵是DirFeedthrough設(shè)為1。因為我們的輸出調(diào)制符號直接由當前輸入的比特對決定不需要等待下一個采樣點更新狀態(tài)。采樣時間ts設(shè)置為[Ts, 0]意味著這個模塊以固定的Ts秒為周期運行。第三步實現(xiàn)輸出計算函數(shù)mdlOutputs這是核心算法所在。注意輸入u是單個比特0或1。但QPSK每次處理2個比特。我們需要一個機制來緩存比特。由于沒有使用狀態(tài)這里有一個關(guān)鍵技巧我們假設(shè)輸入端口的數(shù)據(jù)速率是符號速率的兩倍即每個Ts時間內(nèi)輸入2個比特。但實際上S函數(shù)在這個例子中每個Ts被調(diào)用一次輸入u只有一個值。因此這個簡單版本實際上有缺陷它每次只能處理1個比特無法完成QPSK映射。這引出了我們第一個重要的設(shè)計修正。修正設(shè)計我們需要狀態(tài)來緩存比特即使是無成型濾波QPSK調(diào)制本身也需要記憶2個比特才能映射到一個符號。因此我們必須引入離散狀態(tài)來存儲上一個比特。讓我們修改設(shè)計sizes.NumDiscStates 1;用于存儲上一個比特。并在mdlUpdate中更新它。同時DirFeedthrough應設(shè)為1因為當前輸出由當前輸入u和狀態(tài)x上一個比特共同決定。修正后的初始化sizes.NumDiscStates 1; % 離散狀態(tài)數(shù)量1 (用于緩存上一個比特) sizes.DirFeedthrough 1; % 直接饋通是。輸出依賴于當前輸入和狀態(tài)。修正后的輸出函數(shù)mdlOutputsfunction sysmdlOutputs(t,x,u) % 輸出計算函數(shù) % x: 狀態(tài)這里存儲上一個比特 (0/1) % u: 當前輸入比特 persistent bit_buffer; % 使用持久變量作為另一種緩存方式但更推薦用狀態(tài)x % 但為了教學清晰我們這里用狀態(tài)x作為緩存。 % 將狀態(tài)x(1)解釋為上一個比特 prev_bit x(1); current_bit u(1); % QPSK映射規(guī)則: (I, Q) % 比特對 [prev_bit, current_bit] - (I, Q) % 00 - (-1, -1), 01 - (-1, 1), 10 - (1, -1), 11 - (1, 1) % 映射關(guān)系可以寫成: I 2*prev_bit - 1, Q 2*current_bit - 1 % 因為prev_bit和current_bit是0或12*bit-1得到-1或1。 I_out 2 * prev_bit - 1; Q_out 2 * current_bit - 1; sys [I_out; Q_out]; % 輸出I和Q兩路信號這里我們假設(shè)在仿真開始時狀態(tài)x(1)被初始化為0需要在x0中設(shè)置。這個映射邏輯將連續(xù)的比特流轉(zhuǎn)換為了I/Q符號。第四步實現(xiàn)狀態(tài)更新函數(shù)mdlUpdate我們需要在flag2時處理狀態(tài)更新。在switch語句中添加case 2 % flag 2: 狀態(tài)更新回調(diào) (mdlUpdate) sysmdlUpdate(t,x,u);并實現(xiàn)該函數(shù)function sysmdlUpdate(t,x,u) % 狀態(tài)更新函數(shù)用當前輸入比特更新狀態(tài)即緩存當前比特供下一個符號周期使用 % 注意這里有一個重要的時序邏輯。 % 在當前時刻t我們使用狀態(tài)x(prev_bit)和輸入u(current_bit)計算輸出。 % 然后我們需要為下一個時刻tTs更新狀態(tài)。這個新狀態(tài)應該是當前的輸入比特u。 % 因為下一個符號周期當前的u將變成“上一個比特”。 sys u; % 將當前輸入比特存入狀態(tài)同時修改初始化中的x0x0 [0]; % 初始狀態(tài)假設(shè)第一個“上一個比特”為0。一個完整的、可工作的簡單QPSK調(diào)制器S函數(shù)就完成了。它接收單比特流每兩個比特周期注意由于我們每個Ts處理一個比特對但輸入是每個Ts一個比特所以實際符號周期是2*Ts這里需要仔細思考輸出一個QPSK符號。這里的采樣時間Ts實際是比特周期。符號周期是2*Ts。這是一個需要根據(jù)系統(tǒng)設(shè)計澄清的點。在更正式的設(shè)計中我們可能希望輸入端口接收的是并行的兩路比特I和Q支路或者接收一個整數(shù)輸入0,1,2,3代表2個比特。這體現(xiàn)了S函數(shù)接口設(shè)計的靈活性也反映了實際工程中的權(quán)衡。5. 引入脈沖成型讓S函數(shù)擁有“記憶”矩形脈沖的頻譜太寬在實際通信中不可行?，F(xiàn)在我們引入升余弦滾降濾波器。濾波器是一個有記憶的系統(tǒng)當前輸出是當前輸入和過去多個輸入的加權(quán)和。這意味著我們的S函數(shù)需要更多的狀態(tài)來存儲濾波器的延遲線即過去的一系列輸入符號。我們重新設(shè)計調(diào)制器假設(shè)輸入已經(jīng)是并行的I、Q兩路符號取值±1。這樣輸入端口數(shù)為2每個Ts輸入一個復數(shù)符號由I和Q兩個實數(shù)構(gòu)成。成型濾波器作用在I路和Q路上它們是獨立的、相同的實濾波器。狀態(tài)設(shè)計對于一個長度為L階數(shù)為L-1的FIR濾波器我們需要存儲過去的L-1個輸入符號。我們可以用一個長度為L-1的狀態(tài)向量來實現(xiàn)這個延遲線。每次新的符號到來我們將其插入延遲線頭部并彈出最舊的符號。修改初始化function [sys,x0,str,ts,simStateCompliance]mdlInitializeSizes(Ts, L) % L: 升余弦濾波器的長度抽頭數(shù) sizes simsizes; sizes.NumContStates 0; sizes.NumDiscStates 2*(L-1); % 狀態(tài)數(shù)量為I路和Q路分別存儲(L-1)個過去符號 sizes.NumOutputs 2; % 輸出濾波后的I和Q sizes.NumInputs 2; % 輸入當前的I_symbol和Q_symbol (±1) sizes.DirFeedthrough 1; % 是因為輸出與當前輸入有關(guān)FIR濾波包含當前輸入 sizes.NumSampleTimes 1; sys simsizes(sizes); x0 zeros(2*(L-1), 1); % 初始狀態(tài)設(shè)為0表示初始時延遲線內(nèi)全是0 str []; ts [Ts 0]; simStateCompliance UnknownSimState;修改輸出函數(shù)mdlOutputs這里需要實現(xiàn)FIR濾波運算。假設(shè)我們已經(jīng)有了濾波器的系數(shù)向量b長度為L可以通過S函數(shù)參數(shù)傳入。function sysmdlOutputs(t,x,u,b) % u(1): 當前I路輸入符號, u(2): 當前Q路輸入符號 % x: 狀態(tài)向量前(L-1)個是I路延遲線后(L-1)個是Q路延遲線 % b: FIR濾波器系數(shù)向量 (長度為L) L length(b); % 提取I路延遲線狀態(tài) I_delay_line [u(1); x(1:L-1)]; % 當前輸入 過去狀態(tài) % 提取Q路延遲線狀態(tài) Q_delay_line [u(2); x(L:end)]; % 注意狀態(tài)向量的排列 % FIR濾波計算 I_filtered b(:) * I_delay_line(:); % 向量點乘 Q_filtered b(:) * Q_delay_line(:); sys [I_filtered; Q_filtered];修改狀態(tài)更新函數(shù)mdlUpdate更新延遲線為下一個時刻做準備。function sysmdlUpdate(t,x,u,L) % 更新延遲線狀態(tài)將當前輸入移入最舊的值移出 % x: [I_delay(1:L-1); Q_delay(1:L-1)] % 更新后: [u(1); I_delay(1:end-1); u(2); Q_delay(1:end-1)] sys x; % 先復制原狀態(tài) % 更新I路延遲線整體向下移位最舊的被丟棄新的u(1)放在最前面 sys(1:L-1) [u(1); x(1:L-2)]; % 更新Q路延遲線 sys(L:end) [u(2); x(L:end-1)];現(xiàn)在這個S函數(shù)實現(xiàn)了一個完整的、帶脈沖成型的QPSK調(diào)制器。輸入是符號速率的I/Q符號輸出是經(jīng)過成型濾波的基帶信號。濾波器系數(shù)b和長度L可以作為S函數(shù)的參數(shù)在Simulink模塊對話框中設(shè)置。6. 構(gòu)建完整鏈路在Simulink中集成與調(diào)試有了調(diào)制器和解調(diào)器解調(diào)器是類似的包含匹配濾波和判決判決時需要考慮噪聲影響S函數(shù)后我們可以在Simulink中搭建完整鏈路。創(chuàng)建S-Function模塊在Simulink庫瀏覽器中找到“User-Defined Functions”庫拖入一個“S-Function”模塊。雙擊它在“S-function name”框中填入你的M文件名如qpsk_modulator_pulse_shaping在“S-function parameters”框中填入?yún)?shù)如Ts, L, b。這些參數(shù)會在仿真時傳遞給S函數(shù)。配置輸入源使用“Bernoulli Binary Generator”產(chǎn)生隨機比特流。由于我們的調(diào)制器輸入需要的是I/Q符號±1可能需要一個“Bit to Integer”或自定義的映射模塊將每2個比特映射為一個0-3的整數(shù)再映射為I/Q符號。更直接的方法是使用兩個“Bernoulli Binary Generator”分別產(chǎn)生I路和Q路的比特0/1然后通過2*u-1的運算轉(zhuǎn)換為±1。這可以用一個MATLAB Function模塊或Gain/Bias模塊組合實現(xiàn)。添加噪聲使用“AWGN Channel”模塊。需要設(shè)置Eb/N0、符號周期、每符號能量等參數(shù)。對于QPSK每符號能量Es 2*Eb因為一個符號承載2比特。AWGN模塊的輸入信號功率需要根據(jù)你的信號幅度計算。如果I/Q符號為±1則符號能量Es I^2 Q^2 2。在加噪時需仔細換算。連接解調(diào)器解調(diào)器S函數(shù)輸入是含噪的I/Q兩路信號內(nèi)部先進行匹配濾波系數(shù)與發(fā)射端成型濾波器相同然后采樣在最佳采樣點通常就是符號周期整數(shù)倍時刻最后進行判決。判決邏輯是判斷濾波后的I/Q值最接近哪個星座點(±1, ±1)。觀察與驗證使用“Scope”觀察發(fā)射端I/Q信號、接收端濾波前后信號、以及比特誤差。使用“Error Rate Calculation”模塊計算誤比特率BER并與理論QPSK的BER公式0.5*erfc(sqrt(Eb/N0))進行比較。在較高信噪比下仿真結(jié)果應與理論值接近。使用“Spectrum Analyzer”觀察成型前后信號的頻譜可以看到升余弦濾波器如何有效地限制帶寬。調(diào)試S函數(shù)的實用技巧使用disp或fprintf在S函數(shù)的各個回調(diào)函數(shù)中插入打印語句輸出時間t、輸入u、狀態(tài)x、輸出sys的值。這是最直接的調(diào)試方法。注意在MATLAB Function模塊中也可以這樣做。檢查采樣時間確保所有模塊的采樣時間設(shè)置正確且兼容。特別是S函數(shù)的ts設(shè)置。使用Simulink的“Display - Sample Time - Colors”功能用不同顏色顯示不同采樣時間的信號線。驗證狀態(tài)更新邏輯這是最容易出錯的地方。確保在mdlOutputs中使用的狀態(tài)是“當前時刻”的狀態(tài)而在mdlUpdate中計算的是“下一時刻”的狀態(tài)?？梢允謩佑嬎銕讉€時間步與仿真輸出對比。處理數(shù)據(jù)類型Simulink信號有數(shù)據(jù)類型double, single, int8等。確保你的S函數(shù)內(nèi)部計算能處理輸入的數(shù)據(jù)類型或者使用mdlInitializeSizes中的DWork向量來指定和轉(zhuǎn)換數(shù)據(jù)類型。對于M語言S函數(shù)通常使用double類型最簡單。7. 從M語言到C MEX S函數(shù)追求極致的仿真速度當你用M語言寫完S函數(shù)并驗證功能正確后可能會發(fā)現(xiàn)仿真速度較慢尤其是當濾波器階數(shù)很高、仿真時間很長時。這是因為M語言是解釋執(zhí)行的。對于計算密集的部分我們可以將其升級為C MEX S函數(shù)。C MEX S函數(shù)本質(zhì)是一個用C語言編寫、并被編譯成MEX文件在Windows上是.mexw64文件的S函數(shù)。它被Simulink直接以機器碼調(diào)用速度遠超M語言。編寫C MEX S函數(shù)更復雜需要熟悉C語言和Simulink C S函數(shù)的API。關(guān)鍵步驟獲取模板MATLAB提供了模板文件。你可以使用命令edit sfuntmpl_basic.c來查看一個基礎(chǔ)的C S函數(shù)模板。更簡單的方法是使用Legacy Code Tool傳統(tǒng)代碼工具將已有的C函數(shù)封裝成S函數(shù)或者使用S-Function BuilderGUI工具。實現(xiàn)回調(diào)函數(shù)與M語言類似你需要實現(xiàn)mdlInitializeSizes,mdlInitializeSampleTimes,mdlOutputs,mdlUpdate等函數(shù)但使用的是C語言語法和Simulink的宏如ssGetInputPortSignal,ssGetOutputPortSignal,ssGetDiscStates。編譯在MATLAB命令窗口使用mex命令編譯你的.c文件例如mex qpsk_modulator.c sfun.c。這需要你的系統(tǒng)配置了支持的C編譯器如Microsoft Visual C。在Simulink中使用和M語言S函數(shù)一樣拖入S-Function模塊名稱填寫編譯后的MEX文件名不含擴展名。一個簡單的C MEXmdlOutputs函數(shù)示例概念性static void mdlOutputs(SimStruct *S, int_T tid) { // 獲取輸入、輸出、狀態(tài)指針 real_T *u ssGetInputPortRealSignal(S, 0); // 輸入指針 real_T *y ssGetOutputPortRealSignal(S, 0); // 輸出指針 real_T *x ssGetDiscStates(S); // 離散狀態(tài)指針 // 從S函數(shù)參數(shù)獲取濾波器系數(shù)和長度 real_T *b mxGetPr(ssGetSFcnParam(S, 0)); // 假設(shè)第一個參數(shù)是系數(shù)向量 int_T L mxGetNumberOfElements(ssGetSFcnParam(S, 0)); // 實現(xiàn)FIR濾波計算 (簡化未考慮延遲線索引細節(jié)) real_T I_sum 0.0, Q_sum 0.0; for (int i0; iL; i) { I_sum b[i] * get_I_delay_element(x, i, L); // 需要自定義函數(shù)獲取延遲線元素 Q_sum b[i] * get_Q_delay_element(x, i, L); } y[0] I_sum; // I路輸出 y[1] Q_sum; // Q路輸出 }編寫C MEX S函數(shù)需要細致的內(nèi)存管理和指針操作但帶來的性能提升在大型仿真中是非?？捎^的。對于通信系統(tǒng)仿真將核心算法如濾波、編碼、解碼用C MEX實現(xiàn)是工業(yè)級仿真的常見做法。8. 避坑指南與性能優(yōu)化實戰(zhàn)心得在實現(xiàn)和集成S函數(shù)的過程中我踩過不少坑也總結(jié)了一些優(yōu)化經(jīng)驗1. 采樣時間同步陷阱問題整個鏈路中比特源、調(diào)制器、信道、解調(diào)器可能有不同的采樣時間。如果設(shè)置不當會導致信號維度不匹配或仿真錯誤。對策明確每個模塊的“時間基準”。通常以符號周期為基準最為清晰。比特源的采樣時間應為Ts/2如果每符號2比特調(diào)制器和解調(diào)器為Ts成型濾波器內(nèi)部如果采用多倍采樣上采樣則需要更快的采樣時間。使用“Rate Transition”模塊處理不同速率信號間的轉(zhuǎn)換或者確保S函數(shù)能處理多速率輸入通過設(shè)置ts為[0, 0]繼承輸入端口速率。2. 狀態(tài)初始化與邊界處理問題濾波器在開始階段延遲線是空的初始化為0。這會導致開始的幾個輸出不正確瞬態(tài)效應。在計算誤碼率時通常需要忽略開頭的一段數(shù)據(jù)。對策在S函數(shù)中可以添加一個計數(shù)器狀態(tài)在濾波器充滿之前不進行有效輸出或者直接輸出0。在誤碼率計算模塊中設(shè)置一個“延遲”參數(shù)跳過開始的若干符號。3. 復數(shù)信號處理問題Simulink默認信號是實數(shù)的。處理復基帶信號I jQ時需要拆分成兩路實數(shù)信號I路和Q路進行傳輸和處理這使模型連線翻倍容易混亂。對策可以考慮使用real-imag到complex的轉(zhuǎn)換模塊或者在S函數(shù)內(nèi)部直接處理復數(shù)。對于C MEX S函數(shù)可以創(chuàng)建COMPLEX類型的端口。但在大多數(shù)情況下保持I/Q分離更直觀也便于用Scope觀察。4. 參數(shù)化與模塊封裝問題S函數(shù)參數(shù)如Ts,L,b在模塊對話框中以逗號分隔的字符串形式輸入不易管理和修改尤其是濾波器系數(shù)向量b。對策使用MATLAB工作空間變量。在模塊參數(shù)框中填寫變量名如Ts_symbol,rc_filter_coeff。在仿真腳本或模型回調(diào)函數(shù)如PreLoadFcn中計算這些變量。更好的方法是使用模塊封裝Mask。右鍵點擊S-Function模塊選擇“Mask - Create Mask”。在封裝編輯器中可以創(chuàng)建友好的參數(shù)輸入控件如編輯框、下拉菜單并為每個參數(shù)指定一個變量名。在“Initialization”頁簽你可以寫MATLAB代碼來處理這些輸入例如將字符串轉(zhuǎn)換成數(shù)值向量。封裝后的模塊看起來和Simulink內(nèi)置模塊一樣專業(yè)。5. 性能分析使用Simulink Profiler通過simulink.profiler.start和simulink.profiler.stop來分析仿真時間都花在哪里。如果你發(fā)現(xiàn)你的S函數(shù)占了大部分時間那就是將其轉(zhuǎn)換為C MEX S函數(shù)的最佳信號。向量化運算在M語言S函數(shù)的mdlOutputs中避免使用for循環(huán)處理向量。盡量使用MATLAB的向量化操作。例如FIR濾波用dot(b, delay_line)或卷積conv注意處理邊界會比循環(huán)快很多。通過這個從原理到實現(xiàn)從簡單到復雜從M語言到C語言的QPSK鏈路S函數(shù)開發(fā)旅程你獲得的不僅僅是一個能跑的仿真模型。你獲得的是對Simulink仿真內(nèi)核如何調(diào)度一個自定義組件的深刻理解是對數(shù)字通信基帶處理流程的親手實踐更是將算法思想轉(zhuǎn)化為可執(zhí)行、可集成、可優(yōu)化仿真模塊的工程能力。下次當Simulink庫模塊無法滿足你的需求時你會自信地知道可以自己動手造出最合適的那塊“積木”。 Simulink S函數(shù)自定義模塊開發(fā)：從原理到QPSK通信鏈路實現(xiàn)
1. 從“搭積木”到“造積木”為什么我們需要自定義S函數(shù)模塊如果你用過Simulink肯定熟悉那些琳瑯滿目的庫模塊增益、積分、傳遞函數(shù)、各種信號源和示波器。用這些現(xiàn)成的模塊搭建系統(tǒng)就像玩高級樂高效率很高。但當你需要實現(xiàn)一個標準庫里沒有的算法或者需要與外部硬件、特定數(shù)據(jù)格式、甚至是C/C代碼深度交互時你就會發(fā)現(xiàn)這些“標準積木”不夠用了。比如你想在Simulink里仿真一個自定義的通信調(diào)制解調(diào)算法或者實現(xiàn)一個復雜的非線性控制器標準模塊要么組合起來非常繁瑣要么根本無法實現(xiàn)。這時你就需要自己“造積木”。在Simulink中自定義模塊的終極武器就是S函數(shù)System-Function。它不是一個具體的模塊而是一個編程接口規(guī)范。你可以用MATLAB語言、C、C甚至Fortran來編寫S函數(shù)描述一個動態(tài)系統(tǒng)的輸入、輸出、狀態(tài)以及它們隨時間演化的關(guān)系。然后Simulink就能像調(diào)用內(nèi)置模塊一樣在仿真的每個時間步調(diào)用你的代碼。這相當于你為Simulink擴展了一套專屬的、功能強大的新“積木”。本次我們要“造”的就是一個實現(xiàn)QPSK正交相移鍵控調(diào)制與解調(diào)的通信鏈路“積木”。QPSK是數(shù)字通信中非常經(jīng)典的一種調(diào)制方式它通過改變載波的相位來傳遞信息一個符號可以攜帶2個比特的信息頻譜效率是BPSK的兩倍。在Simulink中雖然有Communications Toolbox提供專業(yè)的通信模塊但通過S函數(shù)親手實現(xiàn)一遍其價值遠超“仿真出一個波形”。你能徹底理解數(shù)據(jù)如何在仿真中流動、狀態(tài)如何更新、離散系統(tǒng)如何與連續(xù)仿真環(huán)境交互這是深入掌握Simulink仿真引擎和數(shù)字通信原理的絕佳路徑。2. 解剖S函數(shù)理解Simulink仿真引擎的“回調(diào)機制”在深入代碼之前我們必須先理解S函數(shù)是如何與Simulink仿真引擎協(xié)同工作的。Simulink的仿真可以看作一個按時間步推進的循環(huán)。在每個時間步引擎需要知道系統(tǒng)當前的狀態(tài)是什么輸入信號的值是多少基于當前狀態(tài)和輸入輸出應該是什么下一個狀態(tài)如果有應該更新為什么仿真初始化、結(jié)束時要做什么S函數(shù)通過一系列預定義好的“回調(diào)函數(shù)”Callback Methods來回答這些問題。你可以把S函數(shù)想象成一個你編寫的“演員”Simulink引擎是“導演”。在仿真的不同階段“導演”會喊出不同的指令調(diào)用對應的回調(diào)函數(shù)“演員”則執(zhí)行相應的動作。主要的回調(diào)函數(shù)包括mdlInitializeSizes: 這是“演員”的第一次亮相。它告訴“導演”我這個系統(tǒng)的基本信息有幾個輸入端口幾個輸出端口是連續(xù)系統(tǒng)、離散系統(tǒng)還是混合系統(tǒng)有沒有內(nèi)部狀態(tài)離散狀態(tài)或連續(xù)狀態(tài)采樣時間是多少這個函數(shù)在仿真一開始只調(diào)用一次用于定義系統(tǒng)的“骨架”。mdlInitializeSampleTimes: 緊接著它告訴“導演”我的“節(jié)奏”。是繼承上一級的采樣時間是固定的離散采樣時間比如0.001秒還是連續(xù)采樣對于我們的QPSK系統(tǒng)通常處理的是離散的數(shù)字符號所以會定義一個固定的離散采樣時間。mdlOutputs: 這是仿真循環(huán)中最重要的回調(diào)之一?！皩а荨眴枴敖o定當前時間、當前輸入和當前狀態(tài)你的輸出是什么”這個函數(shù)負責計算并輸出當前時間步的端口值。對于無狀態(tài)的靜態(tài)系統(tǒng)輸出只依賴于當前輸入對于有狀態(tài)的動態(tài)系統(tǒng)輸出依賴于當前狀態(tài)和/或輸入。mdlUpdate: 對于離散狀態(tài)系統(tǒng)這是另一個核心回調(diào)?！皩а荨眴枴霸谟嬎阃贻敵鲋蠡诋斍拜斎肽愕南乱粋€狀態(tài)在下一個采樣時刻的狀態(tài)應該更新為什么”這個函數(shù)負責更新離散狀態(tài)。注意mdlOutputs使用的是t時刻的狀態(tài)而mdlUpdate計算的是tTs時刻下一個采樣點的狀態(tài)。這個順序保證了仿真的因果性。mdlDerivatives: 如果是連續(xù)系統(tǒng)這個函數(shù)用來計算連續(xù)狀態(tài)的導數(shù)dx/dt供Simulink的數(shù)值積分器如ode45使用。我們的QPSK是離散系統(tǒng)通常用不到這個函數(shù)。mdlTerminate: 仿真結(jié)束時調(diào)用用于執(zhí)行一些清理工作如關(guān)閉文件、釋放內(nèi)存等。對于用MATLAB語言編寫的S函數(shù)即M語言S函數(shù)這些回調(diào)函數(shù)被寫在一個巨大的switch語句里由一個標志flag來區(qū)分。這是M語言S函數(shù)的經(jīng)典結(jié)構(gòu)模板。理解了這個“導演-演員”的交互模型再看代碼就不會覺得是一團亂麻了。注意S函數(shù)有Level-1和Level-2兩種版本。Level-1較為古老使用flag切換Level-2面向?qū)ο蟾鼜姶箪`活支持多實例、更豐富的端口和數(shù)據(jù)類型。對于復雜模塊推薦使用Level-2。但為了清晰理解基本原理我們從經(jīng)典的Level-1 M語言S函數(shù)入手。3. QPSK鏈路核心原理與S函數(shù)設(shè)計映射在動手寫代碼前我們需要把通信原理“翻譯”成S函數(shù)能理解的“系統(tǒng)描述”。一個最簡單的QPSK調(diào)制解調(diào)鏈路可以分解為以下幾個步驟我們將其映射到S函數(shù)的設(shè)計要素上發(fā)射端調(diào)制二進制數(shù)據(jù)源產(chǎn)生隨機的0/1比特流。這可以作為S函數(shù)的輸入也可以直接在S函數(shù)內(nèi)部生成。我們設(shè)計為外部輸入更通用。串并轉(zhuǎn)換將連續(xù)的比特流兩兩一組分成同相I路和正交Q路兩支路。例如比特對[0 0]映射到(I, Q) (-1, -1)。這屬于數(shù)據(jù)處理邏輯將在mdlOutputs中實現(xiàn)。脈沖成型可選但重要為了限制信號帶寬需要將離散的I/Q符號通過一個成型濾波器如升余弦滾降濾波器。這會引入記憶性即當前輸出不僅取決于當前符號還取決于之前的符號。這需要S函數(shù)內(nèi)部有**狀態(tài)State**來存儲濾波器狀態(tài)或之前的符號。狀態(tài)更新在mdlUpdate中完成。載波調(diào)制將基帶I/Q信號分別與載波cos(2πfct)和-sin(2πfct)相乘后相加得到帶通信號。在仿真中我們通常直接處理復基帶信號即I j*Q載波調(diào)制隱含在后續(xù)處理中。因此S函數(shù)的輸出可以是復基帶信號。信道簡化 5.加性高斯白噪聲AWGN在信號中加入噪聲。這可以在另一個S函數(shù)或使用Simulink內(nèi)置的AWGN Channel模塊實現(xiàn)。為了簡化我們可以將噪聲添加集成到調(diào)制S函數(shù)中或者作為獨立的鏈路環(huán)節(jié)。接收端解調(diào)6.匹配濾波/采樣接收端通過匹配濾波器通常與成型濾波器相同最大化信噪比并在最佳時刻采樣。這同樣需要狀態(tài)濾波器狀態(tài)。 7.判決根據(jù)采樣得到的I/Q值判斷它離哪個星座點(±1, ±1)最近并映射回比特對。這屬于無狀態(tài)的決策邏輯在mdlOutputs中實現(xiàn)。我們的S函數(shù)設(shè)計決策為了聚焦S函數(shù)本身我們將實現(xiàn)一個包含成型濾波的QPSK調(diào)制器和一個包含匹配濾波與判決的QPSK解調(diào)器作為兩個獨立的S函數(shù)模塊。它們通過一個簡單的加法器模擬噪聲連接。這樣我們可以清晰地看到調(diào)制器S函數(shù)輸入為二進制比特流輸出為復基帶信號I/Q兩路。內(nèi)部需要狀態(tài)來存儲成型濾波器的延遲線。解調(diào)器S函數(shù)輸入為含噪的復基帶信號輸出為解調(diào)后的二進制比特流。內(nèi)部同樣需要狀態(tài)來存儲匹配濾波器的延遲線。成型濾波器我們選擇簡單的矩形脈沖不整形這樣就沒有碼間串擾也無需狀態(tài)便于入門。在后續(xù)優(yōu)化中我們再引入升余弦濾波器及其狀態(tài)處理。4. 手把手實現(xiàn)QPSK調(diào)制器M語言S函數(shù)讓我們從最簡單的開始一個無脈沖成型即使用矩形脈沖的QPSK調(diào)制器。它沒有記憶性因此不需要內(nèi)部狀態(tài)。第一步創(chuàng)建S函數(shù)模板文件在MATLAB中新建一個.m文件命名為qpsk_modulator_simple.m。開頭是標準的函數(shù)定義function [sys,x0,str,ts,simStateCompliance] qpsk_modulator_simple(t,x,u,flag, Ts, EbN0_dB) % QPSK調(diào)制器S函數(shù) (簡化版矩形脈沖成型) % 輸入?yún)?shù): % t: 當前仿真時間 % x: 當前狀態(tài)向量本例無狀態(tài)為空 % u: 輸入向量本例為單路二進制比特流0/1 % flag: Simulink引擎調(diào)用的回調(diào)標志 % Ts: 符號周期秒通過S函數(shù)參數(shù)傳入 % EbN0_dB: 信噪比Eb/N0 (dB)用于計算噪聲功率本例暫未使用 % 輸出: % sys: 根據(jù)flag不同返回不同的系統(tǒng)信息 % x0: 初始狀態(tài)向量 % str: 保留參數(shù)始終設(shè)為空矩陣 [] % ts: 采樣時間矩陣 [采樣周期, 時間偏移] % simStateCompliance: 保留參數(shù) switch flag case 0 % flag 0: 初始化回調(diào) (mdlInitializeSizes) [sys,x0,str,ts,simStateCompliance]mdlInitializeSizes(Ts); case 3 % flag 3: 輸出計算回調(diào) (mdlOutputs) sysmdlOutputs(t,x,u); case {1, 2, 4, 9} % flag 1,2,4,9: 本例未使用返回空矩陣 sys []; otherwise error([Unhandled flag ,num2str(flag)]); end第二步實現(xiàn)初始化函數(shù)mdlInitializeSizes這個函數(shù)定義了模塊的“規(guī)格說明書”。function [sys,x0,str,ts,simStateCompliance]mdlInitializeSizes(Ts) % 初始化函數(shù) sizes simsizes; % 獲取一個標準的結(jié)構(gòu)體模板 sizes.NumContStates 0; % 連續(xù)狀態(tài)數(shù)量0 sizes.NumDiscStates 0; % 離散狀態(tài)數(shù)量0 (無記憶) sizes.NumOutputs 2; % 輸出數(shù)量2 (I路和Q路我們輸出兩路實數(shù)) sizes.NumInputs 1; % 輸入數(shù)量1 (單路比特流) sizes.DirFeedthrough 1; % 直接饋通是。因為輸出直接由當前輸入決定不依賴狀態(tài)。 sizes.NumSampleTimes 1; % 采樣時間數(shù)量1 sys simsizes(sizes); % 將設(shè)置應用到sys x0 []; % 初始狀態(tài)為空 str []; % 保留 ts [Ts 0]; % 采樣時間[采樣周期, 時間偏移]。偏移為0表示在整數(shù)倍Ts時刻采樣。 simStateCompliance UnknownSimState;這里的關(guān)鍵是DirFeedthrough設(shè)為1。因為我們的輸出調(diào)制符號直接由當前輸入的比特對決定不需要等待下一個采樣點更新狀態(tài)。采樣時間ts設(shè)置為[Ts, 0]意味著這個模塊以固定的Ts秒為周期運行。第三步實現(xiàn)輸出計算函數(shù)mdlOutputs這是核心算法所在。注意輸入u是單個比特0或1。但QPSK每次處理2個比特。我們需要一個機制來緩存比特。由于沒有使用狀態(tài)這里有一個關(guān)鍵技巧我們假設(shè)輸入端口的數(shù)據(jù)速率是符號速率的兩倍即每個Ts時間內(nèi)輸入2個比特。但實際上S函數(shù)在這個例子中每個Ts被調(diào)用一次輸入u只有一個值。因此這個簡單版本實際上有缺陷它每次只能處理1個比特無法完成QPSK映射。這引出了我們第一個重要的設(shè)計修正。修正設(shè)計我們需要狀態(tài)來緩存比特即使是無成型濾波QPSK調(diào)制本身也需要記憶2個比特才能映射到一個符號。因此我們必須引入離散狀態(tài)來存儲上一個比特。讓我們修改設(shè)計sizes.NumDiscStates 1;用于存儲上一個比特。并在mdlUpdate中更新它。同時DirFeedthrough應設(shè)為1因為當前輸出由當前輸入u和狀態(tài)x上一個比特共同決定。修正后的初始化sizes.NumDiscStates 1; % 離散狀態(tài)數(shù)量1 (用于緩存上一個比特) sizes.DirFeedthrough 1; % 直接饋通是。輸出依賴于當前輸入和狀態(tài)。修正后的輸出函數(shù)mdlOutputsfunction sysmdlOutputs(t,x,u) % 輸出計算函數(shù) % x: 狀態(tài)這里存儲上一個比特 (0/1) % u: 當前輸入比特 persistent bit_buffer; % 使用持久變量作為另一種緩存方式但更推薦用狀態(tài)x % 但為了教學清晰我們這里用狀態(tài)x作為緩存。 % 將狀態(tài)x(1)解釋為上一個比特 prev_bit x(1); current_bit u(1); % QPSK映射規(guī)則: (I, Q) % 比特對 [prev_bit, current_bit] - (I, Q) % 00 - (-1, -1), 01 - (-1, 1), 10 - (1, -1), 11 - (1, 1) % 映射關(guān)系可以寫成: I 2*prev_bit - 1, Q 2*current_bit - 1 % 因為prev_bit和current_bit是0或12*bit-1得到-1或1。 I_out 2 * prev_bit - 1; Q_out 2 * current_bit - 1; sys [I_out; Q_out]; % 輸出I和Q兩路信號這里我們假設(shè)在仿真開始時狀態(tài)x(1)被初始化為0需要在x0中設(shè)置。這個映射邏輯將連續(xù)的比特流轉(zhuǎn)換為了I/Q符號。第四步實現(xiàn)狀態(tài)更新函數(shù)mdlUpdate我們需要在flag2時處理狀態(tài)更新。在switch語句中添加case 2 % flag 2: 狀態(tài)更新回調(diào) (mdlUpdate) sysmdlUpdate(t,x,u);并實現(xiàn)該函數(shù)function sysmdlUpdate(t,x,u) % 狀態(tài)更新函數(shù)用當前輸入比特更新狀態(tài)即緩存當前比特供下一個符號周期使用 % 注意這里有一個重要的時序邏輯。 % 在當前時刻t我們使用狀態(tài)x(prev_bit)和輸入u(current_bit)計算輸出。 % 然后我們需要為下一個時刻tTs更新狀態(tài)。這個新狀態(tài)應該是當前的輸入比特u。 % 因為下一個符號周期當前的u將變成“上一個比特”。 sys u; % 將當前輸入比特存入狀態(tài)同時修改初始化中的x0x0 [0]; % 初始狀態(tài)假設(shè)第一個“上一個比特”為0。一個完整的、可工作的簡單QPSK調(diào)制器S函數(shù)就完成了。它接收單比特流每兩個比特周期注意由于我們每個Ts處理一個比特對但輸入是每個Ts一個比特所以實際符號周期是2*Ts這里需要仔細思考輸出一個QPSK符號。這里的采樣時間Ts實際是比特周期。符號周期是2*Ts。這是一個需要根據(jù)系統(tǒng)設(shè)計澄清的點。在更正式的設(shè)計中我們可能希望輸入端口接收的是并行的兩路比特I和Q支路或者接收一個整數(shù)輸入0,1,2,3代表2個比特。這體現(xiàn)了S函數(shù)接口設(shè)計的靈活性也反映了實際工程中的權(quán)衡。5. 引入脈沖成型讓S函數(shù)擁有“記憶”矩形脈沖的頻譜太寬在實際通信中不可行?，F(xiàn)在我們引入升余弦滾降濾波器。濾波器是一個有記憶的系統(tǒng)當前輸出是當前輸入和過去多個輸入的加權(quán)和。這意味著我們的S函數(shù)需要更多的狀態(tài)來存儲濾波器的延遲線即過去的一系列輸入符號。我們重新設(shè)計調(diào)制器假設(shè)輸入已經(jīng)是并行的I、Q兩路符號取值±1。這樣輸入端口數(shù)為2每個Ts輸入一個復數(shù)符號由I和Q兩個實數(shù)構(gòu)成。成型濾波器作用在I路和Q路上它們是獨立的、相同的實濾波器。狀態(tài)設(shè)計對于一個長度為L階數(shù)為L-1的FIR濾波器我們需要存儲過去的L-1個輸入符號。我們可以用一個長度為L-1的狀態(tài)向量來實現(xiàn)這個延遲線。每次新的符號到來我們將其插入延遲線頭部并彈出最舊的符號。修改初始化function [sys,x0,str,ts,simStateCompliance]mdlInitializeSizes(Ts, L) % L: 升余弦濾波器的長度抽頭數(shù) sizes simsizes; sizes.NumContStates 0; sizes.NumDiscStates 2*(L-1); % 狀態(tài)數(shù)量為I路和Q路分別存儲(L-1)個過去符號 sizes.NumOutputs 2; % 輸出濾波后的I和Q sizes.NumInputs 2; % 輸入當前的I_symbol和Q_symbol (±1) sizes.DirFeedthrough 1; % 是因為輸出與當前輸入有關(guān)FIR濾波包含當前輸入 sizes.NumSampleTimes 1; sys simsizes(sizes); x0 zeros(2*(L-1), 1); % 初始狀態(tài)設(shè)為0表示初始時延遲線內(nèi)全是0 str []; ts [Ts 0]; simStateCompliance UnknownSimState;修改輸出函數(shù)mdlOutputs這里需要實現(xiàn)FIR濾波運算。假設(shè)我們已經(jīng)有了濾波器的系數(shù)向量b長度為L可以通過S函數(shù)參數(shù)傳入。function sysmdlOutputs(t,x,u,b) % u(1): 當前I路輸入符號, u(2): 當前Q路輸入符號 % x: 狀態(tài)向量前(L-1)個是I路延遲線后(L-1)個是Q路延遲線 % b: FIR濾波器系數(shù)向量 (長度為L) L length(b); % 提取I路延遲線狀態(tài) I_delay_line [u(1); x(1:L-1)]; % 當前輸入 過去狀態(tài) % 提取Q路延遲線狀態(tài) Q_delay_line [u(2); x(L:end)]; % 注意狀態(tài)向量的排列 % FIR濾波計算 I_filtered b(:) * I_delay_line(:); % 向量點乘 Q_filtered b(:) * Q_delay_line(:); sys [I_filtered; Q_filtered];修改狀態(tài)更新函數(shù)mdlUpdate更新延遲線為下一個時刻做準備。function sysmdlUpdate(t,x,u,L) % 更新延遲線狀態(tài)將當前輸入移入最舊的值移出 % x: [I_delay(1:L-1); Q_delay(1:L-1)] % 更新后: [u(1); I_delay(1:end-1); u(2); Q_delay(1:end-1)] sys x; % 先復制原狀態(tài) % 更新I路延遲線整體向下移位最舊的被丟棄新的u(1)放在最前面 sys(1:L-1) [u(1); x(1:L-2)]; % 更新Q路延遲線 sys(L:end) [u(2); x(L:end-1)];現(xiàn)在這個S函數(shù)實現(xiàn)了一個完整的、帶脈沖成型的QPSK調(diào)制器。輸入是符號速率的I/Q符號輸出是經(jīng)過成型濾波的基帶信號。濾波器系數(shù)b和長度L可以作為S函數(shù)的參數(shù)在Simulink模塊對話框中設(shè)置。6. 構(gòu)建完整鏈路在Simulink中集成與調(diào)試有了調(diào)制器和解調(diào)器解調(diào)器是類似的包含匹配濾波和判決判決時需要考慮噪聲影響S函數(shù)后我們可以在Simulink中搭建完整鏈路。創(chuàng)建S-Function模塊在Simulink庫瀏覽器中找到“User-Defined Functions”庫拖入一個“S-Function”模塊。雙擊它在“S-function name”框中填入你的M文件名如qpsk_modulator_pulse_shaping在“S-function parameters”框中填入?yún)?shù)如Ts, L, b。這些參數(shù)會在仿真時傳遞給S函數(shù)。配置輸入源使用“Bernoulli Binary Generator”產(chǎn)生隨機比特流。由于我們的調(diào)制器輸入需要的是I/Q符號±1可能需要一個“Bit to Integer”或自定義的映射模塊將每2個比特映射為一個0-3的整數(shù)再映射為I/Q符號。更直接的方法是使用兩個“Bernoulli Binary Generator”分別產(chǎn)生I路和Q路的比特0/1然后通過2*u-1的運算轉(zhuǎn)換為±1。這可以用一個MATLAB Function模塊或Gain/Bias模塊組合實現(xiàn)。添加噪聲使用“AWGN Channel”模塊。需要設(shè)置Eb/N0、符號周期、每符號能量等參數(shù)。對于QPSK每符號能量Es 2*Eb因為一個符號承載2比特。AWGN模塊的輸入信號功率需要根據(jù)你的信號幅度計算。如果I/Q符號為±1則符號能量Es I^2 Q^2 2。在加噪時需仔細換算。連接解調(diào)器解調(diào)器S函數(shù)輸入是含噪的I/Q兩路信號內(nèi)部先進行匹配濾波系數(shù)與發(fā)射端成型濾波器相同然后采樣在最佳采樣點通常就是符號周期整數(shù)倍時刻最后進行判決。判決邏輯是判斷濾波后的I/Q值最接近哪個星座點(±1, ±1)。觀察與驗證使用“Scope”觀察發(fā)射端I/Q信號、接收端濾波前后信號、以及比特誤差。使用“Error Rate Calculation”模塊計算誤比特率BER并與理論QPSK的BER公式0.5*erfc(sqrt(Eb/N0))進行比較。在較高信噪比下仿真結(jié)果應與理論值接近。使用“Spectrum Analyzer”觀察成型前后信號的頻譜可以看到升余弦濾波器如何有效地限制帶寬。調(diào)試S函數(shù)的實用技巧使用disp或fprintf在S函數(shù)的各個回調(diào)函數(shù)中插入打印語句輸出時間t、輸入u、狀態(tài)x、輸出sys的值。這是最直接的調(diào)試方法。注意在MATLAB Function模塊中也可以這樣做。檢查采樣時間確保所有模塊的采樣時間設(shè)置正確且兼容。特別是S函數(shù)的ts設(shè)置。使用Simulink的“Display - Sample Time - Colors”功能用不同顏色顯示不同采樣時間的信號線。驗證狀態(tài)更新邏輯這是最容易出錯的地方。確保在mdlOutputs中使用的狀態(tài)是“當前時刻”的狀態(tài)而在mdlUpdate中計算的是“下一時刻”的狀態(tài)?？梢允謩佑嬎銕讉€時間步與仿真輸出對比。處理數(shù)據(jù)類型Simulink信號有數(shù)據(jù)類型double, single, int8等。確保你的S函數(shù)內(nèi)部計算能處理輸入的數(shù)據(jù)類型或者使用mdlInitializeSizes中的DWork向量來指定和轉(zhuǎn)換數(shù)據(jù)類型。對于M語言S函數(shù)通常使用double類型最簡單。7. 從M語言到C MEX S函數(shù)追求極致的仿真速度當你用M語言寫完S函數(shù)并驗證功能正確后可能會發(fā)現(xiàn)仿真速度較慢尤其是當濾波器階數(shù)很高、仿真時間很長時。這是因為M語言是解釋執(zhí)行的。對于計算密集的部分我們可以將其升級為C MEX S函數(shù)。C MEX S函數(shù)本質(zhì)是一個用C語言編寫、并被編譯成MEX文件在Windows上是.mexw64文件的S函數(shù)。它被Simulink直接以機器碼調(diào)用速度遠超M語言。編寫C MEX S函數(shù)更復雜需要熟悉C語言和Simulink C S函數(shù)的API。關(guān)鍵步驟獲取模板MATLAB提供了模板文件。你可以使用命令edit sfuntmpl_basic.c來查看一個基礎(chǔ)的C S函數(shù)模板。更簡單的方法是使用Legacy Code Tool傳統(tǒng)代碼工具將已有的C函數(shù)封裝成S函數(shù)或者使用S-Function BuilderGUI工具。實現(xiàn)回調(diào)函數(shù)與M語言類似你需要實現(xiàn)mdlInitializeSizes,mdlInitializeSampleTimes,mdlOutputs,mdlUpdate等函數(shù)但使用的是C語言語法和Simulink的宏如ssGetInputPortSignal,ssGetOutputPortSignal,ssGetDiscStates。編譯在MATLAB命令窗口使用mex命令編譯你的.c文件例如mex qpsk_modulator.c sfun.c。這需要你的系統(tǒng)配置了支持的C編譯器如Microsoft Visual C。在Simulink中使用和M語言S函數(shù)一樣拖入S-Function模塊名稱填寫編譯后的MEX文件名不含擴展名。一個簡單的C MEXmdlOutputs函數(shù)示例概念性static void mdlOutputs(SimStruct *S, int_T tid) { // 獲取輸入、輸出、狀態(tài)指針 real_T *u ssGetInputPortRealSignal(S, 0); // 輸入指針 real_T *y ssGetOutputPortRealSignal(S, 0); // 輸出指針 real_T *x ssGetDiscStates(S); // 離散狀態(tài)指針 // 從S函數(shù)參數(shù)獲取濾波器系數(shù)和長度 real_T *b mxGetPr(ssGetSFcnParam(S, 0)); // 假設(shè)第一個參數(shù)是系數(shù)向量 int_T L mxGetNumberOfElements(ssGetSFcnParam(S, 0)); // 實現(xiàn)FIR濾波計算 (簡化未考慮延遲線索引細節(jié)) real_T I_sum 0.0, Q_sum 0.0; for (int i0; iL; i) { I_sum b[i] * get_I_delay_element(x, i, L); // 需要自定義函數(shù)獲取延遲線元素 Q_sum b[i] * get_Q_delay_element(x, i, L); } y[0] I_sum; // I路輸出 y[1] Q_sum; // Q路輸出 }編寫C MEX S函數(shù)需要細致的內(nèi)存管理和指針操作但帶來的性能提升在大型仿真中是非?？捎^的。對于通信系統(tǒng)仿真將核心算法如濾波、編碼、解碼用C MEX實現(xiàn)是工業(yè)級仿真的常見做法。8. 避坑指南與性能優(yōu)化實戰(zhàn)心得在實現(xiàn)和集成S函數(shù)的過程中我踩過不少坑也總結(jié)了一些優(yōu)化經(jīng)驗1. 采樣時間同步陷阱問題整個鏈路中比特源、調(diào)制器、信道、解調(diào)器可能有不同的采樣時間。如果設(shè)置不當會導致信號維度不匹配或仿真錯誤。對策明確每個模塊的“時間基準”。通常以符號周期為基準最為清晰。比特源的采樣時間應為Ts/2如果每符號2比特調(diào)制器和解調(diào)器為Ts成型濾波器內(nèi)部如果采用多倍采樣上采樣則需要更快的采樣時間。使用“Rate Transition”模塊處理不同速率信號間的轉(zhuǎn)換或者確保S函數(shù)能處理多速率輸入通過設(shè)置ts為[0, 0]繼承輸入端口速率。2. 狀態(tài)初始化與邊界處理問題濾波器在開始階段延遲線是空的初始化為0。這會導致開始的幾個輸出不正確瞬態(tài)效應。在計算誤碼率時通常需要忽略開頭的一段數(shù)據(jù)。對策在S函數(shù)中可以添加一個計數(shù)器狀態(tài)在濾波器充滿之前不進行有效輸出或者直接輸出0。在誤碼率計算模塊中設(shè)置一個“延遲”參數(shù)跳過開始的若干符號。3. 復數(shù)信號處理問題Simulink默認信號是實數(shù)的。處理復基帶信號I jQ時需要拆分成兩路實數(shù)信號I路和Q路進行傳輸和處理這使模型連線翻倍容易混亂。對策可以考慮使用real-imag到complex的轉(zhuǎn)換模塊或者在S函數(shù)內(nèi)部直接處理復數(shù)。對于C MEX S函數(shù)可以創(chuàng)建COMPLEX類型的端口。但在大多數(shù)情況下保持I/Q分離更直觀也便于用Scope觀察。4. 參數(shù)化與模塊封裝問題S函數(shù)參數(shù)如Ts,L,b在模塊對話框中以逗號分隔的字符串形式輸入不易管理和修改尤其是濾波器系數(shù)向量b。對策使用MATLAB工作空間變量。在模塊參數(shù)框中填寫變量名如Ts_symbol,rc_filter_coeff。在仿真腳本或模型回調(diào)函數(shù)如PreLoadFcn中計算這些變量。更好的方法是使用模塊封裝Mask。右鍵點擊S-Function模塊選擇“Mask - Create Mask”。在封裝編輯器中可以創(chuàng)建友好的參數(shù)輸入控件如編輯框、下拉菜單并為每個參數(shù)指定一個變量名。在“Initialization”頁簽你可以寫MATLAB代碼來處理這些輸入例如將字符串轉(zhuǎn)換成數(shù)值向量。封裝后的模塊看起來和Simulink內(nèi)置模塊一樣專業(yè)。5. 性能分析使用Simulink Profiler通過simulink.profiler.start和simulink.profiler.stop來分析仿真時間都花在哪里。如果你發(fā)現(xiàn)你的S函數(shù)占了大部分時間那就是將其轉(zhuǎn)換為C MEX S函數(shù)的最佳信號。向量化運算在M語言S函數(shù)的mdlOutputs中避免使用for循環(huán)處理向量。盡量使用MATLAB的向量化操作。例如FIR濾波用dot(b, delay_line)或卷積conv注意處理邊界會比循環(huán)快很多。通過這個從原理到實現(xiàn)從簡單到復雜從M語言到C語言的QPSK鏈路S函數(shù)開發(fā)旅程你獲得的不僅僅是一個能跑的仿真模型。你獲得的是對Simulink仿真內(nèi)核如何調(diào)度一個自定義組件的深刻理解是對數(shù)字通信基帶處理流程的親手實踐更是將算法思想轉(zhuǎn)化為可執(zhí)行、可集成、可優(yōu)化仿真模塊的工程能力。下次當Simulink庫模塊無法滿足你的需求時你會自信地知道可以自己動手造出最合適的那塊“積木”。