基于MATLAB的光纤通信仿真

Posted 2021-07-07 fpga&matlab

1.问题描述：

    ·传感器接收到的信息转换为二进制信息（这个我们可以直接模拟出随机的二进制序列来表示传感器的数据，首先设置一组数据，然后进行量化，最后进行转换为二进制。这个部分的仿真结果对应你提供论文的3.2.1部分）

  ·封装，将编码后的数据进行封装，即将数据变为数据包，再把数据包转换为数据帧

     ·码率调整，主要将瞬时数码率不同的多个数码流调整为具有同一较高的数码率的数码流，然后再进行时分复接。

  ·最后由时分复接得到传感复用帧结构。

  ·码分编码封装发送(二进制数据编码，这里，编码方式我们选择性能非常好的LDPC编码)

在发送端，这里根据仿真要求，模拟五个子网

接收端包括：码分解码，解封装，分接，码率恢复，解帧，拆包，译码。。。。。

       最后我们对整个系统进行误码率的分析，对五个子网的数据进行误码率分析，眼图以及Q因子分析。

2.部分程序：

 

clc;
clear;
close all;
warning off;
pack;
addpath 'func\\'

%是否显示中间结果的波形
%fig1:模拟传感器采集到的数据
%fig2:模拟传感器采集到的数据的净荷值
%fig3:五路传感器数据
%fig4:解封装sync位置
%fig5:显示标准眼图
%fig6:显示噪声干扰眼图
%fig7:显示译码之后眼图
is_show = [0 0 0 0 0 0 0];

Ind     = 0;
%参数初始化
EbN0    = [0:0.5:6];
TIMES   = 10;

Err_Rate0  = zeros(1,length(EbN0));
Err_Rate0s = zeros(1,TIMES);
Err_Rate1  = zeros(1,length(EbN0));
Err_Rate1s = zeros(1,TIMES);
Err_Rate2  = zeros(1,length(EbN0));
Err_Rate2s = zeros(1,TIMES);
Err_Rate3  = zeros(1,length(EbN0));
Err_Rate3s = zeros(1,TIMES);
Err_Rate4  = zeros(1,length(EbN0));
Err_Rate4s = zeros(1,TIMES);
Err_Rate5  = zeros(1,length(EbN0));
Err_Rate5s = zeros(1,TIMES);

Q1         = zeros(1,length(EbN0));
Q2         = zeros(1,length(EbN0));
Q3         = zeros(1,length(EbN0));
Q4         = zeros(1,length(EbN0));
Q5         = zeros(1,length(EbN0));
Q0         = zeros(1,length(EbN0));

