dpsk_8_plot.m

4DPSK, 8DPSK在高斯信道下性能仿真. 解调方式不仅包括包括相干解调,还有个人编写的非相干解调,对于急需非相干解调的很有帮助

tic
clear;
SNRindB1=0:1:15;                       %  ready for plotting practical error bit/symbol rario
SNRindB2=0:0.1:15;                     %  ready for plotting theoretical error bit/symbol rario
M=8;                                   %  8 DPSK
k=log2(M);                          
N=10000;                               %  N is the number of octal random integer
N_bit=N*k;                             %  N_bit is the the total numbers binary bit
for SNRindB=1:length(SNRindB1)
    data_bin=rand(k,N)>0.5;            %  generating 3*N 0-1 array
    sym_oct=data_bin(1,:)*4+data_bin(2,:)*2+data_bin(3,:)*1;        %  N octal random integers
    %  Gray mapping
    idx=find(sym_oct==0);
    rad_init(idx)=0/M*2*pi;
    idx=find(sym_oct==1);
    rad_init(idx)=1/M*2*pi;
    idx=find(sym_oct==2);
    rad_init(idx)=3/M*2*pi;
    idx=find(sym_oct==3);
    rad_init(idx)=2/M*2*pi;
    idx=find(sym_oct==4);
    rad_init(idx)=7/M*2*pi;
    idx=find(sym_oct==5);
    rad_init(idx)=6/M*2*pi;
    idx=find(sym_oct==6);
    rad_init(idx)=4/M*2*pi;
    idx=find(sym_oct==7);
    rad_init(idx)=5/M*2*pi;
    %  Differential encoding
    dif_rad(1)=rad_init(1);                                    %  initiating the first differential angle
    for j=2:N
        dif_rad(j)=mod((rad_init(j)+dif_rad(j-1)),2*pi);
    end
    %  I and Q
    data_I=cos(dif_rad);
    data_Q=sin(dif_rad);
    data=data_I+i*data_Q;
    %  generating noise
    snr_in_dec=10.^((SNRindB-1)/10);
    sigma=sqrt(1/(2*k*snr_in_dec));                            %  culculation fomula of sigma
%     [noise_I noise_Q]=gngauss(N,sigma);
%     noise_I_Q=[noise_I;noise_Q];                               %  noise is 2*N gauss random number array by sigma
%     noise=noise_I+i*noise_Q;
%     data_add_noise=data+noise;
data_add_noise=gngauss(N,data,sigma)
    %  data_differential + noise_gauss

    %  Next, we begin decision
    decs_data=zeros(1,N);                                      %  decs_data is for containing octal data
    decs_bin=zeros(3,N);                                       %  decs_bin is for containing binary data
    decs_rad=angle(data_add_noise);                            %  culculating the differential phase
    culcu_rad(1)=decs_rad(1);
    for j=2:N
        culcu_abs_rad(j)=mod((decs_rad(j)-decs_rad(j-1)),2*pi);
    end                                                        %  in this time, the cul_rad is the absolute phase angle added noise
    
    index1=find(culcu_abs_rad<pi/8);                                 %  0 decision. Be careful!!! There are many skills I don't know before.
    index2=find(culcu_abs_rad>15*pi/8);
    decs_bin(:,index1)=repmat([0 0 0]',1,length(index1));
    decs_bin(:,index2)=repmat([0 0 0]',1,length(index2));
    index1=find(culcu_abs_rad>pi/8);                                 %  1 decision
    index2=find(culcu_abs_rad(index1)<3*pi/8);
    decs_bin(:,index1(index2))=repmat([0 0 1]',1,length(index2));
    index1=find(culcu_abs_rad>3*pi/8);                               %  3 decision
    index2=find(culcu_abs_rad(index1)<5*pi/8);
    decs_bin(:,index1(index2))=repmat([0 1 1]',1,length(index2));
    index1=find(culcu_abs_rad>5*pi/8);                               %  2 decision
    index2=find(culcu_abs_rad(index1)<7*pi/8);
    decs_bin(:,index1(index2))=repmat([0 1 0]',1,length(index2));
    index1=find(culcu_abs_rad>7*pi/8);                               %  6 decision
    index2=find(culcu_abs_rad(index1)<9*pi/8);
    decs_bin(:,index1(index2))=repmat([1 1 0]',1,length(index2));
    index1=find(culcu_abs_rad>9*pi/8);                               %  7 decision
    index2=find(culcu_abs_rad(index1)<11*pi/8);
    decs_bin(:,index1(index2))=repmat([1 1 1]',1,length(index2));
    index1=find(culcu_abs_rad>11*pi/8);                              %  5 decision
    index2=find(culcu_abs_rad(index1)<13*pi/8);
    decs_bin(:,index1(index2))=repmat([1 0 1]',1,length(index2));
    index1=find(culcu_abs_rad>13*pi/8);                              %  4 decision
    index2=find(culcu_abs_rad(index1)<15*pi/8);
    decs_bin(:,index1(index2))=repmat([1 0 0]',1,length(index2));
    % END binary decision above
    decs_data=decs_bin(1,:)*4+decs_bin(2,:)*2+decs_bin(3,:)*1;       %  N octal random integers which had been decided
    num_of_err_bit=length(find(xor(data_bin,decs_bin)));             %  culculating the total number of error bit
    num_of_err_symbol=length(find(sym_oct-decs_data));               %  culculating the total number of error symnol
    pb(SNRindB)=num_of_err_bit/N_bit;
    ps(SNRindB)=num_of_err_symbol/N;
end
semilogy(SNRindB1,pb,'k*');
hold on;
semilogy(SNRindB1,ps,'b');
pb_from_ps=ps/3;
semilogy(SNRindB1,pb_from_ps,'r');                                   %  pb_from_ps is probability of error bit by directly culculating pb
grid on;
toc