| 初识OFDM(六)：从零开始的OFDM误码率仿真 | 您所在的位置：网站首页 › 基于matlab的卷积码仿真与分析软件 › 初识OFDM(六)：从零开始的OFDM误码率仿真 |
初识OFDM(六)：从零开始的OFDM误码率仿真
|
文章目录
初识OFDM(六)：从零开始的OFDM误码率仿真零.代码地址一. 加性高斯白噪声对OFDM误码率的影响1. 代码展示2. 代码分析fftshift和ifftshift能量和信噪比问题
二.瑞利信道对OFDM误码率的影响1. 代码展示2. 代码分析瑞利衰落信道是如何通过TDL模型仿真而成的线性卷积，循环卷积和均衡线性卷积输入和输出长度怎么不相等了？
三. 一些还没有思考清楚的问题
初识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窗对准的是 [ − N k 2 , N k 2 ] [-\frac{Nk}{2},\frac{Nk}{2}] [−2Nk,2Nk]这样一段信号，而我们把他放在matlab中进行仿真时，下标则是 [ 1 , N k ] [1,Nk] [1,Nk].因此，这就对应错了啊，因此我们要把前后倒置一下，这才满足理论上的要求。相应的，做完fft后，我们又要把结果恢复到 [ 1 , N k ] [1,Nk] [1,Nk]，这样才方便和输入进行对比。 如何做？注意是顺序是先做ifftshift再做ifft，而接收端是先fft再做fftshift，原因上面已经说的很清楚啦。 能量和信噪比问题 IFFT和FFT前后的能量变化想象一下FFT或是IFFT的公式. IFFT中每一个 x [ n ] x[n] x[n] 是Nk个 1 N k X [ k ] \frac{1}{Nk}X[k] Nk1X[k]组成的,因此能量变为 1 N k 2 × N k = 1 N k \frac{1}{Nk^2}\times Nk=\frac{1}{Nk} Nk21×Nk=Nk1 FFT 中每一个 X [ k ] X[k] X[k]是Nfft个 x [ n ] x[n] x[n]组成的，因此能量变为 N f f t Nfft Nfft倍 能量是会产生变化的，而噪声是在信道时域中加上的，因此我们信噪比使用了时域的能量来计算。 EsN0和EbN0的关系https://blog.csdn.net/sinat_38151275/article/details/79869891可以参考这篇文章 总的来说，当经过数字调制后，一个符号包含了Npsk个比特那么 E s N 0 = E b N 0 = 10 l o g 10 ( N p s k ) EsN0=EbN0=10log_{10}(Npsk) EsN0=EbN0=10log10(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是如何表达瑞利信道的了。 线性卷积，循环卷积和均衡我们知道，在FFT中，循环卷积才是真正的卷积，各种性质都是在循环卷积中证明完成的，这其中当然包括大名鼎鼎的时域的卷积等于频域的点乘，也就是说，时域的循环卷积核频域的点乘结果是相同的。可是信道和信号在时域发生的实际上是线性卷积，这就给我们均衡的时候制造困难了，没有简便的算法啊! 我们在DSP课程中学过用循环卷积去计算线性卷积，方法就是把循环卷积的L扩展到M+N-1时候计算。那么有没有方法可以把线性卷积转化为循环卷积呢？当然有！这就是CP 具体原理可以看https://blog.51cto.com/u_15127585/2669966 这是非常巧妙的
也就是说CP有两个作用 在两个符号之间填充一点东西，这样可以消除掉ISI，当然这个作用用ZP就可以完成实现循环卷积，大大的简化了信道的均衡还有一个就是估计STO和CFO，这可以在前面的文章看到 线性卷积输入和输出长度怎么不相等了？https://tieba.baidu.com/p/5402720283?red_tag=0784237449 这居然也被我查到了 这正是TDL想要的效果，我们多出来的实际上就是TDL产生的最大时延/ T s T_s Ts 三. 一些还没有思考清楚的问题 虚拟子载波放置在什么位置？他和ifftshift的关系是什么？使用ZP时如何解决信道均衡的问题？如何将线性卷积转化为循环卷积？误码率的误差来自哪里？ |
【本文地址】