OFDM通訊系統的MATLAB模擬(2)

羽扇綸巾o0發表於2020-07-23
Matlab

關於OFDM系統的MATLAB模擬實現的第二篇隨筆，在第一篇中，我們討論的是訊號經過AWGN通道的情況，只用新增固定噪聲功率的高斯白噪聲就好了。但在實際無線通道中，通道干擾常常是加性噪聲、多徑衰落的結合。今天我們準備再進一步，讓訊號經過多徑瑞利衰落通道。在這種通道條件下，訊號具體是怎麼怎麼變化的呢？下面將講解系統模擬的各個部分以及實現多徑衰落的方法。

注意：為了整個系統的完整性，第一篇隨筆中的每個步驟這裡也都又寫了一遍，但省略了補充知識部分，在第一篇的基礎上新增了實現多徑衰落的部分。想要看訊雜比計算和噪聲功率計算的同學可以去看第一篇隨筆。

關於OFDM系統我目前參考的是《MIMO-OFDM無線通訊技術及MATLAB實現》這本書，這裡將將作者實現OFDM系統的思路及其程式碼重新理順一遍。注意這篇文章我沒有一來就貼公式，巴拉巴拉講原理，那樣不就和老師上課念PPT一樣了嗎。其實我更喜歡直接學習大佬的模擬程式碼，先對過程有個個大概思路再去推導細節和公式。這裡因為我理解的水平也有限，有不對的地方希望大佬能幫忙指正。如果是沒怎麼接觸過OFDM的萌新，這篇文章可以幫助你對OFDM符號級的模擬有個粗淺的瞭解XD。

首先畫一個我個人認為特別好理解的OFDM符號變化圖來幫助理解程式碼，多徑瑞利衰落在步驟4到步驟5之間，在新增AWGN的前面。接下來我會詳細的介紹每個步驟在幹什麼。


步驟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我們還要進行訊雜比的一些初始化，便於計算噪聲幅度和最後的計算位元誤位元速率。

增加部分：

在步驟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;          % 初始化誤位元速率

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符號的前一半和後一半交換，至於為什麼要做交換，我仍然不是很懂。有大佬知道的話希望能在評論區指導一下，感激不盡！
  • 對交換過後的每個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

增加步驟：

如前面所說的，我們在步驟4和步驟5之間模擬訊號xt經過多徑衰落通道。聽起來一頭霧水，說那麼多有的沒的，其實就是做個卷積啦，就是拿訊號xt與通道衝激響應h做卷積運算就OK了(終於有數字訊號處理內味兒了~)。如何求通道衝激響應呢？這需要小小推導一下。

離散多徑衰落通道的一個簡單數學模型如下：

\[\begin{align} y(n) & = a_1(n)\cdot x(n-\tau_1(n)) + a_2(n)\cdot x(n-\tau_2(n)) + ... + a_N(n)\cdot x(n-\tau_N(n))\notag\\ & = \sum_{i = 1}^{N} a_i(n)x(n-\tau_i(n))\tag{1}\\ \end{align}\]

其中\(x(n)\)表示輸入訊號，\(a_i(n)\)表示第i條路徑上的衰減係數，\(\tau_i(n)\)為第i條路經上的傳播時延。

由於表示的通道是線性通道，故可以用在\(n\)時刻對\(n-\tau\)時刻發射的衝激的響應\(h(\tau,n)\)來表示。我們已知用\(h(\tau,n)\)表示的經過通道的輸入\輸出為卷積關係：

\[y(n) =\sum_\tau h(\tau,n) \cdot x(n-\tau) \tag{2} \]

於是由上述兩個公式我們可以推得多徑衰落通道衝激響應的數學表示式為：

\[h(\tau,n) =a_i(n) \cdot \delta(\tau - \tau_i(n)) \tag{3} \]

瑞利隨機變數產生補充：

在一般的衰落環境中，無線衰落通道可以由復高斯隨機變數\(W1 + jW2\)表示，其中\(W1\)和\(W2\)都是均值為0，方差為\(\delta^2\)的獨立同分布(i.i.d.)高斯隨機變數。

如何產生瑞利隨機變數呢？首先通過MATLAB內建函式randn()產生均值為0，方差為1的兩個i.i.d.高斯隨機變數\(W1\)和\(W2\)。瑞利隨機變數X為：

\[X = \delta \cdot \sqrt{W1^2 + W2^2} \tag{4} \]

所以一旦通過內建函式randn()生成好了\(W1\)和\(W2\)，就可以由公式(4)生成平均功率為\(E(X^2) = 2\delta^2\)的瑞利隨機變數。

在模擬中我們已經提前給出了瑞利通道平均功率\(P_h\)，所以有\(2\delta^2 = p_h\)，推出：

\[\delta = \sqrt{p_h / 2} \tag{5} \]

程式碼：

A_h = (randn(1,NH) + 1i * randn(1,NH)) .* sqrt(P_h / 2); % 由公式(4)(5)生成瑞利隨機變數

h = zeros(1,LH);        % 初始化通道衝激響應模型
h(D_h + 1) = A_h;       % 通道衝激響應(同時體現出衰減係數和通道時延)，公式(3)的程式碼體現
xt1 = conv(xt,h);       % 卷積，輸出通過該通道的訊號，公式(2)的程式碼體現

步驟5>

經過上一步的處理，現在考慮模擬新增高斯白噪聲。由於snr在程式開頭就已經確定好了，所以我們要根據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))); % 根據噪聲標準差新增噪聲

步驟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.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.

相關文章