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教程