初识OFDM(六)：从零开始的OFDM误码率仿真-Toy模板网

这篇具有很好参考价值的文章主要介绍了初识OFDM(六)：从零开始的OFDM误码率仿真。希望对大家有所帮助。如果存在错误或未考虑完全的地方，请大家不吝赐教，您也可以点击"举报违法"按钮提交疑问。

初识OFDM(六)：从零开始的OFDM误码率仿真

零.代码地址

https://github.com/liu-zongxi/OFDM_simulation

代码参考了https://zhuanlan.zhihu.com/p/385096476陈老湿的仿真，但各个函数我都有重新实现，希望写的更规范一些

一. 加性高斯白噪声对OFDM误码率的影响

1. 代码展示

该部分代码来自工程的OFDM_awgn.m

整个OFDM经历

生成一帧 $\rightarrow$ 串并转换 $\rightarrow$ 调制 $\rightarrow$ ifftshift $\rightarrow$ ifft $\rightarrow$ 添加CP $\rightarrow$ 并串转换 $\rightarrow$ 计算能量，信噪比，添加AWGN $\rightarrow$ 接收端串并转换 $\rightarrow$ 去CP $\rightarrow$ FFT $\rightarrow$ FFTshift $\rightarrow$ 解调 $\rightarrow$ 并串转换

%-----------------------仿真OFDM---------------------------%
%-----------------------author:lzx-------------------------%
%% 设置参数
clear;clc;
Nk = 128;           % 子载波个数
Nfft = 128;          % fft长度
Nframe = 6;         % 一帧中有几个OFDM符号
M = 2;              % 调制符号所含比特
SR = 250000;        % 符号速率
BR = SR .* M;       % 比特率
NGI = 32;           % 保护间隔长度
EbN0s = 3:1:10;      % 信噪比
Nsym = Nfft+NGI;    % 系统长度
bers = zeros(1,length(EbN0s));
fprintf('EbN0 \t \t ber\t\t\t per\t\t\t nloop \t\t \n');
%% 函数主体

for kk = 1:length(EbN0s)
    % rng('default')          % 初始化随机种子
    EbN0 = EbN0s(kk);
    nloop = 10000;          % 发送多少帧
    n_biterror = 0;         % 错误的数据
    n_bitdata = 0;          % 一共发送了多少数据
    n_packeterror = 0;      % 有多少错误帧
    n_packetdata = 0;       % 发送了多少帧
    for ii = 1:nloop
        % 生成一帧数据，串并转换，并QPSK，生成一帧
        frame_FDserial = rand(1,Nk*Nframe*M) > 0.5;     % 发送的是bit
        frame_FDparallel = reshape(frame_FDserial,Nk,Nframe*M);% 串并转换
        frame_mod = QPSKMod(frame_FDparallel,Nk,Nframe);     %调制
        % IFFT
        power_FT = sum(sum(abs(frame_mod).^2))/Nk/Nframe;  % 计算下IFFT前的能量，FT表示频域
        frame_mod_shift = ifftshift(frame_mod);         % 频域归零
        frame_ifft = ifft(frame_mod_shift, Nfft);             % ifft
        power_TD = sum(sum(abs(frame_ifft).^2))/Nk/Nframe; % 计算下IFFT前的能量，DT表示时域
        % 添加保护间隔
        frame_withGI = AddGI(frame_ifft, Nfft, NGI, Nframe, "CP");  % 添加保护间隔
        % 并串转换
        frame_TDserial = reshape(frame_withGI,1,Nsym*Nframe);
        x=1:1:160;
        % hold on;
        % plot(x, frame_TDserial(1:160),'b');
        % Channel
        power_TDserial = sum(abs(frame_TDserial).^2)/Nsym/Nframe;   
        EsN0 = EbN0 + 10*log10(M);                                  % 根据信噪比计算噪声能量，幅值，然后加在信号上
        N0 = power_TDserial .* 10.^(-EsN0/10);
        noise_msg = sqrt(N0 / 2) .* (randn(size(frame_TDserial)) + 1j * randn(size(frame_TDserial)));
        frame_recieved = frame_TDserial + noise_msg;
