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);           %十进制转二进制方便误码率计算，每一个元素变四位,但是类型是char
        NRZ_QAM_bit_1 = double(bin2dec(NRZ_QAM_bit(:)));%将char类型的先平铺成一列，再转double
        QPSK(k,:) = Modulation_bit(NRZ,N)/power_QPSK;
%       MQAM(k,:) = qammod(NRZ_QAM,M_QAM);          %0-15范围的十进制数调制，第二个参数是M
        MQAM(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
        % 步骤四：编写子函数，产生直流偏置，把直流偏置叠加到时域信号，并截取负信号，让剩余负信号变成 0
        QPSK_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);
    end
    QPSK_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) == 0
        x((k_r+1)/2) = -1;
    else
        x((k_r+1)/2) = 1;
    end
end
for k_i = 2:2:N             %虚部
    if NRZ(k_i) == 0
        y(k_i/2) = -1;
    else
        y(k_i/2) = 1;
    end
end
z = x+y*sqrt(-1);
s = z(1:end-1);
QAM = s;
end

function [QAM_H] = Hermitain_sym(QAM)
%施密特变换
QAM_H = [0 QAM 0 flip(conj(QAM))];%flip是倒序排列，conj是取共轭
end

function [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
end

function [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) < 0
        QAM_after_DC(k_DC) = 0;
    end
end
end

function [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];
end

function [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)    %添加高斯随机数，第二个参数是标准差因此要开方
end

function [QAM_decrease_CP] = Remove_CP(QAM_after_AWGN,CP)
%删去前缀CP
QAM_decrease_CP = QAM_after_AWGN(CP+1:end);
end

function [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;%保留两位小数用的
end

function [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)) > 0
        QAM_after_demodulation(count5*2-1) = 1;
    else
        QAM_after_demodulation(count5*2-1) = 0;
    end
    if imag(User(count5)) > 0
        QAM_after_demodulation(count5*2) = 1;
    else
        QAM_after_demodulation(count5*2) = 0;
    end
end
end

function [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);
end

function [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基本就是极限了，再低就是一个错误码元都没有了