关于 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.
来源: https://www.cnblogs.com/gjblog/p/13363995.html