% QPSK系統仿真：AWGN信道，BER

clear;
SNRindB1=0:2:8;                                     % 設置SNR范圍
SNRindB2=0:0.1:8;
for ii=1:length(SNRindB1),

counter=0;                                          % 計數器置0
numofsymbolerror=0;
numofbiterror=0;
while (numofbiterror<100)

N=1000;
E=1;                                                % 碼元能量
snr=10^(SNRindB1(ii)/10);                           % 信噪比
sgma=sqrt(E/snr)/2;                                 % 噪聲功率
s00=[1 0]; s01=[0 -1]; s11=[-1 0]; s10=[0 1];       % 雙比特碼元映射

% 信源產生
for i=1:N,
    temp=rand;                                      % 生成取值在0和1之間的服從均勻分布的隨機變量
    if (temp<0.25),                                 % 信源輸出00，10，11，01的概率均為1/4
        dsource1(i)=0; dsource2(i)=0;
        tt(ii,i)=complex(s00(1),s00(2));
    elseif (temp<0.5),                              
        dsource1(i)=0; dsource2(i)=1;
        tt(ii,i)=complex(s01(1),s01(2));
    elseif (temp<0.75),                             
        dsource1(i)=1; dsource2(i)=0;
        tt(ii,i)=complex(s10(1),s10(2));
    else                                            
        dsource1(i)=1; dsource2(i)=1;
        tt(ii,i)=complex(s11(1),s11(2));
    end;
end;

% 差錯檢測及誤碼率計算

for i=1:N,
        mm=raylrnd(0.7); 
    n=sgma*randn(1,2);                              % 生成兩個服從N(0,sgma)的隨機變量
    if ((dsource1(i)==0) & (dsource2(i)==0)), 
        r=mm*s00+n;                                    % 加入噪聲
    elseif ((dsource1(i)==0) & (dsource2(i)==1)), 
        r=mm*s01+n;
    elseif ((dsource1(i)==1) & (dsource2(i)==0)), 
        r=mm*s10+n;
    else 
        r=mm*s11+n;
    end;
    
    % 判決
    rr(ii,i)=complex(r(1),r(2));
    c00=dot(r,s00); c01=dot(r,s01); c10=dot(r,s10); c11=dot(r,s11);
    c_max=max([c00,c01,c10,c11]);
    if (c00==c_max), decis1=0; decis2=0;
    elseif (c01==c_max), decis1=0; decis2=1;
    elseif (c10==c_max), decis1=1; decis2=0;
    else decis1=1; decis2=1;
    end;
    % 計算誤碼率及誤符號率
    
    symbolerror=0;
    if (decis1~=dsource1(i)), numofbiterror=numofbiterror+1; symbolerror=1; 
    end;
    if (decis2~=dsource2(i)), numofbiterror=numofbiterror+1; symbolerror=1; 
    end;
    if (symbolerror==1), numofsymbolerror=numofsymbolerror+1; 
    end;
end;
counter=counter+1;
end;
ps=numofsymbolerror/(N*counter);                    
pb=numofbiterror/(2*N*counter);                     

    smld_bit_err_prb(ii)=pb;
    smld_symbol_err_prb(ii)=ps;
    disp([ps,pb]);
end;

% 理論值
for i=1:length(SNRindB2),
    SNR=10^(SNRindB2(i)/10);                        % 信噪比
    theo_err_prb(i)=0.5*erfc(sqrt(SNR));            % AWGN信道下的誤碼率理論值
    theo_err_ray(i)=0.5*(1-1/sqrt(1+1/SNR));        % Rayleigh衰落信道下的誤碼率理論值
end;

% 畫圖
figure(1)                                           % BER圖
semilogy(SNRindB1,smld_bit_err_prb,'bs');
% legend('AWGN Simulation');
hold on;
semilogy(SNRindB1,smld_symbol_err_prb,'ro');
% legend('Rayleigh Simulation');
hold on;
semilogy(SNRindB2,theo_err_prb,'b');
% legend('AWGN Theory');
hold on;
semilogy(SNRindB2,theo_err_ray,'r');
legend('AWGN Simulation','Rayleigh Simulation','AWGN Theory','Rayleigh Theory');
title('Error performance under the AWGN and Rayleigh fading channels')
grid on;
hold off;