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