OFDM通訊系統的MATLAB模擬(1)
由於是第一篇部落格，想先說點廢話，其實自己早就想把學到的一些東西總結成文章隨筆之類的供自己複習時檢視的了。但是一是覺得自己學的的不夠深入，總結也寫不出什麼很深刻的東西；二是覺得網上也有海量的資料了，需要時查一查根本不需要自己寫。但是恰恰也是網上的資料過於龐大，良莠不齊，導致每次都如海水一樣的知識湧入腦中，最後也如走馬觀花一般了了看下，知識吸收率低的驚人。現在也準備改變一下觀念，儘量把自己學過的東西歸納整理，以隨筆的形式發出來，可能有些地方我還不能理解作者的做法，我也會記錄出來，懂的地方解釋清楚，不懂的地方也標記清楚，同時在之後的不斷學習和總結中補上之前挖的坑，強迫症寫文章真的是太難了~哭。
不過也這樣安慰自己吧，將所學的進行輸出整理其實也是很重要的一步，看似時間浪費了好多，不過讀書人花的時間怎麼能叫浪費呢！整理出來自己看著不爽嗎。如果能同時幫助到其他人的話，那就太好了！希望自己能堅持下去。
關於OFDM系統我目前參考的是《MIMO-OFDM無線通訊技術及MATLAB實現》這本書，我看到網上用這本書做參考的人還挺多的，這裡將將作者實現OFDM系統的思路及其程式碼重新理順一遍。因為我理解的也不是很深入，有不對的地方希望大佬能幫忙指正。如果是沒怎麼接觸過OFDM的萌新，這篇文章可以幫助你對OFDM地模擬有個粗淺的瞭解。
首先畫一個我個人認為特別好理解的OFDM符號變化圖來幫助理解程式碼，接下來我會詳細的介紹每個步驟在幹什麼。
步驟0>
先假裝前景摘要一下。本OFDM系統模擬用到的技術主要有：16QAM調製解調 IFFT與FFT 新增AWGN噪聲，沒用到的有：通道編碼 交織 多徑瑞利通道 通道估計
步驟1>
這個模擬好像暫時沒有時間的概念，單位是按照取樣點來的。假設一幀有三個OFDM符號，每個符號長度為64(剛好在步驟3做IFFT時長度也為64，滿足2的冪次方)。我們首先生成數字基帶訊號，訊號長度為192個取樣點，由於要進行16QAM調製，我們直接隨機生成192個16進位制的數作為基帶訊號X(K)，然後再將X(K)經過16QAM星座圖對映便完成了調製。
另外在步驟1我們還要進行信噪比的一些初始化，便於計算噪聲幅度和最後的計算位元誤位元速率。
程式碼：
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; % 初始化誤位元速率
X = randi([0,15],1,NFFT * NFRAME); % 生成數字基帶訊號
X_mod = qammod(X,M,'gray') / (sqrt(10)); % 16QAM調製，格雷碼星座圖，並歸一化
星座圖補充：
簡單來說，16QAM星座圖對映就是將一個16進位制的數根據一定的規律變成一個複數。具體的數學推導我們這暫時不多加討論，只寫實現方法。比如下圖是一個按格雷碼編碼的16QAM星座圖(為什麼是格雷碼呢，你看它相鄰碼的二進位制都只有1位不同)。看星座圖可知，假設數字訊號X為5，經過qammod()函式後X_mod變為-a+bi;同理假設數字訊號X為13，經過qammod()函式後X_mod變為a+bi,其中a和b都是一個確定的數。注意最後要進行歸一化處理除以\(\sqrt[]{10}\)，如果是QPSK調製，則要除以\(\sqrt[]{2}\)。
信噪比補充：
S：訊號功率
N：噪聲功率
W：傳輸通道頻寬
Eb：訊號單位位元能量
Es：訊號單位符號能量
Rb：訊號位元傳輸速率
Rs：訊號符號傳輸速率
N0：噪聲的功率譜密度
K：通訊系統的調製率，\(\text{K} = log_2 M\)
EbN0 = Eb / N0 ：位元信噪比
EsN0 = Es / N0 ：符號信噪比
snr = S / N；訊號噪聲功率比
其實上面三個是不一樣的，但可以相互推導。初學的時候籠統的叫做信噪比，反而越學越混。在MATLAB模擬中，常常是設出位元信噪比EbN0，然後由以下公式計算出snr。
\[\begin{cases}
snr = \frac{S}{N} = \frac{Eb \cdot Rb}{N0 \cdot W}\\
Rb = Rs \cdot K
\end{cases}\]
得到：
\[snr = \frac{Eb \cdot Rs \cdot K}{N0 \cdot W}
\]
當滾降係數\(\alpha\)為0時，傳輸通道頻寬W等於通道中的符號速率：
\[W = (1+\alpha) \cdot Rs=Rs
\]
帶入snr表示式中，得到：
\[snr = \frac{Eb}{N0} \cdot K
\]
轉換為dB形式，程式碼裡的計算就是按照下面這個公式算的：
\[snr(dB) = EbN0(dB)+10 \cdot log_{10}(K)
\]
步驟2、3、4>
接下來的三個步驟分別如下，注意都是一個符號一個符號處理的，可回去看最開始的符號變化圖：
- 將每個OFDM符號的前一半和後一半交換，至於為什麼要做交換，我也不是很懂。有大佬知道的話希望能在評論區指導一下，感激不盡！
- 對交換過後的每個OFDM符號做IFFT，記錄輸出為x1(n)。
- 對每個OFDM符號新增迴圈字首CP，實際操作很簡單，因為這裡設的CP的長度NCP為16。就是把每個符號的後16個取樣點新增到當前符號的最前面來，每個符號因此就變成了64+16=80個取樣點。
由於這三個步驟都是在一個迴圈裡處理的，所以我也就把步驟2、3、4寫到一起了。
程式碼：
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
步驟5>
經過上一步的處理，現在的訊號已經能夠通過發射機發射出去，本系統只考慮AWGN通道，關於多徑瑞利衰落通道下一篇文章再總結。於是考慮模擬生成高斯白噪聲。由於snr在程式開頭就已經確定好了，所以我們要根據snr計算噪聲幅度。
程式碼：
P_s = xt * xt' / NSYM / NFRAME; % 計算訊號功率
delta1 = sqrt(10 .^ (-snr / 10) * P_s / 2); % 計算噪聲標準差
yr = xt + delta1 * (randn(size(xt)) + 1i * randn(size(xt)));
計算噪聲功率及標準差補充：
知道計算公式後，程式碼只有三行，但噪聲功率和標準差具體是怎麼計算出來的呢？下面是簡單的推導。
首先計算出離散訊號的訊號功率P_s：
\[\begin{align}
P & = \lim_{N\to \infty}(\frac{1}{2N + 1} \sum_{-N}^N |x(n)|^2) \\
& = \frac{1}{N} \sum_{0}^N |x(n)|^2
\end{align}\]
再由snr與P和N的關係有：
\[snr = \frac{S}{N} = \frac{P}{N}
\]
轉換為dB形式：
\[snr(dB) = 10 \cdot log_{10}(\frac{P}{N})
\]
反解得：
\[N = 10^{- \frac{snr(dB)}{10}} \cdot P
\]
由於產生的是均值為零的高斯白噪聲，所以噪聲方差\(\delta^2=N\)。又由於此時我們處理的是複數！！在用randn生成復噪聲的時候方差會變大一倍，所以我們使用\(\delta1^2 =\frac{\delta^2}{2}\)來生成復噪聲，有：
\[\delta1^2 = \frac{\delta^2}{2} = \frac{N}{2} = 10^{- \frac{snr(dB)}{10}} \cdot P/2
\]
最終推導噪聲的標準差如下，程式碼只是實現了一下下面的公式而已
\[\delta1 = \sqrt{10^{- \frac{snr(dB)}{10}} \cdot P/2}
\]
另外一種加噪方式補充：
也可以直接用官方函式，它會根據輸入的訊號向量xt和snr，自動計算出要加的噪聲功率並且加到訊號上輸出，訊號為實數或複數都可以(我特意去看了它的函式內部實現，發現如果輸入訊號是複數的話，噪聲功率也會特意除以2，也應證了我上面的說法)。
yr = awgn(xt,snr,'measured'); % 官方函式直接加噪
步驟6、7>
現在的訊號已經是添加了高斯白噪聲的訊號，我們的模擬已經完成了一半。接下來的兩個步驟與步驟2、3、4是呈映象，倒著實現一遍就行了。步驟分別如下，注意都是一個符號一個符號處理的，可回去看最開始的符號變化圖：
- 對每個OFDM符號去除迴圈字首CP，就是把每個符號的前16個取樣點去掉就好。
- 對每個OFDM符號做FFT，然後將將每個OFDM符號的前一半和後一半交換，記錄輸出為Y(K)。
程式碼：
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)];
len_a = len_a + NFFT;
len_b = len_b + NFFT + NCP;
len_left = len_left + NFFT; len_right = len_right + NFFT;
end
步驟8>
16QAM解調，這裡是直接用的官方自帶函式
程式碼：
Yr = qamdemod(Y * sqrt(10),M,'gray');
步驟9>
16QAM解調完畢後，其實我們已經可以自己在工作區裡對比解調得到的訊號Yr和我們的基帶數字訊號X了。但作為嚴謹的科研小生，怎麼能不進行誤位元速率分析呢？於是當前步驟我們研究一下怎麼分析誤位元速率。其實也很簡單，計算一下Yr和X有幾位元不相同，再計算一下總共有幾位元，把它們相除就得到了我們的位元誤位元速率(BER)。
需要注意的一點是，既然是誤位元率，就要把16進位制的訊號轉換成2進位制，以位元為單位計算錯誤數
程式碼：
Neb = sum(sum(de2bi(Yr,K) ~= de2bi(X,K))); % 轉為2進位制，計算具體有幾bit錯誤
Ntb = NFFT * NFRAME * K; % 模擬的總位元數
BER = Neb / Ntb;
完整程式碼：
最後貼一個完整程式碼，程式碼是參考的《MIMO-OFDM無線通訊技術及MATLAB實現》這本書。我是一行一行自己重新實現了一遍並且加上了詳細的中文註釋，希望能對像我這樣的剛入門的萌新有所啟發。對了，後面有個與理論值相比較的作圖函式有點佔位置，我就暫時不放到這篇文章中了XD。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Version 3.0
%%% 16QAM調製(官方函式)
%%% IFFT(官方函式)
%%% 新增迴圈字首(單徑AWGN中CP好像沒啥用？？)
%%% 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)
EbN0 = 0:1:20; % 設出位元信噪比(dB)
snr = EbN0 + 10 * log10(K); % 由他倆關係推出(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) % 對於每一種位元信噪比，計算該通訊環境下的誤位元速率
%% 16QAM調製(官方函式)
X_mod = qammod(X,M,'gray') / (sqrt(10)); % 16QAM調製，格雷碼星座圖，並歸一化
%% IFFT與迴圈字首新增
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
%% 由snr計算噪聲幅度並加噪
% P_s = xt * xt' / NSYM / NFRAME; % 計算訊號功率
% delta1 = sqrt(10 .^ (-snr(i) / 10) * P_s / 2); % 計算噪聲標準差
% yr = delta1 * (randn(size(xt)) + 1i * randn(size(xt)));
yr = awgn(xt,snr(i),'measured'); % 官方函式直接加噪
%% 去除迴圈字首並且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; % 每一符號分為左右兩邊
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)];
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(i) = Neb / Ntb;
fprintf('EbN0 = %d[dB], BER = %d / %d = %.3f\n', EbN0(i),Neb,Ntb,BER(i))
fprintf(fid, '%d %5.3f\n', EbN0(i),BER(i));
end
%% BER作圖分析
fclose(fid);
disp('Simulation is finished');
plot_ber(file_name,K);
幾個疑問：
有幾個沒搞懂的地方還是記錄一下
- 這個程式中新增的迴圈字首好像並沒有什麼作用，不知道是不是AWGN通道的原因，在多徑瑞利衰落通道中倒是與通道衝激響應有個卷積操作，這裡好像添加了CP也沒用到就去除了。
- 進行IFFT前為什麼要把每個OFDM符號進行左右交換，不是很懂欸。
參考文獻：
[1] 張少侃,呂聰敏,甘浩.數字通訊系統中Eb/N0與SNR轉換方法的研究[J].現代計算機,2019(12):33-36.
[2] Cho Y S, Kim J, Yang W Y, et al. MIMO-OFDM wireless communications with MATLAB[M]. John Wiley & Sons, 2010.