%         attn = sqrt(0.5*power_TDserial*SR/BR*10.^(-EbN0/10));
%         noise_msg = attn .* (randn(size(frame_TDserial)) + 1j * randn(size(frame_TDserial)));
%         frame_recieved = frame_TDserial + noise_msg;
        % plot(x, noise_msg(1:160),'r');
        % hold off;
        % 接收端，串并转换
        frame_recieved_parallel = reshape(frame_recieved,Nsym,Nframe);
        % 去GI
        frame_noGI = RemoveGI(frame_recieved_parallel, Nfft, NGI);
        % FFT
        frame_recieved_FD_shift = fft(frame_noGI, Nfft);
        frame_recieved_FD = fftshift(frame_recieved_FD_shift);
        % QPSK解调
        frame_demod = QPSKDemod(frame_recieved_FD, Nk, Nframe);
        % 并串转换
        frame_output = reshape(frame_demod, 1, Nk*Nframe*M);
       
        % 计算error
        n_biterror_tmp = sum(abs(frame_output-frame_FDserial));
        n_bitdata_tmp = length(frame_FDserial);
        n_biterror = n_biterror + n_biterror_tmp;
        n_bitdata = n_bitdata + n_bitdata_tmp;
        if n_biterror_tmp ~= 0
            n_packeterror = n_packeterror + 1;
        end
        n_packetdata = n_packetdata + 1;
    end
    % 计算在当前信噪比下的误码率
    per = n_packeterror/n_packetdata;
    ber = n_biterror/n_bitdata;
    bers(kk)=ber;
    fprintf('%f\t%e\t%e\t%d\t\n',EbN0,ber,per,nloop);
end
save("BERofdm.mat",'bers');
%% 画图
% bers = load("BERofdm.mat").bers;
awgn_theory = [0.0228784075610853,0.0125008180407376,0.00595386714777866,0.00238829078093281,0.000772674815378444,0.000190907774075993,3.36272284196176e-05,3.87210821552205e-06];
semilogy(EbN0s,awgn_theory,'-*',EbN0s,bers,'-+');
xlabel('比特信噪比');
ylabel('误码率');
title('不同信噪比下误码率仿真曲线');
legend('理论曲线','实验曲线');

2. 代码分析

这个代码主要有这么几点是值得注意的

fftshift和ifftshift

ifftshift和fftshift到底是什么呢？

简单来说，ifftshift是把矩阵做一个向左的长度为矩阵一半的循环位移，fftshift则是把矩阵做一个向右的长度为矩阵一半的循环位移。

作用是什么呢？

当做ifft时，我们希望ifft窗对准的是 $[-\frac{Nk}{2},\frac{Nk}{2}]$ 这样一段信号，而我们把他放在matlab中进行仿真时，下标则是 $[1, N k]$ .因此，这就对应错了啊，因此我们要把前后倒置一下，这才满足理论上的要求。相应的，做完fft后，我们又要把结果恢复到 $[1, N k]$ ，这样才方便和输入进行对比。

如何做？

注意是顺序是先做ifftshift再做ifft，而接收端是先fft再做fftshift，原因上面已经说的很清楚啦。

能量和信噪比问题

IFFT和FFT前后的能量变化

想象一下FFT或是IFFT的公式.

IFFT中每一个 $x [n]$ 是Nk个 $\frac{1}{Nk}X[k]$ 组成的,因此能量变为 $\frac{1}{Nk^2}\times Nk=\frac{1}{Nk}$

FFT 中每一个 $X [k]$ 是Nfft个 $x [n]$ 组成的，因此能量变为 $N f f t$ 倍

能量是会产生变化的，而噪声是在信道时域中加上的，因此我们信噪比使用了时域的能量来计算。

EsN0和EbN0的关系

https://blog.csdn.net/sinat_38151275/article/details/79869891可以参考这篇文章

总的来说，当经过数字调制后，一个符号包含了Npsk个比特那么
$EsN0=EbN0=10log_{10}(Npsk)$

CP对信噪比的影响？

https://zhuanlan.zhihu.com/p/281601152，陈老湿在这里做了很多探讨，不过我更偏向于将CP不看做冗余

因此我的代码是这样的

        power_TDserial = sum(abs(frame_TDserial).^2)/Nsym/Nframe;

二.瑞利信道对OFDM误码率的影响

1. 代码展示

该部分代码来自工程的OFDM_fading.m

整个OFDM经历