for jj1 = 1:length(EbN0)
    jj1
    Ind = Ind + 1;
    Err_Rate2 = zeros(1,TIMES);
    for jj2 = 1:TIMES
        jj2
        Num_NET = 5;%模拟子网的个数
        Lens    = 500;%传感器采集数据长度
        %传感信息的净荷值相关参数
        A       = -40;
        B       = 120;
        Delta   = 0.2;
        Ld      = (ceil(log2((B-A)/Delta)));
        %传感子网数据帧的封装
        Preamble=[0,1,1,1,1,1,1,0];%帧前导码取特殊字符串01111110，表示帧同步，便于传感监控中心判断出帧的起始位置，
        %子网ID号，本案为5个子网，所以ID为三bit数据,但是为了具有扩展性，ID用四个bit表示
        ID      =[0,0,0,0;
                  0,0,0,1;
                  0,0,1,0;
                  0,0,1,1;
                  0,1,0,0];
        %HEC为帧头校验码
        HEC     = [0,1,0,1]; 
        %包头
        Head    = [0,0,1,1,0,0,1,1];
        %address,地址，随着采集数据，地址逐渐加1
        %传感器类型,假设五个子网，每个子网就一种类型
        type1   = [0,0,0,1;
                   0,0,1,0;
                   0,0,1,1;
                   0,1,0,0;
                   0,1,0,1];
        %包尾
        Trail   = [1,1,0,0,1,1,0,0];
        Bao_Size= 10;

        %最后的复用帧的相关参数
        Sync      = [0,0,0,0,0,0,1,1,1,0,1,0,1,1,1,1,1,0,0,0,1,0,1,1,0,0,0,0,0,1,0,1,1,1,0,1,0,1,1,1,1,1,0,1,1,1,0,1,0,0,0,0,0,1,0,1,1,1,0,1,0,1,0,0];
        Separator = [0,1,0,1];
        FCS       = [1,0,0,1,1,0,1,1,1,0,0,1,1,0,0,1];
        N         = 5;
        M         = 5;

        %第一部分，系统发送部分
        %第一部分，系统发送部分

        %本案共有五个子网，所以模拟五组数据
        for nn = 1:Num_NET
            %模拟产生传感器采集到的数据,这里，为了模拟不同的码率，将采集的数据长度设置为不同
            y     = func_sensor_samples(Lens,A,B);
            YY{nn} = y;
            if is_show(1) == 1
               figure(1);
               plot(y,'b');
               xlabel('times');
               ylabel('amplitude')
               title('模拟传感器采集到的数据');
               pause(0.1);
            end
            %光纤光栅传感信号的净荷值的二进制转换
            Lp = func_value2realvalue(y,A,B,Delta);
            if is_show(2) == 1
               figure(2);
               plot(Lp,'b');
               xlabel('times');
               ylabel('amplitude')
               title('模拟传感器采集到的数据的净荷值');
               pause(0.1);
            end
            %计算转换为二进制
            V2 = func_dec2bin(Lp,Ld);
            %将二进制转换为串行流
            Signal{nn} = (reshape(V2',size(V2,1)*size(V2,2),1))';%通过上面的步骤，我们模拟了实际要发送的二进制码流
        end
        if is_show(3) == 1
           figure(3);
           subplot(321);
           plot(YY{1});title('第1个子网传感器采集数据波形');
           subplot(322);
           plot(YY{2});title('第2个子网传感器采集数据波形');
           subplot(323);
           plot(YY{3});title('第3个子网传感器采集数据波形');
           subplot(324);
           plot(YY{4});title('第4个子网传感器采集数据波形');
           subplot(325);
           plot(YY{5});title('第5个子网传感器采集数据波形');
           pause(0.1);
        end

        %Signal即为实际要发送的二进制码流
        %产生数据包数据包,这里，我们假设每个数据包中为10个采样数据，即100bit数据
        address = 0;
        for nn = 1:Num_NET
            tmps2   = [];
            tmps    = [];
            address = 0;
            for i = 1:length(Signal{nn})/(Bao_Size*Ld)
                address  = address + 1;
                address2 = func_dec2bin_addr(address);
                %地址转换为2进制，8bit标识,这样，每个包共128bit
                tmps     = [Head,address2,type1(nn,:),Signal{nn}((i-1)*Bao_Size*Ld+1:i*Bao_Size*Ld),Trail];
                Len_bao  = length(tmps);
                tmps2    = [tmps2 tmps];
            end
            Signal2{nn} = tmps2;
        end
        %由数据帧转换为传感帧
        for nn = 1:Num_NET
            tmps2   = [];
            tmps    = [];
            for i = 1:length(Signal2{nn})/(N*Len_bao)
                tmps     = [Preamble,ID(nn,:),HEC,Signal2{nn}((i-1)*Len_bao*N+1:i*Len_bao*N)];
                Len_zhen = length(tmps);
                tmps2    = [tmps2 tmps];
            end
            Signal3{nn} = tmps2;
        end
        %对五个子网的数据进行时分复接
        %码率调整
        Out1 = func_Rate_sync(Signal3{1});
        Out2 = func_Rate_sync(Signal3{2});
        Out3 = func_Rate_sync(Signal3{3});
        Out4 = func_Rate_sync(Signal3{4});
        Out5 = func_Rate_sync(Signal3{5});
        %时分复接
        Out  = func_tdma(Out1,Out2,Out3,Out4,Out5,Len_zhen);
        %传感复用帧
        tmps2   = [];
        tmps    = [];
        for i = 1:length(Out)/(M*Len_zhen)
            tmps      = [Sync,Separator,Out((i-1)*M*Len_zhen+1:i*M*Len_zhen),Separator,FCS];
            Len_zhen2 = length(tmps);
            tmps2     = [tmps2 tmps];
        end
        Signal4 = tmps2;

        %进行编码，编码矩阵为96*192，每次取80,并补充16个0作为虚拟填充符进行编码，共96个数据进行编码
        load H;
        load G;
        Trans_Signal = [];
        for i = 1:length(Signal4)/80
            Trans_data   = [zeros(1,16),Signal4((i-1)*80+1:(i)*80)];
            %编码
            data_code    = mod(Trans_data*G,2); 
            %BPSK
            Trans_BPSK   = 2*data_code-1; 
            Len_code     = length(Trans_BPSK);
            Trans_Signal = [Trans_Signal,Trans_BPSK];
        end

        %第二部分，信道，信道部分，由于是光纤传输，且您论文中没有详细说明这个部分内容，所以这里暂时仅考虑高斯白噪声
        %第二部分，信道，信道部分，由于是光纤传输，且您论文中没有详细说明这个部分内容，所以这里暂时仅考虑高斯白噪声
        sigma      = sqrt(1./(10^(EbN0(Ind)/10)*0.5)); 
        Rec_Signal = Trans_Signal + sigma*randn(1,length(Trans_Signal));   

        %第三部分，监控中心——接收端
        %第三部分，监控中心——接收端
        %译码
        R_Signal = [];
        for i = 1:length(Signal4)/80
            z_hat       = func_Dec(Rec_Signal((i-1)*Len_code+1:i*Len_code),sigma,H,50);
            x_hat       = z_hat(size(G,2)+1-size(G,1):size(G,2));
            R_Signal    = [R_Signal,x_hat(17:end)'];%去掉16个填充符
        end

        
        if is_show(5) == 1
           st = func_eye(Trans_Signal);
           eyediagram(st,40)
           ylabel('信号幅度');
           title('原始信号眼图');     
           pause(0.1);
        end        
        if is_show(6) == 1
           st = func_eye(Rec_Signal);
           eyediagram(st,40)
           ylabel('信号幅度');
           title('接收带噪声干扰信号眼图');  
           pause(0.1);
        end          
        if is_show(7) == 1
           st = func_eye(R_Signal*2-1);
           eyediagram(st,40)
           ylabel('信号幅度');
           title('译码之后信号眼图'); 
           pause(0.1);
        end         
        
        
        %解封装
        %通过搜索sync来确定每帧的帧头
        [R_Signal,Sync_pos] = func_find_sync(R_Signal,Sync);
        if is_show(4) == 1
           figure(4);
           plot(Sync_pos,'b');
           hold on;
           plot(find(Sync_pos>=40),Sync_pos(find(Sync_pos>=40)),'r*');
           hold off;
           xlabel('times');
           ylabel('amplitude')
           title('解封装sync位置');
           pause(0.1);
        end
        %判断当前帧位置能否检测到帧头
        if_sync = zeros(1,length(Out)/(M*Len_zhen));%这个变量用来存放是否检测到当前帧，如果为0，则该帧未检测到，直接丢弃，检测下一帧
        for i = 1:length(Out)/(M*Len_zhen)
            if Sync_pos(1+(i-1)*Len_zhen2) > 40
               if_sync(i)=1;
            end
        end

        %解封装
        tmps3   = [];
        tmps2   = [];
        tmps    = [];
        for i = 1:length(Out)/(M*Len_zhen)
            if if_sync(i) == 1
               tmps = R_Signal((i-1)*Len_zhen2+1:(i)*Len_zhen2);
            else
               tmps = zeros(1,Len_zhen2);% 如果该帧没有检测到，那么直接赋值0
            end
            %解封装
            tmps2 = tmps(length(Sync)+length(Separator)+1:M*Len_zhen+length(Sync)+length(Separator));
            tmps3 = [tmps3 tmps2];
        end
        Outs = tmps3;

        %数字分接单元
        [I1,I2,I3,I4,I5] = func_tdma2(Outs,Len_zhen,M);
        %对五个子网的数据进行时分复接
        %码率调整
        Outr{1} = func_Rate_sync(I1);
        Outr{2} = func_Rate_sync(I2);
        Outr{3} = func_Rate_sync(I3);
        Outr{4} = func_Rate_sync(I4);
        Outr{5} = func_Rate_sync(I5);

        %将5路传感帧转换为数据帧，进一步拆封
        %由数据帧转换为传感帧
        for nn = 1:Num_NET
            tmps2   = [];
            tmps    = [];
            LL = length(Preamble)+length(ID(nn,:))+length(HEC);
            for i = 1:length(Signal2{nn})/(N*Len_bao)
                tmps     = Outr{nn}((i-1)*Len_zhen+1:i*Len_zhen);
                tmps2    = [tmps2 tmps(LL+1:end)];
            end
            Outr2{nn} = tmps2;
        end
        %将数据包中数据提取
        %Signal即为实际要发送的二进制码流
        %产生数据包数据包,这里，我们假设每个数据包中为10个采样数据，即100bit数据
        address = 0;
        for nn = 1:Num_NET
            tmps2   = [];
            tmps    = [];
            LL      = length(Head)+length(address2)+length(type1(nn,:));
            for i = 1:length(Signal{nn})/(Bao_Size*Ld)
                tmps  = [Outr2{nn}((i-1)*Len_bao+1:i*Len_bao)];
                tmps2 = [tmps2 tmps(LL+1:end-length(Trail))];
            end
            Outr3{nn} = tmps2;
        end

        %误码率统计在转换为十进制之前计算
        [nberr1,rat1] = biterr(Outr3{1},Signal{1});
        [nberr1,rat2] = biterr(Outr3{2},Signal{2});
        [nberr1,rat3] = biterr(Outr3{3},Signal{3});
        [nberr1,rat4] = biterr(Outr3{4},Signal{4});
        [nberr1,rat5] = biterr(Outr3{5},Signal{5});

        Err_Rate0(jj2) = (rat1+rat2+rat3+rat4+rat5)/5;   
        Err_Rate1(jj2) = rat1;   
        Err_Rate2(jj2) = rat2;   
        Err_Rate3(jj2) = rat3;   
        Err_Rate4(jj2) = rat4;   
        Err_Rate5(jj2) = rat5;   
    end
    Err_Rate0s(Ind) = mean(Err_Rate0);
    Err_Rate1s(Ind) = mean(Err_Rate1);
    Err_Rate2s(Ind) = mean(Err_Rate2);
    Err_Rate3s(Ind) = mean(Err_Rate3);
    Err_Rate4s(Ind) = mean(Err_Rate4);
    Err_Rate5s(Ind) = mean(Err_Rate5);
    Q0(Ind)         = sqrt(2)*erfcinv(2*Err_Rate0s(Ind)/10);
    Q1(Ind)         = sqrt(2)*erfcinv(2*Err_Rate1s(Ind)/10);
    Q2(Ind)         = sqrt(2)*erfcinv(2*Err_Rate2s(Ind)/10);
    Q3(Ind)         = sqrt(2)*erfcinv(2*Err_Rate3s(Ind)/10);
    Q4(Ind)         = sqrt(2)*erfcinv(2*Err_Rate4s(Ind)/10);
    Q5(Ind)         = sqrt(2)*erfcinv(2*Err_Rate5s(Ind)/10);
    
    disp('over a cycle');
    pause(5);
end

figure;
semilogy(EbN0,Err_Rate0s,'b-^');
xlabel('EbN0');
ylabel('Ber');
grid on;

figure;
subplot(151);
semilogy(EbN0,Err_Rate1s,'b-^');
xlabel('EbN0');
ylabel('Ber');
title('子网1误码率');
grid on;
subplot(152);
semilogy(EbN0,Err_Rate2s,'b-^');
xlabel('EbN0');
ylabel('Ber');
title('子网2误码率');
grid on;
subplot(153);
semilogy(EbN0,Err_Rate3s,'b-^');
xlabel('EbN0');
ylabel('Ber');
title('子网3误码率');
grid on;
subplot(154);
semilogy(EbN0,Err_Rate4s,'b-^');
xlabel('EbN0');
ylabel('Ber');
title('子网4误码率');
grid on;
subplot(155);
semilogy(EbN0,Err_Rate5s,'b-^');
xlabel('EbN0');
ylabel('Ber');
title('子网5误码率');
grid on;

figure;
plot(EbN0,Q0,'b-o');
xlabel('EbN0');
ylabel('Q');
title('Q因子');
grid on;

3.仿真结论：

A01-68