dpsk_2_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=2;                                   %  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,:)*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;

    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(1,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/2);                                 %  0 decision. Be careful!!! There are many skills I don't know before.
    index2=find(culcu_abs_rad>3*pi/2);
    decs_bin(:,index1)=repmat([0]',1,length(index1));
    decs_bin(:,index2)=repmat([0]',1,length(index2));
    index1=find(culcu_abs_rad>pi/2);                                 %  1 decision
    index2=find(culcu_abs_rad(index1)<3*pi/2);
    decs_bin(:,index1(index2))=repmat([1]',1,length(index2));

    
    % END binary decision above
    decs_data=decs_bin(1,:)*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/1;
semilogy(SNRindB1,pb_from_ps,'r');                                   %  pb_from_ps is probability of error bit by directly culculating pb
grid on;
toc