% 比較了不同調(diào)制情況時的誤比特率
%  BPSK,QPSK,.....
%  QAM,16QAM.....
clear 
n_frame=20;
tblen=32; %回朔深度
 
M=4;    %1--BPSK    2--QPSK or QAM    3---8PSK   4---16PSK or 16QAM
          % psk        psk     qask        psk         psk     qask
Es=10;   % 數(shù)字映射后的符號能量 psk--1
        %                      16QAM--10
L_frame=10000-tblen;% 長度應為（M*a)/2-tblen a是任意正整數(shù)
start=-16;
step=2;
finish=20;
trel=poly2trellis(7,[171,133]); % 定義卷積碼 (2,1,2) 5,7
                                           % (2,1,6) 133,171
                                           % (2,1,7) 
r=0.5; % 編碼效率
Eb=Es/M; % 每比特能量
plot_pe=[];
Q=1;
axis_EbN0=start:step:finish;
%******************************************************************8

for EbN0=start:step:finish
    Liner_EbN0=10^(EbN0/10);
    pe_number=0;
    var=0.5*(Eb/Liner_EbN0)/r;  % 噪聲方差, 來源于a^2/(2*delta^2)=Liner_EbN0,a^2信息幅值的平方即比特能量, 2*delta^2=N0.    N0/r得到編碼后的噪聲功率譜密度
    for i=1:1:n_frame
        x_msg_p=randint(L_frame+tblen,1,2); %產(chǎn)生偽信息
        x_msg=x_msg_p(1:1:L_frame);     %真正的信息
        x_code_msg=convenc(x_msg_p,trel);
        symbol=bi2symbol(x_code_msg,M); % 二進制到十進制轉(zhuǎn)換
        
        modulation_sig=modmap(symbol,1,1,'qask',2^M); % 數(shù)字映射
        
        x_ch=g_ch(modulation_sig,var);  % 經(jīng)過高斯信道，噪聲方差var
        
        de_symbol=demodmap(x_ch,1,1,'qask',2^M); %逆映射
        
        rec=symbol2bit(de_symbol,M); % 十進制轉(zhuǎn)換二進制
        x_decode_msg=vitdec(rec,trel,tblen,'cont','hard');
        x_decode_msg=x_decode_msg(tblen+1:1:L_frame+tblen); % 取出有用信息
        
        pe_number=pe_number+sum(x_msg~=x_decode_msg);
        
        current_time=fix(clock);
       fprintf('i am working %g, total(%g)      %g年  %g月  %g日  %g時  %g分  %g秒\n\n',Q,(((finish-start)/step)+1)*n_frame,current_time(1),current_time(2),current_time(3),current_time(4),current_time(5),current_time(6))
        Q=Q+1;
        fprintf('\n\n')
        
    end
    pe=pe_number/(L_frame*n_frame);
    plot_pe=[plot_pe,pe];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
semilogy(axis_EbN0,plot_pe,'b*-')

xlabel('Eb/N0 信號功率/噪聲功率(db)')
ylabel('BER')
title('誤碼率')
grid on