蒙特卡洛法仿真OOK的NRZ与RZ波形

clear all
close all
q = 1.6e-19;%电子电量
Ib = 202e-6;%平均光电流
N0 = 2*q*Ib;%单边噪声功率谱密度
R = 1;%光电探测器灵敏度
Rb = 1e6;%比特率
Tb = 1/Rb;%码元时间

sig_length = 1e5;%信号长度
nsamp = 10;%每个信号采样点数
Tsamp = Tb/nsamp;%采样间隔
EbN0 = 1:12;%信噪比
SNR = 10.^(EbN0./10);
h = waitbar(0,'计算进度');
%%
%************ simulation of probabbility of errors ****************
for i = 1:length(SNR)
    waitbar(i/length(SNR));%显示计算进度
    p_avg(i) = sqrt(N0*Rb*SNR(i)/(2*R^2));%平均传输光功率
    i_peak(i) = 2*R*p_avg(i);%峰值光电流
    Ep(i) = i_peak(i)^2*Tb;%峰值光能量
    sigma(i) = sqrt(N0/2/Tsamp);%1/Tsamp其实是带宽

    
    
    %%
    %NRZ
    pt = ones(1,nsamp)*i_peak(i);%默认每个信号的采样值都是一样的
    rt = pt;%匹配滤波器
    OOK = round(rand(1,sig_length));%生成ook信号
    Tx_signal = reshape((pt.'*OOK),1,sig_length*nsamp);%发送信号（长度为所有采样点数）
    Rx_signal = Tx_signal+sigma(i)*randn(1,sig_length*nsamp);%匹配滤波器接收信号

    MF_out = conv(Rx_signal,rt)*Tsamp;%卷积，但是为什么要乘采样间隔，或者说除以带宽  匹配滤波器输出
    MF_out_downsamp = MF_out(nsamp:nsamp:end);%nsamp开始，每隔nsamp采样一次直到结束
    MF_out_downsamp = MF_out_downsamp(1:sig_length);%限定长度  截断
    
    Vd_nrz = Ep(i)/2;%判决门限为码元能量

    %%
    %RZ
    pt_rz  = [ones(1,nsamp/2),zeros(1,nsamp/2)]*i_peak(i);%RZ波形
%     rt_rz = pt_rz;%半波长匹配滤波器接收将不再等于原成型滤波
    rt_rz = [zeros(1,nsamp/2),ones(1,nsamp/2)]*i_peak(i);%RZ波形的话，作为匹配滤波器要以T/2为中心对称一下
    Tx_signal_rz = reshape((pt_rz.'*OOK),1,sig_length*nsamp);%发送信号（长度为所有采样点数）
    Rx_signal_rz = Tx_signal_rz+sigma(i)*randn(1,sig_length*nsamp);%匹配滤波器接收信号
    MF_out_rz = conv(Rx_signal_rz,rt_rz)*Tsamp;%卷积， 匹配滤波器输出 详情参考匹配滤波器最大输出值2A^2*TB/nsamp,这里保留了2A^2
    MF_out_downsamp_rz = MF_out_rz(nsamp:nsamp:end);%nsamp开始，每隔nsamp采样一次直到结束
    MF_out_downsamp_rz = MF_out_downsamp_rz(1:sig_length);%限定长度  截断
    Vd_rz = Ep(i)/4;%由表可知，此时的门限是码元能量平均值也就是P/4
    %%
    %蒙特卡洛法判决误码率
    error = 0;error_rz = 0;%归零
    
    for j0 = 1:sig_length
        if MF_out_downsamp_rz(j0) > Vd_rz
            R_rz = 1;
        else
            R_rz = 0;
        end
        if(R_rz~=OOK(j0))
            error_rz = error_rz+1;
        end
    end
    
    for j = 1:sig_length
        if MF_out_downsamp(j) > Vd_nrz
            R0 = 1;
        else
            R0 = 0;
        end
        if(R0~=OOK(j))
            error = error+1;
        end
    end
    ber(i) = error/sig_length;%误码率仿真值
    ber_rz(i) = error_rz/sig_length;
end

%%
%画图
figure;
semilogy(EbN0,ber,'b-*','LineWidth',2);
hold on;
ber_theory = qfunc(sqrt(SNR));%理论值计算误码率
ber_theory_rz = qfunc(sqrt(SNR*0.5));%顶着双倍的带宽了相当于，所以信噪比比别人小很合理吧
semilogy(EbN0,ber_theory,'r-*','LineWidth',2);
hold on;
semilogy(EbN0,ber_rz,'g-*','LineWidth',2);
hold on;
semilogy(EbN0,ber_theory_rz,'y-*','LineWidth',2);
grid on;


legend('NRZ数值解','NRZ解析解','RZ数值解','RZ解析解');
xlabel('Eb/N0  dB');
ylabel('BER');
title('OOK-NRZ调制BER曲线');

哈工深的可见光通信实验二，自己写的