關於OFDM系統的MATLAB仿真實現的第二篇隨筆，在第一篇中，我們討論的是信號經過AWGN信道的情況，只用添加固定噪聲功率的高斯白噪聲就好了。但在實際無線信道中，信道干擾常常是加性噪聲、多徑衰落的結合。今天我們准備再進一步，讓信號經過多徑瑞利衰落信道。在這種信道條件下，信號具體是怎么怎么變化的呢？下面將講解系統仿真的各個部分以及實現多徑衰落的方法。
注意：為了整個系統的完整性，第一篇隨筆中的每個步驟這里也都又寫了一遍，但省略了補充知識部分，在第一篇的基礎上添加了實現多徑衰落的部分。想要看信噪比計算和噪聲功率計算的同學可以去看第一篇隨筆。
關於OFDM系統我目前參考的是《MIMO-OFDM無線通信技術及MATLAB實現》這本書，這里將將作者實現OFDM系統的思路及其代碼重新理順一遍。注意這篇文章我沒有一來就貼公式，巴拉巴拉講原理，那樣不就和老師上課念PPT一樣了嗎。其實我更喜歡直接學習大佬的仿真代碼，先對過程有個個大概思路再去推導細節和公式。這里因為我理解的水平也有限，有不對的地方希望大佬能幫忙指正。如果是沒怎么接觸過OFDM的萌新，這篇文章可以幫助你對OFDM符號級的仿真有個粗淺的了解XD。
首先畫一個我個人認為特別好理解的OFDM符號變化圖來幫助理解代碼，多徑瑞利衰落在步驟4到步驟5之間，在添加AWGN的前面。接下來我會詳細的介紹每個步驟在干什么。
在接觸了工程實現后，發現MATLAB的矩陣化處理對ofdm真是太有幫助了，直接把每個ofdm符號看作是矩陣的一列，對每一列分別進行相同的處理就可以了，代碼可以得到很大的簡化。之后的代碼是按照下圖進行的，是直接按一行數據處理的，其實是比較復雜的。
一：步驟0>
照例假裝前景摘要一下。本OFDM系統仿真用到的技術主要有：16QAM調制解調 IFFT與FFT 多徑瑞利信道 添加AWGN噪聲，沒用到的有：信道編碼 擴頻 交織 信道估計等等，哇，越難的技術越不想學(主要是學不懂)。這些技術的數學理論推導確實很難，但是在MATLAB仿真中往往用幾個自帶的函數就能解決問題，所以要實現一個簡單的OFDM系統還是很容易的，不要被天花亂墜般恐怖的數學公式嚇跑了(所以我最喜歡的就是直接看代碼的運行過程，然后有時間再去研究數學推導23333)。
二：步驟1>
這個仿真好像暫時沒有時間的概念，單位是按照采樣點來的。假設一幀有三個OFDM符號，每個符號長度為64(剛好在步驟3做IFFT時長度也為64，滿足2的冪次方)。我們首先生成數字基帶信號，信號長度為192個采樣點，由於要進行16QAM調制，我們直接隨機生成192個16進制的數作為基帶信號X(K)，然后再將X(K)經過16QAM星座圖映射便完成了調制。注意調制完輸出的X_mod是復信號。
另外在步驟1我們還要進行信噪比的一些初始化，便於計算噪聲幅度和最后的計算比特誤碼率。
2.1 增加部分
在步驟1中，我們增加對信道特征的初始化工作。主要是假設多徑信道個數和信道功率，以及各信道的時延，為之后信號通過多徑信道的計算做准備。
2.2 代碼
clc; clear all; clode all
NFRAME = 3; % 每一幀的OFDM符號數
NFFT = 64; % 每幀FFT長度
NCP = 16; % 循環前綴長度
NSYM = NFFT + NCP; % OFDM符號長度
M = 16; K = 4; % M:調制階數，K:log2(M)
EbN0 = 0; % 設出比特信噪比(dB)
snr = EbN0 + 10 * log10(K); % 由公式推出snr(dB)表達式
BER(1 : length(EbN0)) = 0; % 初始化誤碼率
P_hdB = [0 -8 -17 -21 -25]; % 各信道功率特性(dB)
D_h = [0 3 5 6 8]; % 各信道延遲(采樣點)
P_h = 10 .^ (P_hdB / 10); % 各信道功率特性
NH = length(P_hdB); % 多徑信道個數
LH = D_h(end)+1; % 信道長度(延遲過后)
X = randi([0,15],1,NFFT * NFRAME); % 生成數字基帶信號
X_mod = qammod(X,M,'gray') / (sqrt(10)); % 16QAM調制，格雷碼星座圖，並歸一化
三：步驟2、3、4>
接下來的三個步驟分別如下，注意都是一個符號一個符號處理的，可回去看最開始的符號變化圖：
- 將每個OFDM符號的前一半和后一半交換，至於為什么要做交換，我仍然不是很懂。有大佬知道的話希望能在評論區指導一下，感激不盡！補充：前后交換其實就是iffshift的過程，目的是將頻譜橫坐標中心頻點移動到0Hz處。
- 對交換過后的每個OFDM符號做IFFT，記錄輸出為x1(n)。
- 對每個OFDM符號添加循環前綴CP，實際操作很簡單，因為這里設的CP的長度NCP為16。就是把每個符號的后16個采樣點添加到當前符號的最前面來，每個符號因此就變成了64+16=80個采樣點。
由於這三個步驟都是在一個循環里處理的，所以我也就把步驟2、3、4寫到一起了。
3.1 代碼
x(1 : NFFT * NFRAME) = 0; % 預分配x數組
xt(1 : (NFFT + NCP) * NFRAME) = 0; % 預分配x_t數組
len_a = 1 : NFFT; % 處理的X位置
len_b = 1 : (NFFT + NCP); % 處理的X位置(加上CP的長度)
len_c = 1 : NCP;
len_left = 1 : (NFFT / 2); len_right = (NFFT / 2 + 1) : NFFT; % 每一符號分為左右兩邊
for frame = 1 : NFRAME % 對於每個OFDM符號都要翻轉和IFFT
x(len_a) = ifft([X_mod(len_right), X_mod(len_left)]); % 左右翻轉再ifft
xt(len_b) = [x(len_c + NFFT - NCP), x(len_a)]; % 添加CP后的信號數組
len_a = len_a + NFFT; % 更新找到下一個符號起點位置
len_b = len_b + NFFT + NCP; % 更新找到下一個符號起點位置(CP開頭)
len_c = len_c + NFFT;
len_left = len_left + NFFT; len_right = len_right + NFFT;
end
四：增加步驟：
如前面所說的，我們在步驟4和步驟5之間仿真信號xt經過多徑衰落信道。聽起來一頭霧水，說那么多有的沒的，其實就是做個卷積啦，就是拿信號xt與信道沖激響應h做卷積運算就OK了(終於有數字信號處理內味兒了~)。如何求信道沖激響應呢？這需要小小推導一下。
離散多徑衰落信道的一個簡單數學模型如下：
其中\(x(n)\)表示輸入信號，\(a_i(n)\)表示第i條路徑上的衰減系數，\(\tau_i(n)\)為第i條路經上的傳播時延。
由於表示的信道是線性信道，故可以用在\(n\)時刻對\(n-\tau\)時刻發射的沖激的響應\(h(\tau,n)\)來表示。我們已知用\(h(\tau,n)\)表示的經過信道的輸入\輸出為卷積關系：
於是由上述兩個公式我們可以推得多徑衰落信道沖激響應的數學表達式為：
4.1 瑞利隨機變量產生補充
無線信道傳播環境中的接收信號可以認為是來自無窮多個散射體的信號總和，由中心極限定理，可以用一個高斯隨機變量來表示接收信號。在一般的衰落環境中，無線衰落信道可以由復高斯隨機變量\(W1 + jW2\)表示，其中\(W1\)和\(W2\)都是均值為0，方差為\(\delta^2\)的獨立同分布(i.i.d.)高斯隨機變量。
所以先通過內置函數randn()生成好了\(W1\)和\(W2\)，再生成平均功率為\(E(X^2) = 2\delta^2\)的瑞利隨機變量。
在仿真中我們已經提前給出了瑞利信道平均功率\(P_h\)，所以有\(2\delta^2 = p_h\)，推出：
4.2 代碼
A_h = (randn(1,NH) + 1i * randn(1,NH)) .* sqrt(P_h / 2); % 由公式(4)生成瑞利隨機變量
h = zeros(1,LH); % 初始化信道沖激響應模型
h(D_h + 1) = A_h; % 信道沖激響應(同時體現出衰減系數和信道時延)，公式(3)的代碼體現
xt1 = conv(xt,h); % 卷積，輸出通過該信道的信號，公式(2)的代碼體現
五：步驟5>
經過上一步的處理，現在考慮仿真添加高斯白噪聲。由於snr在程序開頭就已經確定好了，所以我們要根據snr計算噪聲功率(噪聲方差)從而添加噪聲。注意由於卷積過后輸出信號長度會變長，計算信號功率時記得只取原本的長度。
5.1 代碼
xt2 = xt1(1 : NSYM * NFRAME); % 只取卷積過后屬於OFDM符號的部分
P_s = xt2 * xt2' ./ NSYM ./ NFRAME; % 計算信號功率
A_n = sqrt(10 .^ (-snr(i) / 10) * P_s / 2); % 計算噪聲標准差
yr = xt1 + A_n * (randn(size(xt1)) + 1i * randn(size(xt1))); % 根據噪聲標准差添加噪聲
六：步驟6、7>
現在的信號已經是經過多徑瑞利衰落並且添加了高斯白噪聲的信號，不容易啊！我們的仿真已經完成了一半。接下來的兩個步驟與步驟2、3、4是呈鏡像，倒着實現一遍就行了。步驟分別如下，注意都是一個符號一個符號處理的，可回去看最開始的符號變化圖：
- 對每個OFDM符號去除循環前綴CP，就是把每個符號的前16個采樣點去掉就好。
- 對每個OFDM符號做FFT，然后將將每個OFDM符號的前一半和后一半交換，記錄輸出為Y(K)。
注意要進行信道均衡，由經典公式\(Y = HX\)可知，\(X = Y / H\)，知道了信道頻率響應\(H\)就能由接收信號\(Y\)恢復發送信號\(X\)了，注意這里我們直接用理想的信道頻率響應，不通過發送導頻進行信道估計這種方式得到信道頻率響應
6.1 代碼
y(1 : NFFT * NFRAME) = 0; % 預分配y數組
Y(1 : NFFT * NFRAME) = 0; % 預分配Y數組
len_a = 1 : NFFT; % 處理的y位置
len_b = 1 : NFFT; % 處理的y位置
len_left = 1 : (NFFT / 2); len_right = (NFFT / 2 + 1) : NFFT; % 每一符號分為左右兩邊
for frame = 1 : NFRAME % 對於每個OFDM符號先去CP，再FFT再翻轉
y(len_a) = yr(len_b + NCP); % 去掉CP
Y(len_a) = fft(y(len_a)); % 先fft再翻轉
Y(len_a) = [Y(len_right), Y(len_left)] ./ H_shift; % 信道均衡（直接除以理想信道，沒有進行信道估計！）
len_a = len_a + NFFT;
len_b = len_b + NFFT + NCP;
len_left = len_left + NFFT; len_right = len_right + NFFT;
end
七：步驟8>
16QAM解調，這里是直接用的官方自帶函數
7.1 代碼
Yr = qamdemod(Y * sqrt(10),M,'gray');
八：步驟9>
16QAM解調完畢后，其實我們已經可以自己在工作區里對比解調得到的信號Yr和我們的基帶數字信號X了。但作為嚴謹的打工仔，怎么能不進行誤碼率分析呢？於是當前步驟我們研究一下怎么分析誤碼率。其實也很簡單，計算一下Yr和X有幾比特不相同，再計算一下總共有幾比特，把它們相除就得到了我們的比特誤碼率(BER)。
需要注意的一點是，既然是誤比特率，就要把16進制的信號轉換成2進制，以比特為單位計算錯誤數
8.1 代碼
Neb = sum(sum(de2bi(Yr,K) ~= de2bi(X,K))); % 轉為2進制，計算具體有幾bit錯誤
Ntb = NFFT * NFRAME * K; % 仿真的總比特數
BER = Neb / Ntb;
九：完整代碼
最后貼一個主函數代碼，包括三個小函數的完整的程序請點擊這里下載，代碼是參考的《MIMO-OFDM無線通信技術及MATLAB實現》這本書。我是一行一行自己重新實現了一遍並且加上了詳細的中文注釋，希望能對像我這樣的剛入門的萌新有所啟發。對了，后面有個與理論值相比較的作圖函數有點占位置，我就暫時不放到這篇文章中了XD。注意在包含多徑衰落信道的仿真的時候，如果想要仿真不同信噪比時的誤碼率，務必要生成一個狀態種子，保持衰落信道參數在每一次仿真中都不變。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Version 3.1
%%% 16QAM調制(官方函數)
%%% IFFT(官方函數)
%%% 添加循環前綴
%%% 經過多徑瑞利衰減信道
%%% 添加AWGN
%%% 去除循環前綴
%%% FFT(官方函數)
%%% 16QAM解調(官方函數)
%%% BER分析
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;close all;clc;
%% 基帶數字信號及一些初始化
NFRAME = 3; % 每一幀的OFDM符號數
NFFT = 64; % 每幀FFT長度
NCP = 16; % 循環前綴長度
NSYM = NFFT + NCP; % OFDM符號長度
M = 16; K = 4; % M:調制階數，K:log2(M)
P_hdB = [0 -8 -17 -21 -25]; % 各信道功率特性(dB)
D_h = [0 3 5 6 8]; % 各信道延遲(采樣點)
P_h = 10 .^ (P_hdB / 10); % 各信道功率特性
NH = length(P_hdB); % 多徑信道個數
LH = D_h(end)+1; % 信道長度(延遲過后)
EbN0 = 0:1:20; % 設出比特信噪比(dB)
snr = EbN0 + 10 * log10(K); % 由比特信噪比計算出snr(dB)
BER(1 : length(EbN0)) = 0; % 初始化誤碼率
file_name=['OFDM_BER_NCP' num2str(NCP) '.dat'];
fid=fopen(file_name, 'w+');
X = randi([0,15],1,NFFT * NFRAME); % 生成基帶數字信號
%%
for i = 1 : length(EbN0) % 對於每一種比特信噪比，計算該通信環境下的誤碼率
randn('state',0); % 很重要！！保持信道參數在每一次仿真中都不變
rand('state',0);
%% 16QAM調制(官方函數)
X_mod = qammod(X,M,'gray') / (sqrt(10)); % 16QAM調制，格雷碼星座圖，並歸一化
%% IFFT與循環前綴添加
x(1 : NFFT * NFRAME) = 0; % 預分配x數組
xt(1 : (NFFT + NCP) * NFRAME) = 0; % 預分配xt數組
len_a = 1 : NFFT; % 處理的X位置
len_b = 1 : (NFFT + NCP); % 處理的X位置(加上CP的長度)
len_c = 1 : NCP;
len_left = 1 : (NFFT / 2); len_right = (NFFT / 2 + 1) : NFFT; % 每一符號分為左右兩邊？？
for frame = 1 : NFRAME % 對於每個OFDM符號都要翻轉和IFFT
x(len_a) = ifft([X_mod(len_right), X_mod(len_left)]); % 左右翻轉再ifft
xt(len_b) = [x(len_c + NFFT - NCP), x(len_a)]; % 添加CP后的信號數組
len_a = len_a + NFFT; % 更新找到下一個符號起點位置
len_b = len_b + NFFT + NCP; % 更新找到下一個符號起點位置(CP開頭)
len_c = len_c + NFFT;
len_left = len_left + NFFT; len_right = len_right + NFFT;
end
%% 經過多徑瑞利衰減信道
A_h = (randn(1,NH) + 1i * randn(1,NH)) .* sqrt(P_h / 2); % 生成瑞利隨機變量
h = zeros(1,LH); % 初始化信道沖激響應模型
h(D_h + 1) = A_h; % 信道沖激響應(同時體現出衰減系數和信道時延)
xt1 = conv(xt,h); % 卷積，輸出通過該信道的信號
%% 由snr計算噪聲幅度並加噪
xt2 = xt1(1 : NSYM * NFRAME); % 只取卷積過后屬於OFDM符號的部分
P_s = xt2 * xt2' ./ NSYM ./ NFRAME; % 計算信號功率
A_n = sqrt(10 .^ (-snr(i) / 10) * P_s / 2); % 計算噪聲標准差
yr = xt1 + A_n * (randn(size(xt1)) + 1i * randn(size(xt1))); % 根據噪聲標准差添加噪聲
%% 去除循環前綴並且FFT
y(1 : NFFT * NFRAME) = 0; % 預分配y數組
Y(1 : NFFT * NFRAME) = 0; % 預分配Y數組
len_a = 1 : NFFT; % 處理的y位置
len_b = 1 : NFFT; % 處理的y位置
len_left = 1 : (NFFT / 2); len_right = (NFFT / 2 + 1) : NFFT; % 每一符號分為左右兩邊
H= fft([h zeros(1,NFFT-LH)]); % 信道頻率響應
H_shift(len_a)= [H(len_right) H(len_left)];
for frame = 1 : NFRAME % 對於每個OFDM符號先去CP，再FFT再翻轉
y(len_a) = yr(len_b + NCP); % 去掉CP
Y(len_a) = fft(y(len_a)); % 先fft再翻轉
Y(len_a) = [Y(len_right), Y(len_left)] ./ H_shift; % 信道均衡（直接除以理想信道，沒有進行信道估計！）
len_a = len_a + NFFT;
len_b = len_b + NFFT + NCP;
len_left = len_left + NFFT; len_right = len_right + NFFT;
end
%% 16QAM解調(官方函數)
Yr = qamdemod(Y * sqrt(10),M,'gray');
%% BER計算(多次迭代算均值會更准確)
Neb = sum(sum(de2bi(Yr,K) ~= de2bi(X,K))); % 轉為2進制，計算具體有幾bit錯誤
Ntb = NFFT * NFRAME * K; %[Ber,Neb,Ntb]=ber(bit_Rx,bit,Nbps);
BER(i) = Neb / Ntb;
fprintf('EbN0 = %3d[dB], BER = %4d / %8d = %11.3e\n', EbN0(i),Neb,Ntb,BER(i))
fprintf(fid, '%d %11.3e\n', EbN0(i),BER(i));
end
%% BER作圖分析
fclose(fid);
disp('Simulation is finished');
plot_ber(file_name,K);
十：參考文獻
[1] Tse D, Viswanath P. Fundamentals of wireless communication[M]. Cambridge university press, 2005.
[2] Cho Y S, Kim J, Yang W Y, et al. MIMO-OFDM wireless communications with MATLAB[M]. John Wiley & Sons, 2010.
[3] Goldsmith A. Wireless communications[M]. Cambridge university press, 2005.