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

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

因爲是第一篇博客，想先說點廢話，其實本身早就想把學到的一些東西總結成文章隨筆之類的供本身複習時查看的了。可是一是以爲本身學的的不夠深刻，總結也寫不出什麼很深入的東西；二是以爲網上也有海量的資料了，須要時查一查根本不須要本身寫。可是偏偏也是網上的資料過於龐大，參差不齊，致使每次都如海水同樣的知識涌入腦中，最後也如蜻蜓點水通常了了看下，知識吸取率低的驚人。如今也準備改變一下觀念，儘可能把本身學過的東西概括整理，以隨筆的形式發出來，可能有些地方我還不能理解做者的作法，我也會記錄出來，懂的地方解釋清楚，不懂的地方也標記清楚，同時在以後的不斷學習和總結中補上以前挖的坑，強迫症寫文章真的是太難了~哭。
數組

不過也這樣安慰本身吧，將所學的進行輸出整理其實也是很重要的一步，看似時間浪費了好多，不過讀書人花的時間怎麼能叫浪費呢！整理出來本身看着不爽嗎。若是能同時幫助到其餘人的話，那就太好了！但願本身能堅持下去。
less

關於OFDM系統我目前參考的是《MIMO-OFDM無線通訊技術及MATLAB實現》這本書，我看到網上用這本書作參考的人還挺多的，這裏將將做者實現OFDM系統的思路及其代碼從新理順一遍。由於我理解的也不是很深刻，有不對的地方但願大佬能幫忙指正。若是是沒怎麼接觸過OFDM的萌新，這篇文章能夠幫助你對OFDM地仿真有個粗淺的瞭解。


函數

首先畫一個我我的認爲特別好理解的OFDM符號變化圖來幫助理解代碼，接下來我會詳細的介紹每一個步驟在幹什麼。
學習


步驟0>

先僞裝前景摘要一下。本OFDM系統仿真用到的技術主要有：16QAM調製解調 IFFT與FFT 添加AWGN噪聲，沒用到的有：信道編碼 交織 多徑瑞利信道 信道估計等等，哇，沒用到的技術還有這麼多，慚愧慚愧。這些技術的數學理論推導確實很難，可是在MATLAB仿真中每每用幾個自帶的函數就能解決問題，因此要實現一個簡單的OFDM系統仍是很容易的，不要被天花亂墜般恐怖的數學公式嚇跑了(因此我最喜歡的就是直接看代碼的運行過程，而後有時間再去研究數學推導2333)。


編碼

步驟1>

這個仿真好像暫時沒有時間的概念，單位是按照採樣點來的。假設一幀有三個OFDM符號，每一個符號長度爲64(恰好在步驟3作IFFT時長度也爲64，知足2的冪次方)。咱們首先生成數字基帶信號，信號長度爲192個採樣點，因爲要進行16QAM調製，咱們直接隨機生成192個16進制的數做爲基帶信號X(K)，而後再將X(K)通過16QAM星座圖映射便完成了調製。spa

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

代碼：

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}\)。code

信噪比補充：

S：信號功率

N：噪聲功率

W：傳輸信道帶寬

Eb：信號單位比特能量

Es：信號單位符號能量

Rb：信號比特傳輸速率

Rs：信號符號傳輸速率

N0：噪聲的功率譜密度

K：通訊系統的調製率，\(\text{K} = log_2 M\)
blog

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) \]


步驟二、三、4>

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

  • 將每一個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

步驟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'); % 官方函數直接加噪

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

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