生成一帧 $\rightarrow$ 串并转换 $\rightarrow$ 调制 $\rightarrow$ ifftshift $\rightarrow$ ifft $\rightarrow$ 添加CP $\rightarrow$ 并串转换 $\rightarrow$ 经历信道衰落，和信道卷积 $\rightarrow$ 计算能量，信噪比，添加AWGN $\rightarrow$ 接收端串并转换 $\rightarrow$ 去CP $\rightarrow$ FFT $\rightarrow$ FFTshift $\rightarrow$ $\rightarrow$ FFT信道矩阵h $\rightarrow$ 信道均衡 $\rightarrow$ 解调 $\rightarrow$ 并串转换

%-----------------------用TDL仿真多径衰落-------------------%
%-----------------------author:lzx-------------------------%
%-----------------------date:2022年3月25日-----------------%
%% 设置参数
clear;clc;
Nk = 128;           % 子载波个数
Nfft = 128;          % fft长度
Nframe = 6;         % 一帧中有几个OFDM符号
Npsk = 2;              % 调制符号所含比特
M = 2^Npsk;          % 调制数
SR = 250000;        % 符号速率
BR = SR .* Npsk;       % 比特率
NGI = 32;           % 保护间隔长度
EbN0s = 0:2:20;      % 信噪比
Nsym = Nfft+NGI;    % 系统长度
bers = zeros(1,length(EbN0s));  % 误码率储存数组
PowerTDL_dB = [0 -8 -17 -21 -25];   % TDL中信道抽头的功率,dB为单位
Delay = [0 3 5 6 8];                % TDL中信道时延
PowerTDL = 10.^(PowerTDL_dB/10);    % TDL中信道抽头的功率
Nchannel=length(PowerTDL_dB);       % 信道抽头数
Tau_maxTDL = Delay(end)+1;          % 最大时延除以帧长,就是归一化的最大时延
fprintf('EbN0 \t \t ber\t\t\t per\t\t\t nloop \t\t \n');
%% 函数主体

for kk = 1:length(EbN0s)
    % rng('default')          % 初始化随机种子
    EbN0 = EbN0s(kk);
    nloop = 10000;          % 发送多少帧
    n_biterror = 0;         % 错误的数据
    n_bitdata = 0;          % 一共发送了多少数据
    n_packeterror = 0;      % 有多少错误帧
    n_packetdata = 0;       % 发送了多少帧
    for ii = 1:nloop
%--------------------------发射端-------------------------------%
        % 生成一帧数据，串并转换，并QPSK，生成一帧
        frame_FDserial = rand(1,Nk*Nframe*Npsk) > 0.5;     % 发送的是bit
        frame_FDparallel = reshape(frame_FDserial,Nk,Nframe*Npsk);% 串并转换
        frame_mod = QPSKMod(frame_FDparallel,Nk,Nframe);     %调制
        % IFFT
        power_FT = sum(sum(abs(frame_mod).^2))/Nk/Nframe;  % 计算下IFFT前的能量，FT表示频域
        frame_mod_shift = ifftshift(frame_mod);         % 频域归零
        frame_ifft = ifft(frame_mod_shift, Nfft);             % ifft
        % frame_ifft = ifft(frame_mod, Nfft);
        power_TD = sum(sum(abs(frame_ifft).^2))/Nk/Nframe; % 计算下IFFT前的能量，DT表示时域
        % 添加保护间隔
        frame_withGI = AddGI(frame_ifft, Nfft, NGI, Nframe, "CP");  % 添加保护间隔
        % 并串转换
        frame_TDserial = reshape(frame_withGI,1,Nsym*Nframe);
            % x=1:1:160;
            % hold on;
            % plot(x, frame_TDserial(1:160),'b');
%--------------------------Channel-------------------------------%
        % 信号先经历衰落
        channel = Rayleigh_model(Nchannel, PowerTDL);
        h = zeros(1, Tau_maxTDL);
        h(Delay+1) = channel;
        frame_conv = conv(frame_TDserial, h);
        frame_fading = frame_conv(:,1:length(frame_TDserial));        % 看似是线性卷积，实际上由于CP变成了循环卷积
        % 添加高斯白噪声
        power_TDserial = sum(abs(frame_TDserial).^2)/Nk/Nframe;     % 计算出的能量和理论不符啊，没发现问题在哪
        EsN0 = EbN0 + 10*log10(Npsk);                                  % 根据信噪比计算噪声能量，幅值，然后加在信号上
        N0 = power_TDserial .* 10.^(-EsN0/10);
        noise_msg = sqrt(N0 / 2) .* (randn(size(frame_TDserial)) + 1j * randn(size(frame_TDserial)));
        frame_recieved = frame_fading + noise_msg;
        % 陈老湿方法添加高斯白噪声，本质上是一样的
