OFDM通訊系統的MATLAB仿真(2)

時間  2020-07-23
標籤 ofdm 通訊系統 matlab 仿真 欄目 MATLAB 简体版
原文   原文鏈接

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

注意：爲了整個系統的完整性，第一篇隨筆中的每一個步驟這裏也都又寫了一遍，但省略了補充知識部分，在第一篇的基礎上添加了實現多徑衰落的部分。想要看信噪比計算和噪聲功率計算的同窗能夠去看第一篇隨筆。less

關於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是覆信號。spa

另外在步驟1咱們還要進行信噪比的一些初始化，便於計算噪聲幅度和最後的計算比特誤碼率。3d

增長部分：

在步驟1中，咱們增長對信道特徵的初始化工做。主要是假設多徑信道個數和信道功率，以及各信道的時延，爲以後信號經過多徑信道的計算作準備。code

代碼：

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調製，格雷碼星座圖，並歸一化

步驟二、三、4>

接下來的三個步驟分別以下，注意都是一個符號一個符號處理的，可回去看最開始的符號變化圖：blog

  • 將每一個OFDM符號的前一半和後一半交換，至於爲何要作交換，我仍然不是很懂。有大佬知道的話但願能在評論區指導一下，感激涕零！
  • 對交換事後的每一個OFDM符號作IFFT，記錄輸出爲x1(n)。
  • 對每一個OFDM符號添加循環前綴CP，實際操做很簡單，由於這裏設的CP的長度NCP爲16。就是把每一個符號的後16個採樣點添加到當前符號的最前面來，每一個符號所以就變成了64+16=80個採樣點。

因爲這三個步驟都是在一個循環裏處理的，因此我也就把步驟二、三、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))); % 根據噪聲標準差添加噪聲

步驟六、7>

如今的信號已是通過多徑瑞利衰落而且添加了高斯白噪聲的信號，不容易啊！咱們的仿真已經完成了一半。接下來的兩個步驟與步驟二、三、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.

相關文章
相關標籤/搜索
每日一句
    每一个你不满意的现在，都有一个你没有努力的曾经。
本站公眾號
   歡迎關注本站公眾號,獲取更多信息
聯繫我們 最近搜索 最新文章 沪ICP备13005482号-10 MyBatis教程 SQL 教程 MySQL教程 Java 教程 Thymeleaf 教程 Hibernate教程 Spring教程 Redis教程