DCO-OFDM信道采用QPSK与16QAM的蒙特卡洛仿真与误码率分析
可见光通信关于16QAM和QPSK在OFDM下的仿真与误码率分析
clear all;
close all;
clc;
N = 64; %子载波数
N_data_symbol = N/2-1; %调制码元数
CP = N/4; %循环前缀大小
M = 4; %进制
BitperSymbol = log2(M); %每个码元的信息量
N_Iteration = 1000; %迭代次数(理解成实验次数吧大概)
SNR_dB = 0:1:15; %信噪比
SNR_pow = 10.^(SNR_dB/10); %信噪比倍率形式
Nsym = 256; %OFDM块数/帧
Es = 1; %码元能量归一化
Eb = Es/BitperSymbol; %每比特的能量
K = 3.2; %直流偏置
DC_QAM = K*2;%%
%16QAM
M_QAM = 16; %调制阶数
k_QAM = log2(M_QAM);
QAM_16_rand_number = (N/2-1)*k_QAM;
power_QAM = sqrt((2+10+18+10)/4); %16QAM的平均符号功率，用于后续能量归一化
% power_QPSK = sqrt((1*4+0+2*4)/9); %QPSK的平均符号功率
power_QPSK = sqrt(2); %QPSK的平均符号功率
Eb_MQAM = Es/k_QAM; %16QAM的码元能量for snr=0:1:15 Eb_N0=10.^(snr/10);N0=Eb./Eb_N0; %QPSK的噪声（能量归一化）N0_MQAM = Eb_MQAM/Eb_N0; %16QPSK的噪声（能量归一化）%参数初始化QPSK = zeros(N_Iteration,N_data_symbol); %QPSK调制后QPSK_H = zeros(N_Iteration,N); %施密特变换后QPSK_after_IFFT = zeros(N_Iteration,N); %IFFT后QPSK_after_DC = zeros(N_Iteration,N); %直流偏置后QPSK_add_CP = zeros(N_Iteration,N+CP); %添加CP后QPSK_after_AWGN = zeros(N_Iteration,N+CP); %经过AWGN后QPSK_decrease_CP = zeros(N_Iteration,N); %去掉CP后QPSK_after_FFT = zeros(N_Iteration,N); %FFT变换后QPSK_after_demodulation = zeros(N_Iteration,N-2); %解调后QPSK_BER = zeros(1,N_Iteration); %计算误码率MQAM = zeros(N_Iteration,N_data_symbol);MQAM_H = zeros(N_Iteration,N);MQAM_after_IFFT = zeros(N_Iteration,N);MQAM_after_DC = zeros(N_Iteration,N);MQAM_add_CP = zeros(N_Iteration,N+CP);MQAM_after_AWGN = zeros(N_Iteration,N+CP);MQAM_decrease_CP = zeros(N_Iteration,N);MQAM_after_FFT = zeros(N_Iteration,N);MQAM_after_demodulation = zeros(N_Iteration,N_data_symbol);MQAM_BER = zeros(1,N_Iteration);for k = 1:N_Iteration%% Input bit streams for modulation % 步骤一：编写子函数，将产生的二进制比特数调制成 M-QAM 码NRZ = round(rand(1,N)); %QPSK用的N位随机二进制数，后面只取N-2位，这里懒得改了NRZ_QAM = round(15*rand(1,N_data_symbol)); %QAM的31个10进制数范围0-15，因为qammod函数是0-15之间的数调制NRZ_QAM_bit = dec2bin(NRZ_QAM,4); %十进制转二进制方便误码率计算，每一个元素变四位,但是类型是charNRZ_QAM_bit_1 = double(bin2dec(NRZ_QAM_bit(:)));%将char类型的先平铺成一列，再转doubleQPSK(k,:) = Modulation_bit(NRZ,N)/power_QPSK;
% MQAM(k,:) = qammod(NRZ_QAM,M_QAM); %0-15范围的十进制数调制，第二个参数是MMQAM(k,:) = qammod(NRZ_QAM,M_QAM)/power_QAM; %0-15范围的十进制数调制，第二个参数是M%% Hermitian Symmetry % 步骤二：编写子函数，将QAM信号变换为具有赫密特对称的形式QPSK_H(k,:) = Hermitain_sym(QPSK(k,:)); MQAM_H(k,:) = Hermitain_sym(MQAM(k,:));%% IFFT% 步骤三：编写子函数，将具有赫密特对称的复数信号进行 IFFT 后，产生双极性时域信号QPSK_after_IFFT(k,:) = do_IFFT(QPSK_H(k,:),N); MQAM_after_IFFT(k,:) = do_IFFT(MQAM_H(k,:),N); %% add DC and clipping% 步骤四：编写子函数，产生直流偏置，把直流偏置叠加到时域信号，并截取负信号，让剩余负信号变成 0QPSK_after_DC(k,:) = DC_Clip(QPSK_after_IFFT(k,:),K);MQAM_after_DC(k,:) = DC_Clip(MQAM_after_IFFT(k,:),DC_QAM);%% Guard Interval insertion and the CP addition% 步骤五：编写子函数，时域信号增加循环前缀QPSK_add_CP(k,:) = add_cp(QPSK_after_DC(k,:),CP);MQAM_add_CP(k,:) = add_cp(MQAM_after_DC(k,:),CP);%% Received signal with AWGN% 步骤六：编写子函数，考虑 AWGN 信道，接收信号只受到噪声影响QPSK_after_AWGN(k,:) = Channel(QPSK_add_CP(k,:),N0,N+CP);MQAM_after_AWGN(k,:) = Channel(MQAM_add_CP(k,:),N0_MQAM,N+CP);%% Remove CP % 步骤七：编写子函数，去掉接收信号循环前缀QPSK_decrease_CP(k,:) = Remove_CP(QPSK_after_AWGN(k,:),CP);MQAM_decrease_CP(k,:) = Remove_CP(MQAM_after_AWGN(k,:),CP);%% FFT % 步骤八：编写子函数，接收信号经过 FFT 后变成频域信号QPSK_after_FFT(k,:) = do_FFT(QPSK_decrease_CP(k,:),N)*power_QPSK;%把归一化的功率乘回来解调用
% MQAM_after_FFT(k,:) = do_FFT(MQAM_decrease_CP(k,:),N);MQAM_after_FFT(k,:) = do_FFT(MQAM_decrease_CP(k,:),N)*power_QAM;%% demodulation% 步骤九：编写子函数，把频域信号在每一个载波上进行解调QPSK_after_demodulation(k,:) = demodulation(QPSK_after_FFT(k,:),N);MQAM_after_demodulation(k,:) = qamdemod(MQAM_after_FFT(k,2:N/2),M_QAM);% 取2-32位进行解调%% BER calculation% 步骤十：编写子函数，计算系统误码率QPSK_BER(k) = BER_calculation(QPSK_after_demodulation(k,:),NRZ);MQAM_BER(k) = BER_calculation_QAM(MQAM_after_demodulation(k,:),NRZ_QAM_bit_1,k_QAM);endQPSK_BER_per_snr(snr+1) = mean(QPSK_BER);MQAM_BER_per_snr(snr+1) = mean(MQAM_BER);
end
%%
%理论误码率
parameter1 = 4*(sqrt(M)-1)/(sqrt(M)*log2(M));
parameter2 = qfunc(sqrt(3*log2(M)*SNR_pow/(M-1)));
parameter3 = 4*(sqrt(M)-2)/(sqrt(M)*log2(M));
parameter4 = qfunc(3*sqrt(3*log2(M)*SNR_pow/(M-1)));
BER_theory_QPSK = parameter1*parameter2 + parameter3*parameter4;parameter_QAM1 = 4*(sqrt(M_QAM)-1)/(sqrt(M_QAM)*log2(M_QAM));
parameter_QAM2 = qfunc(sqrt(3*log2(M_QAM)*SNR_pow/(M_QAM-1)));
parameter_QAM3 = 4*(sqrt(M_QAM)-2)/(sqrt(M_QAM)*log2(M_QAM));
parameter_QAM4 = qfunc(3*sqrt(3*log2(M_QAM)*SNR_pow/(M_QAM-1)));
BER_theory_MQAM = parameter_QAM1*parameter_QAM2 + parameter_QAM3*parameter_QAM4;%%
% 画图
figure;
semilogy(SNR_dB,QPSK_BER_per_snr,'r-*'); % 此处添加自己定义的 BER 变量
hold on
grid on
semilogy(SNR_dB,BER_theory_QPSK,'b--'); % 此处绘制解析 BER
hold on
semilogy(SNR_dB,MQAM_BER_per_snr,'g-*'); % 此处添加自己定义的 BER 变量
hold on
semilogy(SNR_dB,BER_theory_MQAM,'m--'); % 此处绘制解析 BER
xlabel('EbNo(dB)');
ylabel('BER')
title('BER for DCO-OFDM');
legend('QPSK仿真值','QPSK理论值','MQAM仿真值','MQAM理论值');%%
%自定义的函数
function [QAM] = Modulation_bit(NRZ,N)
%QAM调制
QAM = zeros(1,N/2-1); %初始化
for k_r = 1:2:N %实部if NRZ(k_r) == 0x((k_r+1)/2) = -1;elsex((k_r+1)/2) = 1;end
end
for k_i = 2:2:N %虚部if NRZ(k_i) == 0y(k_i/2) = -1;elsey(k_i/2) = 1;end
end
z = x+y*sqrt(-1);
s = z(1:end-1);
QAM = s;
endfunction [QAM_H] = Hermitain_sym(QAM)
%施密特变换
QAM_H = [0 QAM 0 flip(conj(QAM))];%flip是倒序排列，conj是取共轭
endfunction [QAM_after_IFFT] = do_IFFT(QAM_H,N);
%傅里叶逆变换
QAM_after_IFFT = ifft(QAM_H,2^nextpow2(N))*sqrt(N);
% QAM_after_IFFT = ifft(QAM_H,2^nextpow2(N));
%第二个参数是长度
%scaled是sqrtN，coeff是N，实数默认N，复数默认1
endfunction [QAM_after_DC] = DC_Clip(QAM_after_IFFT,DC)
%直流偏置
% QAM_after_DC = QAM_after_IFFT+DC*sqrt(mean(QAM_after_IFFT.^2));
QAM_after_DC = QAM_after_IFFT+DC
for k_DC = 1:1:length(QAM_after_DC)if QAM_after_DC(k_DC) < 0QAM_after_DC(k_DC) = 0;end
end
endfunction [QAM_add_CP] = add_cp(QAM_after_DC,CP)
%添加循环前缀
QAM_add_CP = [QAM_after_DC(length(QAM_after_DC)-CP+1:length(QAM_after_DC)),QAM_after_DC];
endfunction [QAM_after_AWGN] = Channel(QAM_add_CP,N0,N_CP)
%AWGN信道
%QAM_after_AWGN = awgn(QAM_add_CP,SNR,'measured');
% QAM_after_AWGN = awgn(QAM_add_CP,sqrt(SNR));
QAM_after_AWGN = QAM_add_CP+normrnd(0,sqrt(N0),1,N_CP) %添加高斯随机数，第二个参数是标准差因此要开方
endfunction [QAM_decrease_CP] = Remove_CP(QAM_after_AWGN,CP)
%删去前缀CP
QAM_decrease_CP = QAM_after_AWGN(CP+1:end);
endfunction [QAM_after_FFT] = do_FFT(QAM_decrease_CP,N)
%傅里叶变换
QAM_after_FFT = fft(QAM_decrease_CP,2^nextpow2(N))/sqrt(N);%除以根号N归一化
% QAM_after_FFT = fft(QAM_decrease_CP,2^nextpow2(N));
QAM_after_FFT = round(QAM_after_FFT.*100)/100;%保留两位小数用的
endfunction [QAM_after_demodulation] = demodulation(QAM_after_FFT,N)
%对每个子载波解调
User = QAM_after_FFT(2:N/2);%取2~32位准备解调
QAM_after_demodulation = zeros(1,N-2);
%解调
for count5 = 1:1:length(User)if real(User(count5)) > 0QAM_after_demodulation(count5*2-1) = 1;elseQAM_after_demodulation(count5*2-1) = 0;endif imag(User(count5)) > 0QAM_after_demodulation(count5*2) = 1;elseQAM_after_demodulation(count5*2) = 0;end
end
endfunction [QAM_BER_0] = BER_calculation(QAM_after_demodulation,NRZ)
%计算误码率
NRZ_real = NRZ(1:end-2);
error = sum(abs(NRZ_real-QAM_after_demodulation));%原始01数据和解调数据相减，去绝对值，求和为错误总和
QAM_BER_0 = error/length(NRZ_real);
endfunction [QAM_BER_0] = BER_calculation_QAM(QAM_after_demodulation,NRZ,k_QAM)
%计算误码率QAM_after_demodulation_2 = dec2bin(QAM_after_demodulation,4); %把每个十进制数转为4个二进制数
demodulation_bit = double(bin2dec(QAM_after_demodulation_2(:)));%转浮点数后
error = sum(abs(NRZ - demodulation_bit));%原始01数据和解调数据相减，去绝对值，求和为错误总和
QAM_BER_0 = error/length(NRZ);
end结果图如下，基于运算精度和速度的问题，数据量量级为10e6左右，因此误码率在10e-5基本就是极限了，再低就是一个错误码元都没有了
本文来自互联网用户投稿，文章观点仅代表作者本人，不代表本站立场，不承担相关法律责任。如若转载，请注明出处。 如若内容造成侵权/违法违规/事实不符，请点击【内容举报】进行投诉反馈！