%       attn = sqrt(0.5*power_TDserial*SR/BR*10.^(-EbN0/10));
%       noise_msg = attn .* (randn(size(frame_TDserial)) + 1j * randn(size(frame_TDserial)));
%       frame_recieved = frame_TDserial + noise_msg;
            % plot(x, noise_msg(1:160),'r');
            % hold off;
%--------------------------接收端-------------------------------%
        % 接收端，串并转换
        frame_recieved_parallel = reshape(frame_recieved,Nsym,Nframe);
        % 去GI
        frame_noGI = RemoveGI(frame_recieved_parallel, Nfft, NGI);
        % FFT
        frame_recieved_FD_shift = fft(frame_noGI, Nfft);
        frame_recieved_FD = fftshift(frame_recieved_FD_shift);
        % frame_recieved_FD = fft(frame_noGI, Nfft);
        % 信道均衡
        H = fftshift(fft([h zeros(1, Nfft-Tau_maxTDL)].', Nfft));
        frame_equalization = frame_recieved_FD ./ repmat(H, 1, Nframe);
        % QPSK解调
 
        frame_demod = QPSKDemod(frame_equalization, Nk, Nframe);
        % 并串转换
        frame_output = reshape(frame_demod, 1, Nk*Nframe*Npsk);
       
        % 计算error
        n_biterror_tmp = sum(abs(frame_output-frame_FDserial));
        n_bitdata_tmp = length(frame_FDserial);
        n_biterror = n_biterror + n_biterror_tmp;
        n_bitdata = n_bitdata + n_bitdata_tmp;
        if n_biterror_tmp ~= 0
            n_packeterror = n_packeterror + 1;
        end
        n_packetdata = n_packetdata + 1;
    end
    % 计算在当前信噪比下的误码率
    per = n_packeterror/n_packetdata;
    ber = n_biterror/n_bitdata;
    bers(kk)=ber;
    fprintf('%f\t%e\t%e\t%d\t\n',EbN0,ber,per,nloop);
end
save("BERofdm_rayleigh.mat",'bers');
%% 画图
% bers = load(!"BERofdm_rayleigh.mat").bers;
rayleigh_theory = 0.5.*(1-(1-(1./(10.^(EbN0s./10).*(48/80)+1))).^0.5);
semilogy(EbN0s,rayleigh_theory,'-*',EbN0s,bers,'-+');
xlabel('比特信噪比');
ylabel('误码率');
title('不同信噪比下误码率仿真曲线');
legend('理论曲线','实验曲线');

2. 代码分析

我们把信道卷积的代码拿出来看一下

PowerTDL_dB = [0 -8 -17 -21 -25];   % TDL中信道抽头的功率,dB为单位
Delay = [0 3 5 6 8];                % TDL中信道时延
PowerTDL = 10.^(PowerTDL_dB/10);    % TDL中信道抽头的功率
Nchannel=length(PowerTDL_dB);       % 信道抽头数
Tau_maxTDL = Delay(end)+1;          % 最大时延除以帧长,就是归一化的最大时延


% 信号先经历衰落
channel = Rayleigh_model(Nchannel, PowerTDL);
h = zeros(1, Tau_maxTDL);
h(Delay+1) = channel;
frame_conv = conv(frame_TDserial, h);
frame_fading = frame_conv(:,1:length(frame_TDserial));        % 看似是线性卷积，实际上由于CP变成了循环卷积

% 信道均衡
H = fftshift(fft([h zeros(1, Nfft-Tau_maxTDL)].', Nfft));
frame_equalization = frame_recieved_FD ./ repmat(H, 1, Nframe);

这里其实就涵盖了所有第二个仿真和第一个仿真的区别了，这个代码主要有这么几点是值得注意的

瑞利衰落信道是如何通过TDL模型仿真而成的

TDL模型里只包含两个内容，第一个是信道对功率的变化，这是用瑞利分布乘以功率来完成的。第二个是时延，这是用信道的长度和线性卷积来完成的。我们看到代码中有一项Tau_maxTDL = Delay(end)+1; 这一般被称作信道长度。就是最大时延在TDL这个单位时延下的长度嘛。这样就能理解TDL是如何表达瑞利信道的了。